An integrated genomic and epigenomic approach predicts therapeutic response to zebularine in human liver cancer

Abstract

Epigenomic changes such as aberrant hypermethylation and subsequent atypical gene silencing are characteristic features of human cancer. Here, we report a comprehensive characterization of epigenomic modulation caused by zebularine, an effective DNA methylation inhibitor, in human liver cancer. Using transcriptomic and epigenomic profiling, we identified a zebularine response signature that classified liver cancer cell lines into two major subtypes with different drug responses. In drug-sensitive cell lines, zebularine caused inhibition of proliferation coupled with increased apoptosis, whereas drug-resistant cell lines showed up-regulation of oncogenic networks (for example, E2F1, MYC, and TNF) that drive liver cancer growth in vitro and in preclinical mouse models. Assessment of zebularine-based therapy in xenograft mouse models demonstrated potent therapeutic effects against tumors established from zebularine-sensitive but not zebularine-resistant liver cancer cells, leading to increased survival and decreased pulmonary metastasis. Integration of the zebularine gene expression and demethylation response signatures allowed differentiation of patients with hepatocellular carcinoma according to their survival and disease recurrence. This integrated signature identified a subclass of patients within the poor-survivor group that is likely to benefit from therapeutic agents that target the cancer epigenome.

Figure 1

Development of the zebularine response gene signature. (A) Unsupervised hierarchical clustering of gene expression data using the 323-gene zebularine response signature. The zebularine signature was derived from the gene lists of differential expressed genes computed at P≤.001 for each individual cell line using a corresponding control sample as a reference. The signature included a total of 631 genes universally changed by zebularine in HCC and 451 genes in CC cell lines. Comparison of the gene lists with those in the corresponding non-tumorigenic cell lines identified gene expression changes which were unique to either HCC (238 genes) or CC cell lines (85 genes). Integration and clustering of the signature grouped the cell lines into two clusters, resistant (red, non-responder) and sensitive (green, responder) depending on the zebularine-sensitivity. (B) Methylation profiling of cell lines treated with zebularine. A list of 133 demethylated genes (P≤.01, Δβ≥0.17) was identified by computing signatures for each cell line using a corresponding untreated sample as a reference (table S2). Note, that the cell lines clustered according to zebularine-sensitivity and corresponding xenographs response to zebularine therapy.

Figure 2

Evaluation of the efficacy of the zebularine therapy in vivo. (A) In vivo bioluminescence imaging of Huh7^luc+-xenografts. Tumor bearing mice were randomly assigned either to control (raloxifene, n=9) or treatment (zebularine and raloxifene, n=9) based on the intensity of bioluminescence before zebularine and raloxifene therapy (day 9). Representative images at day 7, 18, and 27 are shown. (A) Quantification of bioluminescence. The total flux is plotted as photon/s (C) Liver/body weight. (D) Pulmonary metastases. Treated mice (red bar) developed less lung metastases as compared to controls (blue bar). The isolated lungs were reimaged ex vivo and total flux was plotted as photon/s (log₁₀ scale). (E) Brain metastasis did not respond to therapy. The isolated brains were reimaged ex vivo and the total photon flux was reported as photon/s (log₁₀ scale). (F) Survival analysis for Huh7^luc+-recipients representing a separate experiment. Kaplan-Meier plot and Mantel Cox statistics demonstrate that the reduction in tumor growth was associated with increased overall survival. (G–I) The overall survival of Huh1^luc+- (G), KMCH^luc+- (H) and WRL68^luc+-recipients (I), following zebularine therapy. Note, either beneficial (G, H) or inferior (I) effects of zebularine therapy on survival of tumor-bearing mice. Bars represent the mean±s.d. Statistical significance was determined by Mann Whitney test (two-tailed).

Figure 3

Gene expression changes in response to zebularine therapy in vivo. (A) Hierarchical cluster analysis of liver tumors derived from the zebularine sensitive Huh7^luc+-xenografts identifies a list of 484 differentially expressed genes by Bootstrap-t (P≤.01, 5,000 repetitions). Control tumors, raloxifene alone (n=6), and the treatment group zebularine and raloxifene (n=7). (B, C) Ingenuity Pathway Analysis shows zebularine and raloxifene treatment caused a preferential deregulation of genes controlled by the oncogene MYC (B) and up-regulation of initiator caspases (e.g. CASP2 and CASP9) (C).

Figure 4

Prognostic survival genes. (A) Gene Set Enrichment Analysis (GESA) of the zebularine classifier. Analysis of the classifier revealed an Enrichment Score (ES) of the rank-ordered genes which demonstrated a significant positive correlation with the poor survival subtype A HCC. (B) Hierarchical cluster analysis of the validation data set. Genes which were significantly (P≤.01) associated with the disease outcome were identified by applying a Cox proportional hazards model and Wald statistics. 20 genes independently demonstrated a prognostic ability, based on the improved overall survival (C) and time to recurrence (D). Kaplan-Meier and Mantel-Cox statistics were used to determine levels of significance. (E) Meta-analysis using the zebularine classifier. The zebularine classifier is significantly associated to a study with a clinical outcome if the odds ratio >2, P<.0001 (Oncomine). Forest plot represent the odds ratio ± 95% CI.

Figure 5

Prediction of drug response. (A) Unsupervised clustering of subtype A HCCs based on the demethylation signature. The 133-gene list was significantly reduced to 16 genes which correctly predicted the classification using 7 different algorithms. (B) Sensitivity of the gene list to correctly predict the classification of patients within subtype A. The specificity is represented by the area under the receiver operating characteristic (ROC) curve (AUC, 95% CI) using Bayesian Compound Covariate prediction modeling during leave-one-out cross-validation. (C) Kaplan-Meier curves for survival risk prediction within subtype A HCC. Permutation p-value of the log-rank test statistic between risk groups was based on 100 permutations. (D) Meta-analysis using the demethylation signature. The signature is significantly associated to a study if the odds ratio>2, P<.0001 (Oncomine). In each incidence >25% of genes in the signature is represented in the top 10% under-expressed genes (responders) within the gene set analyzed. Forest plot represent the odds ratio ± 95% CI.