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

DOI: https://doi.org/10.21203/rs.3.rs-33112/v1

Abstract

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.

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 final pathology, the prognosis of HCC patients remains dismal [46]. It has been well known that vascular invasion is a significant factor linked to HCC aggressive progression and poor prognosis [79]. 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.

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 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

 

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 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 [2225]. 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.

Conclusion

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.

Materials And Methods

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.

Declaration

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.

References

  1. Villanueva A. Hepatocellular Carcinoma. N Engl J Med. 2019;380:1450–62.
  2. Massarweh NN, El-Serag HB. Epidemiology of Hepatocellular Carcinoma and Intrahepatic Cholangiocarcinoma. Cancer Control. 2017;24:1073274817729245.
  3. Zhang B, Zhang B, Zhang Z, Huang Z, Chen Y, Chen M, Bie P, Peng B, Wu L, Wang Z, Li B, Fan J, Qin L, et al. 42,573 cases of hepatectomy in China: a multicenter retrospective investigation. Sci China Life Sci. 2018;61:660–70.
  4. Kanwal F, Singal AG. Surveillance for Hepatocellular Carcinoma: Current Best Practice and Future Direction. Gastroenterology. 2019;157:54–64.
  5. Yang JD, Hainaut P, Gores GJ, Amadou A, Plymoth A, Roberts LR. A global view of hepatocellular carcinoma: trends, risk, prevention and management. Nat Rev Gastroenterol Hepatol. 2019;16:589–604.
  6. Doycheva I, Thuluvath PJ. Systemic Therapy for Advanced Hepatocellular Carcinoma: An Update of a Rapidly Evolving Field. J Clin Exp Hepatol. 2019;9:588–96.
  7. Finn RS, Zhu AX, Farah W, Almasri J, Zaiem F, Prokop LJ, Murad MH, Mohammed K. Therapies for advanced stage hepatocellular carcinoma with macrovascular invasion or metastatic disease: A systematic review and meta-analysis. Hepatology. 2018;67:422–35.
  8. Goh BK, Chow PK, Teo JY, Wong JS, Chan CY, Cheow PC, Chung AY, Ooi LL. Number of nodules, Child-Pugh status, margin positivity, and microvascular invasion, but not tumor size, are prognostic factors of survival after liver resection for multifocal hepatocellular carcinoma. J Gastrointest Surg. 2014;18:1477–85.
  9. Koizumi S, Yamashita S, Matsumura S, Takeda K, Minagawa T, Kobayashi S, Hibi T, Shinoda M, Endo I, Tanabe M, Yamamoto M, Otsubo T. Significance of a preoperative tumor marker gradient for predicting microvascular invasion in cases of hepatocellular carcinoma. Mol Clin Oncol. 2020;12:290–4.
  10. Banerjee S, Wang DS, Kim HJ, Sirlin CB, Chan MG, Korn RL, Rutman AM, Siripongsakun S, Lu D, Imanbayev G, Kuo MD. A computed tomography radiogenomic biomarker predicts microvascular invasion and clinical outcomes in hepatocellular carcinoma. Hepatology. 2015;62:792–800.
  11. Yin J, Wang L, Zhu JM, Yu Q, Xue RY, Fang Y, Zhang YA, Chen YJ, Liu TT, Dong L, Shen XZ. Prp19 facilitates invasion of hepatocellular carcinoma via p38 mitogen-activated protein kinase/twist1 pathway. Oncotarget. 2016;7:21939–51.
  12. Yao Z, Luo J, Hu K, Lin J, Huang H, Wang Q, Zhang P, Xiong Z, He C, Huang Z, Liu B, Yang Y. ZKSCAN1 gene and its related circular RNA (circZKSCAN1) both inhibit hepatocellular carcinoma cell growth, migration, and invasion but through different signaling pathways. Mol Oncol. 2017;11:422–37.
  13. Zhang Y, Wei Y, Li X, Liang X, Wang L, Song J, Zhang X, Zhang C, Niu J, Zhang P, Ren Z, Tang B. microRNA-874 suppresses tumor proliferation and metastasis in hepatocellular carcinoma by targeting the DOR/EGFR/ERK pathway. Cell Death Dis. 2018;9:130.
  14. Zhang Z, Zhang Y, Sun XX, Ma X, Chen ZN. microRNA-146a inhibits cancer metastasis by downregulating VEGF through dual pathways in hepatocellular carcinoma. Mol Cancer. 2015;14:5.
  15. Lou W, Chen J, Ding B, Chen D, Zheng H, Jiang D, Xu L, Bao C, Cao G, Fan W. Identification of invasion-metastasis-associated microRNAs in hepatocellular carcinoma based on bioinformatic analysis and experimental validation. J Transl Med. 2018;16:266.
  16. Zhu GF, Xu YW, Li J, Niu HL, Ma WX, Xu J, Zhou PR, Liu X, Ye DL, Liu XR, Yan T, Zhai WK, Xu ZJ, et al. Mir20a/106a-WTX axis regulates RhoGDIa/CDC42 signaling and colon cancer progression. Nat Commun. 2019;10:112.
  17. Zhang K, Dong C, Chen M, Yang T, Wang X, Gao Y, Wang L, Wen Y, Chen G, Wang X, Yu X, Zhang Y, Wang P, et al. Extracellular vesicle-mediated delivery of miR-101 inhibits lung metastasis in osteosarcoma. Theranostics. 2020;10:411–25.
  18. Ben Khelifa H, Soyah N, Ben-Abdallah-Bouhjar I, Gritly R, Sanlaville D, Elghezal H, Saad A, Mougou-Zerelli S. Xp22.3 interstitial deletion: a recognizable chromosomal abnormality encompassing VCX3A and STS genes in a patient with X-linked ichthyosis and mental retardation. Gene. 2013;527:578–83.
  19. Gonzalez-Huerta L, Mendiola-Jimenez J, Del Moral-Stevenel M, Rivera-Vega M, Cuevas-Covarrubias S. Atypical X-linked ichthyosis in a patient with a large deletion involving the steroid sulfatase (STS) gene. Int J Dermatol. 2009;48:142–4.
  20. Deng H, Zeng J, Zhang T, Gong L, Zhang H, Cheung E, Jones C, Li G. Histone H3.3K27M Mobilizes Multiple Cancer/Testis (CT) Antigens in Pediatric Glioma. Mol Cancer Res. 2018;16:623–33.
  21. Taguchi A, Taylor AD, Rodriguez J, Celiktas M, Liu H, Ma X, Zhang Q, Wong CH, Chin A, Girard L, Behrens C, Lam WL, Lam S, et al. A search for novel cancer/testis antigens in lung cancer identifies VCX/Y genes, expanding the repertoire of potential immunotherapeutic targets. Cancer Res. 2014;74:4694–705.
  22. Wang L, Li B, Zhang L, Li Q, He Z, Zhang X, Huang X, Xu Z, Xia Y, Zhang Q, Li Q, Xu J, Sun G, et al. miR-664a-3p functions as an oncogene by targeting Hippo pathway in the development of gastric cancer. Cell Prolif. 2019;52:e12567.
  23. Bao Y, Chen B, Wu Q, Hu K, Xi X, Zhu W, Zhong X, Chen J. Overexpression of miR-664 is associated with enhanced osteosarcoma cell migration and invasion ability via targeting SOX7. Clin Exp Med. 2017;17:51–8.
  24. Wu L, Li Y, Li J, Ma D. MicroRNA-664 Targets Insulin Receptor Substrate 1 to Suppress Cell Proliferation and Invasion in Breast Cancer. Oncol Res. 2019;27:459–67.
  25. Yang Y, Liu H, Wang X, Chen L. Up-regulation of microRNA-664 inhibits cell growth and increases cisplatin sensitivity in cervical cancer. Int J Clin Exp Med. 2015;8:18123–9.
  26. He JH, Han ZP, Liu JM, Zhou JB, Zou MX, Lv YB, Li YG, Cao MR. Overexpression of Long Non-Coding RNA MEG3 Inhibits Proliferation of Hepatocellular Carcinoma Huh7 Cells via Negative Modulation of miRNA-664. J Cell Biochem. 2017;118:3713–21.
  27. Wang X, Zhou Z, Zhang T, Wang M, Xu R, Qin S, Zhang S. Overexpression of miR-664 is associated with poor overall survival and accelerates cell proliferation, migration and invasion in hepatocellular carcinoma. Onco Targets Ther. 2019;12:2373–81.
  28. Yang H, Cho ME, Li TW, Peng H, Ko KS, Mato JM, Lu SC. MicroRNAs regulate methionine adenosyltransferase 1A expression in hepatocellular carcinoma. J Clin Invest. 2013;123:285–98.
  29. Liang W, Liao Y, Li Z, Wang Y, Zheng S, Xu X, Ran F, Tang B, Wang Z. MicroRNA-644a promotes apoptosis of hepatocellular carcinoma cells by downregulating the expression of heat shock factor 1. Cell Commun Signal. 2018;16:30.
  30. Abuzenadah A, Al-Saedi S, Karim S, Al-Qahtani M. Role of Overexpressed Transcription Factor FOXO1 in Fatal Cardiovascular Septal Defects in Patau Syndrome: Molecular and Therapeutic Strategies. Int J Mol Sci. 2018; 19.
  31. Kuo HC, Li SC, Guo MM, Huang YH, Yu HR, Huang FC, Jiao F, Kuo HC, Andrade J, Chan WC. Genome-Wide Association Study Identifies Novel Susceptibility Genes Associated with Coronary Artery Aneurysm Formation in Kawasaki Disease. PLoS One. 2016;11:e0154943.