Discovery of Potential Biomarkers Related to Aggressive Progression and Vascular Invasion of Hepatocellular Carcinoma

Background Aggressive progression together with vascular invasion limit the treatment of hepatocellular carcinoma (HCC). In this study, we aim to identify novel candidate targets related to aggressive stage progression and vascular invasion in HCC and explore the underlying molecular mechanism. Results We rst found that high expression of VCX3A was signicantly associated with aggressive progression and vascular invasion in HCC, and indicated poor prognosis. Subsequently, hsa-miR-664a-3p and TUBA3C was found to be the most potential upstream miRNA and downstream targets of VCX3A, respectively. Conclusions Our research demonstrates that hsa-miR-664a-3p-VCX3A-TUBA3C axis may be a novel pathway linked to stage progression and vascular invasion of HCC, which provides promising prognostic biomarkers and effective therapeutic approaches for HCC. Methods Firstly, RNA sequencing data were obtained from The Cancer Genome Atlas. Differentially expressed genes with consistency in three comparison sets (normal VS stage I/II; stage I/II VS stage III/IV; none invasion VS invasion) were identied, after which survival analysis was performed. Next, hub gene expression was validated using clinical samples. Then, upstream miRNA was found according to miRNA predication, clinical pathological bioinformatic analysis and expression analysis. Besides, downstream molecular was found in terms of co-expressed analysis, survival analysis and expression analysis. ZKSCAN1 gene and its related circular RNA (circZKSCAN1) both inhibit hepatocellular carcinoma cell


Introduction
Hepatocellular carcinoma (HCC) is the most common type of primary liver cancer and ranks the fourth leading cause of cancer-associated deaths worldwide [1,2]. Estimated by the World Health Organization, more than one million liver cancer patients will die in 2030 [1]. What's more, China is the biggest country that owns liver diseases including liver cancer [3]. Although the early-diagnosed HCC patients are increasing today and those patients can be treated by surgical resection, local ablation or liver transplantation, while about 16% hepatocellular carcinoma patients have distant metastases at the diagnosis and 85% accompany with lymphovascular invasion on nal pathology, the prognosis of HCC patients remains dismal [4][5][6]. It has been well known that vascular invasion is a signi cant factor linked to HCC aggressive progression and poor prognosis [7][8][9]. Besides, microvascular invasion is reported to be an independent and powerful predictor for poor overall survival and early recurrence of HCC patients after surgical treatment [10]. Thus, it is urgent need for us to identify potential targets to manage and prevent the vascular invasion or aggressive progression for HCC patients, thereby improving hepatocellular carcinoma patient's outcome.
Generally, RNA consists of coding RNA (mRNA) and non-coding RNA (ncRNA). Mounting evidences suggest that ncRNA can dysregulate of speci c candidate genes and contribute to in uence the prognosis of HCC through modulating vascular invasion and tumor progression [11,12]. Nowadays, cancer research focus on ncRNA, especially miRNA (22 nucleotides small ncRNA), make it become a hotspot with the cumulatively deepened understanding of crucial functions in HCC through regulating of gene expression in the post-transcriptional manner [13,14].
In this research, we aimed to excavate novel candidate targets involved in stage progression and vascular invasion and insight into molecular mechanism in HCC by performing a series of bioinformatic analyses and experimental validation as shown in Fig. 1. Firstly, we screened out differentially expressed genes (DEGs) corelated with stage progression and vascular invasion of HCC using TCGA database. Then, through validating their prognostic values in HCC using Kaplan-Meier plotter, performing survival analysis in vascular invasion-HCC patients and employing expression analysis, we obtained a novel candidate gene, VCX3A. Additionally, based on miRNA prediction, tumor stage and grade analysis in OncomiR and expression validation, the most potential upstream miRNA (hsa-miR-664a-3p) was identi ed.
Subsequently, co-expression genes of VCX3A were found. Next, by using overall survival analysis in Kaplan-Meier plotter, correlation analysis in R software and expression validation in clinical samples, the most reliable downstream regulatory target (TUBA3C) was explored. The hsa-miR-664a-3p-VCX3A-TUBA3C network may provide novel promising therapeutic strategies for HCC patients.

Results
Screening key differentially expressed genes related with poor progression and vascular invasion of HCC in TCGA HCC patients in TCGA only contain both gene expression data and enough clinical information like pathological stage and vascular invasion were chosen in our study. Firstly, we divided 256 HCC patients into early stage HCC groups (stage I/II) and 90 HCC patients into advanced stage HCC groups (stage III/IV). However, 6 stage I/II HCC cases, 1 stage III/IV HCC case that lack of gene expression values and 24 HCC cases that lack of pathologic stage clinical information were excluded. What's more, according to different vascular invasion clinical information, HCC patients were classi ed into none invasion-HCC groups (n=206) and vascular invasion-HCC groups (n=109). Vascular invasion-HCC groups containing 93 micro-vascular invasion-HCC patients and 16 macro-vascular invasion-HCC patients. 4 none invasion-HCC cases from a total of 210 none invasion-HCC, 1 micro-vascular invasion-HCC case from a total of 94 micro-vascular invasion-HCC, 1 macro-vascular invasion-HCC case from a total of 17 macro-invasion HCC that lose gene expression values and 56 HCC cases that lack of vascular invasion clinical information were also be discarded. Next, by using edgeR package in R software and set the |log 2 FC| > 1.5 and p values < 0.05 as restricted condition, we found 2279 DEGs were upregulated and 492 DEGs were downregulated in early stage HCC compared with normal liver tissues ( Figure 2A and Table S1). Besides, 430 genes were found upregulated and 103 genes were found downregulated in advanced stage HCC in comparison with early stage HCC ( Figure 2B and Table S2). In addition, 375 differentially expressed genes between none invasion-HCC and vascular invasion-HCC were discovered, including 285 upregulated and 90 downregulated genes ( Figure 2C and Table S3). In this research, we aimed to analyze the hub genes related with the regulatory mechanism of the progression and the metastasis of HCC. Hence, 29 upregulated DEGs with consistency in three comparison sets were obtained, while no downregulated DEGs were found. Consequences were shown in Figure 2D-F. Finally, 29 hub genes were selected as candidate genes for further analysis.
Further prognosis and expression analysis found a hub gene: VCX3A in HCC Subsequently, we using Kaplan-Meier plotter to investigate the prognostic values of 29 hub gens in HCC.
Discovering TUBA3C as the most potential downstream target of VCX3A In order to explore the mechanism, further analyses were performed to obtain downstream regulatory molecule of VCX3A. 94 correlated genes with correlation values ≥ 0.5 from cBioPortal database were chose as potential regulatory partners ( Table 2). After that, starBase database was used to verify the correlations of those VCX3A-correlated genes pairs (Table 3). Only genes with correlation values also ≥ 0.5 in starBase were selected for the next research. 83 correlated genes were found have correlation values in starBase, while 11 correlated genes could not nd correlation values. Among 83 correlated genes, we discovered 57 key correlated genes with correlation values more than 0.5 ( Figure 6A). Next, we performed survival analyses of those 57 key genes through Kaplan-Meier plotter, the consequences indicated that only 21 genes were closely related with poor prognosis in HCC in accordance with VCX3A.  Figure 6B). Thus, 21 genes were selected for the further research. In the next step, correlation analysis of 21 key co-expressed genes in HCC patients with vascular invasion were performed through R software. As presented in Figure 6C, TUBA3C was found to be the most closely co-expression gene that interacted with VCX3A with r=0.66. Expression validations in TCGA and clinical samples demonstrated that TUBA3C was upregulated in tumor than in normal tissues, as well as in vascular invasion-HCC compared with none invasion-HCC tissues ( Figure 6D-E). Altogether, based on all the analyses, VCX3A-TUBA3C was found to be the most reliable interaction network in HCC (Figure 7).

Discussion
Aggressive progression along with vascular invasion have been reported to be related with worse prognosis and limit therapies of HCC patients. Though growing studies have identi ed the correlation of hub genes with the progression or metastasis of HCC, respectively, while we rst combined with stage progression and vascular invasion to screen key genes in HCC. It is valuable and important for us to explore the molecular mechanism of HCC progression and vascular invasion.
In this study, we rst discovered a hub gene VCX3A, which was increased in HCC tissues and associated with aggressive progression and higher vascular invasiveness in HCC through comprehensive expression analysis and survival analysis in TCGA. VCX3A belongs to the VCX/Y gene family. Although, many researches have revealed that VCX3A was closely involved in X-linked ichthyosis [18,19]. However, recent evidences demonstrated that VCX3A may works in multiple human cancers. For example, Deng, H et al.
suggested that knockdown the expression of VCX3A signi cantly inhibited the proliferation ability of pediatric glioma cells carrying the K27M mutation. The results indicated that VCX3A might sever as an oncogenic gene in pediatric high-grade gliomas [20]. Besides, Taguchi, A et al. found that the protein expression level of VCX3A could be detected in approximately 35% lung squamous cell carcinomas tissues and 20% lung adenocarcinomas tissues, while no expression could be detected in normal lung tissues. Thus, VCX3A might act as a cancer/testis (CT) antigen and a potential immunotherapeutic target in lung cancer [21]. However, research about VCX3A involved in cancer, especially in HCC is still extremely limited. Thus, additional studies of VCX3A are necessary to investigate the function and molecular mechanism in HCC. Importantly, we further explore the upstream and downstream mechanisms of VCX3A in HCC and discovered an important upstream molecule (hsa-miR-664a-3p) and a signi cant compared with corresponding peri-cancerous tissues and overexpression it could increase apoptosis of HCC in vitro and in vivo [29]. Study in our research also indicated that miR-664a-3p was a suppressor of HCC and closely associated with stage progression and vascular invasion. For TUBA3C, previous articles have indicated that TUBA3C was an alpha tubulin gene, and associated with Patau Syndrome and Kawasaki disease [30,31]. However, the role of TUBA3C in cancer is still unclear, more studies should be launched in the future.
Our data indicated that hsa-miR-664a-3p-VCX3A-TUBA3C axis was a novel pathway associated with HCC aggressive prognosis and vascular invasion, may provide three independent useful prognostic biomarkers and new insights into the underlying molecular mechanism of hepatocellular carcinoma.

Conclusion
In conclusion, by performing comprehensive bioinformatic analysis and experiment validation, we for the rst time suggest that hsa-miR-664a-3p-VCX3A-TUBA3C axis may be a critical pathway associated with stage progression and vascular invasion of HCC. The components in the pathway may be utilized as promising prognostic biomarkers and therapeutic targets of HCC in the future. More experimental trials should be launched.

Survival analysis
Kaplan-Meier plotter (http://kmplot.com/analysis/index.php?p=background) was utilized to evaluate the prognostic in uences of those hub genes and microRNAs in hepatocellular carcinoma. The result was presented as a forest plot using GraphPad prism software (version 6). Additionally, all of the HCC patients with vascular invasion in TCGA were classi ed into high-and low-expression groups in terms of the median expression values of hub genes and miRNAs, then survival curves were predicted using GraphPad prism software (version 6). Only logrank P-value < 0.05 was considered as statistically signi cant.

Construction of regulatory network of VCX3A
Three predicting programs including TargetScan (http://www.targetscan.org/vert 71/), miRWalk (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk/) and miRDB (http://www.mirdb.org/) were utilized to predict potential miRNAs binding to VCX3A. miRNAs appeared twice or more than twice in predicting programs were selected for the further analysis. In the next step, we evaluated the prognostic values of potential miRNAs in HCC through Kaplan-Meier plotter. Then "Tumor Stage and Grade" analyses on OncomiR (http://www.oncomir.org/) were performed. In summary, expression analysis in clinical samples demonstrated that hsa-miR-664a-3p was higher expressed in none invasion-HCC tissues than in vascular invasion-HCC tissues, on the contrary, hsa-miR-31-5p was higher expressed in vascular invasion-HCC than in none invasion-HCC. Thus, hsa-miR-664a-3p was the most potential upstream miRNA binding to VCX3A.
In order to obtain downstream regulatory molecule, cBioPortal (http://www.cbioportal.org/) database was performed to gain the co-expressed genes of VCX3A and starBase (www.sysu.edu.cn/) database was used to verify the correlations between co-expressed genes and VCX3A. Finally, a total of 67 co-expressed genes with correlation value ≥ 0.5 in two databases were selected for the next survival analysis using Kaplan-Meier plotter. As a result, 21 co-expressed genes with poor prognosis in HCC were selected for the further analysis. Besides, correlation analyses of those hub co-expressed genes expression in HCC patients with vascular invasion were performed in R software and TUBA3C was found to be the most reliable downstream target of VCX3A.

Clinical samples and sample preparation
The present study was approved by the Ethics Committee of the First A liated Hospital of Zhejiang University and obtained the informed consent from each patient. A total of 10 liver tissues of HCC patients without vascular invasion and 10 HCC patients with vascular invasion which diagnosed during 2017 to 2019 were gained from the First A liated Hospital of Zhejiang University (Hangzhou, China).
Total RNA was harvested using RNAiso plus Reagent (TaKaRa, Kusatsu, Japan). Besides, PrimeScript RT Reagent Kit (TaKaRa, RR0037A) and SYBR Premix Ex Taq (TaKaRa, RR420A) were used to perform quanti cational real-time polymerase chain reaction (qRT-PCR) according to the manufacturer's instructions. All the primers were designed and purchased from RiboBio Co. Ltd (Guangzhou, China) and BGI.