Comprehensive molecular and immunological characterization of hepatocellular carcinoma.

BACKGROUND: Hepatocellular carcinoma (HCC) is a heterogeneous disease with various etiological factors, and ranks as the second leading cause of cancer-related mortality worldwide due to multi-focal recurrence. We herein identified three major subtypes of HCC by performing integrative analysis of two omics data sets, and clarified that this classification was closely correlated with clinicopathological factors, immune profiles and recurrence patterns.
METHODS: In the test study, 183 tumor specimens surgically resected from HCC patients were collected for unsupervised clustering analysis of gene expression signatures and comparative analysis of gene mutations. These results were validated by using genome, methylome and transcriptome data of 373 HCC patients provided from the Cancer Genome Atlas Network. In addition, omics data were obtained from pairs of primary and recurrent HCC.
FINDINGS: Comprehensive molecular evaluation of HCC by multi-platform analysis defined three major subtypes: (1) mitogenic and stem cell-like tumors with chromosomal instability; (2) CTNNB1-mutated tumors displaying immune suppression; and (3) metabolic disease-associated tumors, which included an immunogenic subgroup characterized by macrophage infiltration and favorable prognosis. Although genomic and epigenomic analysis explicitly distinguished between HCC with intrahepatic metastasis (IM) and multi-centric HCC (MC), the phenotypic similarity between the primary and recurrent tumors was not correlated to the IM/MC origin, but to the classification.
INTERPRETATION: Identification of these HCC subtypes provides further insights into patient stratification as well as presents opportunities for therapeutic development. FUND: Ministry of Education, Culture, Sports, Science and Technology of Japan (16H02670 and 18K19575), Japan Agency for Medical Research and Development (JP15cm0106064, JP17cm0106518, JP18cm0106540 and JP18fk0210040).

Evidence before study

Caner stratification will be essential for improvement in prognosis following therapy. Several studies have challenged to categorize HCC by mutation, DNA methylation and expression profiles, but the links between the molecular and clinicopathological traits have not been fully unveiled.

Added value of this study

We have classified HCC into three subtypes characterized by distinctive molecular and clinical features, of mitogenic, CTNNB1-mutated and metabolic disease-associated tumors, through the integration of genome, epigenome and transcriptome analysis. Cell-type deconvolution of gene expression data further identified the second subtype as an immunosuppressive phenotype, and the immune-related subgroup with encouraging outcome in the third subtype. Comparative analysis between primary and recurrent tumors clarified that they resembled each other in subtype whether they shared cells-of-origin or not.

Implication of all the available evidence

This phenotypic classification of HCC has important implications for clinical practice and drug development of targeted therapies including immunotherapy, particularly in predicting post-operative recurrence of HCC.

Alt-text: Unlabelled Box

Introduction

Liver cancer is the sixth common cancer with an estimated global incidence of 782,000 cases in 2012, and ranks as the second leading cause of mortality in patients with malignancy worldwide [1]. Hepatocellular carcinoma (HCC) accounts for 90% of all primary liver cancers, and is a heterogeneous disease with a variety of etiological factors. Chronic infection with hepatitis B virus (HBV) and hepatitis C virus (HCV), alcohol consumption and metabolic syndrome which includes obesity and diabetes trigger liver injury, and progressive destruction and regeneration of liver contribute to inflammation, fibrosis and then carcinogenesis. Understanding HCC diversity to develop targeted therapies will require unraveling the molecular events underlying the process [2,3].

Genome-wide analysis of gene mutation [[4], [5], [6]], DNA methylation [7] and mRNA expression profiles [[8], [9], [10], [11]] has been devoted to this purpose over the past two decades. Schulze et al. identified three major clusters of associated alterations, centered on TP53, CTNNB1 and AXIN1, by analyzing the whole coding sequences of 243 liver tumors surgically resected [6]. Villanueva et al. quantified the array-based DNA methylation levels in 221 and 83 HCC samples from the training and validation sets, respectively, and clarified that methylation signature of 36 unique probes divided patient survival [7]. Gene expression patterns of HCC have been examined by several laboratories [[8], [9], [10], [11]], and mainly classified into two subclasses [[12], [13], [14]], “proliferative and non-proliferative” or “aggressive and less aggressive”. Thus, although clinical samples are stratified in each genome, methylome or transcriptome study, these omics data have not been integrated sufficiently, or the relationships between the molecular and clinicopathological features have not been elucidated completely.

A serious conundrum in HCC treatment is represented by frequent recurrence, which is categorized into two types, intrahepatic metastasis arising from a primary tumor (IM) and multi-centric occurrence (MC). Although discrimination and characterization of the recurrence types provide great benefits for determination of post-operative adjuvant therapy and care strategy, there is no consensus on accurate diagnosis from clinical or pathological traits. Whole-genome sequencing of pairs of primary and recurrent HCC clearly distinguished between the IM and MC types with high (15.53 to 93.02%) and low (0 to 0.28%) shared mutation call rates, respectively, but was unable to discover genes specifically mutated in the IM or MC types [15].

Hence, this study performed integrative analysis of two HCC datasets, and established a molecular classification that resolved into three distinct subtypes. The molecular subtypes of paired primary and recurrent tumors closely matched each other, although there was no significant correlation between the molecular classification and IM/MC origin designation.

Materials and methods

Human tissue samples

A total of 183 patients who underwent curative hepatic resection for HCC between 2006 and 2013 at Tokyo Medical and Dental University Hospital participated in integrative analysis of classification (Classification study). Fifteen adjacent liver tissues obtained from patients with metastasis of colorectal cancer who had not received chemotherapy were used as control. In addition, primary and recurrent HCC samples surgically resected from 18 patients at the hospital were utilized for integrative analysis of recurrence (Recurrence study). With Institutional Review Board approval, written informed consent was obtained from all patients (permission No. G2017-018). Patients were anonymously coded in accordance with ethical guidelines, as instructed by the Declaration of Helsinki.

DNA and RNA extraction

Genomic DNA and total RNA were extracted from tissue specimens and cell lines by using QIAamp DNA Mini Kit (QIAGEN, Tokyo, Japan) and RNeasy Protect Mini Kit (QIAGEN), respectively. Contaminating DNA was removed by digestion with RNase-Free DNase Set (QIAGEN). The integrity of the obtained DNA and RNA was confirmed by using 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA).

Microarray analysis

Microarray data were analyzed as follows in accordance with the manufacturer's instructions. Complementary RNA was prepared from 100 ng of total RNA from each sample with 3′ IVT Express Kit (Affymetrix, Santa Clara, CA). Hybridization and signal detection of the GeneChip Human Genome U133 Plus 2.0 Array (Affymetrix) were carried out.

Integrative analysis for the TMDU study

The microarray data set of 183 tumor tissues from HCC patients and 15 normal liver controls (GSE112790, GEO) was normalized by using the robust multiarray average method in R statistical software (version 3.0.3) and the affy Bioconductor package. Then, 938 genes were selected according to two criteria; genes significantly upregulated or downregulated between tumor tissues and normal liver tissues by using Mann-Whitney U test and the Benjamini-Hochberg procedure for calculating P-values and false discovery rates (FDRs), respectively (FDR < 0.05); genes differentially expressed by estimating interquartile ranges (IQRs) among tumor samples (IQR > 1.5). Consensus clustering analysis was performed by using the k-means algorithm in the ConsensusClusterPlus Bioconductor package.

Genome analysis

For the Classification and Recurrence studies, genomic DNA were obtained from 33 pairs of HCC and adjacent liver controls and from primary and recurrent tumors and adjacent liver tissues of 18 HCC patients, respectively (hum0041.v1, NBDC). The genomic DNA was fragmented with Covaris S220 (Covaris Inc., Woburn, MA), and DNA library was prepared by using the SureSelect XT Human All Exon V4 kit (Agilent Technologies, Santa Clara, CA) according to the manufacturer's protocol. After high throughput sequencing for each DNA library was performed by using HiSeq 2000 platform (Illumina, San Diego, CA) with 100-bp paired-end sequencing, potential PCR duplicates with the same sequences in both pairs of reads were removed, and then sequence reads were aligned to the reference genome (hg19) by using the Burrows-Wheeler Aligner (BWA, version 0.6.1) [16]. Genome Analysis Toolkit (GATK, version 1.5.30) was utilized for base calibration and local realignment by following the GATK best practices [17]. Single nucleotide variants and small insertions/deletions were identified by MuTect (version 1.0.27783) [18] and SomaticIndelDetector (version 1.5.30, included in the GATK software) with the default parameters, respectively. The somatic variants were annotated by ANNOVAR (version 2012/03/08) [19]. The whole-exome sequencing data (hum0041.v1) are deposited in National Bioscience Database Center (NBDC).

Methylation analysis

The microarray data set of eight human HCC cell lines (Hep3B, HepG2, HLE, HLF, Huh1, Huh7, JHH2 and JHH4) exposed and unexposed to 5 μM 5-azacitidine, a demethylating agent, for three days (GSE112788, GEO) was normalized by using the robust multiarray average method in R statistical software and the affy Bioconductor package. For the classification and recurrence studies, genomic DNA were obtained from 29 pairs of HCC and adjacent liver controls (GSE113017, GEO) and from primary and recurrent tumors and adjacent liver tissues of 18 HCC patients (GSE113019, GEO), respectively. Bisulfite conversion, amplification, fragmentation, hybridization, labeling and scanning were performed according to the manufacturer's protocols. DNA methylation profiles were analyzed by using Infinium HumanMethylation450 BeadChip (Illumina) according to the manufacturer's instruction. The methylation levels of CpG sites were reported as β-values by GenomeStudio software (Illumina). In the Classification study, since methylation levels are increased around the promoter regions but not necessarily contribute to gene silencing, 112 genes were extracted following four criteria; genes hypermethylated in HCC by evaluating the difference of median β-values of cancerous and non-cancerous specimens (βcancer - βnon-cancer > 0.1); genes differentially methylated between paired tumor tissues and adjacent liver tissues by using Wilcoxon signed-rank test (p < .01); genes differentially methylated among the molecular subtypes by using Kruskal-Wallis test (p < .01); genes upregulated in more than one HCC cell line among the eight HCC cell lines under the treatment of 5-azacitidine (fold change >1.5). In the Recurrence study, kernel density estimation for β-values of the same 112 genes in the 18 pairs of primary and recurrent HCC was plotted in the ks package of R.

Integrative analysis for the TCGA study

Public omics and clinical data of 373 HCC patients were provided from the Cancer Genome Atlas Network (TCGA) and National Cancer Institute, and downloaded from the cBioPortal and Genomic Data Commons (GDC) Data Portal sites, respectively. All the mRNA expression data (RNASeq Version 2 RSEM) were transformed into log 2 base, and 726 genes were matched to the 938 genes determined in the TMDU study. Significant focal copy number alterations were identified from segment data by using GISTIC 2 (version 2.0.23) [20]. DNA methylation profiles of 373 cancerous and 46 non-cancerous tissues by using the same 112 CpG probe set of Infinium HumanMethylation450 BeadChip defined in the TMDU study were collected from the GDC Data Portal. The tables of the TCGA sample ID and the subtype of the iCluster classification [21] and the immune classification [22] have been published in the papers.

Gene set enrichment analysis and aggregate score

The gene set enrichment analysis (GSEA) was performed with the MSigDB gene sets (H: hallmark gene sets, C2: chemical and genetic perturbations and C5: GO biological process; version 5.0) [23]. Aggregate scores for gene sets were calculated as previously explained [24]. Briefly, the expression of each gene in the pathway was transformed into percentiles and the activity of each pathway was calculated as the average percentile score of all genes in a pathway minus 50, that is, the expected median activity of a pathway. Gene expression signatures were determined by aggregate scores of representative gene sets to prevent a mosaic pattern of gene sets including similar components.

Immune analysis

Immune-related gene sets were previously listed [25]. In accordance with the manufacturer's instructions, CIBERSORT analysis of mRNA expression data of all the 373 HCC samples provided from TCGA was performed in the Absolute mode with LM22 signature matrix on the web (https://cibersort.stanford.edu/) [26].

Statistical analysis

Statistical analyses were performed by using R statistical software. Mann-Whitney U test was used to analyze for differences between values of two independent groups. Kruskal-Wallis test with Steel-Dwass post hoc test was also done for those of more than three groups. Hierarchical clustering analysis was performed by using Ward's method in the hclust package of R. Fisher's exact test and χ2 test were applied to analyze categorical variables. Survival curves were constructed by using the Kaplan-Meier method and compared by using the log-rank test with univariate and multivariate Cox proportional hazards analysis. p < .05 was considered statistically significant.

Results

Classification based on transcriptome analysis

To define molecular subgroups of HCC, we first extracted 938 genes upregulated or downregulated between tumor tissues and normal tissues (FDR < 0.05, Mann-Whitney U test and the Benjamini-Hochberg procedure) and differentially expressed among tumors (IQR > 1.5) by using microarray data from 183 HCC samples surgically excised in Tokyo Medical and Dental University (TMDU) Hospital (Supplementary Table S1a). We next performed unsupervised clustering analysis of them, and yielded two major groups termed as A and B (Fig. 1a) on the basis of consensus clustering analysis inferring that the optimal number of groups was two (Fig. 1a and Supplementary Fig. S1a). We then examined biological meaning ascribed to the two subgroups. The Gene Set Enrichment Analysis (GSEA) demonstrated that gene sets involved in mitogenic activity and cellular metabolism were enriched in Group A and B, respectively (Supplementary Fig. S1b and Supplementary Table S1b). Group A was remarkably characterized by the increased expression levels of genes playing essential roles in G2/M cell cycle phase (FoxM1, aurora kinases and polo-like kinases), and stem/progenitor markers including epithelial cell adhesion molecule (EPCAM), and CD133 (PROM1) cytokeratin 19 (CK19) as shown in Supplementary Fig. S1c. In contrast, genes specifically upregulated in Group B included GLUL encoding glutamine synthetase as well as FABP4 and SPARCL1 [27,28].

Fig. 1

Comprehensive molecular classification of HCC in the TMDU test study identifying three molecular subtypes. (a) Transcriptomic classification of HCC. Unsupervised hierarchical clustering analysis of gene expression identified two groups; Group A (red, n = 94) and B (black, n = 89). The panel of molecular features is a heatmap displaying the relative expression levels of genes specifically upregulated in Group A or B. (b) Kaplan-Meier analysis of patients stratified by group. (c) Genomic and transcriptomic classification of HCC. Genome analysis divided Group B into HCC samples with (MS2) or without (MS3) CTNNB1 mutations which were detected only in Group B. MS1 (red, n = 17), MS2 (green, n = 6) and MS3 (blue, n = 10). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

We then compared the clinical characteristics and observed significant differences of not liver conditions but etiological features between the HCC samples (Supplementary Table S1c). Group A was mainly composed of carriers of hepatitis virus (p = .002), whereas Group B was associated with patients with obesity and diabetes (p = .045 and 0.012). Alcohol abuse had no relationship with this classification. The serum levels of tumor markers for HCC, α-fetoprotein (AFP) and protein induced by vitamin K absence or antagonists II (PIVKA-II), were significantly elevated in Group A (p < .001 and p = .012). Pathological findings indicated that Group A was composed of tumors with poor differentiation and portal vein (PV) invasion of HCC (both p < .001), resulting in advanced TNM and BCLC stages (p = .001 and p < .001). These molecular and histological traits conferred worse disease-free survival (DFS) and overall survival (OS) in the patients of Group A (Fig. 1b). Molecular classification of Group A and B was an independent prognostic factor of DFS (HR = 0.62, 95% CI = 0.43 to 0.89, p = .009; Supplementary Table S1d).

Integrative analysis of exome, methylome and transcriptome profiles

Since the large-scale transcriptome analysis identified the two distinct subgroups, we randomly selected 17 and 16 HCC samples from Group A and B, respectively, and further conducted integrative analysis. They were subjected to whole-exome sequencing, detecting 24,878 somatic mutations comprised of 23,050 single nucleotide variants and 1828 small insertions/deletions. Comparative analysis of the genomic data discovered nine genes mutated specifically in each group (Supplementary Table S1e). We highlighted somatic mutations of CTNNB1 observed only in Group B for two reasons; CTNNB1 ranked as one of the top genes differentially mutated between Group A and B; active mutation of CTNNB1 is frequently detected and well-known as a driver in HCC [[4], [5], [6]]. We re-categorized the 33 HCC samples into three molecular subtypes (MS) as below; the MS1 was equal to Group A, and the MS2 and MS3 were Group B with or without CTNNB1 mutations, respectively (Fig. 1c). The classification could not additionally provide any other information about genes frequently mutated in HCC such as TERT, TP53, chromatin remodelers, and mediators of the Wnt/β-catenin and RTK/RAS/PI(3)K pathways.

Infinium Human Methylation450 arrays were used for quantifying the methylation levels of the tumors and normal pairs, and gene expression profiles of eight HCC cell lines treated with and without 5-azacitidine were examined for exploring gene silenced in HCC. We extracted a set of 112 probes targeting CpG sites highly methylated in cancerous tissues relative to non-cancerous tissues (βcancer - βnon-cancer > 0.1; p < .01, Wilcoxon signed-rank test) and differentially modified among the three molecular subtypes (p < .01, Kruskal-Wallis test) considering the results from the comparative expression analysis of HCC cell lines (genes >1.5-fold upregulated under the treatment of 5-azacitidine; Supplementary Fig. S1d and Supplementary Table S1f). DNA methylation profiles with this probe set revealed that the MS2 was hypermethylated followed in order by the MS1 and MS3, possibly resulting from a protein interaction between β-catenin and a DNA methyltransferase DNMT1 [29].

We compared gene expression patterns of the three subtypes with each other by introducing the GSEA with the Hallmark gene sets of the MSigDB collections, and elucidated that gene sets implicated in mitosis, the Wnt/β-catenin pathway, and inflammatory response via the IL6/STAT3 and TNFα/NFκB pathways were specifically enriched in the MS1, MS2 and MS3, respectively. We defined G2M_CHECKPOINT, WNT_BETA_CATENIN_SIGNALING and INFLAMMATORY_RESPONSE as the representatives of gene expression signatures for each subtype, which was statistically validated by computing aggregate scores for determining the fittest in each specimen (p = 3.49 × 10−7). As illustrated in Fig. 1a, the MS1 exhibited distinct gene expression patterns of proliferative and stem/progenitor cell-like phenotypes. The expression levels of GLUL and SPARCL1 were selectively upregulated in the MS2 bearing active mutations of CTNNB1, since they are direct downstream genes of the Wnt/β-catenin signaling pathway [10,30].

Clinicopathological properties in the integrative study were determined from those in the transcriptome study (Supplementary Table S1c and S1g). It was worth noting that the MS1 was composed of patients with hepatitis virus infection (p < .001), and that the MS3, not the MS2, was tightly associated with metabolic risk factors of obesity, diabetes and dyslipidemia (p = .066, 0.024 and 0.060).

Validation study for molecular classification by using public omics data

We used public data of the Cancer Genome Atlas Network (TCGA) study [21] provided from the cBioPortal and Genomic Data Commons (GDC) Data Portal to evaluate the accuracy of our classification and further characterize each subtype. We first performed hierarchical clustering analysis of mRNA expression data from 373 HCC samples in the TCGA study, similarly to the TMDU study. The 938 genes differentially expressed between Group A and B (Supplementary Table S1a) could also divide them into two major clusters, and the left cluster was enriched with genes contributing to mitogenic and stem cell-like activity, which indicated that the left and right clusters might correspond to Group A and B, respectively. Whole-genome sequencing uncovered 58 tumors with CTNNB1 mutations in Group B, especially in the left node of Group B. Considering the molecular effects of CTNNB1 active mutations on gene expression profiles, we placed the three nodes to the MS1, MS2 and MS3 subtypes respectively as shown in Fig. 2.

Fig. 2

Comprehensive molecular classification of HCC in the TCGA validation study. Hierarchical clustering analysis with the gene set used in Fig. 1 could separate 373 HCCs into the MS1 (red, n = 114), MS2 (green, n = 74) and MS3 (blue, n = 185) as defined in Fig. 1c. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Genomic analysis discovered 49 genes differentially mutated among the three subtypes (Supplementary Table S2a), and that TP53 and CTNNB1 mutations were substantially augmented in the MS1 and MS2, respectively (p = 1.18 × 10−4 and 5.53 × 10−27). Somatic copy number alteration (CNA) analysis by using the GISTIC identified 40 amplifications and 34 deletions (Supplementary Fig. S2a), and clarified that the MS1 harbored extensive copy number aberrations, which could imply chromosomal instability, compared with the MS2 and MS3 subtypes. While focal amplifications of TERT and focal deletions of CDKN2A were commonly detected in the three subtypes, focal deletions of tumor suppressor genes MACROD2 and ACVR2A were characteristic of the MS2 and MS3, respectively (Supplementary Table S2b). Methylome analysis of the 373 tumor tissues and 46 adjacent liver tissues by using the same probe set of the Human Methylation450 array (Supplementary Table S1f) distinguished cancerous and non-cancerous tissues from each other, and revealed hypermethylation phenotype of the MS2 samples (Supplementary Fig. S2b). An aggregate score for each gene expression signature defined in the TMDU test study was calculated in each cases, and confirmed the correspondence relationship between the signatures and subtypes (Supplementary Table S2c; p = 1.14 × 10−25). The expression patterns of the nine genes selected for representing each subtype in the TMDU test study resembled those in the TCGA validation study.

Clinicopathological features in the TCGA data set were highly consistent with those in the TMDU data set (Supplementary Table S2d); the MS1 was associated with increased serum AFP level, vascular invasion scores and TNM stages (all, p < .001). The MS3 was composed of overweight individuals (p = .003). Although patients infected with hepatitis virus were accumulated in the MS1 in the TMDU study, this tendency was not clearly observed in the TCGA study (p = .042) due to racial and ethnic gaps in the enrolled population. In DFS and OS, Group A showed worse prognosis compared with Group B, and the MS1 more than the MS3, but not more than the MS2 (Supplementary Fig. S2c and S2d), which was confirmed as an independent predictor of patient outcome in HCC by univariate and multivariate Cox proportional hazards model analysis (Supplementary Table S2e).

Comparative analysis of molecular classifications

Several groups have proposed HCC molecular classifications mainly based on transcriptome data during the past two decades, and Lee, Hoshida, Boyault and Chiang's classifications [[8], [9], [10], [11]] are widely accepted. We then compared our classification based on multi-platform data from the TMDU and TCGA studies (Fig. 3a, Supplementary Fig. S3 and Supplementary Table S3a). In the TMDU test study, Group A was linked to the subclasses of SURVIVAL_DN in Lee's classification, S1 and S2 in Hoshida's classification, G1/2 and G3 in Boyault's classification, and PROLIFERATION in Chiang's classification. All the subclasses have been well-documented as the “proliferative” and “aggressive” phenotypes characterized by high serum AFP levels and poor patient prognosis [[12], [13], [14]]. As predicted, these malignant subclasses were also enriched in the MS1 of the TMDU test study and the TCGA validation study. However, as previous investigations have described that the “non-proliferative” and “less aggressive” phenotypes are heterogeneous in nature, only the MS2, the G5/6 in Boyault's classification, and the CTNNB1 in Chiang's classification revealed patients with CTNNB1 somatic mutations in Group B as shown in Supplementary Table S3a (p = 5.53 × 10−25, 5.28 × 10−21, and 7.82 × 10−36). Conversely, the MS3, the SURVIVAL_UP in Lee's classification, and the S3 in Hoshida's classification were associated with obesity, suggesting that only our subtyping could distinctively separate the “non-proliferative” and “less aggressive” phenotypes into the MS2 (CTNNB1-mutated subtype) and MS3 (metabolic disease-related subtype). We finally compared our molecular classification with the iCluster classification proposed by the Cancer Genome Atlas Network [21], and observed that the iCluster 1, 2 and 3 were corresponded to the MS1, 3 and 2, respectively (Supplementary Table S3a). Our molecular classification showed much stronger correlation with the four other classification than the iCluster classification. Thus, our classification may have important advantages for determining the relationship of the molecular characteristics to clinicopathological features such as viral infection, obesity and vascular invasion (Supplementary Table S3b).

Fig. 3

Summary of molecular classification of HCC. (a) Comparison of aggregate scores with gene sets associated with the previously defined molecular classifications of HCC in the TMDU test study (upper) and TCGA validation study (lower). (b) Schematic representation of molecular subtypes.

Taken together, the unique traits of each molecular subtype were summarized in Fig. 3b; the MS1 harboring TP53 mutations displayed extensive mitotic activity eliciting chromosomal instability and stem cell-like property, and has significant correlation with high serum AFP levels, aggressive vascular invasion, and unfavorable prognosis; CTNNB1 mutations activate the Wnt/β-catenin pathway in the MS2 and induced hypermethylation; the MS3 is composed of patients with metabolic risk factors, such as obesity and diabetes, and inflammatory response may be enhanced via the IL6/STAT3 and TNFα/NFκB pathways.

Classification based on immunological expression profiles

Given the potential effects of infiltrating immune cells on carcinogenesis and the significant improvement of therapies targeting immune checkpoint pathways, immune phenotyping of cancer has recently been addressed by several studies, which is accompanied by the development of computational approaches that allow virtual dissection of tumors [21,25,31]. We first calculated aggregate scores of each tissue specimen surveyed in the TMDU test study for immune-related gene sets, and performed hierarchical clustering analysis (Fig. 4a and Supplementary Fig. S4a). Among the distinct expression patterns, one cluster was attributed to the presence of inflammatory conditions and immune cells, and was designated as the Immune class. The Immune class accounted for 42% of all the cases, and contained 90% of the MS3 samples (p = 4.22 × 10−4), consistent with the activation of inflammatory response pathways in the MS3. In the TCGA validation study, the cluster corresponding to the Immune class was also identified (38%), and marked by the enrichment of the MS3 cases (46%, p = 5.61 × 10−3). We then divided each of the three molecular subtypes into the Immune and Non-immune classes for the evaluation of individual outcomes, and DFS in the Immune class of the MS3 was better than that in the Non-immune class of the MS3 (p = .003), whereas there were no significant differences in survival time between the other Immune and Non-immune classes (Supplementary Fig. S4b and Supplementary Table S4a). Considering this discrepancy between the Immune and Non-immune classes of the MS3 (Supplementary Fig. S4c), we re-defined them as the MS3i and MS3n, respectively. The MS3i had longer DFS than the other three subtype (Fig. 4b), and showed weak positive correlation with alcohol consumption (Supplementary Table S4b and S4c).

Fig. 4

Immunological evaluation of primary HCC. (a) Hierarchical clustering analysis of aggregate scores with immune-related gene sets in the TCGA validation study. (b) Kaplan-Meier analysis of patients stratified by molecular subtype and immune class. (c) Cumulative CIBERSORT score for various types of immune cells in each sample. Horizontal lines show the median values. (d) CIBERSORT score for immune cells in each molecular subtype. Boxes in violin plots represent the interquartile range (range from the 25th to the 75th percentile), and horizontal lines show the median values.

Next, we investigated correlation between our molecular classification and the immune molecular subtypes [22], and observed that the C1/C2 (wound healing/IFNγ dominant), C3 (inflammatory) and C4 (lymphocyte depleted) were closely corresponded to the MS1, MS3i and MS2/MS3n, respectively (p = 1.28 × 10−25; Supplementary Fig. S4d and Supplementary Table S4d). These findings were intriguing in two ways; the C1 had mitogenic activity and chromosomal instability, similarly to the MS1; the C3 showed better prognosis, which was characteristic of the MS3i subclass.

To further explore which types of immune cells resided in each subtype, we utilized the CIBERSORT for the deconvolution of tumor samples, and revealed much higher cumulative scores of immune cells in the MS3i (Fig. 4c). Activated M1 and M2 macrophages infiltrated specifically in tumor tissues of the MS3i, and CD8+ T cells were significantly depleted in the MS2 and MS3n (Fig. 4d), which was accurately predicted by the immune molecular classification described above. Notably, such features in cumulative scores of immune cells were also observed when the HCC samples were simply classified into the three molecular subtypes (Supplementary Fig. S4b, S4e and S4f).

Inheritance of the phenotypic subtypes during tumor recurrence

To investigate phenotypic alterations during tumor recurrence, whole-exome sequencing was performed for 18 primary and 19 recurrent tumors surgically resected from 18 individuals. The mutation patterns of single nucleotide variants and frameshifts discriminated them into two types on the basis of the number of common mutations (Supplementary Table S5); one shared much more mutations than the other, indicating genetically determined HCC with intrahepatic metastasis (gIM) and multi-centric HCC (gMC), respectively (Fig. 5a). Unfortunately, gIM- or gMC-specific mutations were not detected in this study [15], and frequently mutated genes such as TERT, TP53 and CTNNB1 were not enriched in the gIM or gMC. We next compared DNA methylation statuses of the 112 probes listed in Supplementary Table S1f (Fig. 5b), and identified that the gIM pairs had higher correlation coefficient of β-values than the gMC (p = 2.26 × 10−7). These results supported the hypothesis that tumors derived from the same origin could inherit the DNA methylation patterns (R2 = 0.520), and also suggested that even multi-focal tumors developing in the same background could resemble each other with respect to epigenetic signature (R2 = 0.331) because of field cancerization [32]. Finally, we categorized the primary and recurrent HCC into the three molecular subtypes according to the optimized criteria which were established from the data set of the TMDU test study by using logistic regression analysis (Supplementary Table S6). Surprisingly, the recurrent HCC exhibited the same phenotype of the primary HCC, but the gIM/gMC diagnosis was not correlated to the subtype transition (Fig. 5c). These findings implied that the molecular subtype of HCC could be controlled by extrinsic factors such as tumor microenvironment.

Fig. 5

Molecular evaluation of recurrent HCC. (a) Genomic landscape of the gIM (orange, n = 8) and gMC pairs (lime, n = 10). Genes commonly mutated in more than two pairs are shown. White, gray and black bars represent primary, recurrent and re-recurrent HCC, respectively. (b) Density plots of genes which were frequently methylated in HCC and differentially done among the molecular subtypes in the gIM and gMC pairs. (c) Subtype transition between primary and recurrent HCC. The criteria for the determination of subtypes are described in Supplementary Table S6. Upper (light) and lower (dark) bars shows the gIM and gMC, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Discussion

Molecular classifications with gene expression profiles have achieved clinical success in the prediction of patient outcome and the decision for cancer therapy, and advances in DNA sequencing and methylation array technologies have led to a precise categorization and a better understanding of carcinogenesis. In HCC, the four transcriptome-based classifications (Fig. 3a) are broadly accepted [[8], [9], [10], [11]], but lacks genomic and epigenomic analysis. Meta-analysis of these classifications suggests that HCC can be divided into two major groups [[12], [13], [14]], “proliferative and non-proliferative” or “progressive and less progressive”. The former is well-characterized by mitogenic and stem cell-like properties, HBV infection, malignant phenotype such as high AFP, poor differentiation and frequent vascular invasion, and then unfavorable prognosis, whereas the latter is heterogeneous and the clinical and pathological features are ambiguous.

Integrated with genomic, epigenomic and clinical data, the 938 differentially expressed genes which were generated from the TMDU dataset of 183 HCC samples (Supplementary Table S1a) discriminated the non-proliferative/less progressive group into the MS2 and MS3 subtypes, that is, subtypes harboring CTNNB1 mutations and complicated with metabolic disease, respectively. We reached similar conclusion by using 1157 differentially expressed genes between 123 pairs of HCC tissues and adjacent normal liver tissues of the TMDU test study (Supplementary Fig. S5). In contrast, 881 differentially expressed genes from the TCGA data set could only separate the subtype characterized by active mutation of CTNNB1 from the others (Supplementary Fig. S6). In the MS2 tumors, the Wnt/β-catenin signaling pathway was aberrantly activated, which resulted in the elevated expression levels of downstream genes like GLUL and SPARCL1 [10,30]. Consistent with previous findings in methylome analysis of HCC [25], DNA methylation was highly enriched at the CpG sites in CTNNB1-mutated HCC, presumably because constitutively active β-catenin protein recruits a DNA methyltransferase DNMT1 to the binding sites [29]. The MS3 tumors were usually developed in individuals with metabolic syndrome including obesity and diabetes, and could be derived from non-alcoholic steatohepatitis (NASH). Gene sets involved in inflammation were remarkably augmented in the MS3, and these findings are similar to the results previously described. Lipid accumulation stimulates inflammatory pathways through endoplasmic reticulum dysfunction and reactive oxygen species production [33]. Combinatory analysis of genomics, transcriptomics and metabolomics identifies the C2 subtype linked to excess body weight and severe inflammation among 199 Asian liver cancer [31].

The MS2 harboring CTNNB1 mutations demonstrated immunosuppressive phenotype, which was confirmed by the CIBERSORT analysis as well as the immune classification. Three recent papers have reported the close correlation between functional mutation of β-catenin and potential exhaustion of immune cells including T cells [22,25,34] (Spranger et al., 2015; Sia et al., 2017; Thorsson et al., 2018). Regulatory T cells (Treg), suppressing immune responses, were rather decreased in the MS2 subtype than in the stem cell-like MS1. Although Spranger et al. described that the Wnt/β-catenin signaling pathway upregulated cytokine secretion which inhibited tumor immune microenvironment (TIME) in melanoma cells, further investigation is necessary for understanding the molecular mechanism in HCC.

Chronic inflammation induced by viral hepatitis and NASH promotes hepatocarcinogenesis. Immunological expression profiles [25] elucidated the enrichment of immunogenic HCC in the MS3 subtype, and the immune subclass of the MS3 (MS3i) had two hallmarks, better patient prognosis and increased immune cell infiltration, named as “hot”. As previously reported, patients with “hot” tumors of the Immune class [25] and the C3 subtype [22] exhibited better survival. We most recently addressed that fatty acid binding protein 4 (FABP4) was overexpressed in intratumoral hepatic stellate cells of metabolic disease-associated HCC, and that FABP4 contributed to the secretion of inflammatory cytokines through NF-κB nuclear translocation, resulting in the recruitment of macrophages [28]. M2 Macrophages are attractive targets for cancer treatment [35], and pexidartinib, an inhibitor of the CSF1/CSF1-R signaling which plays essential roles in proliferation and differentiation of M2 macrophages, may improve the prognosis of patients with the MS3 subtype. Shalapour et al. have recently reported that the PD-L1/PD-1 signaling contributes to liver carcinogenesis in individuals with NASH [36]. These results are consistent with the immune profiles of the MS3 subtype (Fig. 4 and Supplementary Fig. S4), and then imply that PD-1 blockade may inhibit the MS3 HCC growth and development. Since anti-CSF1R and anti-PD-1 antibodies can synergistically attenuate BRAF-driven melanoma cells in mice [37], the combination is also a promising therapy for the MS3 subtype.

It is an important issue which of intrinsic or extrinsic signaling pathways mainly cause cancer initiation and progression. If intrinsic factors are dominant, the IM pairs of primary and recurrent HCCs, which are originated from the same parental cells and share gene mutations and epigenetic patterns, could resemble each other in the phenotype, but not the MC pairs. Our molecular classification however revealed that the subtypes, not the IM/MC classification, influenced the phenotypic similarity between the paired primary and recurrent tumors, supporting the superior role of extrinsic factors such as inflammation and vascularization in hepatocarcinogenesis. This observation provides evidence for the use of immunotherapy. This idea is particularly attractive in the MS3 due to the underlying systemic diseases such as overweight and hyperglycemia. Therefore, optimal weight and glycemic management could be useful for prevention of recurrence in patients with the MS3 tumors.

This study presented three critical findings in molecular classification of HCC. First, combinational gene expression profiles of the mitosis and immune axes identified unique subtypes. The gene set defined in the TMDU test study was closely connected to mitotic process and clearly classified HCC samples into three molecular subtypes, namely the proliferative subtype (MS1) and the non-proliferative subtypes with and without CTNNB1 mutation (MS2 and MS3). Immune signature discovered significant accumulation of HCC with enhanced inflammatory response in the MS3, and further divided this subtype into immunogenic and non-immunogenic subclasses (MS3i and MS3n), resulting in favorable prognosis of the MS3i. Second, immunological characteristics were essentially different among the molecular subtypes. Immune cell-type deconvolution, recently emerging and providing profound insights into TIME, elucidated that the MS2 and MS3i exhibited immunosuppressive and immunogenic traits (decreased CD8+ T cells and increased macrophages, respectively). There is a limitation that molecular mechanisms causing the immunological difference remain obscure. The immunological classification could guide immunotherapy for HCC, such as 4-1BB agonist antibody for the immunosuppressive subtype and PD-1 blocking antibody for the immunoactive subtype [38]. Third, a cluster of metabolic disease-associated HCC was determined. We clarified that the MS3, the non-proliferative subtype without CTNNB1 mutation, was intimately linked with metabolic risk factors such as diabetes and obesity in the TMDU test study. These observations are important since HCC in patients suffering from metabolic syndrome, not infected with hepatitis virus, is gradually rising worldwide in the past decade, and it would be difficult to formulate the concept of this subtype were it not for the rich clinical information of the TMDU test study. Taken together, these subtypes demonstrated distinct features including activated signaling pathways, infiltrating immune cells and clinicopathological factors, and the classification system could therefore be used to design more intelligent clinical trials in the future so as to develop strategies in a subtype specific manner [2,3].

Funding sources

This work was supported by Grant-in-Aid for Scientific Research A (16H02670) and Challenging Research Exploratory (18K19575) from the Ministry of Education, Culture, Sports, Science and Technology of Japan; Research Grant from the Princess Takamatsu Cancer Research Fund; P-DIRECT (JP15cm0106064), P-CREATE (JP17cm0106518 and JP18cm0106540) and Division of Infectious Disease Research, Research Program on Hepatitis (JP18fk0210040) from AMED (Japan Agency for Medical Research and Development).

Conflicts of interest

The authors declare no competing interests.

Author contributions

S.T. conceptualized, designed and supervised the study. F.T., S.W., T. O., K.O., H.O., Y.M., D.B., A.K., S.A., M.T., and S.T. prepared surgical tissue samples of HCC. S.S., Y.A., F.T., and S.W. obtained data. S.S. and K.M. analyzed data. S.S. and K.M. wrote the first draft of manuscript. J. R. W. and S.T. consolidated the data and refined the manuscript.