Network-based assessment of HDAC6 activity predicts preclinical and clinical responses to the HDAC6 inhibitor ricolinostat in breast cancer

Abstract

Inhibiting individual histone deacetylase (HDAC) is emerging as well-tolerated anticancer strategy compared with pan-HDAC inhibitors. Through preclinical studies, we demonstrated that the sensitivity to the leading HDAC6 inhibitor (HDAC6i) ricolinstat can be predicted by a computational network-based algorithm (HDAC6 score). Analysis of ~3,000 human breast cancers (BCs) showed that ~30% of them could benefice from HDAC6i therapy. Thus, we designed a phase 1b dose-escalation clinical trial to evaluate the activity of ricolinostat plus nab-paclitaxel in patients with metastatic BC (MBC) (NCT02632071). Study results showed that the two agents can be safely combined, that clinical activity is identified in patients with HR⁺/HER2^- disease and that the HDAC6 score has potential as predictive biomarker. Analysis of other tumor types also identified multiple cohorts with predicted sensitivity to HDAC6i's. Mechanistically, we have linked the anticancer activity of HDAC6i's to their ability to induce c-Myc hyperacetylation (ac-K148) promoting its proteasome-mediated degradation in sensitive cancer cells.

Extended Data Fig. 1. HDAC6 score in BC

(a) The number of samples and data set of origin used to evaluate the HDAC6 regulon. (b) Overlap of the original and updated HDAC6 regulons. P value was estimated using two-tailed Fisher’s exact test. N was determined by the total number of genes for network inference. (c) The network plot of HDAC6 and updated HDAC6 regulon. Edge width is corresponding to the correlation strength measured by mutual information. Red and blue edges indicate positive and negative correlations between HDAC6 and each of the regulon genes. The table below shows the pathway enrichment of the genes in the regulon, showing its association with unfolded protein response. P value was estimated using two-tailed Fisher’s exact test. (d) New HDAC6 score comparing IBCs vs non-IBCs. (e) HDAC6 scores of all BCs from TCGA and METABRIC are divided into molecular subtypes. (f) HDAC6 scores in 45 ductal metastatic breast cancer samples from the MBC Project divided into histological molecular subtypes. In d, e and f, the center line indicates the median value. The lower and upper hinges represent the 25th and 75th percentiles, respectively, and whiskers denote 1.5x interquartile range; the red line represents the median of the HDAC6 scores in IBC samples and the numbers over each whisker plot indicate the percentage of samples over this value in each clinical subtype. Sample size (n=number of samples) of each group was indicated in the axis labels. P value was estimated using two-tailed t test.

Extended Data Fig. 2. Anticancer activity of HDAC6 inhibitors

(a) The western blot shows the titration of ricolinostat to identify an effective dose (accumulation of Ac-α-Tubulin) without off-target effects in class-I HDACs (accumulation of Ac-H3K27). SAHA is used as a control Pan-HDAC inhibitor. WT-blot results were reproduced n=3 times from independent experiments. (b) Chemical structure of the different HDAC6 inhibitors used in Fig. 1E. (c) Normalized cell number 6 days after transfection with individualized siRNAs targeting HDAC6 or non-targeting control (NTC) (n=3 independent independent experiments per siRNA). The WT-blots show the silencing efficiency. WT-blot results were reproduced n=3 times from independent experiments. All error bars represent Mean±SD. P value was estimated by two-tailed t test. (d) The graphic shows the lack of synergistic activity between ricolinostat and commonly used chemotherapy (paclitaxel and doxorubicin) in MDA-MB-436 cells. In contrast, cells sensitive to ricolinostat MDA-MB-453, SK-BR-3 and MDA-MB-474 show synergistic activity between ricolinostat and chemotherapy. R and S indicate ricolinostat resistance and sensitivity respectively. N=3 independent replicate experiments per drug combination and concentration. (e) Histological intratumoral evaluation of H&E, Caspase-3, and Ki-67 in tumor samples from Fig. 1B. Quantification is also shown in bar graphs. Notice that the combo treatment (Pac+Ric) is not shown because all tumors regressed with this treatment. The white asterisks indicate necrotic areas and the white arrows indicate Caspase-3 positive stained cells. All error bars represent Mean±SD. P value was estimated by two-tailed t test. N= 6 samples from individual tumors. (f) The list shows all the transgenic mouse models evaluated by the HDAC6 score in Fig. 2c and d and indicates their correlation with human BCs. (g) Kaplan – Meier graphic showing the survival of the MMTV tumors in Fig. 2e. Control n=7; ricolinostat (Ric) n=8; paclitaxel (Pac) n=8; Ric+Pac n=8. P value was estimated using two tailed Log-Rank test.

Extended Data Fig. 3

Characteristics of the evaluable patients enrolled in the clinical trial

Extended Data Fig. 4. Biomarker evaluation of HDAC6 score calculated by using the NetBID and VIPER algorithms in the phase Ib trial

The Supplementary Extended Data Fig. shows the similarities between the HDAC6 scores inferred by NetBID (a) and VIPER (b) in responders and non-responders, and (c) ROC curve plot of HDAC6 score inferred by VIPER (similar plot for NetBID is in Fig. 3d). In a and b, one dot represents one patient sample and the center line indicates the median value. The lower and upper hinges represent the 25th and 75th percentiles, respectively, and whiskers denote 1.5x interquartile range. P value was estimated by two-tailed t test. For all panels n=3 non-responders and n=7 responders.

Extended Data Fig. 5. The HDAC6 score does not correlate with the response in paclitaxel-only treated patients

The Supplementary Extended Data Fig. shows the HDAC6 scores of all BCs from Hatzis et al. (JAMA, 2011) divided into clinical (a) and molecular (b) subtypes. (c) HDAC6 scores of patients divided by response to paclitaxel, Residual Disease (CD) and Pathological Complete Response (pCR), showing the lack of correlation between the HDAC6 score and the response to paclitaxel. In a, b and c the number of patient samples is indicate as (n). d, Kaplan – Meier graphic showing the survival of breast cancer patients treated exclusively with paclitaxel in the neoadjuvant setting separated by HDAC6 score (high/low=higher and lower −0.36, based on the ROC analysis in our clinical trial). N = 35 patient samples for HDAC6 score-low and N = 71 patient samples for HDAC6 score-high group. P value was estimated using two-tailed Log-Rank test. The 95% confidence interval of the regression lines were displayed. In a, b and c, the center line indicates the median value. The lower and upper hinges represent the 25th and 75th percentiles, respectively, and whiskers denote 1.5x interquartile range; the red line represents the median of the HDAC6 scores in IBC samples and the numbers under each whisker plot indicate the percentage of samples over this value in each clinical subtype. Sample size of each group was indicated in the axis labels. P value was estimated using two-tailed t test.

Extended Data Fig. 6. HDAC6 scores in other human cancers

(a) Correlation between HDAC6 regulon in different tumor types. The tumors sets are the same as in figure 4c. The number of patient sample for each set is indicated there. LAML: Acute Myeloid Leukemia; LIHC: Liver Hepatocellular Carcinoma; UVM: Uveal Melanoma; KIRP: Kidney Renal Papillary Cell Carcinoma; PCPG: Pheochromocytoma and Paraganglioma; DLBC: Diffuse Large B-Cell Lymphoma; ACC: Adrenocortical Carcinoma; UCS: Uterine Carcinosarcoma; KIRC: Kidney Renal Clear Cell Carcinoma; THYM: Thymoma; CHOL: Cholangiocarcinoma; SKCM: Skin Cutaneous Melanoma; CRC: Colorectal Carcinoma; LGG: Brain Lower Grade Glioma; UCEC: Uterine Corpus Endometrial Carcinoma; KICH: Kidney Chromophobe; MESO: Mesothelioma; TGCT: Testicular Germ Cell Tumors; GBM: Glioblastoma; PRAD: Prostate Adenocarcinoma; CESC: Cervical Squamous Cell Carcinoma and Endocervical Adenocarcinoma; BLCA: Bladder Urothelial Carcinoma; SARC: Sarcoma; OV: Ovarian Serous Cystadenocarcinoma; THCA: Thyroid Carcinoma; ESCA: Esophageal Carcinoma; BRCA: BC; LUAD: Lung Adenocarcinoma; HNSC: Head and Neck Squamous Cell Carcinoma; PAAD: Pancreatic Adenocarcinoma; STAD: Stomach Adenocarcinoma; LUSC: Lung Squamous Cell Carcinoma. P value was estimated by two-tailed Fisher’s exact test. (b) List with all the cell lines evaluated by dose-response to ricolinostat and HDAC6 scores. The correlation (R) and P value between the response to ricolinostat and HDAC6 scores were estimated by two-tailed Spearman correlation test. (c) Graphic showing the correlation between the HDAC6 score and the response to ricolinostat in individual cell types (only tumor types with more than 6 cell lines are shown). N=8 individual independent experiments for each ricolinostat dose. The curve was fitted by stat_smooth algorithsm using lm smoothing method and y~log2(x) formula.

Extended Data Fig. 7. The correlation of HDAC6 score with immune infiltrates

(a) Scatter plot showing the correlation between HDAC6 score and immune score by ESTIMATE in TCGA-BRCA primary patient samples (n=1109). The correlation coefficient (R) and P value were estimated using Spearman correlation test. (b), Violin plot showing the distribution of immune score across IHC-based breast cancer subtypes. Sample size of each group was indicated in the axis labels. P value was estimated using two-tailed t test. The center line indicates the median value. The lower and upper hinges represent the 25th and 75th percentiles, respectively, and whiskers denote 1.5x interquartile range.

Extended Data Fig. 8. Ricolinostat treatment reduces the expression of c-MYC in sensitive cell lines

(a) Heatmap representing GSEA analysis of hallmark signatures during ricolinostat exposure in sensitive breast cancer cell lines and TgMMTV-Neu model. P value was estimated by two-tailed t test. The Z-scores were transformed from these P values and further combined using Stouffer’s method. Only significant (combined Z > 1.96 or < −1.96) sets are shown. (b) The graphic shows summarized z-scores in cell lines and TgMMTV-Neu sensitive to ricolinostat during a time curse treatment (6, 12, 24 hours and 4 weeks). For a and b N = 2 individual independent experiments for cell lines and N=3 individual tumors for TgMMTV-Neu. P value was estimated using two-tailed t test. (c) Bubble plot representing GSEA analysis of multiple MYC-associated signatures after ricolinostat in sensitive and resistant cells. N = 3 independent experiments per cell line. P value was estimated by two-tailed t test. The Z-score was transformed from the P values and further combined by Stouffer’s method. (d) QRT-PCR of c-Myc mRNA expression after 6 hours of exposure to ricolinostat. N = 3 independent experiments for each time point All effort bars represent Mean±SD. P value was estimated by two-tailed t test. (e) The WT-blot shows the changes in the protein expression of c-Myc and ac-Tubulin in SUM-149 cells after ricolinostat is added to the culture media. WT-blot results were reproduced n=3 times from independent experiments. (f) Graphic showing an efficient reduction in the HDAC6 score after treatment with ricolinostat in multiple cell lines, n≥2 independent experiments per time point. The center line indicates the median value. The lower and upper hinges represent the 25th and 75th percentiles, respectively, and whiskers denote 1.5x interquartile range.

Extended Data Fig. 9. Acetylation of c-Myc in Lys148 after inhibition of HDAC6

a) MA plots showing the peptides upregulated upon HDAC6 knockout (above) and HDAC6 catalytic domain 2 mutants (below). In the MA-plots, each dot represents a peptide. The significantly upregulated peptides were identified by fold change > 1.5 and p-value < 0.05 and highlighted in red. P value was estimated by two-tailed t test. (b) Scatter plot showing the correlation of the Z-score of the comparison between HDAC6 KO and wild type with that of the comparison between HDAC6 mutant and wild type. The curve was fitted by stat_smooth algorithm using lm smoothing method and y~x formula. The correlation coefficient (R) and P value were estimated using two-tailed Spearman correlation test. Each dot represents a peptide. For a and b the N = 2 independent proteomic replica studies per cell line. c) The western blot shows the accumulation of ac-K148-c-Myc after HDAC6 is inhibited by ricolinostat in MDA-MB-453 and SK-BR-3 BC lines. WT-blot results were reproduced n=3 times from independent experiments.

Extended Data Fig. 10. Acetylation of c-Myc in Lys148 after inhibition of HDAC6 in sensitive and resistant BC cancer cells

The western blot shows the accumulation of ac-K148-c-Myc after HDAC6 is inhibited by ricolinostat in MDA-MB-453 (sensitive) and MDA-MB-436 (resistant) lines. WT-blot results were reproduced n=3 times from independent experiments.

Fig. 1.. The HDAC6 score identifies breast cancers sensitive to the HDAC6 inhibitor ricolinostat.

HDAC6 scores of BC patients from TCGA and METABRIC cohorts (a) and BC cell lines from CCLE database (b) are divided into subtypes. The red line represents the median of the HDAC6 scores in IBC samples and the numbers under each whisker plot indicate the percentage of samples over this value in each clinical subtype. Sample size of each group was indicated in the axis labels. P value was estimated using two-tailed t test. For (a) and (b) the center line indicates the median value. The lower and upper hinges represent the 25th and 75th percentiles, respectively, and whiskers denote 1.5x interquartile range. (c) Graphical representation showing the correlation between HDAC6 score and response to ricolinostat among 14 breast cancer cell lines (n=8 independent experiments/per ricolinostat concentration). The curve was fitted by stat_smooth algorithm using lm smoothing method and y~log₂(x) formula. The correlation coefficient (R) and P value were estimated using the Spearman correlation test. (d) Annexin-V/PE staining comparing the apoptotic response after Ricolinostat treatment of sensitive (MDA-MB-453) vs resistant (MDA-MB-436) BC cells. Quantification of alive, apoptotic, and dead cells is provided and visualized in stacked bar plot at the bottom of the panel. (e) Growth curves of sensitive (S) vs resistant (R) BC cells treated with four different HDAC6 inhibitors (n=3 independent experiments/per drug concentration). (f) Normalized cell number 6 days after transfection with siRNAs targeting HDAC6 or non-targeting control (NTC) (n=3 independent experiments/siRNA). The WT-blots show the silencing efficiency (The WT-blot results were reproduced n=3 times from independent experiments/siRNA). All error bars represent Mean±SD. P value was estimated by two-tailed t test.

Fig. 2.. Anticancer activity of ricolinostat in vivo.

(a) Treatment of ricolinostat sensitive (MDA-MB-453) and resistant (MDA-MB-436) growing as xenografts in SCID mice. The cartoon illustrates the treatment regimen. The combinatorial effect with paclitaxel was also investigated in sensitive cells. On resistant cells, only ricolinostat was used because no effect was observed with the combo in vitro. The western blots show the accumulation of acetylated tubulin in tumors treated with ricolinostat. Additionally, the absence of off-target effect in class-I HDACs is shown by the minimum changes seen in the levels of acetylated Histone-3-K27 (two independent tumor samples are shown), (WT-blot results were reproduced n=3 times from independent experiments). For the growth curve, n=number of tumors for each of the treatment cohorts. (b) HDAC6 scores in tumors emerging in transgenic mouse models that recapitulate the molecular characteristics of human BCs (the number of samples is shown for each group in the figure). (c) HDAC6 scores in tumors emerging in 27 different transgenic mouse models of BC (with and without molecular characteristics of human BCs), (the number of samples is shown for each group in the figure). In C and D the red line represents the mean of the HDAC6 scores in IBC samples. In c and d, the red line represents the mean of the HDAC6 scores in IBC samples; The center line indicates the median value. The lower and upper hinges represent the 25th and 75th percentiles, respectively, and whiskers denote 1.5x interquartile range. (d) Treatment of BC tumors emerging in the MMTV_Neu transgenic mouse model. The beginning and end of treatment are indicated by the red arrows. The cartoon illustrates the treatment regimen. The combinatorial effect of ricolinostat plus paclitaxel was also investigated. The western blots indicate the same that in panel B (WT-blot results were reproduced n=3 times from independent experiment). For the growth curve, n=number of tumors for each of the treatment cohorts. (e) Histological intratumor evaluation of H&E, Caspase-3, and Ki-67 in tumor samples from panel E. Quantification is also shown in bar graphs. The white asterisks indicate necrotic areas and the white arrows indicate Caspase-3 positive stained cells. All error bars represent mean±SD and P value was estimated by two-tailed t test.

Fig. 3.. Phase 1b trial of ricolinostat combined with nab-paclitaxel in metastatic breast cancer.

(a) Graphical description of the clinical study. (b) Waterfall plot showing the tumor best response for patients with measurable disease. Of the 16 patients, 3 were TNBCs (1 showed stable disease (SD) and 2 showed progressive disease (PD). The rest 13 patients were of the HR+/HER2− subtype (2 showed a partial response (PR), 9 showed SD and 2 showed PD. Note that one evaluable patient with stable disease did not have measurable disease and is not included in the waterfall plot (n=15). (c) The bar graph shows the HDAC6 scores in the patients in the trial (labeled in red) together with all the BC samples evaluated and separated by subtype. Labeled in blue are the IBC and the matched non-IBC series. The center line indicates the median value. The lower and upper hinges represent the 25th and 75th percentiles, respectively, and whiskers denote 1.5x interquartile range. The red line represents the median of the HDAC6 scores in IBC samples and the numbers above each whisker plot indicate the percentage of samples over this value in each clinical subtype. Sample size of each group was indicated in the axis labels. The full list of p values can be found summarized in the Source Data and was estimated using two-tailed t test. (d) Receiver operating characteristic (ROC) curve plot for evaluation of HDAC6 score to predict the response of patients with breast cancer to ricilinostat from the clinical trial. The recommended cutoff of the HDAC6 score and corresponding sensitivity, specificity, and accuracy were inside the box (n=2 independent HDAC6 score replicates per patient) (e) Kaplan – Meier graphic showing the survival of the patients in the study separated by HDAC6 score (high/low= higher and lower than −0.36, the cutoff ofHDAC6 score based on the ROC analysis in the study. In this study 10 out of 16 evaluable patients had tissue available for translational analyses. P value was estimated using two-tailed Log-Rank test.

Fig. 4.. The HDAC6 score correlates with the response to ricolinostat in other cancer types.

(a) Graphic showing the correlation between HDAC6 score and response to ricolinostat in 72 human cancer cell lines. The curve was fitted by stat_smooth algorithm using lm smoothing method and y~log₂(x) formula. The correlation coefficient (R) and P value were estimated using two-tailed Spearman correlation test. N=6 independent experiments per cell line and ricolinostat dose. (b) HDAC6 scores were calculated for 1,156 different cancer cell lines available in the CCLE database representing 20 different human cancers and (c) for over 10,000 molecular profiles representing 32 different types of human cancer that have been collected in the TCGA database. For b and c the number of cell line/patient samples is shown for each group on the figure. The center line indicates the median value. The lower and upper hinges represent the 25th and 75th percentiles, respectively, and whiskers denote 1.5x interquartile range. As N we used all independent replicates available in the CCLE and TCGA data sets, precise information for each sample is available in CCLE and TCGA. LAML: Acute Myeloid Leukemia; LIHC: Liver Hepatocellular Carcinoma; UVM: Uveal Melanoma; KIRP: Kidney Renal Papillary Cell Carcinoma; PCPG: Pheochromocytoma and Paraganglioma; DLBC: Diffuse Large B-Cell Lymphoma; ACC: Adrenocortical Carcinoma; UCS: Uterine Carcinosarcoma; KIRC: Kidney Renal Clear Cell Carcinoma; THYM: Thymoma; CHOL: Cholangiocarcinoma; SKCM: Skin Cutaneous Melanoma; CRC: Colorectal Carcinoma; LGG: Brain Lower Grade Glioma; UCEC: Uterine Corpus Endometrial Carcinoma; KICH: Kidney Chromophobe; MESO: Mesothelioma; TGCT: Testicular Germ Cell Tumors; GBM: Glioblastoma; PRAD: Prostate Adenocarcinoma; CESC: Cervical Squamous Cell Carcinoma and Endocervical Adenocarcinoma; BLCA: Bladder Urothelial Carcinoma; SARC: Sarcoma; OV: Ovarian Serous Cystadenocarcinoma; THCA: Thyroid Carcinoma; ESCA: Esophageal Carcinoma; BRCA: BC; LUAD: Lung Adenocarcinoma; HNSC: Head and Neck Squamous Cell Carcinoma; PAAD: Pancreatic Adenocarcinoma; STAD: Stomach Adenocarcinoma; LUSC: Lung Squamous Cell Carcinoma. (d) Correlation of ranks of HDAC6 score among 22 cancer types between cancer cell lines from CCLE and patient cancer samples from TCGA. The cell lines used for this comparison were the same used in panels b and c (n=patients and cancer lines used in those panels). The curve was fitted by stat_smooth algorithsm using lm smoothing method and y~x formula, and the 95% confidence interval of the regression line is displayed as shaded. The correlation coefficient (R) and P value were estimated using two-tailed Spearman correlation test.

Fig. 5.. Treatment with ricolinostat induces a robust reduction of MYC expression and activity.

(a) Bubble plot representing GSEA analysis of MYC signatures during ricolinostat exposure in sensitive breast cancer cells. The size and intensity of the bubble indicate the statistical significance. N = 2 independent experiments for each time point of each BC cell line or mouse model. P value was estimated by two-tailed t test. The Z-scores were transformed from these P values and further combined using Stouffer’s method. (b) The upper panel (WT-blots) shows the reduction of MYC protein expression after ricolinostat treatment (2.5 μM) in sensitive but not resistant cancer cells. WT-blot results were reproduced n=3 times from independent experiments. The lower panel shows the summarized z-scores, comparing Ricolinostat-treated (6hrs) vs. untreated in the same cell lines shown in the upper panel. The Z-score was transformed from the P value estimated by two-tailed t test. (c) WT-blots showing the reduction of MYC in sensitive cell lines when HDAC6 was inhibited by the small molecule inhibitor CAY10603 (1μM, upper panel) or by RNAi (100nM, lower panel) (c=non-targeting siRNA control, H6=siRNA targeting HDAC6, M= siRNA targeting MYC). WT-blot results were reproduced n=3 times from independent experiments. (d) Cell viability (cell number) after MYC is silenced by siRNA in ricolinostat sensitive cell lines (notice that the western blot showing efficient silencing of MYC is shown in Fig. 5c). N=3 independent experiments per siRNA and cell line. Data are presented as mean values +/−SD. P values was estimated by two-tailed t test.

Fig. 6.. HDAC6 modulates the acetylation level of MYC.

(a) Co-IP of MYC and HDAC6 in HEK-293T cells. The asterisk and arrow indicate the endogenous and the transduced c-MYC respectively (transduced construct expresses a slightly larger form). IP results were reproduced n=3 times from independent experiments. (b) Schematic description of the proteomic study described in the text. The WT-blot shows the accumulation of acetylated α-Tubulin in HDAC6 deficient HAP1 cells. The MS/MS spectrum shows an example assigned to the peptide containing MYC K148 acetylation site with b- and y- ions corresponding to the N- and C-terminal fragments, respectively. Peaks that match to theoretically calculated fragmented ions of the lysine-acetylated peptide are indicated. Modifications on specific residuals are indicated for TMT (+229 Da) and acetylation (+42 Da), respectively. WT-blot results were reproduced n=3 times from independent experiments. The dot plot shows the top differentially acetylated proteins from the proteomic study. The vertical and horizontal dash lines indicated the cutoffs of fold change > 1.5 and P < 0.05, respectively. N = 2 for each group. P value was estimated using two-tailed t test. (c) Degradation of MYC correlates with the accumulation of the ac-K148 form during treatment of MDA-MB-453 with ricolinostat (2.5μM). WT-blot results were reproduced n=3 times per time point from independent experiments. (d) Proteasome inhibition by Bortezomib (2.5nM) blocks the reduction of c-MYC protein induced by ricolinostat (2.5μM). WT-blot results were reproduced n=3 times per drug combination from experimental replicates. (e) Mechanistic model of the effect of HDAC6 inhibition on MYC expression. C-MYC stability is influenced by posttranslational modification. Acetylation of K148 promotes degradation by the proteasome and it is prevented by HDAC6. Thus, HDAC6 inhibition leads to hyperacetylation of MYC leading to its degradation.