Abstract
Understanding the intratumoral heterogeneity of hepatocellular carcinoma is instructive for developing personalized therapy and identifying molecular biomarkers. Here we applied whole-exome sequencing to 69 samples from 11 patients to resolve the genetic architecture of subclonal diversification. Spatial genomic diversity was found in all 11 hepatocellular carcinoma cases, with 29% of driver mutations being heterogeneous, including TERT, ARID1A, NOTCH2, and STAG2. Similar with other cancer types, TP53 mutations were always shared between all tumor regions, that is, located on the “trunk” of the evolutionary tree. In addition, we found that variants within several drug targets such as KIT, SYK, and PIK3CA were mutated in a fully clonal manner, indicating their therapeutic potentials for hepatocellular carcinoma. Temporal dissection of mutational signatures suggested that mutagenic processes associated with exposure to aristolochic acid and aflatoxin might play a more important role in early, as opposed to late, stages of hepatocellular carcinoma development. Moreover, we observed extensive intratumoral epigenetic heterogeneity in hepatocellular carcinoma based on multiple independent analytical methods and showed that intratumoral methylation heterogeneity might play important roles in the biology of hepatocellular carcinoma cells. Our results also demonstrated prominent heterogeneity of intratumoral methylation even in a stable hepatocellular carcinoma genome. Together, these findings highlight widespread intratumoral heterogeneity at both the genomic and epigenomic levels in hepatocellular carcinoma and provide an important molecular foundation for better understanding the pathogenesis of this malignancy. Cancer Res; 77(9); 2255–65. ©2017 AACR.
Introduction
As the second leading cause of cancer mortality, the incidence of liver cancer is increasing in almost all countries, representing a major health problem, particularly in Asia (1). Hepatocellular carcinomas contribute to over 90% of all liver cancers, with well-studied risk factors including infection of hepatitis B virus (HBV) and hepatitis C virus (HCV), alcohol consumption, aflatoxin exposure, and metabolic abnormalities, among others (2, 3). Previous Hepatocellular carcinoma genomic analysis using single sampling per tumor have identified a number of driver mutations in this malignancy (4–11); however, much remains unclear regarding the extent of genomic diversity and clonal evolution of primary hepatocellular carcinomas. In addition, despite the few exogenous mutational processes (e.g., aflatoxin ingestion) associated with particular mutational signatures (6, 9, 12), the temporal dynamics of these processes and their contribution to subclonal diversification during hepatocellular carcinomas evolution remain unclear. Moreover, the intratumoral epigenetic heterogeneity of hepatocellular carcinomas is unexplored and its biological significance remains unknown. Recently, our group addressed these important questions in esophageal tumors and observed profound intratumoral heterogeneity at both genetic and epigenetic levels (13).
Precise understanding of both the genomic and epigenomic architecture of primary hepatocellular carcinomas tumors is crucial for developing personalizing patient treatment and molecular-based biomarkers (14). In this study, we address these critical issues through integrative molecular approaches, including multiregion whole-exome sequencing (M-WES), phylogenetic tree construction, and multiregion methylation profiling. Furthermore, we analyzed hepatocellular carcinomas genomic data from The Cancer Genome Atlas (TCGA) to dissect temporally the mutational processes and clonality of driver mutations.
Materials and Methods
Patients and sample processing
All of the samples including tumor tissues (n = 52), morphologically nonmalignant liver tissues (n = 6), and matched peripheral blood (germline, n = 11) from 11 hepatocellular carcinoma patients were collected after diagnosis from Sun Yat-Sen Memorial Hospital during the years of 2015 to 2016 (Supplementary Table S1). Samples from the case HCC8031 were collected after transarterial chemoembolization; all other samples were obtained before treatment. All patients presented evidence of HBV infection except HCC8392, and no patient reported a history of alcohol consumption. The median patient age was 48 years (range 32–71). Tumor stages ranged from II to IIIC. Detailed clinical and pathological characteristics of these hepatocellular carcinoma cases are provided in Supplementary Table S2. This study has been approved by the Ethics Committee at Sun Yat-Sen Memorial Hospital. To study intratumoral heterogeneity, five spatially-isolated tumor specimens were obtained per patient (except for HCC8392, HCC8716, and HCC5647, where only four different tumor regions from each individual were available). Each specimen was at least 3 to 4 cm away from the others. In six out of 11 patients, we also collected morphologically nonmalignant liver tissues, which were at least 3 to 4 cm away from the nearest tumor region. We carefully reviewed available hematoxylin and eosin (H&E) slides of each tumor region to verify that the tumor cell content of the selected regions were comparable and were at least greater than 70%. All of the tissue samples were snap frozen in liquid nitrogen and preserved at −80°C before analysis.
Whole exome sequencing
DNA was extracted using Qiagen DNeasy Blood & Tissue Kit. Whole exome capture was performed on 2 to 3 μg of genomic DNA with the Agilent SureSelect Human All Exon v4 (51 Mb) Kit according to the manufacturer's instructions, and captured nucleotides were subjected to massively parallel sequencing using Illumina HiSeq4000 platform at Beijing Genome Institute (BGI), China.
Construction of phylogenetic tree
To increase the sensitivity of detecting mutations that had been called from one tumor section in other tumor sections from the same case at a low VAF, we adopted the method published by Stachler and colleagues (15). Specifically, read counts for all somatic variants across all tumor sections were extracted using Bam-Readcount (URL). A variant with VAF less than 0.02 were considered as absent. A binary table was then generated across all tumor regions for each hepatocellular carcinomas case. The binary tables were used to construct phylogenetic trees, using the method of Discrete Character Parsimony, implemented in PHYLogeny Inference Package (URL). Germline sample was set as an outgroup root. Inferred trees were manually recolored, with branch/trunk lengths proportional to the number of mutations.
DNA methylation analysis
DNA methylation profiles of 22 tumor regions and four morphologically nonmalignant liver tissues (three out of four were matched nonmalignant livers) from five hepatocellular carcinoma individuals were performed using the Illumina Infinium HumanMethylation450K (HM450) platform (Illumina) at the University of Southern California Norris Comprehensive Cancer Center Genomics Core. The identical DNA extracted from these 22 tumor regions were also subjected to M-WES as described earlier. For data analysis, preprocessing steps including background correction and normalization to get corrected β values were performed using the TCGA pipeline [methylumi R package as described in Triche and colleagues (16)]. Probes with a detection P-value over 0.01 were removed, as were probes overlapping with dbSNP SNPs, and probes on the X chromosome. We defined low methylated probes in nonmalignant tissues as those for which β values of all nonmalignant samples were less than or equal to 0.3; highly methylated probes in nonmalignant tissues were defined as those for which all nonmalignant samples were larger than 0.6. For intratumoral analysis, computational methods were performed as we have described recently (13). Specifically, we defined the variable probes as those where the difference in β values of at least 1 pair of tumor regions was at least 0.3, and invariable probes as those where the differences in β values of all the pairs were less than 0.1.
Accession codes
Digital sequencing and HM450 bead array files have been deposited into Sequence Read Archive (SRP076602) and Gene Expression Omnibus (GSE83691), respectively. Other Materials and Methods are described in the Supplementary Information.
Results
Hepatocellular carcinoma genomic architecture and evolution
To investigate hepatocellular carcinoma genomic evolution, we first performed M-WES on a total of 52 tumor regions and 11 matched germline DNA (peripheral blood), collected from 11 hepatocellular carcinoma cases (Supplementary Tables S1 and S2). In six out of these 11 patients, we also sequenced whole exome of the DNA from morphologically nonmalignant liver tissues, which were 3 to 4 cm away from cancerous regions, to investigate whether cancer “field effect” exists in hepatocellular carcinoma. Field effect refers to the presence of clonal expansions in nonmalignant tissues that provide a background for malignancy development, which was observed by us and others in the tumors from colon (17), prostate (18), and head and neck (19).
M-WES achieved an average target depth of 158X (range, 72X–232X; Supplementary Table S1), and identified an average of 129 somatic variants per tumor region (Supplementary Tables S3–S5), a mutational rate comparable to recent sequencing reports on hepatocellular carcinoma (4–6, 9, 20). Notably, the case-wise mutational rate (159 mutations/case) was significantly higher (P < 0.001; Supplementary Table S3), underscoring the stronger power and improved resolution of our multibiopsy approach for genomic analysis. To evaluate intratumoral heterogeneity and to illustrate hepatocellular carcinoma genome evolution, we constructed phylogenetic trees with the trunk representing mutations present in all tumor regions, branches representing mutations present in some but not all regions. Spatial intratumoral heterogeneity was evident in all 11 hepatocellular carcinoma cases, showing a median of 38.9% of variants being heterogeneous (range, 5.5–91.7%; Fig. 1, Supplementary Fig. S1, Supplementary Table S3). The phylogenetic tree structure varied considerably between hepatocellular carcinoma cases (Fig. 1). For example, HCC8392 had a longer branch than its trunk, whereas HCC8716 displayed a very homogenous mutational pattern. Importantly, we validated the intratumoral heterogeneity at the protein level by immunohistochemistry staining of the tumor regions that had available matched formalin-fixed paraffin-embedded (FFPE) samples. In HCC5647, FAT4 mutation was only observed in T4 region, and correspondingly FAT4 protein expressed at markedly lower level compared with other tumor regions (Fig. 2A). In contrast, in HCC6690, which harbored a truncal TP53 mutation, we observed strong staining of nuclear p53 protein in all five tumor regions (Fig. 2A).
Phylogenetic trees of 11 hepatocellular carcinomas constructed on the basis of M-WES. Phylogenetic trees were constructed from all somatic variants by Wagner parsimony method using PHYLIP (see Materials and Method). Blue, green, and orange lines represent trunk, shared branch, and private branch, respectively. Lengths of trunk and branch are proportional to their number of mutations. Putative driver events are mapped along the trees as indicated. Heat maps indicate the presence (red) or absence (gray) of a mutation in each tumor region (T) or matched nonmalignant liver tissue (N). For case HCC5647, trunk length was reduced for display purpose.
Spatial heterogeneity and CCF of putative driver mutations in hepatocellular carcinoma. A, IHC staining of FAT4 and p53 in different tumor regions that have been profiled by M-WES. Mutational status is indicated on top of each region. Scale bars, 100 μm. WT, wild type. B, Heat map of the CCF of putative driver mutations. Numeric number in each square shows the CCF. Columns, tumor regions; rows, genes.
Spatial heterogeneity and clonal status of driver mutations
Deciphering intratumoral heterogeneity of driver mutations is instructive for both better understanding of the cancer genome evolution and improving precision medicine strategies. We, therefore, identified potential driver mutations in this hepatocellular carcinoma cohort by comparing against COSMIC gene census (21), Pan-cancer analysis (22), and recent large-scale hepatocellular carcinoma sequencing results (4–7, 20, 23) on the basis of a set of criteria evaluating the functional impact of the variant (see Materials and Methods). These candidate driver events were then mapped on the phylogenetic tree structures. Similarly as observed in breast (24) and esophageal adenocarcinoma (25), 29% (10/34) of putative driver mutations were present in the branches, including TERT, ARID1A, NOTCH2, STAG2, and so on. Consistent with previous reports in other cancer types, TP53 mutations were always truncal. These results suggest that single-biopsy genome analysis may not be sufficient to evaluate comprehensively hepatocellular carcinoma driver events.
We next calculated cancer cell fraction (CCF) of all mutations to assess their clonal status through integrative analysis of tumor cellularity, variant copy number, and variant allele frequency (VAF) as described previously (26, 27). Notably, several driver mutations were fully clonal within some tumor regions, yet could not be detected in other areas. For example, TERT promoter hotspot mutation, which has frequently been identified in hepatocellular carcinoma, was absent in T1 region in case HCC8031 but appeared to be clonally dominant in all of the rest tumor regions (Fig. 2B). Likewise, ARID1A mutation was assessed to be present in 100% tumor cells in T3 and T4 regions, yet was undetectable in T1 and T2 regions in HCC8257 (Fig. 2B). Overall, although a number of mutations targeting driver genes were estimated to be fully clonal in all tumor regions sequenced (e.g., TERT and CTNNB1 in HCC6046, TP53 and SYK in HCC6952), it is evident that many of the driver mutations can be subclonal, indicating that they occur as relatively late events during the development of hepatocellular carcinoma.
To corroborate further these findings, we calculated the relative timing and clonal status of driver mutations using TCGA hepatocellular carcinoma datasets (see Materials and Methods; ref. 23). We observed that 45.5% (87 out of 191) TCGA hepatocellular carcinoma patients had more than 10% subclonal mutations (Supplementary Fig. S2). Importantly, in line with our earlier results, we found that in TCGA hepatocellular carcinoma samples a number of driver genes [defined by MutSig (28) in TCGA cohort] harbored a significantly higher proportion of subclonal mutations compared with that of total driver variants, such as CTNNB1, PTEN NOTCH2, and PTPRB (Fig. 3B). Overall, compared with passenger genes, driver genes still had significantly more clonal mutations in the TCGA cohort, in accordance with the observations in other types of cancers (ref. 23 and Fig. 3A).
Clonal status of hepatocellular carcinoma driver mutations sequenced by TCGA. A, The proportion of clonal and subclonal mutations in both driver and passenger events. B, CCF analysis of mutations in representative driver genes from TCGA hepatocellular carcinoma dataset. Each line represents an individual mutation. Round dot, upper and lower end of each line represents CCF, upper and lower bound of confidence interval, respectively. Clonal and subclonal CCF are shown as dark blue and orange, respectively. Driver genes were identified using MutSig algorithm (28) on the basis of the TCGA data (FDR q < 0.1). P values were derived by hypergeometric tests comparing the frequency of subclonal mutations in each gene against that in all driver genes. In those genes without any subclonal mutations, P value was not calculated.
However, we did not find compelling evidence of a “field effect” in our cohort of hepatocellular carcinomas. Only an average of two nonsilent mutations per sample (range 0–5) occurred in nonmalignant liver tissues from these hepatocellular carcinoma cases; and none of them seemed to be functionally relevant to cancer biology (Fig. 1). CCF analysis of these variants also suggested that many of them exhibited lack of clonal expansions within nonmalignant liver tissue (Supplementary Table S5).
We next examined hepatocellular carcinoma intratumoral heterogeneity at the copy number level, and compared our somatic copy number alteration (SCNA) calls with genomic identification of significant targets in cancer (GISTIC) results from the TCGA hepatocellular carcinoma project. Consistent with TCGA and other previous hepatocellular carcinoma SCNA studies (4, 9), a number of important copy number changes were noted in the present cohort, including copy number gain of 1q21, 5p15 (encompassing TERT), 8q24 (harboring MYC), and 11q13 (harboring CCND1) as well as loss of 1p36, 10q26, and 14q32. Overall, SCNA profiles were more similar within the same individual than between different ones. Nevertheless, extensive copy number intratumoral heterogeneity was observed, such that the majority of recurrent SCNAs were spatially heterogeneous in one or more cases (Supplementary Fig. S3). For instance, gain of 8q24 (harboring MYC) was ubiquitous in some cases but occurred as a heterogeneous aberration in other cases. Gain of 5p15 (encompassing TERT) was detected to be often spatially heterogeneous. Interestingly, our deep sequencing results also found a branch hotspot mutation in TERT promoter region (Fig. 1). These data suggest that the hyperactivation of telomerase resulting from both genetic mutations and copy number gains are likely late events during hepatocellular carcinoma pathogenesis. We also noticed that two significant SCNAs, copy number gain of 11q13 (harboring CCND1) and loss of 11q14 were consistently ubiquitous across intratumoral regions, underscoring the importance of these SCNAs as initiating forces during hepatocellular carcinoma transformation. Together, these results demonstrate that similar to DNA single nucleotide mutations, SCNAs also exhibit significant intratumoral heterogeneity in hepatocellular carcinoma, congruent with the observations in other types of cancers (25, 29, 30).
Temporal dissecting mutational spectra and signatures
We next analyzed the mutational spectra of both truncal and branch mutations to determine the dynamics of the mutagenic processes operative in hepatocellular carcinoma. Although we only observed moderate changes in the six mutational classes between truncal and late branch mutations (Fig. 4A and Supplementary Fig. S4), conspicuous differences in the 96 trinucleotide mutational signatures were noted (Fig. 4B). We next calculated the contributions of individual mutational signatures to each tumors, and identified several dominant signatures in these tumors, including Signature 1 (associated with age), Signatures 22 (associated with aristolochic acid exposure) and 24 (associated with aflatoxin exposure), and Signatures 6 and 15 (associated with DNA mismatch repair), which were in agreement with previous results from other groups (6, 9, 12). The robustness and stability of the analysis was confirmed by extensive permutation tests (Supplementary Fig. S5). Interestingly, a number of tumors exhibited prominent decreases of the contribution of Signatures 22, 24, and 25 (unknown etiology) in the branch compared with truncal mutations, albeit without obtaining statistical significance due to the relatively small number of tumors analyzed. These results indicate that exogenous factors, such as exposure to aristolochic acid and aflatoxin, contribute significantly to the mutagenic processes in the early stage of hepatocellular carcinoma development, whereas other endogenous factors might play more important roles in shaping the mutational spectrum after the most recent common tumor ancestor is established.
Temporal dissection of hepatocellular carcinoma mutational spectra and signatures. A, Pairwise fraction analysis of truncal and branch variants on the basis of the six mutation classes. P values were derived by paired Student t test upon verification of normality and variance within each group, with those over 0.1 not shown. B, Mutational signatures of all truncal and branch variants was inferred by deconstructSigs. Signatures are displayed according to the 96-substitution classification defined by the substitution class and sequence context (12). C, Dot plots display the contributions of individual mutational signatures to individual cases, with each dot representing one case. Signatures 1–30 were based on the Wellcome Trust Sanger Institute Mutational Signature Framework (12).
Intratumoral heterogeneity of DNA methylation in hepatocellular carcinoma
As in other common cancers, epigenetic alterations contribute significantly to hepatocellular carcinoma pathogenesis; and DNA methylation has been found to alter the expression of critical hepatocellular carcinoma cancer genes (31–33). To investigate intratumoral heterogeneity at the epigenetic level, we profiled the global DNA methylation of a total of 22 tumor tissues and four nonmalignant liver tissues from five hepatocellular carcinoma patients (all of which had matched WES results) using Illumina Human Methylation 450k (HM450) Bead array. First, we identified variably-methylated CpG sites (see Materials and Methods) across different tumor regions within each individual and found intratumoral heterogeneity in all individuals (Supplementary Fig. S6), a phenomenon that was recently observed in glioma, prostate, and esophageal cancers (34–37). All of the tumor regions had numerous differences from the methylation pattern of matched nonmalignant liver tissues (Supplementary Fig. S6; in two cases where matched nonmalignant samples were unavailable, tumor samples were compared against nonmalignant tissues from other individuals). Importantly, although different tumor regions showed variable methylation patterns, nonmalignant liver tissues were very homogenous, as demonstrated by the tight clustering; this high degree of interindividual homogeneity of nonmalignant liver tissues is consistent with previous epigenetic studies in hepatocellular carcinoma (31–33), and was further confirmed by our analysis of TCGA nonmalignant liver samples (Supplementary Fig. S7).
To explore potential biological relevance of the intratumoral heterogeneity of DNA methylation in hepatocellular carcinoma, we examined how the differentially-methylated probes were distributed across different genomic contexts. To this end, we defined both hyper- and hypomethylated CpG sites in these five individuals compared with matched nonmalignant liver tissues, and further separated them into ones with variable methylation between tumor regions (Fig. 5A) and ones with invariable methylation (Supplementary Fig. S8). This classification revealed a wide range in the degree of intratumoral methylation heterogeneity as measured by the number of variable probes in these five hepatocellular carcinomas (Fig. 5A and Supplementary Table S6). For example, only 1.8% of total hypermethylated probes exhibited heterogeneity in HCC6952; in sharp contrast, the majority (63.4%) of hypermethylated CpG sites in HCC8257 were variable among different tumor regions (Supplementary Table S6). To rule out and mitigate the confounding effects from nontumor DNA contamination, we performed approaches described previously by us and others to infer and correct for noncancer cells (37–39), using genome-wide Infinium methylation profiles of various flow-sorted immune cell types (40) to adjust β values of each probe accordingly. Importantly, the extent of intratumoral methylation heterogeneity in these hepatocellular carcinomas was almost identical after this correction (Supplementary Table S6), suggesting that the intratumoral methylation heterogeneity was largely driven by the cancer cells themselves. Moreover, we validated the relative degree of methylation heterogeneity using a recently-described metric based on Shannon entropy (41), after carefully controlling for the number of regions per case. This independent analysis agreed strongly with the degree of methylation heterogeneity inferred from the number of variably methylated probes (Supplementary Fig. S9), and confirmed intratumoral methylation heterogeneity is a highly robust feature of these samples.
Epigenetic intratumoral heterogeneity in hepatocellular carcinoma. A, Unsupervised hierarchical clustering of intratumoral methylation profiling of five hepatocellular carcinoma cases. Rows of the heat maps denote the methylation levels of variably hypermethylated (top) or hypomethylated (bottom) CpG sites across different tumor regions. Columns represent samples. B, Enrichment plots showing the distribution of both variably and invariably hypermethylated (top) or hypomethylated (bottom) CpG sites across a variety of functional genomic domains. nCGI-Prom, non-CG island promoter; CGI-nProm, CG islands not in promoter regions; CGI-Prom, CG island promoters; nCGI-PMD, PMD excluding CGI probes. All P values of hypergeometric enrichment test (comparing the frequency of each variable and invariable probe set category to that of array background) are shown in the Supplementary Table S6. C, Dot plots displaying enriched GO biological processes for the genes associated with variably and invariably hypermethylated promoters, with each dot representing one individual case. Red line, P value of 0.01.
We next compared variable and invariable loci to different genomic contexts, including relevant regulatory features such as promoters, CpG Islands (CGIs), and enhancers (Fig. 5B). Invariably-methylated CpG sites showed known patterns of cancer-specific methylation changes (42), including hypermethylated probes that were strongly enriched in CGI promoter regions, and depleted in both partially methylated domains (PMD) and enhancer regions after removing CGIs. Conversely, invariably-hypomethylated probes were markedly depleted in CGI promoters whereas enriched in PMDs as well as enhancer regions (Fig. 5B and Supplementary Table S7). Strikingly, variably-methylated CpG sites strongly resembled the distribution pattern of their invariable counterpart (Fig. 5B and Supplementary Table S7) in all five cases. Given the rich body of evidence showing the direct contribution of cancer-specific methylation to tumorigenesis and progression (e.g., suppression of tumor suppressors through promoter hypermethylation, activation of oncogenes through long-range hypomethylation; refs. 43, 44), our results suggest that heterogeneity of intratumoral methylation is likely to play important roles in the biology of hepatocellular carcinoma cells. In support of this postulate, gene ontology (GO) functional analysis of the genes with variably-hypermethylated promoters showed that they were significantly enriched in a number of cancer-related processes, including cell proliferation, differentiation, death, migration, adhesion, and transcriptional regulation (Fig. 5C). As expected, invariably-hypermethylated promoters also exhibited stronger enrichment in most processes (Fig. 5C). Close examination of the methylation targeted genes showed that a number of critical cancer genes (e.g., APC, IRF6, RUNX3, CDKN2A) could be either invariably or variably hypermethylated in their promoter regions, indicating that their expression might be differentially suppressed in different tumor regions (Supplementary Tables S8 and S9).
Prognostic value of genomic and epigenomic intratumoral heterogeneity
Recent studies in other cancer types have found a high degree of genomic intratumoral heterogeneity was associated with worse clinical outcome (45). Although clinical outcome information was unavailable for the cases studied here, we attempted to gain insights into this question using the TCGA hepatocellular carcinoma dataset. TCGA cohort contains only single-region samples, and is thus not ideal for studying intratumoral heterogeneity. Nevertheless, a recent mathematical algorithm has been introduced to estimate genomic intratumoral diversity from bulk-tumor WES data, using the Mutant-Allele Tumor Heterogeneity (MATH) approach (46, 47).
We first determined the MATH score of each individual (see Materials and Methods), then performed Log-rank analysis to test the association between MATH scores and patients' survival probability. Hepatocellular carcinoma individuals with higher MATH scores (>median) had worse outcome compared with those with lower scores (<median), albeit the result did not achieve statistical significance (Supplementary Fig. S10B, P = 0.19). We reasoned that more extreme difference in MATH scores (corresponding to larger variation in the degree of genomic heterogeneity) might reflect more difference in cancer malignancy, and therefore we compared patients with the highest 10% MATH scores to those with the lowest 10%. Here, cases with low MATH values exhibited significantly better survival rate (Supplementary Fig. S10A, P = 0.05).
Inspired by this result, we defined a methylation intratumoral heterogeneity (mITH) score similar in form as the MATH metric but using beta values from the HM450 array. Briefly, the mITH measures the median absolute deviation of probe beta values (see Materials and Methods for details). We restricted our analysis to hypermethylated probes, because the process of hypermethylation (which occurs primarily at CpG island promoters) is a more constrained and well-understood process in cancer than hypomethylation. When TCGA cases were divided into two equal groups based on mITH score, we observed the same trend as for the MATH comparison—high-mITH cases had worse survival rates than low-mITH ones; like the MATH comparison, this trend did not reach statistical significance (Supplementary Fig. S10D, P = 0.17). However, when we did the analysis of the mITH extremes (10% highest mITH vs. 10% lowest mITH), the result was highly significant (Supplementary Fig. S10C, P = 0.001). Although we believe that multiregion sampling will ultimately be required to validate fully these inference-based results, they do suggest that intratumoral heterogeneity may have prognostic value in hepatocellular carcinoma, and reinforce the importance of fully capturing intratumoral heterogeneity using a multiregion approach.
Discussion
Hepatocellular carcinoma is a worldwide leading malignancy with surgical resection and liver transplantation being the current mainstay treatment options. Morphologic and immunophenotypic intratumoral heterogeneity of hepatocellular carcinoma have been observed previously (48–50). By focusing recently on TP53 and CTNNB1 mutations in hepatocellular carcinoma, Friemel and colleagues (50) also reported that those variants could be either clonal or subclonal, which is in line with our findings. In contrast, hepatocellular adenoma, which is another form of rare, benign tumor of the liver, was found to exhibit homogenous histology and morphology (51, 52). Accordingly, genotypic and phenotypic classifications of hepatocellular adenoma has been established, and begin to have implications for clinical management (53). However, unresolved intratumoral heterogeneity impedes the development of both molecular classifications and targeted therapies in hepatocellular carcinoma. Sorafenib, a small-molecular inhibitor of several tyrosine protein kinases, is the only targeted therapeutic drug available for advanced hepatocellular carcinoma patients, yet it only offers limited survival benefit (54). Therapeutic merit of targeting Wnt/β-catenin pathway has not been established for hepatocellular carcinoma treatment, which could partially be attributed to the heterogeneous alterations of the activation of this signaling pathway. Clearly, an urgent need exists to decipher hepatocellular carcinoma genomic architecture, clonal evolution as well as subclonal diversification.
This study represents the first comprehensive genomic and epigenomic investigations of hepatocellular carcinoma intratumoral heterogeneity. Notably, we found evidence of spatial diversity with respect to DNA mutations in all 11 hepatocellular carcinoma cases. In support of this finding, we noted that 45.5% (87/191) TCGA-sequenced hepatocellular carcinoma patient harbored at least 10% subclonal mutations. The extent of subclonal diversification of hepatocellular carcinoma was comparable with esophageal and lung adenocarcinoma (25, 29, 55), but lower than clear cell renal cell carcinoma (30).
We demonstrated that many driver mutations (10/34) were mapped in the branches, suggesting that they occurred after the founding clone was established and contributed to the growth advantage of subpopulations. Further CCF analysis of each variant in each tumor region showed an illusion of clonal dominance of a number of driver events, where they appeared fully clonal in some tumor regions but were undetectable in others. This drastic CCF variation was not caused by tumor cellularity as skewed clonality was not observed in any of hepatocellular carcinoma samples, and overall the cancer cell content was comparable among intratumoral regions (Supplementary Fig. S11). This result challenges the accuracy of the approaches that estimate the clonal status of driver events from a single biopsy for developing biomarkers and therapeutic targets. Although independently-acquired somatic variants in distinct subclones were observed in some other tumor types (29, 30), this phenomenon was not present in this hepatocellular carcinoma cohort.
Conceivably, therapeutic targeting of truncal drivers might provide more anti-neoplastic efficacy than targeting branch divers. Indeed, it has been shown that inhibiting subclonal drivers in multiple myeloma only achieved partial treatment efficacy at best (26). In some cases targeting subclonal drivers even resulted in growth acceleration of non-mutated subpopulations, underscoring the importance of analysis of intratumoral heterogeneity for treatment designing and stratification (26). Our findings that TP53 mutations are always truncal might have clinical relevance for hepatocellular carcinoma treatment, albeit p53-based cancer therapies have yet to achieve clinical success. In addition, variants in several potential drug targets, including c-KIT and SYK, also ubiquitously occurred in all tumor regions in this cohort. Analysis of TCGA data also revealed that most of KIT mutations in hepatocellular carcinoma were clonal (Fig. 3), suggesting its potential therapeutic merit for this malignancy. In addition to the implications for targeted therapy, our analysis of the single-biopsy WES data from the TCGA hepatocellular carcinoma cohort revealed the potential prognostic value of intratumoral genomic heterogeneity, which is in agreement with the findings in other cancer types suggesting that higher magnitude of intratumoral heterogeneity is associated with more aggressive malignancy (45).
Although M-WES is able to resolve more subclonal diversity than single-biopsy strategy, to determine the true extent of intratumoral heterogeneity in hepatocellular carcinoma is still difficult. In addition, the prevailing view that Darwinian selection drives cancer evolution was challenged by recent reports (56, 57). Future studies through multiregional whole-genome sequencing of more tumor regions with deeper coverage are required to decipher the molecular mechanisms associated with the clonal expansion of hepatocellular carcinoma (56, 58). In the present cohort, 10 of 11 individuals were HBV positive (Supplementary Table S2), and future work with virus free hepatocellular carcinomas will shed light on the intratumoral heterogeneity resulting from viral integrations. Moreover, longitudinal multiregional genomic studies will help unravel whether and how intratumoral heterogeneity contribute to drug resistance and distal metastasis of hepatocellular carcinoma cells.
Hepatocellular carcinoma methylomes have been characterized using a single-sampling approach; however, little is known regarding intratumoral diversity at the epigenetic level or its relationship to changes in the genome. Here, we have shown profound epigenetic intratumoral heterogeneity in hepatocellular carcinoma. To minimize and alleviate the confounding effects from nontumor DNA contamination, additional approaches were performed to account for the potential influence from immune cells, and similar results of intratumoral methylation heterogeneity were observed (Supplementary Table S6, see Materials and Methods). Nevertheless, one should keep in mind that while intratumoral methylation heterogeneity possibly reflects complex heterogeneous populations of cancer cells, it may partially result from “phenotypic plasticity,” where biological processes (e.g., epigenetic modifications) of the same cells can change depending on their own differentiation state as well as the microenvironment. Close inspection of heterogeneously altered methylation sites revealed that a number of known cancer-related genes, such as APC, IRF6, RUNX3, CDKN2A, were differentially methylated between different tumor regions. These changes may have an impact on the cellular biology between different subpopulations. Importantly, albeit the number of differentially-methylated probes varied substantially between different individuals, our enrichment analysis (Fig. 5B) and GO pathway annotation (Fig. 5C) strongly suggest that intratumoral methylation heterogeneity likely plays important roles in the biology of hepatocellular carcinoma. Interestingly, although a clear interplay and codependence was found between phylogenetic and phyloepigenetic evolution in our investigations on esophageal tumors (13), we did not observe such a relationship in hepatocellular carcinoma in the present study (data not shown). It is also worthwhile noting that even tumors with “intratumoral-quiet” genome could still exhibit profound epigenetic heterogeneity (as demonstrated by the cases HCC5647 and HCC8010), suggesting that complex intratumoral heterogeneity can occur on multiple layers. Finally, we inferred methylation intratumoral heterogeneity based on the TCGA hepatocellular carcinoma cohort using a new analytic metric, mITH, and found that a higher degree of methylation heterogeneity was associated with a worse clinical outcome in hepatocellular carcinoma. These results, along with GO functional classes associated with variably methylated genes, suggest that epigenetic diversity could have functional relevance for hepatocellular carcinoma biology.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Web Resources
GATK best practices, https://www.broadinstitute.org/gatk/guide/bp_step.php?p=1.
PHYLIP, http://evolution.genetics.washington.edu/phylip.html.
BWA-MEM, http://arxiv.org/abs/1303.3997v2.
Bam-readcount, (https://github.com/genome/bam-readcount).
TCGA data portal, https://tcga-data.nci.nih.gov/tcga/.
Authors' Contributions
Conception and design: D.-C. Lin, P. Huang, J. Wang, B.P. Berman, E.-W. Song
Development of methodology: H.Q. Dinh, L.-w. Ding, B.P. Berman
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): L.-w. Ding
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): D.-C. Lin, A. Mayakonda, H.Q. Dinh, L.-w. Ding, B.P. Berman, D. Yin
Writing, review, and/or revision of the manuscript: D.-C. Lin, H.Q. Dinh, B.P. Berman, E.-W. Song, H.P. Koeffler
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): H.Q. Dinh, P. Huang, L. Lin, X. Liu, J. Wang
Study supervision: D.-C. Lin, B.P. Berman, E.-W. Song, D. Yin, H.P. Koeffler
Grant Support
This research was supported by the National Research Foundation Singapore under its Singapore Translational Research Investigator Award (NMRC/STaR/0021/2014) and administered by the Singapore Ministry of Health's National Medical Research Council (NMRC), the NMRC Centre Grant awarded to National University Cancer Institute of Singapore, the National Research Foundation Singapore, and the Singapore Ministry of Education under its Research Centres of Excellence initiatives to H.P. Koeffler. D-C. Lin was supported by National Center for Advancing Translational Sciences UCLA CTSI Grant UL1TR000124 and NSFC (81672786). D. Yin was supported by grants from the NSFC (81572484, 81420108026, 81621004) and Guangdong Science and Technology Department (2014A050503026, 2015B050501004). B.P. Berman and H.Q. Dinh were partially supported by Cedars-Sinai Department of Biomedical Sciences and Samuel Oschin Comprehensive Cancer Institute.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
Acknowledgments
We thank Dr. Hui Shen (Van Andel Institute) for providing a R script to analyze immune cell fraction and Dr. Dan Weisenberger (University of Southern California) for the help on methylation analysis. We are grateful for the generous donation from Blance and Steven Korgler as well as the Melamed family.
Footnotes
Note: Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).
- Received October 19, 2016.
- Revision received November 14, 2016.
- Accepted February 4, 2017.
- ©2017 American Association for Cancer Research.
References
Volume 77, Issue 9
