Identification of VRK1 as a Novel Potential Biomarker for Prognosis and Immunotherapy in Hepatocellular Carcinoma

Introduction

A member of the vaccine-related kinase (VRK) family of serine/threonine protein kinases, VRK1 exhibits unusual activity similar to that of nucleosome kinase, chromatin kinase Mars, or Ser-Thr kinase.¹ Through stable and direct interactions with various chromatin proteins implicated in carcinogenesis, VRK1 can take part in biological processes on the cell cycle, apoptosis, and autophagy regulation.² Research has demonstrated that VRK1 was significantly expressed in a large number of human malignancies and could influence patient prognosis. High VRK1 expression was closely associated with poor prognosis of bladder cancer patients, while down-regulation of VRK1 significantly inhibited bladder cancer cell proliferation.³ The lung squamous cell carcinoma cells’ motility, apoptosis, proliferation, and cell cycle progression were all greatly impacted by the downregulation of VRK1, which implies that VRK1 is a good therapy target.⁴ By initiating the p53 signaling pathway, luteolin could decrease the expression of VRK1, causing cell proliferation to be inhibited and apoptosis and cell cycle arrest at G2/M in ovarian cancer.⁵ In HCC, by regulating the amounts of cell cycle regulators linked to the G1/S transition, VRK1 encouraged tumor development.⁶ Unfortunately, no research has documented a connection between VRK1, immune cell infiltration, and the response to immunotherapy in HCC, which deserves further investigation.

In this study, we explored the expression profiles and potential prognostic value of VRK1 in multiple publicly available datasets and further validated them in clinical samples using qRT-PCR and IHC to systematically investigate the relationship between VRK1 and pathological features and clinical prognosis of HCC patients. In addition, the study evaluated the predictive value of VRK1 expression on immunotherapy response in HCC patients. Finally, the study also explored the correlation between VRK1 expression and immune infiltration and m6A-related genes. These findings provide important insights into the potential significance of VRK1 in HCC.

Methods

Samples Collection

We downloaded the normalized expression data as well as relevant clinical features of HCC and normal samples from The Cancer Genome Atlas-LIHC cohort (TCGA, https://xenabrowser.net/), International Cancer Genome Consortium-LIRI-JP cohort (ICGC, https://dcc.icgc.org/) and Gene Expression Omnibus-GSE14520 (GEO, https://www.ncbi.nlm.nih.gov/geo) databases, respectively. The log2 (FPKM+1) transformation was used to standardize the transcriptome data. The clinical features of these patients are shown in Supplementary Table S1.

2.2 qRT-PCR and Immunohistochemistry (IHC) Assay

The VRK1 mRNA expression was determined in tumor tissues and paracancerous tissues from 50 HCC patients using the qRT-PCR technique, as was previously done.⁷ The primer synthesized sequences of VRK1 were as follows: forward: 5ʹ-CCTCGTGTAAAAGCAGCTCAA-3ʹ, reverse: 5ʹ- GCCAATGGGTAATCCTACTTTCC-3ʹ. Paraffin-embedded tissues from HCC patients were subjected to independent IHC assays using anti-VRK1 (Proteintech, China). After secondary antibodies were added to the slides, a Leica DM 2500 microscope was used to capture images. The necessary proteins’ immunostaining intensity was assessed and ranked independently by two different observers. A final score was obtained by adding the extent of expression score (no positive cells = 0, <10% = 1, 10–50% = 2, positive staining of >50% = 3) and intensity score (negative = 0, weak = 1, moderate = 2, strong = 3) was used to distinguish between negative (≤4) and positive (>4) expression of VRK1 protein in HCC and paracancerous tissues.

Immune Cell Infiltration Analysis

The infiltration of different kinds of immune cells in tumor microenvironment (TME) employing the CIBERSORT,⁸ xCELL,⁹ MCPcounter,¹⁰ TIMER databases,¹¹ and single-sample genome enrichment analysis (ssGSEA) algorithm was done to better investigate the variability in HCC.

Immunotherapy Analysis

We examined the variations in the expression of numerous immunomodulatory genes between the VRK1-positive and VRK1-negative subgroups. As there are fewer clinical cohorts of immunotherapy for HCC patients, we further determined the predictive role of VRK1 for immunotherapy using an anti-PD1/PD-L1 inhibitor cohort (IMVigor 210) with considerable richness of transcriptomic data and information on response to immunotherapy. Based on their response to immunotherapy, individuals in the IMvigor 210 cohort were divided into four groups: stable disease (SD), partial response (PR), disease progression (PD), and complete response (CR). We examined the disparities between patients with low and high expression of VRK1 in terms of survival outcomes and responses to anti-PD-L1 immunotherapy. Finally, we also compared the difference in neoantigen burden between the two groups.

Potential Tumor-Sensitive Drug Prediction

In particular, we examined the relationships between VRK1 expression and medications included in the CellMiner database.¹² If a medication’s adjusted P-value is less than one in a thousand and its Pearson correlation coefficient is greater than four-tenths, the medication is deemed tumor-sensitive. To show therapeutic sensitivity, the half-maximal inhibitor dose (IC50) of the targeted drug was then estimated using gene expression levels.

The Pathway Richness Analysis

To assess the variation in biological processes between the VRK1-negative and VRK1-positive subgroups in the TGCA cohort, we used Gene Set Variation Analysis (GSVA) and Gene Set Enrichment Analysis (GSEA) with the R package “GSVA” and “ClusterProfiler”, respectively.

Statistical Analysis

All statistical analyses were carried out by R-4.2.1. In the context of quantitative data, variables with non-normal distributions were analyzed for statistical significance using the Wilcoxon rank sum test, while variables with regular distributions were checked for statistical significance using the Student’s t-test. The Benjamini-Hochberg approach was used to adjust the false discovery rate (FDR) for multiple hypothesis testing in all two-sided comparisons with an alpha level of 0.05.

Results

VRK1 Was Highly Elevated in HCC and Linked to Poor Progression

We performed gene expression analysis using RNA sequencing data obtained from the three publicly available datasets and showed that VRK1 was significantly higher in HCC than in normal tissues (Figure 1A). We also evaluated the ability of VRK1 to distinguish cancer from normal samples using ROC curves and found that VRK1 was able to distinguish HCC tissues from the neighboring normal tissues well, with AUC values of 0.878, 0.865, and 0.878, respectively (Figure 1B). In addition, a significant difference in the protein expression of VRK1 in normal tissues and HCC tissues was observed in the Human Protein Atlas database (https://www.proteinatlas.org/) (Figure 1C). The results of the qRT-PCR (Figure 1D) and IHC (Figure 1E) assay in clinical samples were further validated by their consistency with the bioinformatics analysis. Furthermore, the expression of VRK1 was found to be mainly localized in the nucleus by immunofluorescence staining analysis (Figure 1F). Finally, we explored the relationship between VRK1 expression and clinicopathological characteristics of HCC patients and found that VRK1 expression was closely associated with tumor recurrence (Figure 2A), Grade stage (Figure 2B), TNM stage (Figure 2C), T stage (Figure 2D), and Alpha-fetoprotein (AFP) value (Figure 2E), while it was not associated with gender (Figure 2F), cirrhosis (Figure 2G), and microvascular invasion status (Figure 2H). According to the results above, VRK1 was overexpressed in HCC and has a favorable correlation with tumor progression.

Figure 1 VRK1 was highly elevated in HCC. (A) VRK1 was significantly higher in HCC than in normal tissues in the three publicly available datasets. (B) VRK1 was able to distinguish HCC tissues from the neighboring normal tissues well. (C) A significant difference in the protein expression of VRK1 in normal tissues and HCC tissues was observed in the Human Protein Atlas database. (D–E) The results of the qRT-PCR and IHC assay in clinical samples. (F) Immunofluorescence staining analysis. ***p < 0.001.

Figure 2 VRK1 was linked to poor progression in HCC. (A) VRK1 expression was closely associated with tumor recurrence, (B) Grade stage, (C) TNM stage, (D) T stage, and (E) AFP value, while it was not associated with (F) gender, (G) cirrhosis, and (H) microvascular invasion status. NS, not significant; **p <0.01; ***p < 0.001.

VRK1 Was a Novel Potential Prognosis Biomarker for HCC

Kaplan-Meier survival analysis demonstrated that high VRK1 expression was strongly associated with poor survival outcomes in HCC patients, and the higher the VRK1 expression, the shorter the survival time of the patients (Figure 3A). Moreover, VRK1 expression was also a good predictor of 1-, 2- and 3-year survival in HCC patients (Figure 3B). Finally, due to the univariate Cox regression analysis (Figure 3C) and multivariate Cox regression analysis (Figure 3D), increased VRK1 was an independent effector of poor prognosis in HCC patients. According to the results above, increased VRK1 was a notable effector in HCC patients who had a poor survival outcome.

Figure 3 VRK1 was a novel potential prognosis biomarker for HCC. (A) Kaplan-Meier survival analysis. (B) VRK1 expression was also a good predictor of 1-, 2- and 3-year survival in HCC patients. (C–D) Univariate Cox regression analysis and multivariate Cox regression analysis.

Relationship Between VRK1 Expression and TME

The results of TIMER (Figure 4A), MCPcounter (Figure 4B), CIBERSORT (Figure 4C), Xcell (Figure 4D), and ssGSEA (Figure 4E) analyses showed that the infiltration levels of multiple immune cells were also higher in the TME of samples with high VRK1 expression. Furthermore, samples with higher VRK1 tended to have higher major histocompatibility complexes (MHCs) expression than those with lower VRK1 (Figure 4F). Considering that CD8 T cells differed in all five methods, we evaluated their infiltration in HCC clinical samples using IHC and found that the level of CD8 T cell infiltration was significantly higher in VRK1-positive samples than in VRK1-negative samples (Figure 4G). According to the results above, compared to samples with lower VRK1, those with higher VRK1 often exhibited a stronger immune infiltration in TME.

Figure 4 Relationship between VRK1 expression and TME. The results of TIMER (A), MCPcounter (B), CIBERSORT (C), Xcell (D), and ssGSEA (E). (F) Samples with higher VRK1 tended to have higher MHCs expression than those with lower VRK1. (G) The level of CD8 T cell infiltration was significantly higher in VRK1-positive samples than in VRK1-negative samples. ns, not significant; *p < 0.05; **p <0.01; ***p < 0.001.

Relationship Between VRK1 Expression and Immunotherapy Sensitivity

To explore potential sensitivity to immunotherapy, we compared the expression of several immunomodulatory genes, including PD1-related genes, CTLA4-related genes, and other antagonists or agonists of T cell activation-related genes, in samples with high or low VRK1 expression. The results showed that the expression of these genes was upregulated in VRK1 high-expressing samples (Figure 5A–C), suggesting that HCC patients with high VRK1 expression may be more sensitive to immunotherapy. We subsequently validated this speculation in the IMvigor210 cohort. After receiving anti-PD1/PD-L1 therapy, patients with high VRK1 expression had a better prognosis (Figure 5D) and were more likely to respond favorably to anti-PD1/PD-L1 therapy (Figure 5E). In addition, we found that neoantigen burden was also higher in samples with high VRK1 expression (Figure 5F). Finally, we examined the expression levels of PD1 and PD-L1 in 50 clinical HCC samples and found that although there was no correlation between VRK1 and PD1, there was a strong positive correlation between VRK1 and PD-L1 (Figure 5G), and PD1 and PD-L1 expression was higher in VRK1-positive samples compared to VRK1-negative samples (Figure 5H). According to the results above, HCC patients with high VRK1 expression may have better survival outcomes with anti-PD-L1 immunotherapy.

Figure 5 Relationship between VRK1 expression and immunotherapy sensitivity. The expression of (A) PD1 and CTLA4-related genes, and other antagonists (B) or agonists of T cell activation-related genes (C) were upregulated in VRK1 high-expressing samples. (D–E) After receiving anti-PD1/PD-L1 therapy, patients with high VRK1 expression had a better prognosis and were more likely to respond favorably to anti-PD1/PD-L1 therapy. (F) Neoantigen burden was also higher in samples with high VRK1 expression. (G) There was a strong positive correlation between VRK1 and PD-L1. (H) PD1 and PD-L1 expression was higher in VRK1-positive samples compared to VRK1-negative samples. *p < 0.05; **p <0.01; ***p < 0.001.

Prediction of Tumor-Sensitive Drugs Targeting VRK1

Based on the results of CellMiner analysis, we identified six tumor-sensitive drugs (Figure 6A), in which Chelerythrine and Nelarabine were significantly positively correlated with the expression of VRK1, while INK-128, GDC-0349, LY-3023414, and PQR-620 were significantly negatively correlated with the expression of VRK1. Among these six drugs, only the IC50 values of LY-3023414 and PQR-620 were differentiated between the samples with high or low VRK1 expression (Figure 6B). All these suggest that LY-3023414 and PQR-620 may be more effective in antitumor therapy for patients with high VRK1 expression.

Figure 6 Prediction of tumor-sensitive drugs targeting VRK1. (A) Six tumor-sensitive drugs targeting VRK1. (B) Among these six drugs, only the IC50 values of LY-3023414 and PQR-620 were differentiated between the samples with high or low VRK1 expression. ns, not significant; *p < 0.05.

Relationship Between VRK1 Expression and m6A-Related Genes

Considering that m6A modification plays a crucial role in HCC progression, we analyzed the relationship between the expression of 20 m6A-related genes and VRK1 expression.¹³ It was found that the expression of these m6A-related genes was not only significantly and positively correlated with VRK1 expression (Figure 7A), but their expression was also significantly elevated in samples with high VRK1 expression (Figure 7B). All these suggest that VRK1 may be involved in the progression of HCC by mediating m6A modifications.

Figure 7 Relationship between VRK1 expression and m6A-related genes. (A) The expression of these m6A-related genes was not only significantly and positively correlated with VRK1 expression, (B) but their expression was also significantly elevated in samples with high VRK1 expression. *p < 0.05; ***p < 0.001.

Prognostic Nomogram Construction

Considering that the TNM stage was also an independent factor influencing the poor prognosis of HCC patients in all three publicly available datasets, we constructed a nomogram including the TNM stage and VRK1 expression to predict the overall survival of HCC patients separately (Figure 8A). Subsequently, we further validated the predictive performance of the nomogram in the three datasets and found that the nomogram predicted well (Figure 8B) and its prediction was better than either the TNM stage or VRK1 (Figure 8C). All these suggest that the nomogram could independently predict survival outcomes of HCC patients very well.

Figure 8 Prognostic nomogram construction. (A) We constructed a nomogram including the TNM stage and VRK1 expression to predict the overall survival of HCC patients separately. (B) The nomogram predicted well in the three datasets. (C) The nomogram was better than either the TNM stage or VRK1.

Function Enrichment Analyses

To explore the potential functions of VRK1, we performed enrichment analysis of DEGs between HCC samples with higher or lower VRK1 expression. A total of 512 up-regulated and 12 down-regulated DEGs were identified at logFC >1 and p < 0.05 (Figure 9A). GSEA analysis revealed that these DEGs were mainly associated with DNA replication, cell cycle control, and fatty acid metabolism (Figure 9B), which was consistent with the results of GSVA analysis (Figure 9C).

Figure 9 Function enrichment analyses. (A) A total of 512 up-regulated and 12 down-regulated DEGs were identified. (B) GSEA analysis. (C) GSVA analysis.

Discussion

HCC is a main pathological form of liver cancer with a high death rate, poor clinical prognosis, and a propensity for metastasis and recurrence. Immunosuppressive TME is caused by massive interaction between antitumor, stromal, tumor parenchymal, and immunosuppressive cells. This reduces the frequency and function of antitumor cells while increasing that of related pro-tumorigenic cells, which ultimately leads to the spread of malignant tumors.¹⁴ Immunotherapy is becoming more and more popular as a cutting-edge treatment option for HCC because of the shortcomings of traditional medicines and the absence of reliable early detection techniques.¹⁵ To identify additional important targets and particular biomarkers and to create more efficient early diagnostic and tailored treatment strategies for hepatocellular carcinoma, it is imperative to clarify and comprehend the signaling pathways and molecular processes of cellular crosstalk in TME. This study revealed, for the first time, that there were variations in VRK1 expression not only between tumor and paracancerous tissues described in publicly accessible databases but also in the clinical tissue samples we collected. HCC patients who expressed high levels of VRK1 molecules had a worse prognosis and exhibited a stronger immune infiltration in TME. VRK1 expression was an independent risk factor affecting the prognosis of HCC. Additionally, patients with high expression levels of VRK1 were better candidates for anti-PD-L1 immunotherapy. These findings consistently show how important it may be to use VRK1 to determine the prognosis of HCC.

Immune cells can infiltrate tumor tissues and change their metabolism and function, which can improve immunosuppression and immune evasion.¹⁶ Researchers have shown that the existence of many immune cells, including macrophages, T cells, and B cells, within tumor tissues is linked to a poor prognosis in HCC patients.^17–21 Our findings demonstrated that VRK1 expression may be associated with the microenvironment. Interestingly, CD8+ T cells often have anti-tumor effects, whereas patients with high VRK1 expression have a poorer prognosis, suggesting that CD8+ T cells infiltrated in our VRK1 high-expression samples are likely to be in a state of exhaustion, and these deserve further experimental validation in future work.

Growing interest has been shown in immunotherapy to stop tumor growth by targeting the TME. The TME expresses several targets linked to immunotherapy, including PDCD1, PD-L1, LAG3, and CTLA4, which makes it easier to treat HCC with immunotherapy. When treating HCC, a range of immune checkpoint drugs that target PD1 and PDL1 are more effective than sorafenib and have a promising future.^22,23 At the same time, immunotherapy resistance in HCC is a result of immune system complexity, and one of the biggest unmet challenges in immunotherapy research is the absence of specific biomarkers to predict the impact of immunotherapy in HCC. Therefore, there is an urgent need to explore new biomarkers for accurately predicting the efficacy of immunotherapy in HCC patients. The expression of several immunomodulatory genes was upregulated in VRK1-overexpressing samples, including PD1-related genes, CTLA4-related genes, and other antagonists or agonists of T-cell-activation-related genes, and patients with VRK1-overexpressing patients had a better prognosis in the IMvigor210 cohort, suggesting that patients with VRK1-overexpressing HCC may have a better survival outcome after receiving anti-PD-L1 immunotherapy. Moreover, since m6A-modified genes typically have an oncogenic function in cancer, m6A-related therapies, such as m6A modification modulation or inhibition, may provide promising cancer treatment approaches.²⁴ In this investigation, we discovered a favorable correlation between m6A modification genes and VRK1 expression. Moreover, there was a substantial increase in these m6A modifier genes in HCC samples that had high VRK1 expression. These findings imply that VRK1 may be involved in the prediction of m6A-related treatments for the management of HCC.

Ly3023414, also known as Samotolisib, is a PI3K/mTOR pathway inhibitor that has been used in a variety of clinical antitumor therapies, such as prostate cancer,^25,26 endometrial cancer,²⁷ pancreatic adenocarcinoma,²⁸ and skin squamous cell carcinoma cell.²⁹ Unfortunately, Samotolisib has not been reported for antitumor therapy in patients with HCC. Similarly, PQR620 is an mTOR kinase inhibitor that has shown promising anti-tumor activity in non-small cell lung cancer³⁰ and lymphomas.³¹ In this study, we found that VRK1 was strongly and positively correlated with Samotolisib and PQR620 in the CellMiner database, and the IC50 of Samotolisib and PQR620 was lower in samples with high VRK1 expression, suggesting that HCC patients with high VRK1 expression are more sensitive to Samotolisib or PQR620 treatment, which, of course, needs to be verified by recruiting a certain number of HCC patients for clinical drug trials in the future.

Conclusion

In conclusion, abnormally elevated VRK1 in HCC is associated with the immune microenvironment and may serve as a novel biomarker for predicting the prognosis of HCC patients as well as their response to immunotherapy.