Identification of portal vein tumor thrombus with an independent clonal origin in hepatocellular carcinoma via multi-omics data analysis

Shupeng Liu, Zaixin Zhou, Yin Jia, Jie Xue, Zhiyong Liu, Kai Cheng, Shuqun Cheng, Shanrong Liu

1Clinical Research Center, Changhai Hospital, Second Military Medical University, Shanghai 200433, China; 2Department of Laboratory Diagnostics, Changhai Hospital, Second Military Medical University, Shanghai 200433, China; 3Department of Hepatic Surgery VI, Eastern Hepatobiliary Surgery Hospital, Second Military Medical University, Shanghai 200433, China

ABSTRACT Objective:Multiple mechanisms underlying the development of portal vein tumor thrombus (PVTT) in hepatocellular carcinoma(HCC) have been reported recently. However, the origins of PVTT remain unknown. Increasing multi-omics data on PVTTs in HCCs have made it possible to investigate whether PVTTs originate from the corresponding primary tumors (Ts).Methods:The clonal relationship between PVTTs and their corresponding primary Ts was investigated using datasets deposited in public databases. One DNA copy number variations dataset and three gene expression datasets were downloaded for the analyses.Clonality analysis was performed to investigate the clonal relationship between PVTTs and Ts from an individual patient.Differential gene expression analysis was applied to investigate the gene expression profiles of PVTTs and Ts.Results:One out of 19 PVTTs had no clonal relationship with its corresponding T, whereas the others did. The PVTTs with independent clonal origin showed different gene expression and enrichment in biological processes from the primary Ts. Based on the unique gene expression profiles, a gene signature including 24 genes was used to identify pairs of PVTTs and primary Ts without any clonal relationship. Validation in three datasets showed that these types of pairs of PVTTs and Ts can be identified by the 24-gene signature.Conclusions:Our findings show a direct evidence for PVTT origin and consolidate the heterogeneity of PVTTs observed in clinic.The results suggest that PVTT investigation at a molecular level is clinically necessary for diagnosis and treatment.

KEYWORDS Hepatocellular carcinoma; portal vein tumor thrombus; clonal origin; copy number variation; bioinformatics

Introduction

Portal vein tumor thrombus (PVTT) is one of the main hepatocellular carcinoma (HCC) complications, found in approximately 50% of HCC patients1. The prognosis of HCC patients with PVTT is poor, and the median overall survival(OS) time of patients with PVTT was primarily less than 1 year2,3. Although clinicians in Europe and America recommend sorafenib alone, and multidisciplinary therapy has been adapted in Asia4,5, no worldwide consensus or guidelines have shown effectiveness for specific PVTT treatment, which might be due to a lack of knowledge about the nature of PVTT.

Although currently, hypoxia, non-coding RNAs, and cancer stem cells have been found to contribute to PVTT development6-9, the clonal origins of PVTT still remain unknown. The identical gene expression profiles between PVTT and the corresponding primary HCCs demonstrated that PVTT could originate from HCC primary nodules by metastasis10,11. Contradictorily, in small HCC patients, some PVTTs are clinically defined, that are distant from the liver parenchyma tumor (PT) nodules, and some PVTTs are found without liver parenchyma tumor nodules12,13.Moreover, some PVTTs showed significantly different gene expression profiles compared with those of the paired primary tumors (Ts)11. These findings indicate that PVTTs display relatively high heterogeneity and may have different clonal origins from the paired primary Ts.

Recently, high-throughput technology, such as microarray and sequencing, has been widely used to investigate the mechanisms underlying the development of Ts, including HCC. Among them, 3 published studies involved the genomic or transcriptomic datasets of HCC as well as PVTT9,11,14. In the present study, previously published datasets (GSE77509, GSE69164, GSE77275, and GSE74656),and freshly obtained samples were used for analysis, and we found a type of PVTT that has a different clonal origin from that of the corresponding primary Ts in HCC. Additionally,this type of PVTT was found to express a specific set of genes,which varied from that of the corresponding Ts.

Materials and methods

Data resources

Four previously published datasets (Cytoscan HD array,GSE77275; RNA-seq, GSE77509 and GSE69164; and cDNA microarray, GSE74656) used in this study were downloaded from the NCBI Gene Expression Omnibus (GEO,http://www.ncbi.nlm.nih.gov/geo/). GSE77275 included the copy number variation (CNV) profiles of 60 matched samples (PVTT/T/PT from 20 patients) assessed by Cytoscan HD array. GSE77509 included the gene expression profiles of these 60 samples assessed by RNA-seq. DNA and total RNA from the 60 samples were used for CNV and gene expression analysis, respectively, as previously reported11.

The samples of patient 14 (P14), in these two datasets,were not included in the present study because they were reported to be contaminated11. GSE69164 included the gene expression of 33 matched samples from 11 patients assessed by RNA-seq14. GSE74656 included the gene expression data of 15 matched samples from 5 patients assessed by cDNA analysis.

Clonality testing

To investigate the clonal relationship between PVTTs and Ts from an individual patient, the CNV profiles of 19 pairs of PVTTs and Ts were extracted from GSE77275 using R package Affy2sv and loaded to R package Clonality for further analysis. The CNV of the chromosome region in each patient was estimated using the circular binary segmentation(CBS) analysis based on the normalized weighted Log2ratio.A significant CNV in chromosome regions was determined according to the threshold value (median absolute deviations = 1.25). The clonal relationship between PVTTs and Ts from an individual patient was investigated according to the likelihood ratio 2 (LR2)15,16.

Differential gene expression analysis

The gene expression profiles of PVTTs and Ts from the datasets GSE77509, GSE69164, and GSE74656 were extracted using the appropriate methods11. EBSeq from NovelBrain BioCloud (https://cloud.novelbrain.com) was applied to analyze differentially expressed genes between PVTTs and Ts from individual patients in GSE77509. Genes with fold changes ＞ 3 and FDR ＜ 0.05 were identified as significantly differentially expressed.

Gene ontology and KEGG pathway enrichment analyses

The differentially expressed genes between PVTT13 and T13 were loaded to Metascape for Gene ontology (GO) and KEGG pathway enrichment analyses17. Protein-protein interaction (PPI) network analysis was also performed in Metascape, and the PPI network was exported using Cytoscape 3.4.

Gene signature enrichment analysis

Gene expression signature enrichment was evaluated using single-sample gene set enrichment analysis (ssGSEA) and nearest template prediction (NTP) analyses (Gene Pattern modules)18. Molecular Signature Database gene sets,(MSigDB, www.broadinstitute.org/msigdb ) including hematopoietic stem cells and liver cancer stem cells, and previously published gene-expression signatures representing different cell types including bile duct cells (BDC),hepatocytes, hepatic progenitor cells (HPC), and hepatic stellate cells (HSC)19(Supplementary Table S1), were tested using ssGSEA. Previously, published HCC molecular classifications were also analyzed using NTP analysis20,21.

Tumor tissues

Eight pairs of human PVTTs and Ts were obtained from HCC patients who underwent curative resection in the Eastern Hepatobiliary Surgery Hospital (Shanghai, China)after obtaining informed consent. All 8 pairs of PVTTs and Ts were used to perform real-time polymerase chain reaction(PCR) analysis. The study was approved by the Ethics Committee of Eastern Hepatobiliary Surgery Hospital.

Quantitative real-time PCR

Total RNA was isolated using TRIzol reagent (Invitrogen,CA, USA), and reverse-transcription (RT) was performed using the PrimeScript RT Reagent Kit (Takara, Kusatsu,Japan), according to manufacturer's instructions. Real-time PCR was performed with SYBR Premix Ex Taq (Takara,Kusatsu, Japan) using the StepOnePlus Real-Time PCR System (Applied Biosystems, CA, USA). The primers used are listed in Supplementary Table S2. The gene expression levels were calculated relative to the expression of β-actin.

Cluster analysis

All pairs of PVTT and Ts were classified using one minus Pearson correlation hierarchical cluster analysis by the 24-gene signatures (Morpheus). The 24 genes are listed in Supplementary Table S3.

Statistical analysis

All P values and FDR values were obtained by the appropriate methods used in the different gene expression and enrichment analyses. Statistical analysis of the gene signature enrichment score from ssGSEA was not performed due to only one PVTT-T pair identified to have no clonal relationship.

Results

Analysis of the clonal relationship between PVTTs and corresponding primary Ts

To investigate the clonal relationship of PVTTs with corresponding Ts, DNA copy number variation (CNV)profiles of 19 paired PVTTs and Ts were analyzed using the R package Clonality22. The clonal-relatedness of the paired tissues was evaluated according to the LR2 value. Overall, 18 out of 19 PVTT-T pairs showed a significantly higher LR2 value (P ＜ 0.001), whereas only one pair (patient 13, P13)showed a much lower LR2 value (P ＞ 0.05) (Table 1). As a low LR2 value indicates no clonal relatedness, PVTT13 was identified to have an independent clonal origin from T13,whereas the other PVTTs had a clonal origin from the corresponding Ts. For further verification, the CNV of the chromosome regions and genes was investigated using the R package DNA copy23. The DNA CNV patterns of chromosome regions were estimated based on the threshold values (threshold = 0.3). Aberrant chromosome regions showed very different CNV patterns between PVTT13 and T13, whereas these differences were not observed in other pairs (Figure 1A). The region of 2q24.1-q31.1 was deleted in PVTT13 but amplified in T13. Regions of 5q13.2-q35.2 and 15q11.2-q21.1 were amplified in PVTT13 but deleted in T13.These three chromosome regions displayed similar CNV patterns between PVTTs and Ts from other patients (Figure 1A). CNV analysis of genes showed that 10 HCC driver genes have different copy numbers between the PVTTs and Ts from individual patients (Figure 1B). Less than 3 different mutant genes were observed in most PVTT-T pairs (n = 15),whereas 3 or more different mutant genes were observed in pairs from patient 6 (P6), patient 11 (P11), P13, and patient 22 (P22) (Figure 1B). Among these 4 PVTT-T pairs with more than 2 mutant genes with different copy numbers, only the pair from P13 had 2 mutant genes, with 2 and 3 different copies. Five copies of ADCY2 and 3 copies of CTNNB1 were observed in PVTT13, whereas 2 copies of ADCY2 and 5 copies of CTNNB1 were observed in T13 (Figure 1B).Additionally, Hoshida S1-S2-S3 classes of the PVTTs and Ts were evaluated by NTP analysis using their gene expression profiles20. All PVTTs and Ts were classified into one of the Hoshida S1-S2-S3 classes (Bonferroni P = 0.019, FDR ＜ 0.01,Figure 1C). Among the 19 PVTT-T pairs, 6 pairs including PVTT13-T13 had their PVTTs and Ts in different subclasses(Figure 1C). PVTT13 was classified into the S3 class, which is associated with hepatocyte differentiation, and T13 was classified into the S1 class, which reflects an aberrant activation of the Wnt-signaling pathway (Figure 1C).Collectively, these data indicated that PVTT13 and T13 had different origins and that some PVTTs may not originate from the primary Ts.

Table 1 Clonality analysis of PVTT and T (n = 19) according to copy number variations at the probe level

Analysis of aberrantly expressed genes in PVTT13

Figure 1 Analysis of the clonal relationship between PVTTs and corresponding primary Ts. (A) Copy number variations (CNVs) of chromosome regions in PVTTs and corresponding primary Ts. The specific CNV regions in PVTT13 and T13 were labeled by rectangles. The CNV threshold = 0.3; red, gain of CNV; blue, loss of CNV. (B) Copy numbers of HCC driver genes in PVTTs and Ts. (C) Molecular subtype of PVTTs and Ts according to Hoshida S1-S2-S3 classes.

To investigate aberrantly expressed genes in PVTTs and Ts,EBSeq was used to analyze the differently expressed genes(DEGs) between PVTTs, Ts, and parenchyma tumor tissue(PT) from individual patients. Gene expression profiles of these tissues (RNA-seq, GSE77509) were downloaded from the GEO database. Genes with fold change ＞ 3 and FDR＜ 0.05 were identified as having significantly different expression. More than 2,000 DEGs were observed in 5 PVTTT pairs including P13, whereas less than 1,000 DEGs were observed in the other pairs (Figure 2A). The Venn diagram in Supplementary Figure S1A shows that 512 and 3,257 DEGs were identified in PVTT13 and T13 compared with PT13, respectively. Among these 1,795 DEGs, 1,006 upregulated and 789 downregulated genes were identified as genes with specifically aberrant expression in PVTT13, and 1,540 DEGs with 595 upregulated and 945 downregulated genes were identified as genes with specifically aberrant expression in T13 (Supplementary Figure S1B). GO analysis showed that upregulated DEGs in PVTT13 were enriched in 15 biological processes (BPs) including the glutamate receptor signaling pathway, cell adhesion, and cell cycle(Figure 2B, Supplementary Table S4). Downregulated DEGs in PVTT13 were enriched in 16 BPs, including positive regulation of immune response, cytokine production, and peptidyl-tyrosine phosphorylation (Figure 2B,Supplementary Table S4). Upregulated DEGs in T13 were enriched in 15 BPs, including nuclear division, cell cycle checkpoints, and regulation of protein kinase activity (Figure 2C, Supplementary Table S4). Downregulated DEGs in T13 were enriched in 17 BPs, including inflammatory response,cell-cell adhesion, and negative regulation of cell proliferation (Figure 2C, Supplementary Table S4). KEGG enrichment analysis showed that upregulated DEGs in PVTT13 were enriched in pathways for nicotine addiction and type II diabetes mellitus, and The downregulated DEGs were enriched in signaling pathways, mainly PI3K-Akt, TGF-β, and p53 signaling (Figure 2D, Supplementary Table S4).The upregulated DEGs in T13 were enriched in signaling pathways, such as IL-17, p53 and cell cycle signaling, and the downregulated DEGs in T13 were enriched in signaling pathways including cGMP-PKG, Hippo and Rap1 (Figure 2E, Supplementary Table S4). The GO and KEGG analyses showed that DEGs in PVTT13 and T13 were enriched in different BPs and pathways. These data suggest that PVTT13 and T13 were different tumor types and that different mechanisms were involved in PVTT13 and T13 development.

Enrichment of hepatocyte gene signature in PVTT13

To further investigate the independent origin of PVTT13,enrichment of the gene signatures of the main cell types in liver tissues were analyzed using ssGSEA. The gene signatures of bile duct cells (BDCs), hepatocytes, hepatic progenitor cells (HPCs), and hepatic stellate cells (HSCs) were included due to their potential role in liver cancer24.

The hepatocyte gene signature showed higher enrichment value in PVTT13 (7223.62) than that in T13 (5985.81), and the HPC gene signature showed lower enrichment value in PVTT13 (2608.81) than that in T13 (3212.07) (Figure 3A).Differences in the BDC and HSC gene signature enrichment values were not observed between PVTT13 and T13 (Figure 3A), whereas other PVTT-T pairs showed similar gene signature enrichment (Figure 3A). Further gene expression analysis showed a higher expression of hepatocyte-specific genes and lower expression of HSC-specific genes in PVTT13 than in T13 (Figure 3B). In addition, previous reports suggest that cancer stem cells (CSCs) and hematopoietic stem cells may be involved in PVTT development8,12. Our findings show that gene markers of these cells, such as NANOG and C-kit, were the nodes of the gene network or aberrantly expressed in PVTT13, further consolidating their potential roles in PVTT development (Supplementary Figure S2).Therefore, enrichment of gene signatures of hematopoietic stem cells and liver cancer stem cells were also investigated. A higher enrichment value of hematopoietic stem cell gene signatures was observed in PVTT13 (3114.07) than in T13(1031.52), whereas a lower enrichment value of the EpCAM+liver cancer cell gene signature was observed in PVTT13(-315.38) than in T13 (1401.61) (Figure 3C). Most PVTTs and Ts, except for PVTT6, PVTT11, PVTT19, T22, and pairs of P10, P16, and P25, had no enrichment of these signatures(Figure 3C). Gene expression analysis showed a higher expression of hematopoietic stem cell-specific genes, such as CD34 and KIT, and a lower expression of EpCAM+liver cancer cell-specific genes, such as EpCAM and CD44 in PVTT13 (Figure 3D). These data suggested that PVTT13 expressed high levels of hepatocyte- and hematopoietic stem cell-specific genes, and that these two cell types may be involved in PVTT13 development.

Gene signatures identifying PVTTs with no clonal relatedness with Ts

Figure 2 Analysis of aberrantly expressed genes in PVTT13 and T13. (A) Number of differently expressed genes between PVTTs and corresponding primary Ts. Top enriched biological processes in PVTT13 (B) and T13 (C) from GO enrichment analysis based on aberrantly expressed genes. Top enriched KEGG pathways in PVTT13 (D) and T13 (E) from KEGG pathway enrichment analysis based on aberrantly expressed genes.

Having found that some PVTTs had no clonal relatedness with the corresponding Ts and that they have unique gene expression profiles, gene signatures identifying whether PVTTs have clonal relatedness with Ts may be very helpful for PVTT treatment. An integration analysis of CNV and gene expression profiles showed that 160 genes have differentially expressed copy number variations between PVTT13 and T13 (Supplementary Figure S3). The following PPI analysis found 24 genes to be the nodes in the network,suggesting that they may have crucial roles in PVTT13 development (Supplementary Figure S3). The ability of the 24-gene signature to identify PVTTs and Ts with no clonal relatedness was validated using three publicly available datasets with unsupervised clustering. In the dataset GSE77509, which includes 19 pairs of PVTT and Ts, the 24-gene signature classified PVTT13 and T13 into two different groups and classified other PVTTs and Ts from the same patient into the same group (Figure 4A). In the dataset GSE69164, including 11 pairs of PVTTs and Ts, PVTT9 and T9 were classified into two different groups, whereas other PVTTs and Ts from an individual patient were classified into the same group (Figure 4B). In the dataset GSE74656, which includes 5 pairs of PVTTs and Ts, all the PVTTs and Ts from an individual patient were classified into the same group(Figure 4C). In addition, the ability of the 24-gene signature to identify clonal relatedness was also validated using newly collected samples (n = 8). PVTT and T from patient 2 were,however, classified into two different groups (Figure 4D). All these data indicated that the 24-gene signature may be used to identify PVTTs and Ts with no clonal relatedness.

Figure 3 Enrichment analysis of gene signatures in PVTTs and Ts. (A) Enrichment of gene signatures of bile duct cells (BDC), hepatocytes,hepatic progenitor cells (HPC), and hepatic stellate cells (HSC) in PVTTs and Ts by ssGSEA based on their gene expression profiles. (B)Expression of the genes included in gene signatures of different liver cells in PVTTs and Ts from the dataset GSE77509. (C) Enrichment of gene signatures of hematopoietic stem cells and cancer stem cells in PVTTs and Ts by ssGSEA. (D) Expression of the genes included in gene signatures of stem cells in PVTTs and Ts from the dataset GSE77509.

Discussion

Figure 4 Validation of identification ability of 24-gene signature cluster analysis of pairs of PVTTs and Ts from GSE77509 (A), GSE69164 (B),GSE74656 (C), and newly collected samples (D) based on expression of 24 genes assessed by RNA-seq, cDNA microarray, or real time PCR.Orange line indicates the two subgroups classified based on the 24-gene signature.

In the present study, we found that PVTTs may have a different clonal origin than that of the corresponding HCC Ts. These types of PVTTs have different gene expression profiles from those of the corresponding Ts, whereas PVTTs originating from primary Ts have similar gene expression profiles to those of their origin Ts tissues. Previous studies reported the involvement of several factors, such as hypoxia and non-coding RNA in PVTT development, and similar gene expression profiles of PVTTs with Ts7-10, suggesting that PVTT may be a special type of HCC intrahepatic metastasis.Here, the clonality analysis, which was established to identify paired Ts with clonal origin, was applied to investigate the clonal relationship between PVTTs and Ts. Based on the clonality analysis and 24-gene signature identification, we found that most of the PVTTs (18 out of 19 in GSE77509, 10 out of 11 in GSE69164, and all 5 of GSE74656) possibly originate from the paired Ts (Table 1, Figure 4).Additionally, we also found another type of PVTTs, which received a very low LR2 value in the clonality analysis and were classified into a different group with the paired Ts in the cluster analysis (Table 1, Figure 4). These types of PVTTs were considered to represent tumor tissues having a different clonal origin than that of the paired Ts. These findings suggest the existence of two types of PVTTs according to their clonal origin, wherein one type has the same origin as the paired Ts and the other does not. Moreover, our findings can explain why PVTTs were observed in some patients without liver parenchyma tumor nodules13. The results can also explain why PVTTs from different patients showed very different gene expression profiles11. Of note is the observation that different CNV patterns were observed between PVTT13 and T13. These CNVs may be caused by gene translocation or genomic instability in the tumor cells.Genes transfer onto other chromosomes or extrachromosomal DNA (ecDNA) is amplified in tumor cells25. Genomic instability is a well-known characteristic of tumor cells, inducing alterations in the genomes26. Hepatitis virus infection is another important factor that causes genomic instability in HCC27. Hence, the mechanism underlying the different CNV patterns between PVTTs and Ts needs more investigation.

Our study suggests that integrated analysis of the published datasets provides a better understanding of Ts, and that such analyses may lead to improved therapeutic strategies. With the development of high-throughput technology and multi-omics analyses, increasingly more omics data from Ts are uploaded into the public databases such as GEO and The Cancer Genome Atlas (TCGA)28. An integrated analysis of these datasets makes it possible to obtain a more comprehensive understanding of Ts at multidimensional levels. The comprehensive integrated analysis of HCC using six distinct data platforms determined subtypes of HCC associated with poorer prognosis, and the gene expression signatures correlated with poor survival and potential therapeutic targets for HCC29. Another multiplatform integrative analysis based on TCGA found that molecular signatures provide a higher accuracy for the clinical outcome prediction than that of the currently used tissue-of-origin-based classification in some cases30. By integrated analysis of CNV and gene expression profiles previously used to investigate aberrantly expressed lncRNAs11, we analyzed the origin of PVTTs and developed a gene signature to identify the type of PVTTs. It has been reported that the “cell-of-origin”-based features mainly dominate the molecular taxonomy of tumor types30.Therefore, the gene signatures developed in our study may provide a platform for a high-accuracy identification of PVTTs with independent origins. However, follow-up studies and additional samples are necessary to validate the findings reported here.

Acknowledgements

This work was supported by grants of China National Funds for Distinguished Young Scientists (Grant No. 81425019),National Natural Science Foundation of China (Grant No.81672899), the State Key Program of National Natural Science Foundation of China (Grant No. 81730076), Shanghai Science and Technology Committee Program (Grant No.18XD1405300) and Specially-Appointed Professor Fund of Shanghai (Grant No. GZ2015009).

Conflict of interest statement

No potential conflicts of interest are disclosed.