DOI: https://doi.org/10.21203/rs.3.rs-33112/v1
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 first found that high expression of VCX3A was significantly 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 identified, 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.
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 final pathology, the prognosis of HCC patients remains dismal [4–6]. It has been well known that vascular invasion is a significant factor linked to HCC aggressive progression and poor prognosis [7–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 specific candidate genes and contribute to influence 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 identified. 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.
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 classified 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 |log2FC| > 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. As shown in Figure 3A-B, among 29 candidate DEGs, only 5 genes (VCX3A, PNCK, SLC6A3, MMP1, CTSE) show positively associations between high expression and significantly unfavorable overall survival (OS) in HCC patients. Whereas overexpression of HHATL, SOHLH1, CDH10, STAC2, NKX2-5, RNF151, CDK5R2, PHOX2A, HRK, TRIM63, KLK15, GAD2, ISL1, ATOH1, HOXD13, GPR50 show better OS in HCC. Furthermore, no significant relationships between the expression of PSCA, CALCA, DUOX2, CHST5 and their prognostic values in HCC were found. In addition, no prognostic values of CXorf67, NAT16, HNRNPCL3 and FOXL2NB were found in Kaplan-Meier plotter. Thus, we only selected VCX3A, PNCK, SLC6A3, MMP1, CTSE for the next analysis. Subsequently, for improving accuracy of our analysis, we intended to further evaluated the prognostic values of VCX3A, PNCK, SLC6A3, MMP1 and CTSE only in HCC patients with vascular invasion in TCGA. All the gene expression information and survival clinical information of HCC patients with vascular invasion in TCGA were extracted alone. Patients were divided into down-expression groups (n=54) and up-expression groups (n=55) according the median values of hub gene expressions. Finally, survival analyses were performed through GraphPad prism software (version 6). As the results presented in Figure 3C–G, only increased mRNA expressions of VCX3A and SLC6A3 contribute to worse OS in vascular invasion-HCC patients. However, no significant survival values of PNCK, MMP1 and CTSE were found. Lastly, we validated the expression values of VCX3A and SLC6A3 in 10 none-invasion liver tissues and 10 vascular invasion-HCC tumor tissues in clinical samples (Figure 3H-I). Notably, we found that VCX3A was high-expressed in tumor tissues than normal liver tissues, while for SLC6A3, no significant difference was discovered. In conclusion, VCX3A is supposed to be the most significant gene which plays function related with HCC progression and metastasis. Besides, clinical characteristics among VCX3A and HCC in TCGA also showed that high expression of VCX3A was more likely to cause vascular invasion in HCC patients (Table 1).
Table 1: Correlations the clinical characteristics among VCX3A and HCC in TCGA. (The significant P value is marked with Bold type. NA=Not Applicable.)
Variables |
|
HCC |
|
VCX3A |
|||
low/high expression case (n) |
P value |
||
Gender |
N |
|
|
Male |
250 |
231/19 |
0.2079
|
Female |
121 |
107/14 |
|
Age at diagnosis |
|
|
|
>=60 |
201 |
179/22 |
0.1359
|
<60 |
169 |
158/11 |
|
NA |
1 |
1/0 |
|
T stage |
|
|
|
T1/T2 |
275 |
252/23 |
0.4858
|
T3/T4 |
93 |
83/10 |
|
TX |
1 |
1/0 |
|
NA |
2 |
2/0 |
|
N stage |
|
|
|
N0 |
252 |
233/19 |
0.2792
|
N1 |
4 |
3/1 |
|
NX |
114 |
101/13 |
|
NA |
1 |
1/0 |
|
M stage |
|
|
|
M0 |
266 |
245/21 |
1.0000
|
M1 |
4 |
4/0 |
|
MX |
101 |
89/12 |
|
Pathologic stage |
|
|
|
I/II |
257 |
235/22 |
0.3081
|
III/IV |
90 |
79/11 |
|
NA |
24 |
24/0 |
|
Vascular invasion |
|
|
|
None |
206 |
194/12 |
0.0313
|
Micro/Macro |
109 |
95/14 |
|
NA |
56 |
49/7 |
Identifying hsa-miR-664a-3p-VCX3A axis play a key role in HCC
First, we predicted upstream miRNAs of VCX3A through TargetScan, miRWalk and miRDB three databases, since increasing studies have indicated that miRNAs take part in the invasion and metastasis of multiple tumors and play significant functions in cancer development [15-17]. As presented in Figure 4A-B, a total of 15 miRNAs (hsa-miR-2861, hsa-miR-5739, hsa-miR-4505, hsa-miR-6770-3p, hsa-miR-7855-5p, hsa-miR-4299, hsa-miR-6871-3p, hsa-miR-31-5p, hsa-miR-181a-5p, hsa-miR-450b-5p, hsa-miR-3615, hsa-miR-5696, hsa-miR-664b-3p, hsa-miR-579-3p, hsa-miR-664a-3p) appeared twice or more than twice in three predicting databases were selected for the further analysis. After that, we evaluated the prognostic values of 15 miRNAs in HCC using Kaplan-Meier plotter. The analytic results showed in Figure 4C-D indicated that high expression of 11 miRNAs including hsa-miR-2861, hsa-miR-5739, hsa-miR-4505, hsa-miR-6770-3p, hsa-miR-7855-5p, hsa-miR-4299, hsa-miR-6871-3p, hsa-miR-31-5p, hsa-miR-5696, hsa-miR-664b-3p, hsa-miR-664a-3p were found to be positively correlated with longer OS in HCC. On the contrary, upregulated the expression of hsa-miR-579-3p was found to be negatively associated with better OS in HCC. Additionally, no significant correlations between the expression values of hsa-miR-181a-5p, hsa-miR-450b-5p, hsa-miR-3615 and their prognostic values in HCC were found. Therefore, we chose hsa-miR-2861, hsa-miR-5739, hsa-miR-4505, hsa-miR-6770-3p, hsa-miR-7855-5p, hsa-miR-4299, hsa-miR-6871-3p, hsa-miR-31-5p, hsa-miR-5696, hsa-miR-664b-3p and hsa-miR-664a-3p for subsequent research. We particularly focused on miRNAs which were greatly related with HCC progression. As clinical parameters like histologic grade, pathologic stage, pathologic T status, pathologic N status, pathologic M status are indicators of tumor progression. Thus, we used the “Tumor Stage and Grade” analysis on OncomiR database to discover key miRNAs in HCC, and the results indicated that there were 0 clinical parameters in HCC associated with hsa-miR-2861, hsa-miR-5739, hsa-miR-4505, hsa-miR-6770-3p, hsa-miR-7855-5p, hsa-miR-4299, hsa-miR-6871-3p, hsa-miR-5696, hsa-miR-664b-3p (Figure 5A). Importantly, hsa-miR-31-5p was found to be closely related with histologic grade, pathologic stage, pathologic T status, pathologic N status, pathologic M status and sex clinical parameters in HCC and hsa-miR-664a-3p was significantly correlated with histologic grade, pathologic T status, pathologic N status and pathologic M status in HCC (Figure 5B-C). Then, we evaluated the expression values of hsa-miR-31-5p and hsa-miR-664a-3p in 10 none invasion-HCC tissues and 10 vascular invasion-HCC tissues (Figure 5D-E). We found that only hsa-miR-664a-3p was significantly downregulated in vascular invasion-HCC tissues compared with none invasion-HCC tissues, as expected. However, the expression of hsa-miR-31-5p was found to be significantly higher in vascular invasion-HCC than in none invasion-HCC tissues. It’s obvious that target gene negatively regulated by the miRNA, since overexpression the expression value of VCX3A was closely related with poor prognosis in HCC patients and VCX3A was importantly downregulated in none invasion-HCC tissues than in vascular invasion-HCC tissues, the upstream miRNA should downregulated in vascular invasion-HCC tissues than in none invasion-HCC tissues, like hsa-miR-664a-3p. Thus, hsa-miR-664a-3p serves as the most potential upstream miRNA binding to VCX3A.
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 find 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. The genes list as following: MAGEA6, MKRN3, LIN28B, CTAG2, DCAF8L1, MAGEA12, RFPL4B, DDX53, MAGEA3, GABRA3, SYT1, SSX1, CSMD1, CSAG1, FAM133A, C12ORF56, TUBA3C, GTSF1, MAGEB6, CLEC2L, ZNF257. However, 28 genes including VCX, VCX3B, NAA11, CARD18, RPL10L, DCAF4L2, DSCR4, MAGEB2, GNGT1, OR56A3, CNTNAP4, TGIF2LX, DCAF8L2, CT47B1, SLCO6A1, MAGEB1, SSX5, MAGEC1, MDGA2, OR8A1, XAGE5, NPFFR2, MAGEA11, DSCR8, MAGEA8, PAGE2, PAGE1 and GLB1L3 were found to be associated with better prognosis in HCC, not as expected. In addition, we found that there were no significant correlations between the expression values of COX7B2, MAGEC2, FSTL5, PAGE5 and prognosis in HCC. Moreover, 4 genes (XAGE1B, PPP4R3C, LINC01139, CT83) lose prognostic values information in Kaplan-Meier plotter (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).
Table 2 : Co-expression genes of VCX3A in cBioPotral.
correlated gene |
Spearman's Correlation |
p-Value |
q-Value |
|
VCX VCX3B NAA11 XAGE1B COX7B2 CARD18 RPL10L DCAF4L2 MAGEA6 MKRN3 SSX4 LIN28B CTAG2 DCAF8L1 CTAG1B DSCR4 MAGEB2 MAGEC2 MAGEA12 RFPL4B DDX53 PPP4R3C MAGEA3 GABRA3 SYT1 SSX1 SSX3 GNGT1 OR56A3 SSX2 MAGEA2 CNTNAP4 GAGE4 TGIF2LX CSMD1 DCAF8L2 CSAG1 FAM133A SSX6P LINC01139 BAGE FSTL5 CT47B1 CDH9 SLCO6A1 RHOXF2B LOC440173 MAGEB1 FLJ36000 SSX5 MAGEC1 C12ORF56 GAGE12J GABRG2 TUBA3C EZHIP SLC44A5 MDGA2 TP53TG3B GTSF1 MAGEB6 OR8A1 SOHLH2 XAGE5 NPFFR2 MAGEA11 DSCR8 CLEC2L MAGEA8 PAGE2 HERC2P4 CSAG3 MUC15 KCNK9 ZNF257 PAGE5 CNGB3 ANKRD30BP2 NXF2 ZNF716 GAGE2D PAGE1 GAGE12D COL24A1 PHF2P1 CSAG2 GAGE1 CT83 ADAMTS20 NAP1L6 SSX8P MAGEB16 HTR2C GLB1L3 |
0.807 0.701 0.694 0.675 0.671 0.67 0.66 0.653 0.653 0.651 0.645 0.645 0.64 0.639 0.636 0.636 0.636 0.632 0.631 0.63 0.63 0.629 0.629 0.628 0.628 0.628 0.628 0.624 0.622 0.618 0.618 0.617 0.617 0.615 0.612 0.61 0.605 0.603 0.599 0.592 0.589 0.589 0.589 0.588 0.582 0.581 0.581 0.577 0.576 0.575 0.574 0.573 0.572 0.57 0.568 0.568 0.568 0.566 0.562 0.559 0.558 0.556 0.554 0.552 0.549 0.548 0.547 0.544 0.543 0.538 0.537 0.536 0.535 0.535 0.535 0.532 0.532 0.526 0.52 0.52 0.52 0.515 0.512 0.507 0.506 0.505 0.503 0.503 0.501 0.501 0.5 0.5 0.5 0.5 |
3.67E-81 8.80E-53 2.63E-51 1.26E-47 6.02E-47 8.97E-47 5.86E-45 9.57E-44 1.11E-43 2.94E-43 2.98E-42 2.99E-42 1.53E-41 2.21E-41 7.56E-41 8.34E-41 8.62E-41 3.70E-40 4.51E-40 5.97E-40 6.54E-40 1.04E-39 1.11E-39 1.26E-39 1.28E-39 1.50E-39 1.57E-39 6.27E-39 1.19E-38 4.66E-38 5.70E-38 5.95E-38 7.18E-38 1.54E-37 3.84E-37 7.47E-37 4.59E-36 7.66E-36 2.47E-35 2.66E-34 6.11E-34 6.92E-34 7.98E-34 9.78E-34 6.19E-33 7.87E-33 9.29E-33 2.87E-32 3.39E-32 5.72E-32 6.03E-32 8.88E-32 1.39E-31 2.22E-31 3.60E-31 3.68E-31 4.53E-31 8.30E-31 2.50E-30 4.72E-30 7.36E-30 1.15E-29 2.12E-29 3.75E-29 8.06E-29 1.08E-28 1.42E-28 3.17E-28 4.24E-28 1.58E-27 1.97E-27 2.64E-27 3.40E-27 3.56E-27 3.59E-27 7.38E-27 8.87E-27 3.65E-26 1.55E-25 1.57E-25 1.88E-25 5.21E-25 1.27E-24 3.58E-24 5.42E-24 6.35E-24 1.03E-23 1.15E-23 1.48E-23 1.66E-23 1.85E-23 2.06E-23 2.08E-23 2.25E-23 |
7.37E-77 8.83E-49 1.76E-47 6.35E-44 2.42E-43 3.00E-43 1.68E-41 2.40E-40 2.48E-40 5.91E-40 5.00E-39 5.00E-39 2.36E-38 3.17E-38 1.01E-37 1.02E-37 1.02E-37 4.13E-37 4.77E-37 6.00E-37 6.26E-37 9.51E-37 9.72E-37 1.03E-36 1.03E-36 1.16E-36 1.17E-36 4.50E-36 8.23E-36 3.12E-35 3.69E-35 3.73E-35 4.37E-35 9.07E-35 2.21E-34 4.17E-34 2.49E-33 4.05E-33 1.27E-32 1.34E-31 2.99E-31 3.31E-31 3.73E-31 4.46E-31 2.76E-30 3.44E-30 3.97E-30 1.20E-29 1.39E-29 2.30E-29 2.37E-29 3.43E-29 5.26E-29 8.25E-29 1.31E-28 1.32E-28 1.60E-28 2.88E-28 8.53E-28 1.58E-27 2.42E-27 3.73E-27 6.77E-27 1.18E-26 2.49E-26 3.30E-26 4.27E-26 9.36E-26 1.23E-25 4.54E-25 5.56E-25 7.36E-25 9.35E-25 9.61E-25 9.61E-25 1.95E-24 2.31E-24 9.40E-24 3.94E-23 3.94E-23 4.66E-23 1.28E-22 3.07E-22 8.55E-22 1.28E-21 1.48E-21 2.37E-21 2.61E-21 3.34E-21 3.70E-21 4.08E-21 4.50E-21 4.50E-21 4.81E-21 |
|
Table 3: Co-expression analysis of VCX3A in starBase. (NA=Not Applicable)
VCX3A co-expression gene |
R |
P |
VCX VCX3B NAA11 XAGE1B COX7B2 CARD18 RPL10L DCAF4L2 MAGEA6 MKRN3 SSX4 LIN28B CTAG2 DCAF8L1 CTAG1B DSCR4 MAGEB2 MAGEC2 MAGEA12 RFPL4B DDX53 PPP4R3C MAGEA3 GABRA3 SYT1 SSX1 SSX3 GNGT1 OR56A3 SSX2 MAGEA2 CNTNAP4 GAGE4 TGIF2LX CSMD1 DCAF8L2 CSAG1 FAM133A SSX6P LINC01139 BAGE FSTL5 CT47B1 CDH9 SLCO6A1 RHOXF2B LOC440173 MAGEB1 FLJ36000 SSX5 MAGEC1 C12ORF56 GAGE12J GABRG2 TUBA3C EZHIP SLC44A5 MDGA2 TP53TG3B GTSF1 MAGEB6 OR8A1 SOHLH2 XAGE5 NPFFR2 MAGEA11 DSCR8 CLEC2L MAGEA8 PAGE2 HERC2P4 CSAG3 MUC15 KCNK9 ZNF257 PAGE5 CNGB3 ANKRD30BP2 NXF2 ZNF716 GAGE2D PAGE1 GAGE12D COL24A1 PHF2P1 CSAG2 GAGE1 CT83 ADAMTS20 NAP1L6 SSX8P MAGEB16 HTR2C GLB1L3 |
r=0.771 r=0.813 r=0.632 r=0.615 r=0.612 r=0.606 r=0.631 r=0.568 r=0.637 r=0.504 r=0.383 r=0.623 r=0.577 r=0.561 r=0.224 r=0.714 r=0.626 r=0.589 r=0.622 r=0.622 r=0.534 r=0.558 r=0.620 r=0.635 r=0.518 r=0.586 r=0.480 r=0.625 r=0.539 r=0.481 r=0.447 r=0.575 NA r=0.596 r=0.581 r=0.535 r=0.604 r=0.531 NA r=0.606 NA r=0.546 r=0.513 r=0.497 r=0.573 r=0.487 NA r=0.565 NA r=0.545 r=0.548 r=0.608 r=0.308 r=0.445 r=0.687 NA r=0.467 r=0.540 r=0.000 r=0.655 r=0.513 r=0.570 r=0.487 r=0.593 r=0.622 r=0.536 r=0.664 r=0.559 r=0.672 r=0.598 r=0.437 NA r=0.490 r=0.487 r=0.513 r=0.560 r=0.497 r=0.386 r=0.225 r=0.448 NA r=0.560 r=0.000 r=0.479 NA NA r=0.461 r=0.522 r=0.474 r=0.485 NA r=0.467 r=0.449 r=0.505 |
p=8.44e-75 p=2.91e-89 p=3.72e-43 p=2.90e-40 p=8.01e-40 p=8.52e-39 p=6.29e-43 p=2.38e-33 p=6.97e-44 p=1.59e-25 p=1.65e-14 p=1.31e-41 p=1.21e-34 p=1.95e-32 p=1.17e-05 p=1.69e-59 p=3.85e-42 p=2.71e-36 p=2.04e-41 p=2.29e-41 p=5.41e-29 p=5.74e-32 p=4.31e-41 p=1.56e-43 p=4.27e-27 p=6.63e-36 p=5.60e-23 p=6.56e-42 p=1.34e-29 p=5.31e-23 p=9.43e-20 p=2.91e-34 NA p=2.29e-37 p=4.16e-35 p=4.48e-29 p=1.51e-38 p=1.24e-28 NA p=7.88e-39 NA p=1.71e-30 p=1.88e-26 p=9.66e-25 p=4.89e-34 p=1.22e-23 NA p=5.90e-33 NA p=2.77e-30 p=1.09e-30 p=3.21e-39 p=1.19e-09 p=1.36e-19 p=1.56e-53 NA p=1.13e-21 p=1.10e-29 p=1.00e+00 p=4.18e-47 p=1.69e-26 p=1.40e-33 p=1.11e-23 p=5.97e-37 p=2.39e-41 p=3.68e-29 p=6.34e-39 p=4.56e-32 p=1.81e-50 p=1.19e-37 p=6.81e-19 NA p=5.11e-24 p=1.02e-23 p=1.81e-26 p=2.92e-32 p=1.01e-24 p=9.19e-15 p=1.08e-05 p=7.87e-20 NA p=2.59e-32 p=1.00e+00 p=7.67e-23 NA NA p=4.90e-21 p=1.58e-27 p=2.49e-22 p=2.06e-23 NA p=1.23e-21 p=6.47e-20 p=1.26e-25 |
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 identified the correlation of hub genes with the progression or metastasis of HCC, respectively, while we first 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 first 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 significantly 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 significant downstream molecule (TUBA3C) interact with VCX3A for the first time. Lots of researches have been launched to investigated the functions of miR-664a-3p in multiple tumors and different analyses found different roles of miR-664a-3p in cancer [22–25]. In hepatocellular carcinoma, there have been many disputes arising out of the roles of miR-664a-3p. For example, He, J. H, et al. demonstrated that long non-coding RNA MEG3 might act as a suppressor of HCC cells through sponging miR-664 that mediated the regulation of ADH4 [26]. Wang et al. indicated that overexpression of miR-664 could promote the growth, invasion and migration of HCC cells [27]. In addition, Yang et al. suggested that knockdown of miR-664 could reduced proliferation and invreased apoptosis of HCC cells through inducing MAT1A expression [28]. All of the evidences above demonstrated miR-664 might function as a promoter in the development of HCC. However, research of Liang et al. supported that miR-644a was lower expression in HCC tissues 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.
In conclusion, by performing comprehensive bioinformatic analysis and experiment validation, we for the first 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.
Differential expression validation in TCGA
The edgeR package (http://bioconductor.org/packages/edgeR/) in R software was employed to gain the differentially expressed genes (DEGs) and differentially expressed miRNAs (DEmiRNAs) in The Cancer Genome Atlas (TCGA) database (https://cancergenome.nih.gov/). |Fold change (FC)| > 1.5 and p values < 0.05 were set as cut-off criteria. HCC patients were subdivided into different groups, such as early stage HCC (n=256), advanced stage HCC (n=90), none invasion-HCC (n=206) or vascular invasion-HCC (n=109), according to the clinical pathological information downloaded from TCGA. Finally, upregulated genes and downregulated genes were inputted into VENNY 2.1.0 (http://bioinfogp.cnb.csic.es/tools/venny/index.html) to gain the Veen diagrams, respectively.
Survival analysis
Kaplan-Meier plotter (http://kmplot.com/analysis/index.php?p=background) was utilized to evaluate the prognostic influences 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 classified 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 significant.
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 Affiliated 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 Affiliated 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 quantificational 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.
Statistical analysis
The paired Student’s t-test and Chi-square test methods were employed to analysis the data using GraphPad prism software (version 6), the results were shown as mean ± SD. Value of p < 0.05 was considered as a statistically different result.
Acknowledgements
Not applicable
Funding
This research was supported by the National Natural Science Foundation of China (81874225).
Availability of data and materials
Not applicable
Authors’ Contributions
Bisha Ding, Weiyang Lou: main text and illustration; Weiyang Lou and Weimin Fan: concept, discussion.
Ethics approval and consent to participate
This study was approved by the ethics committee of the First Affiliated Hospital of Zhejiang University, College of Medicine.
Patient consent for publication
Not applicable.
Competing interests
The authors declare that there is no conflict of interests regarding the publication of this paper.