Introduction

Hepatocellular carcinoma (HCC) is the fifth most frequent cancer worldwide, occupying approximately 90% of primary liver cancer,1 while it is also a leading cause of death related to major cancers globally.2 In contrast to the observed reduction in mortality from all other common cancers (eg, breast, lung, and prostate), the death rate from HCC continually increases by 2–3% each year.3 According to current studies, chronic hepatitis B virus (HBV) infection, chronic hepatitis C virus (HCV) infection, and non-alcoholic fatty liver disease (NAFLD) are the most prominent risk factors for HCC.4 In recent years, tumor markers have displayed superior clinical value in assisting the diagnosis of malignant tumors, and among them, serum alpha-fetoprotein (AFP) has been the most prevalent and essential indicator for detecting HCC and monitoring the therapeutic efficacy. However, despite its overall specificity for predicting HCC as high as 80% to 94%, the sensitivity is only 25% to 65%, and there still exists a controversy about the specific threshold of AFP for survival or recurrence.5 Also, as a typical tumor related to inflammation, HCC is characterized by immune evasion during the process of development and evolution,6 which is mediated by the accumulation of immunosuppressive cells in the tumor and the activation of multiple inhibitory receptor-ligand pathways.7 Hence, it is imperative to explore alternative biomarkers for diagnosing and predicting the effectiveness of HCC treatment, while also impeding immune evasion.

The Parvalbumin (PVALB) gene is positioned on the long arm of chromosome 22 (22q13),8 which encodes a protein belongs to the EF-hand protein superfamily, and depicted as a Ca2+ buffering protein.9,10 Proteins of this family exhibit a conserved helix-loop-helix structural motif, and their combined functions encompass the modulation of numerous vital cellular processes such as gene transcription, protein phosphorylation, nucleotide metabolism, and ion transport.11–13 Under normal conditions, as a calcium-binding protein, PVALB has a high affinity for Ca2+ and functions in a variety of organs. For example, in the central nervous system, PVALB exhibits high expression level in inhibitory GABAergic interneurons, thereby managing the excitability of pyramidal cells.14 In the kidney, PVALB also modulates the expression of human sodium-chloride cotransporter and responds to extracellular ATP signaling in distal tubule cells by regulating intracellular Ca2+ in epithelial cells lining tubular subpopulations of distal nephron.15 Furthermore, as reported in previous studies, PVALB may potentially contribute to the progression of malignant tumors. In thyroid cancer, overexpression of PVALB can reduce Ca2+ influx into mitochondria, leading to modification in mitochondrial morphology, increasing the amount of mitochondria and changing subcellular localization, along with the possibility of suppressing cell proliferation and triggering cell death through the AKT / GSK-3β pathway.16 Besides, PVALB has been identified to be down-regulated in glioma and glioma patients with high expression of PVALB have a favorable progression-free survival and overall survival.17 Nevertheless, the precise function and mechanism of PVALB in HCC, as well as its connection with the prognosis of HCC patients are still unclear, which appeal to us to further study it.

In this study, we have screened PVALB by Weighted Gene Co-expression Network Analysis (WGCNA) combined with Least Absolute Shrinkage and Selection Operator (LASSO) prognostic regression and assessed the expression of PVALB in normal and HCC tissue, the impact of PVALB itself on the existence of HCC sufferers, as well as the level of PVALB promoter methylation and its correlation with the survival of HCC sufferers by several databases and samples. What’s more, we conducted functional enrichment analyses of PVALB and probed the underlying mechanism of its role in HCC from varying perspectives such as tumor immune infiltration, Competing endogenous RNA (ceRNA) regulatory network and drug sensitivity.

Materials and Methods Data Collection and Processing

We obtained PVALB expression data from the The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov), including 374 LIHC samples and 50 normal samples. In addition, we also selected the GSE54236 dataset from the Gene Expression Omnibus database (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) which includes 81 tumor samples and 80 normal samples for our study.

Cell Culture

Human HCC cell line LM3 was transfected in 37°C DMEM (HyClone, Germany) and 10% fetal bovine serum (Gibco, USA). HCCLM3 was derived from the Cell Bank of Type Culture Collection of the Chinese Academy of Sciences and the Shanghai Institute of Cell Biology in China.

WGCNA

WGCNA performs weighted gene co-expression network analysis based on R packages, looking for interconnection among modules to facilitate the identification of candidate biomarkers or therapeutic targets.18 We have first prepared packages related to TCGA and GSE54236 cohorts. Then, β = 18 and β = 12 were selected as appropriate soft thresholds to merge similar modular genes by means of the hierarchical clustering dendrogram. Afterwards, from the result of the module-trait maps in TCGA and GSE54236 cohorts, we filtered out the purple and pink modules.

UALCAN Analysis

UALCAN (http://ualcan.path.uab.edu/) is an easy-to-use web that provides insightful analysis into the gene expression data of TCGA, which is instrumental in accelerating cancer research.19 In this work, we have explored the expression of PVALB in HCC tissue and normal tissue. Besides, UALCAN database was applied to assess the relevance of PVALB expression with clinicopathological features, as well as the association of PVALB promoter methylation level to various clinical features.

Database of Hepatocellular Carcinoma Expression Atlas (HCCDB) Analysis

The HCCDB (http://lifeome.net/database/hccdb) establishes a global differential gene expression atlas for HCC by analyzing genes that are differentially expressed in multiple datasets 20. With the help of the “Expression pattern” module in HCCDB, we performed a final validation of PVALB expression in HCC samples and adjacent samples.

MethSurv Analysis

MethSurv (https://biit.cs.ut.ee/MethSurv/) provides assistance to subscribers lacking coding ability, and dedicated to the preliminary evaluation of methylation-based cancer biomarkers.20 In our research, the methylation level of sites related to PVALB in HCC was visualized by “Gene visualization” module, and “Single CpG” module was also been adopted to examine the correlation of methylation level at the cg25625146 locus with HCC patients survival.

The Shiny Methylation Analysis Resource Tool (SMART) Analysis

SMART (http://www.bioinfo-zs.com/smartapp/) is a amiable and intuitive application according to the data from TCGA, and it could be applied to assist clients in the multidimensional exploration of DNA methylation in 33 cancer types.21 With SMART, in the “Home” module, we not only determined the distribution of methylation probes linked to PVALB on chromosomes, but also acquired more detailed information on the location of CpG genomes related to PVALB. Moreover, via the “Methylation DIY” and “Survival” modules, we evaluated the methylation level of PVALB at the cg25625146 locus in normal and HCC tissue, along with the association between the methylation level at this locus and the survivorship of HCC patients.

Kaplan Meier Plotter Analysis

Kaplan-Meier plotter database (https://kmplot.com/analysis/) introduces a web-based survival analysis tool.22 Here we initially inquired about the relation of PVALB expression with the survival of HCC sufferers through mRNA RNA-seq data, displaying survival conditions on the basis of OS (Overall Survival), PFS (Progression Free Survival), RFS (Recurrence free survival), and DSS (Disease-specific survival). More in depth, we also compared the prognostic impact of PVALB in HCC patients under varied enrichment of immune infiltrating cells.

LinkedOmics Analysis

LinkedOmics (http://www.linkedomics.org/ login.php) collects multi-omics data and clinical data from 11,158 patients with 32 cancer types, all from the TCGA project.23 With this tool, we identified the co-expressed genes of PVALB and carried out Gene Ontology (GO) pathway analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) route analysis in which pathway they are involved. Beyond that, LinkedOmics was also been utilized to create scatter plots for exploring the relationship between PVALB and CGR3A (CD16), NCAM1 (CD56) and FASLG (FasL).

TIMER Analysis

TIMER (https://cistrome.shinyapps.io/timer/) is a web server that allows a thorough study of the molecular characterization of tumor-immune interactions.24 In our study, the TIMER database was applied for investigating the relationship of PVALB expression with multiple immune cells in HCC. In addition, its upgraded TIMER2 database was further explored to learn the relevance between PVALB expression and distinct immune cell infiltration in pan-cancer under different algorithms.

TISIDB Analysis

TISIDB (http://cis.hku.hk/TISIDB/) leverages multi-omics data from TCGA to aid us in investigating the relevance of specific genes to immunity in the context of tumorigenesis.25 Building on this database, we examined the relation of PVALB expression with lymphocytes, MHCs, chemokines and chemokine receptors. In addition, the variation of PVALB expression among different subgroups was studied by “Subtype” module, as well as drugs targeting PVALB were analyzed by “Drug” module.

Constructing the ceRNA Network

ceRNA is a regulatory mechanism for genes. By means of the Targetscan database (http://www.targetscan.org/)26 and the DIANA (http://diana.imis.athena-innovation.gr/DianaTools/index.php?r=site/pageandview=software)27 database, we jointly predicted and filtered the target microRNA (miRNA) that was negatively correlated with the mRNA expression of PVALB, hsa-miR-6735-5p. Then, upstream of the miRNA-mRNA silencing mechanism, we further anticipated the long noncoding RNA (lncRNA) (lnc-YY1AP1-3) that might interact with hsa-miR-6735-5p using the LncRNASNP2 database (http://bioinfo.life.hust.edu.cn/LncRNASNP/),28 analyzing and modelling the HCC-associated comprehensive lncRNA-miRNA-mRNA (PVALB) ceRNA network.

GeneMANIA Analysis and Docking Analysis

Using GeneMANIA (http://www.genemania.org),29 a robust website for gene function analysis, we mapped the PPI functional network of PVALB and analyzed the proteins that physically interact with PVALB, and SREBF1 was selected through this method. The part of the “Mutations” in cBioPortal website (http://cbioportal.org)30 was used to characterize secondary structures of PVALB and SREBF1 (study ID, LIHC-TCGA-Firehouse Legacy). From the SWISS-MODEL Database (https://swissmodel.expasy.org/)31 (PDB ID: 1AM9), we retrieved SREBF1 advanced structures.31–33 Moreover, the sophisticated structure of PVALB was predicted by the AlphaFold Database (https://alphafold.com/)34 (AlphaFold ID: P20472). In the end, the HDOCK (http://hdock.phys.hust.edu.cn/) server predicted the interactive docking model of PVALB and SREBF1, and displayed them using PyMOL software.35

CTD (Comparative Toxicogenomics Database) Analysis

CTD (http://ctdbase.org/) creates an initiative digital ecosystem by linking toxicological information on chemicals, genes and diseases.36 Our study deployed cancer-related drugs interacting with PVALB from the part of the “Chemical Interactions” in CTD.

Gene Set Cancer Analysis (GSCA)

GSCA (http://bioinfo.life.hust.edu.cn/GSCA/#/) combines data from 4 public databases for 33 cancer types and carries out hypothesis checking at the genomic, pharmacogenomic and immunogenomic level.37 Employing the “Drug” module, we compared the drug sensitivity of PVALB and SREBF1 utilizing both Genomics of Drug Sensitivity in Cancer (GDSC) and The Cancer Therapeutics Response Portal (CTRP) methods.

Colony Formation Assay

After culturing the transfected cells for 48h, approximately 0.6×103 cells were carefully retrieved and seeded into 6-well plates until colonies of suitable dimensions were formed. Then, these cells were fixed with 4% paraformaldehyde for 30 minutes, followed by staining with 1.0% crystal violet for an equivalent period. Once visible colonies had formed, the number of colonies can be counted and tallied across 10 distinct fields.

Quantitative RT-PCR Analysis

Total mRNA was extracted using a standard Trizol-based protocol (Invitrogen, USA) and reverse transcription reactions were performed by the PrimeScript RT reagent kit (Invitrogen, USA). qPCR was then carried out with SYBR Premix Ex Taq (TaKaRa, China). Lastly, we conducted semi-quantitative analysis, and employed this method to detect mRNA expression of a range of cytokines.

Western Blot

Detailed instructions for the comprehensive preparation of total protein extraction for Western blot can be located in reference.38 To summarize, proteins were meticulously isolated on ice with RIPA buffer (Beyotime, Shanghai, China) from a mixture of protease and inhibitor compounds (Thermo Fisher Scientific, New York, USA). Following centrifugation, the protein concentration was closely monitored through employment of the BCA Protein Assay Kit (Thermo Scientific, Waltham, MA, USA), and the proteins in the sodium dodecyl sulfate (SDS) polyacrylamide gel electrophoresis were transferred to a polyvinylidene fluoride (PVDF) membrane via electroblotting. Thereafter, the primary antibody was applied and left for an overnight incubation at a temperature of 4°C, and then the closed membrane was washed three times with TBST while incubated with the secondary antibody under room temperature conditions for a duration of 1 hour. In the final step, the expression of the proteins was detected with electrochemiluminescence (ECL).

RNA Binding Protein Immunoprecipitation Assay

RIP assays were performed using the Magna RNA immune-precipitation kit (Millipore). In brief, following the use of RIP lysates for HCC cells, the lysates were incubated with anti-AGO2 (Cat#ab186733, Abcam) or IgG (Cat#ab172730, Abcam) antibody at 4°C until overnight. The above magnetic bead-antibody complexes were resuspended with 150ul Proteinase K Buffer. In the end, RNA was extracted by q-PCR.

Dual-Luciferase Reporter Gene System

The dual-luciferase reporter gene assay utilized firefly luciferase (Genepharma Company (Shanghai, China)) as the reporter gene, while the plasmid carrying the firefly luciferase gene (Genepharma Company (Shanghai, China)) was used as the control plasmid to co-transfect cells with the reporter gene. We cloned the binding site containing hsa-miR-6735-5p into pmiRGLO Vector (Genepharma Company (Shanghai, China)) (vector name). The plasmid was positively transfected into cells using Liposome 3000 and PLUS reagent (Genepharma Company (Shanghai, China)) (name of transfection reagent), and the gene expression was obtained by observing the fluorescein luminescence.

Cell Counting Kit-8(CCK-8) Assay

We determined HCC cell viability using the Cell Reagent Counting Kit-8 (manufacturer). Briefly, after incubating the cells in 96-well plates for 24 h, 10 ul of CCK8 solution was added to each well and incubated for 2 h. Subsequently, the optical density (OD) index at 450 nm was determined using an enzyme meter.

Statistical Analysis

This article utilizes R for analyzing all of the data. Initially, based on the R package “MCPcounter”, Our study computed the level of immune cell abundance and stromal cell infiltration in the TCGA-LIHC and GSE54236 cohorts. Next, the C-index was set as a filter, and we constructed the LASSO prognostic regression model by performing an analysis of the cleaned data to obtain the values of the variable lambda with the package “glmnet”. Except for the correlation chord plot, whose result was visualized with the “circlize” package, the “ggplot2” package served for all the other results. Meanwhile, the packages “stats” and “car” were designed to count the expression of PVALB in normal and HCC tissue, “pROC” was deployed for ROC analysis of the data, while the “clusterProfilerer” was employed for GO and KEGG analysis, and “msigdbr” package was additionally applied for gene set enrichment analysis (GSEA) analysis to refer the datasets. Apart from that, the correlation of PVALB with the three scores was done by the package “ESTIMATE”, while the immune infiltration of PVALB was derived from the single-sample gene set enrichment analysis (ssGSEA) algorithm provided in the R package “GSVA”. Eventually, the relationship between PVALB and microsatellite instability (MSI) was realized by the “ggstatsplot” package, and the immunotherapy efficacy was calculated by the tumor immune dysfunction and exclusion (TIDE) algorithm through the “ggpubr” package.

Results PVALB Was Identified as a Key Gene Associated with NK Cells in HCC Patients

In an effort to screen for critical genes capable of influencing the progression of HCC, we have performed a series of studies. Primarily, on the basis of all gene expression data in the TCGA and GSE54236 datasets related to HCC, we undertook differential expression analyses, and as presented in volcano Figure 1A–B, there were 12,471 and 9395 DEGs in the two cohorts respectively. Subsequently, with the MCP-counter algorithm, we calculated the abundance of all DEGs infiltrated in immune cells from both cohorts. Interestingly, it was revealed that most immune cells specifically NK cells were infiltrating at high abundance in the tumor microenvironment (Supplementary Figure 1A and D), Which invited us to exploit the WGCNA method, incorporating adaptive immune response and NK cells into the filter conditions, so that we can further hunt for molecular targets that can potentially inhibit the progression of HCC.

Figure 1 Profiling to mine potential biomarker for HCC. Volcano plots exhibited all differentially expressed genes in (A) TCGA and (B) GSE54236 cohorts. (C) Purple module in TCGA cohort and (D) pink module in GSE54236 cohort were filtered by the relationship analysis between modules and features. (E) Venn diagram displayed the 43 pivotal genes intersected by DEGs and modules. Screening of (F) LASSO regression coefficients and plotting of (G) prognostic LASSO variable trajectories based on TCGA-LIHC data.

In the WGCNA, depending on the scale-free R2 (R2 = 0.8), we set the soft thresholds in the TCGA and GSE54236 cohorts to 18 and 12 (Supplementary Figure 1B and E), and 29 and 20 modules were acquired by average hierarchical clustering and dynamic tree clipping (Supplementary Figure 1C and F). Afterwards, by examining the correlation of each module with HCC, NK cells and adaptive immune response, we believed that the purple module in the TCGA cohort and the pink module in the GSE54236 cohort represented the most significant modules, since both modules possessed the highest negative correlation with the “HCC” column as well as the high positive correlation with the “NK” and “GOBP_ADAPTIVE_IMMUNE_RESPONSE” columns (Figures 1C and D). Moreover, at the convergence of DEGs, purple module and pink module, we have found 43 hub genes (Figure 1E). What’s more, with the importation of these 43 pivotal genes into the LASSO prognostic regression model of TCGA-LIHC, we identified three key genes affecting the prognosis of HCC patients, which were DNASE1L3, LCAT and PVALB (Figure 1F and G). Among them, the roles of DNASE1L3 and LCAT have been widely reported in HCC, so we chose PVALB as the target of this paper.

PVALB Expression Was Down-Regulated in HCC

Target at investigating the differential expression of PVALB in HCC, we have employed a variety of tools. Initially, we assessed the mRNA expression of PVALB in an entire assortment of 26 different types of tumor tissue and normal tissue, and the outcome revealed that compared to normal tissue, PVALB expression was low in 17 different types of tumor tissue, in particular it was dramatically down-regulated in HCC (Supplementary Figure 2A). Next, 374 HCC samples and 50 normal samples were retrieved via the TCGA database, and then we used R software to draw differentiated expression plot and paired differentiated expression plot. The findings indicated that the HCC tissue showed markedly lower expression of PVALB than the normal tissue (Supplementary Figure 2B–C). what’s more, as Supplementary Figure 2D–E displayed, we analyzed PVALB expression using the GSE54236 dataset downloaded from the GEO database, as well as the data from TCGA and Genotype-Tissue Expression (GTEx) databases, which both obtained the same result consistently. In addition, making use of UALCAN online website, we further determined that the mRNA and protein expression of PVALB were notably down-regulated in HCC tissue compared to normal tissue (Supplementary Figure 2F–G). Finally, we also applied the HCCDB for verification, and the result also confirm the above (Supplementary Figure 2H). In conclusion, What we found in our study was that PVALB showed lower expression level in HCC tissue compared to normal tissue.

Relationship Between the Expression of PVALB and Clinicopathologic Features in HCC Patients

Via the UALCAN website, we geared to take the association of PVALB expression with clinicopathologic features into account for further elucidation of what impact PVALB has in HCC. According to the findings, PVALB expression varied significantly depending on gender and race (Supplementary Figures 3A and B). Additionally, we observed a link between PVALB and the cancer stages of HCC patients, with stage 2 showing an increased expression level than other stages (Supplementary Figure 3C). Furthermore, there were statistically significant differences among subgroups in several histological subtypes (Supplementary Figure 3D). Overall, the findings suggest that there exists a strong association between several clinicopathologic characteristics and PVALB expression.

Low Expression of PVALB Was Associated with a Worse Prognosis for HCC Patients

To discover the link between PVALB expression and prognosis in patients with HCC, we undertook a survival analysis according to the Kaplan Meier plotter database. The results indicated that OS, PFS, RFS and DSS were shorter in patients with low PVALB expression compared to those with high PVALB expression (Figure 2A–D, P<0.05). Subsequently, we ran a rate-of-occupancy curve (ROC) analysis to define the diagnostic worth of PVALB expression. It was found that the PVALB expression level could be precisely differentiated between HCC tissue and normal tissue with an AUC as high as 0.964 (Figure 2E). We examined PVALB expression in different cell lines and finally selected HCCLM3 cell line with the lowest expression for subsequent experiments (Supplementary Figure 7A). We then generated Flag-PVALB cell lines and examined them (Supplementary Figure 7B). Findings from the colony formation assay indicated that up-regulation of PVALB expression diminished the colony-forming ability of HCC cells (Figure 2F). CCK8 assay also demonstrated that overexpression of PVALB inhibited the proliferation of HCC cells (Supplementary Figure 7C). In summary, we discovered that low expression of PVALB portends a worse outcome for HCC patients and PVALB can be a valid separate predictor of survival.

Figure 2 Relationship between PVALB differential expression and prognosis of HCC patients. Kaplan-Meier plotter was employed to draw (A) OS (overall survival), (B) PFS (progression-free survival), (C) RFS (relapse-free survival), and (D) DSS (disease-specific survival) survival curves. (E) Diagnostic ROC curves was used to distinguish the diagnostic value of PVALB.(F) Colonies formed by HCC cells transfected with shRNA targeting PVALB. Statistical chart shows the quantitative analysis of the result of the colony formation assay (**p < 0.01).

DNA Methylation of PVALB in HCC

It is known that the changes in the DNA methylation level of a gene in HCC could affect its own expression level,39,40 so we started with the effect of DNA methylation on the PVALB expression, investigating the methylation level of PVALB. As a first step, we drew a heatmap using the MethSurv database and found that six sites were hypomethylated when PVALB was lowly expressed (Figure 3A). Then, through the SMART database, we explored the sites where the PVALB gene undergoes methylation and visualized the distribution of these sites on the chromosome (Figure 3B). The detailed distribution patterns of the sites were additionally displayed, with one of them on the N_shelf (cg16293983), five of them on the island (cg25625146, cg23520948, cg21789136, cg14336003, and cg03263929), six on the S_shelf (cg21385983, cg19136704, cg26244661, cg02978737, cg07856761, cg04688476) and cg06066601 was on the N_shore (Figure 3C). Also, the relationship between promoter methylation level of PVALB and various clinical features was assessed. Through the UALCAN website, we discovered that there was less PVALB promoter methylation in HCC tissue in comparison with normal tissue (Supplementary Figure 4A). Besides, the degree of PVALB promoter methylation was related to patient age, gender, presence of TP53 mutation and cancer stage according to further data. Specifically speaking, the methylation level of the PVALB promoter was strongly inversely linked with age (Supplementary Figure 4B), and it was lower in males than in females (Supplementary Figure 4C). Among them, it was captivating that TP53 mutant patients possessed a lower degree of methylation than TP53 non-mutant patients (Supplementary Figure 4D) and HCC patients at varying stages had distinct PVALB methylation (Supplementary Figure 4E). Apart from that, we once more examined the methylation of the sites and its effect on survival applying the SMART online tool, and MethSurv database. Figures 3A and E of the results indicated that cg25625146 was less methylated in HCC tissue than in normal tissue, with demethylation of this locus associated with adverse outcome in HCC patients (Figure 3D–F). All in all, hypomethylation at the cg25625146 locus might affect PVALB expression and induced poor survival in HCC patients.

Figure 3 DNA methylation of PVALB in HCC. (A) Heatmap indicated the extent to which CpG sites were methylated. (B) Overall distribution for PVALB methylation sites on chromosome. (C) Localization of CpG sites linked to PVALB. (D) HCC patients had a worse prognosis when the cg25625146 locus was hypermethylated. (E) comparison of methylation level at PVALB cg25625146 site in normal and HCC tissue. (F) SMART database verified that HCC patients had a worse prognosis when hypomethylated at the cg25625146 locus.

Enrichment Analysis of PVALB Gene Co-Expression Network in HCC

Striving to have a better understanding of the biological significance of PVALB, we attempted to pinpoint co-expressed genes with PVALB as well as analyzed the functional pathways in which they were involved. The LinkedOmics website was firstly used for the study, with volcano plot exhibiting a large number of genes significantly positively (red dots) and negatively (green dots) associated with PVALB (Figure 4A). Meanwhile two lollipop charts revealed some co-expressed genes showing the strongest plus and minus correlation with PVALB (Figure 4B and C). Following that, we performed GO (biological process) pathway profiling on these commonly expressed genes via the LinkedOmics website again, which suggested that they were majorly enriched in immune-related pathways, like type 2 immune response and adaptive immune response (Figure 4D), and the further KEGG analysis also revealed that PVALB and its related genes were engaged in immune-related pathways, such as primary immunodeficiency (Figure 4E). Within the next step, we fetched the intersection of 4519 prognostic factors from HCC with 100 genes most related to PVALB as illustrated in Figure 4B and C, and ultimately obtained 8 genes (Figure 4F). As displayed by Figure 5G, we further verified the relationship between the 8 genes at the intersection and PVALB expression, finding that when PVALB was highly expressed, the expression of genes TYROBP and CARD9 upregulated, and the remaining 6 genes were all down-regulated. Finally, we performed GO and KEGG analyses on these 8 HCC-related genes that were most correlated with PVALB as well as capable of influencing HCC patients survival, and the network diagram displayed that these PVALB co-expressed genes were functionally most relevant to the positive regulation of tumor necrosis factor superfamily cytokine production and myeloid leukocyte mediated immunity (Figure 4H). Taken together, our study revealed that genes co-expressed in PVALB were majorly enriched in pathways related to immunity.

Figure 4 Enrichment analysis of PVALB gene co-expression network in HCC. (A) Volcano plot showed co-expressed genes associated with PVALB expression. (B–C) Lollipop charts showed the top 50 co-expressed genes each with the strongest positive and strongest negative correlation with PVALB expression. (D) Enrichment analysis of Gene Ontology (GO) terms for PVALB co-expressed genes (GO-BP). (E) Enrichment analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) terms for PVALB co-expressed genes. (F) Eight genes were identified at the intersection of 4519 HCC prognosis factors and 100 genes with the strongest positive and negative correlations with PVALB. (G) Co-expression heatmap displayed the relationship between PVALB and 8 genes. (H) Significant enrichment of the GO and KEGG pathways of the 8 genes.

Figure 5 Relationship between PVALB and immune infiltration in HCC. (A–C) Correlation of PVALB expression level with StromalScore, ImmuneScore and ESTIMATEScore was demonstrated by scatter plots based on TCGA-LIHC data. Correlation of PVALB expression with multiple immune cell infiltration based on (D) TCGA-LIHC data and (E) TIMER database. (F) Grouped comparison plot demonstrated the enrichment scores of different immune cells in high and low PVALB expression groups. (***p < 0.001, **p < 0.01, *p < 0.05, ns: non-significant).

The Link Between PVALB and the Immune Microenvironment of HCC

According to previous findings, the tumor-infiltrating lymphocytes and the cytokines it generates have been shown to remodel the tumor immune microenvironment, thus exerting an effect on the existence of tumor sufferers.41 Besides, based on the results of enrichment analysis, it was observed that PVALB and its co-expressed genes were mostly engaged in pathways related to immune (Figure 4), and this encouraged us to further probe the connection between PVALB and the immune microenvironment in HCC. First, on the basis of TCGA-LIHC data, we computed the association of PVALB expression with StromalScore, ImmuneScore and ESTIMATEScore in HCC by applying the ESTIMATE algorithm, and it pointed a strong positive relationship of PVALB expression with three types of scores, which were 0.35, 0.40 and 0.41, respectively (Figures 5A–C). Then, we independently examined the correlation of PVALB expression with the different immune cells infiltration using the data from TCGA and TIMER databases. In HCC, both findings indicated that PVALB was most highly correlated with macrophages infiltration, while PVALB expression also exhibited a significantly positive correlation with the infiltration of various immune cells such as NK cells, B cells, CD8+ T cells, neutrophils, dendritic cells, TFH etc. (Figure 5D and E). Additionally, the infiltration of specific immune cell in each tumor on the basis of various algorithms was demonstrated in Supplementary Figure 5. At the same time, classifying PVALB into high and low expression groups, an R-based ssGSEA analysis on HCC revealed significantly higher immune cell infiltration in the PVALB high expression group across a wide range of immune cells (Figure 5F).

Next, to validate the above outcomes, we assessed the relationship between tumor-infiltrating lymphocytes (TILs) and the expression of PVALB by utilizing the TISIDB database, and the finding displayed that almost all lymphocytes were positively correlated with PVALB expression in HCC, of which macrophages and MDSC cells had the highest relevance at 0.432 and 0.441, respectively (Figure 6A). Alternatively, after further characterizing the association of PVALB expression with major histocompatibility complex (MHC), we similarly discovered that PVALB maintained a positive association with most of the MHCs, especially HLA-DPB1 (rho=0.399, p=1.68e-16) and HLA-DQB1 (rho=0.396, p=1.02e-15) (Figure 6B). Lastly, we comprehensively explored the association of PVALB expression with all cytokines, particularly chemokines and chemokine receptors capable of recruiting circulating NK cells to the site of tumorigenesis in HCC, and the outcomes showed a considerable correlation between PVALB expression and the chemokines CXCL12 and CX3CL1, as well as the chemokine receptors CXCR3 (rho=0.252, p=8.72e-07), CXCR4 (rho=0.26, p=3.86e-07), CCR7 (rho=0.198, p=0.000126), CX3CR1 (rho=0.136, p=0.00858) (Figure 6C and D). Apart from that, it was noted that PVALB expression was highest in iCluster1 and C2 (IFN-gamma dominant) subtypes when resolving the expression of PVALB in diverse molecular and immune subtypes (Figures 6E and F). These findings suggested that PVALB was closely associated with the immune microenvironment in HCC and may play a role in regulating tumor immunity.

Figure 6 Relationship of PVALB with lymphocytes, MHCs and cytokines as well as their approximate expression level. Correlation of PVALB expression level with (A) lymphocytes (B) MHCs (C) chemokines (D) chemokine receptors in HCC. Expression of PVALB in different (E) molecular subtypes and (F) immunosubtypes in HCC.(**p < 0.01, *p < 0.05).

PVALB Expression May Impact the Survival of HCC Sufferers Through Immune Cell Infiltration

Given the confirmed association between PVALB expression and the prognosis of HCC patients, as well as immune cell infiltration, it is plausible to suggest that PVALB may potentially impact HCC prognosis through specific immune cells. To test this conjecture, we examined the impact of PVALB on HCC patients survival utilizing the Kaplan-Meier plotter database, comparing outcomes for patients with enriched or decreased immune-infiltrating lymphocytes in the tumor immune microenvironment. Strikingly, when Natural killer T-cells and Type 1 T-helper cells were decreased, the prognosis of HCC patients with low PVALB expression was worse, whereas when Natural killer T-cells and Type 1 T-helper cells were enriched, the prognosis of HCC patients were not affected, regardless of the high or low expression of PVALB (Figure 7D and E). Meanwhile, the survival situation of the HCC patients would not be influenced when the infiltration of other immune cell change, no matter how PVALB expressed (Figure 7A–7C, 7F–7I). Therefore, this suggests us that PVALB may influence the survival of HCC sufferers by Natural killer T-cells and Type 1 T-helper cells.

Figure 7 Relationship between PVALB expression and prognosis of HCC patients based on different immune cell subgroups. (A) B cell; (B) CD4+ memory; (C) CD8+ T-cell; (D) Natural killer T-cell; (E) Type 1 T-helper cell; (F) Type 2 T-helper cell; (G) Regulatory T cell; (H) Mesenchymal stem cell; (I) Macrophage.

PVALB May Affect NK Cell Infiltration Through Specific Pathways in HCC

Considering the effect of NK cell enrichment on the survival of HCC patients, we would like to investigate the relation between PVALB and the Fas pathway from a bioinformatic perspective in depth. We initially ran a GSEA study, and the mountain range plot displayed that PVALB was enriched in five related pathways including the Fas pathway (Supplementary Figure 6A), with more detailed enrichment analyses as shown in Supplementary Figure 6B–F. These plots suggested PVALB was enrichment in the Tnfrelated Weak Inducer of Apoptosis Tweak Signaling Pathway, Pi3k akt Signaling Pathway, Fas Pathway, WNT Pathway and Adaptive Immune System, which may exert vital biological functions. Following this, we created scatter plots to investigate the correlation between PVALB and NK cell surface markers as well as ligands for Fas, discovering that when the expression of the PVALB was down-regulated, FCGR3A (CD16), NCAM1 (CD56) and FASLG (FasL) were all decreased (Supplementary Figure 6G–6I). Combined with the analysis of our previous results, we hypothesized that low expression of PVALB would affect the expansion content of some key proteins in the Fas pathway by influencing NK cell infiltration, and then activate or inhibit this pathway in order to reduce apoptosis, which ultimately promotes the immune escape of HCC cells to boost their proliferation.

The Lnc-YY1AP1-3/Hsa-miR-6735-5p/PVALB Axis May Modulate the Progression of HCC

Recently, there has been a growing amount of researches demonstrating that the lncRNA-miRNA-mRNA regulatory network had a pivotal role in driving the progression and pathogenesis of HCC.42,43 Consequently, we constructed a ceRNA regulatory network focused on PVALB in HCC and attempted to elucidate the underlying mechanism of PVALB from another perspective. Via TargetScan and DIANA databases, we jointly predicted 16 miRNAs potentially binding to PVALB (Figure 8A), while only one of them (hsa-miR-6735-5p) was identified to satisfy a negative correlation with PVALB (Figure 8B). Utilizing the TargetScan database once again, we also pinpointed promising sites on PVALB mRNA that bind to hsa-miR-6735-5p and subsequently introduced mutation at these sites (Figure 8C). We performed RIP-qPCR and demonstrated that PVALB could bind to hsa-miR-6735-5p (Figure 8D). Based on the discovery of a dual-luciferase reporter gene system, we ascertained that the hsa-miR-6735-5p mimic could reduce the luciferase activity of cells that transfected with wild-type PVALB (Figure 8E). Additionally, as seen in the comparison with the mock group, the hsa-miR-6735-5p mimic was capable of inhibiting the expression of PVALB in HCCLM3 cells, whereas the hsa-miR-6735-5p inhibitor has conversely increased the expression level of PVALB (Figures 8F–G). Afterwards, exploiting the LncRNASNP2 database, we delved deeper into the potential lncRNAs that might bind to hsa-miR-6735-5p, and identified the only one of these lncRNAs (lnc-YY1AP1-3) that conformed to negative correlation with hsa-miR-6735-5p and positive correlation with PVALB according to the experimental results. Similarly, the prospective binding sites of hsa-miR-6735-5p to lnc-YY1AP1-3 were displayed in Figure 8H, and the binding sites of hsa-miR-6735-5p were further mutated. Alternatively, the result of the dual-luciferase reporter gene system revealed that the lnc-YY1AP1-3 mimic reduced luciferase activity of cells that transfected with the hsa-miR-6735-5p wild type, whereas this did not occur of cells that transfected with the mutant-type hsa-miR-6735-5p (Figure 8I). Besides, it was not difficult to notice from Figure 8J that lnc-YY1AP1-3 mimic can suppress the expression of hsa-miR-6735-5p, and the overexpression of lnc-YY1AP1-3 can also elevate the expression level of PVALB in HCCLM3 cells (Figure 8K). Thereafter, relying on a rescue experiment, we also observed that lnc-YY1AP1-3 overexpression was able to rescue the reduction in PVALB expression level caused by hsa-miR-6735-5p overexpression (Figure 8L). Ultimately, as indicated in the colony formation assay, hsa-miR-6735-5p could facilitate the proliferation of HCCLM3 cells, in direct contrast to the effect of lnc-YY1AP1-3. Notably, hsa-miR-6735-5p overexpression had the ability to rescue the attenuation of the cellular proliferative capacity caused by lnc-YY1AP1-3 overexpression (Figure 8M). Altogether, our findings suggest that the lnc-YY1AP1-3/hsa-miR-6735-5p/PVALB axis may be involvement in modulating PVALB expression, with potential implication for influencing the proliferation of HCC cells.

Figure 8 Construction of the ceRNA regulatory network of PVALB in HCC. (A) Venn diagram showed miRNAs jointly predicted by TargetScan and DIANA databases. (B) Correlation of PVALB with hsa-miR-6735-5p. (C) Potential binding sites of PVALB to hsa-miR-6735-5p predicted by TargetScan database. (D) RIP-qPCR showed the enrichment of hsa-miR-6735-5p detected by the PVALB. (E) Potential binding sites of PVALB to hsa-miR-6735-5p validated by the dual-luciferase reporter gene system. (F) mRNA expression level of PVALB under the influence of hsa-miR-6735-5p mimic and (G) hsa-miR-6735-5p inhibitor. (H) Potential binding sites of lnc-YY1AP1-3 to hsa-miR-6735-5p predicted by the LncRNASNP2 database and (I) validated by the dual-luciferase reporter gene system. (J) Expression level of hsa-miR-6735-5p under the influence of lnc-YY1AP1-3 mimic. (K) Protein expression level of PVALB in mock and lnc-YY1AP1-3 groups were explored by Western blot. (L) Changes in protein expression level of PVALB under four different conditions. (M) Changes in the number of colonies in four different groups.(***p < 0.001, **p < 0.01, *p < 0.05, ns: non-significant).

PVALB Expression Could Serve as a Prognostic Marker for Immunotherapy Response in HCC Patients

Over the past few years, immunotherapy has gradually come into the limelight and is rapidly growing as a robust clinical strategy for the healing of cancer.44 In this article, we carried out a lot of works with the aim of further elucidating the role of PVALB in immunotherapy (mainly ICB). Initially, employing TCGA-LIHC data, an association between PVALB expression and immune checkpoint genes was investigated, and we observed significant correlation between PVALB expression and a majority of these factors (Figure 9A). Specifically, it was meaningfully and positively linked to 9 immune checkpoint genes except TGFBR1 (Figure 9B). In addition, the 3 most correlated genes were presented in scatter plots, which were CD244 (R=0.403, p<0.001), TIGIT (R=0.400, p<0.001), CD96 (R=0.375, p<0.001) (Figure 9C–E). As a critical tumor feature, MSI may be connected with the reaction to immunotherapy,45 therefore, we probed the link between PVALB expression and MSI score in HCC. It was revealed that the two maintained a negative correlation, suggesting that the probability of mutation was higher when PVALB was lowly expressed in HCC, at which time a better result was achieved with immunotherapy (Figure 9F). Simultaneously, the TIDE score also implied that lower score may be received when PVALB was lowly expressed, which hinted us that immune checkpoint blockade therapy was more therapeutic when the expression of PVALB was low in HCC (Figure 9G). All of these findings indicated that PVALB may be instrumental in forecasting the response to immunotherapy for HCC patients, and better efficacy may be harvested with immunotherapy in the presence of low PVALB expression.

Figure 9 Role of PVALB in immunotherapeutic response. (A) Chordal plot demonstrated the overall connection between PVALB and immune checkpoint genes (B) Heatmap displayed the specific relationship between PVALB expression and 10 immune checkpoint genes (C–E) Scatter plots showed the correlation between PVALB expression and CD244, TIGIT, CD96 (F) Association between PVALB expression and MSI score (G) Box plot demonstrates the difference in TIDE score between high and low PVALB expression groups.(****p < 0.0001, *p < 0.05).

Network and Drug Sensitivity Analysis Across PVALB and Its Interacting Genes

In search of other related proteins that possibly interact with PVALB, we constructed the interaction network with GeneMANIA online tool for PVALB, in which PVALB had an explicit physical interaction with SREBF1 (Figure 10A). Then, we have further downloaded the secondary structure of some proteins, which were encoded by these two genes, and analyzed their structural domains. The analysis gave the results that PVALB had an EF-hand_7 structural domain and SREBF1 possessed an HLH structural domain (Figure 10B). Subsequently, we obtained the spatial structure of SREBF1 from the SWISS-MODEL database and made a prediction of the spatial structure of PVALB with the AlphaFold database, then utilizing HDOCK Server to predict the docking model of its protein-protein interaction (Figure 10C). Besides, by modeling the PVALB-SREBF1 molecular docking, we also observed that it contains 11 hydrogen bonds (Figure 10D). Moreover, in order to better probe the clinical value of PVALB for HCC patients, we thus examined the relevance of PVALB expression to multiple drugs related to cancer via the CTD, finding that 8 drugs could affect PVALB expression, of which copper could boost its expression level (Figure 10E). Afterwards, we mapped using the TISIDB and recognized PVALB as one of the targets of the drug DB01942, also known as formic acid (Figure 10F). On a final note, the GSCA database was employed to perform drug sensitivity analysis, as shown in Figure 10G. To our surprise, with low PVALB expression, cells were only resistant to Navitoclax, while sensitive to the other 16 drugs. Whereas, in the CTRP drug sensitivity analysis, cells with low PVALB expression were sensitive to CAL-101, TGX-221, dasatinib and lovastatin (Figure 10H). In conclusion, our findings provide new and alternative treatment options for HCC patients accompanied by low PVALB expression.

Figure 10 Molecular docking analysis of PVALB and its drug sensitivity analysis. (A) Construction of the interaction network of PVALB with other genes by GeneMANIA database. (B) Protein secondary structures underlying PVALB and SREBF1. (C–D) Molecular docking model of PVALB-SREBF1 (E) Cancer-related drugs from CTD that interact with PVALB. (F) Drugs targeting PVALB studied by TISIDB. Pharmacosensitivity analysis of PVALB and SREBF1 based on the (G) GDSC and (H) CTRP databases (The large column represent PVALB and the small columns represent SREBF1).

Discussion

HCC is the fifth largest cancer in the world,46 with high malignancy, poor long-term prognosis and easy recurrence. The five-year survival rate of HCC is only 18%. At present, the treatment of HCC is still mainly based on surgical treatment, supplemented by radiotherapy, chemotherapy and ablation.47–49 In addition, there are other novel approaches, such as TACE, minimally invasive interventional therapy, molecular targeted therapy, immune checkpoint inhibition, etc,50–53 which are also the research hotspots of HCC treatment. Some existing molecular biomarkers, such as AFP, are used for the auxiliary diagnosis of HCC in clinical practice, but their sensitivity and accuracy are still insufficient. Finding a more reliable biomarker is therefore essential for the diagnosis, therapy and prognostic assessment of HCC patients. In this research, an array of rigorous bioinformatic analyses and experimental validation identified PVALB as an independent prognostic factor and a prospective therapeutic candidate that can regulate the HCC immune microenvironment and affect HCC progression.

Immunological escape is a key factor in the refractoriness of HCC, which possesses a complicated immunological microenvironment.54,55 The liver has long been recognized as the site of immune tolerance. This property results from the interaction between para-hepatic cells and peripheral leukocytes, involving high expression of immune checkpoints, IL-10 and TGFβ mediated immune suppression microenvironment, etc. NK cells in the liver can directly kill stress cells and tumor cells, and mediate antibody-dependent cytotoxicity (ADCC) with the participation of CD16.56 Once carcinogenesis occurs, these regulatory mechanisms are destroyed, the immune microenvironment of the liver will also be disordered, the infiltration of immune cells will be changed, the infiltration of dendritic cells, tumor-associated macrophages, regulatory CD4+T cells and other cells will be increased,57–59 CD8+T cells will be functionally depleted.60,61 NK cells will also show a depleted phenotype, and express high level of programmed death ligand 1 (PD-L1).62,63 The secretion of immunosuppressive cytokines such as IL-10 causes immune escape of tumor cells, and HCC develops in the direction of malignancy. Therefore, it is important to explore the immune microenvironment of HCC and the mechanism of immunosuppression. In our study, MCP-counter findings also suggested the importance of exploring the immune microenvironment in HCC.

In order to find the key factors that can prevent the development of immunosuppressive microenvironment and immune escape of HCC, we used WGCNA, based on TCGA and GEO databases, to screen out the tumor suppressor modules closely related to NK cells and adaptive immunity. Through LASSO regression analysis, combined with the background report of the gene, we finally selected PVALB. Our results firstly indicated that HCC patients with elevated PVALB expression were inclined to exhibit better prognostic outcomes, and low expression of PVALB was a biomarker indicating a unfavorable outcomes in HCC patients. We also found that PVALB was hypomethylated, and cg25625146 was hyp