Expression of HDACs in pan-cancer

A comprehensive visual representation of the research process can be found in Fig. 1. Leveraging data from the TCGA repository, we investigated the expression levels of the HDAC family in tumor and normal samples across a variety of malignancies. Firstly, our results determined that most of the 11 members of HDACs are highly expressed in pan-cancer, especially HDAC1, while low expression was observed for HDAC9 (Fig. 2A). And, expression of each member in the 33 cancer tissues illustrated separately and ranked them from higher- lowest (Additional file 1: Figure S1). Then, results shown across 18 distinct cancer forms, there had been significant variations among HDAC family members. For example, HDAC7 was the most highly expressed in CHOL and HDAC11 was the least expressed in GBM (Fig. 2B). The level of expression of 11 HDAC family members was then found to be highly associated with one another in all types of cancers using Spearman’s correlation analysis. We learned that HDAC3 and HDAC1, HDAC6 and HDAC5 were the two most significantly positively correlated gene pairs (correlation coefficient = 0.33); HDAC5 and HDAC1 were the two most negatively correlated genes (correlation coefficient = -0.21) (Fig. 2C). Weak correlations between them can be seen, and these results suggest intrinsic differences in HDAC expression among different HDAC family members or different cancer types. Thus, a specific HDAC gene may exhibit either oncogenic or anti-oncogenic properties in different cancer types, highlighting the necessity for individual in-depth study of each HDAC family member.

Fig. 1

Diagram showing the study’s flow regarding HDAC family genes and human cancer

Fig. 2

Differential analysis of gene expression in the HDAC family. (A) A box plot illustrating the levels of transcriptional expression for HDACs. (B) Based on TCGA data, this heatmap illustrates the variations in HDAC activity between nearby healthy tissues and cancerous cells across a variety of cancer types. (C) Correlation coefficients by co-expression analysis between each two HDAC family genes

The levels of mRNA expression of HDACs are compared between cancer tissues and typical tissues (Additional file 2: Figure S2). Most HDAC-related genes have been discovered to be significantly expressed in a variety of tumor types, for instance HDAC1 expression was only decreased in KICH, KIRC, and KIRP; HDAC3 showed low expression in KICH, THCA; HDAC8 was down-regulated only in KIRC; and HDAC10 exhibited low expression only in KICH. However, we also found a certain degree of intra- and inter-cancer heterogeneity in the expression levels of certain HDACs, such as HDAC4 was overexpressed only in CHOL and LIHC, HDAC2 was down-regulated only in KICH and KIRC, which was in contrast to other tumor types.

Analysis of therapeutic correlation in pan-cancer

Using univariate Cox regression analysis, we studied the prognostic value of HDAC levels across different forms of cancer. As indicated in Fig. 3, the expression of some HDAC genes was related with OS, DFS, PFS and DSS in various types of cancer. The prognostic impact of HDACs exhibited some degree of variability in different cancer types. The patients were divided into groups with low or high expression following the average expression of gene values. Subsequently, we performed Kaplan–Meier analysis for the HDAC1-11 genes to assess OS, DFS, DSS, and PFS (Additional file 3: Table S1; Additional file 4: Additional Data1). Among these patients, high HDAC1 activity has been related to an unfavorable outcome in KIRC sufferers (Fig. 4A) and LGG (Fig. 4B). Conversely THYM and KIRC individuals showed less favorable results when their HDAC2 expression was decreased (Fig. 4C, D). In contrast, those having significant levels of HDAC3 in SKCM and THCA had a negative prediction (Fig. 4E and F, respectively). Besides, HDAC4 (Fig. 4G) and HDAC10 (Fig. 4K) are high-risk genes in patients with ACC; HDAC9 plays a high-risk role in THYM (Fig. 4J). Contrarily, increased expression was observed for HDAC7 in BLCA (Fig. 4H) and SKCM (P < 0.001; Fig. 4I) patients had a favorable prognosis. Similarly, the higher levels of HDAC11 was related to better prognostic results in KIRP (Fig. 4L). Significant correlation was observed between the TNM staging and the expression of HDAC family. Specifically, in KICH and TGCT, the expression level of HDAC3 showed an increase as the stage progressed (Fig. 5A). Conversely, the levels of HDAC11 expression in stage I versus stage II/III/IV were found to be the highest in KIRC and KIRP, gradually decreasing as the stage progressed (Fig. 5B). The supplementary findings of other family genes are represented in Additional file 5: Figure S3.

Fig. 3

The forest graphs show the association between the activity of HDAC and OS (A), PFS (B), DFS (C), and DSS (D) by the univariate COX method in pan-cancer

Fig. 4

K-M survival curves highlight the link between overall survival and the HDAC expression variances of both low and high expression cohorts. HDAC1 in KIRC (A) and in LGG (B), HDAC2 in KIRC (C) and in THYM (D), HDAC3 in THCA (E) and in SKCM (F), HDAC4 in ACC (G), HDAC7 in KIRC (H) and in SKCM (I), HDAC9 in THYM (J), HDAC10 in ACC (K), HDAC11 in KIRP (L)

Fig. 5

The association of HDACs expression and clinicopathological stage in various cancer types. (A, B) Boxplot graph indicate the characteristics of HDAC11 and HDAC3 in relation to clinicopathological staging, respectively

Genetic alterations of HDACs in pan-cancer

We examined the situation of genetic variations in HDAC using the cBioPortal database. The HDAC family exhibited changes in 27.32% of the total 2683 cases across various tumor types, as observed from the ICGC-TCGA dataset (Fig. 6A). According to Fig. 6B, the occurrence rates of genetic mutations in HDAC varied from 1.5% to 10%. Notably, HDAC6 and HDAC9 exhibited higher mutation frequencies, with 8% and 10%, respectively, compared to other members. Amplification was the most frequent mutation type across most members in pan-cancer patients. With a frequency of mutations over 40%, the variant rate of HDACs was comparatively high in the following tumors: endometrial, head and neck, soft tissue sarcoma, lung, ovarian, esophagogastric, and bladder. In terms of amplification mutations, the incidence was higher in endometrial cancer, breast cancer, esophagogastric cancer, lung cancer, melanoma, cervical cancer (Fig. 6C). Furthermore, the more intuitive percentage of genetic changes of each HDAC family member in pan-cancer and the degree of genetic changes for each member in various cancers are detailed in Additional file 6: Figure S4. With a p-value of 0.082, the KM plotter results revealed that there was no significant difference in the OS between the groups with or without the HDAC familial gene substitution (Fig. 6D).

Fig. 6

Alterations of the HDAC family in pan-cancer. (A) The HDACs altered in 27.32% of 2683 cases. (B, C) HDACs variant types and frequency in pan-cancer. (D) Analysis of overall survival in patients with HDACs mutation status using the cBioPortal database

In terms of CNVs, shallow deletions and gain were more prominent than deep deletion or amplification. Moreover, a statistically robust rise in the expression levels of most HDACs was observed with the rise in copy number variation (CNV). This finding aligns with the comprehensive analysis of whole genome data conducted by the ICGC/TCGA pan-cancer consortium in 2020 (Fig. 7). The study of the pan-cancer findings retrieved from the ICGC/TCGA, we also determined that the non-mutated regions accounted for a larger proportion compared to those with missense, truncating, multiple, inframe, or splice mutations. Additionally, significant differences in expression levels were observed for HDAC3/9 genes among different genomic variant types, while no statistical differences were found for other members (Additional file 7: Figure S5). These data sets demonstrated that the upregulation of HDAC mRNA in pan-cancer may mostly be attributed to an increase in the number of HDAC genomic copies.

Fig. 7

Correlation analysis between the expression of each HDAC gene and its copy number (A-K)

Interrelation between HDAC expression and immune subtypes, tumor microenvironment, and stemness indices

To examine the correlation involving tumor immune cell infiltration and HDACs, we delved into the association between immunological cells and HDACs. A previous analysis of human malignancies identified six distinct immune subtypes: C1 (wound healing), C2 (INF-r dominant), C3 (inflammation), C4 (lymphocyte failure), C5 (Immune quiet), and C6 (TGF-β predominance), and influencing tumor progression and development [27]. It has been demonstrated that patients belonging to immune subtypes C3 and C5 exhibit significantly higher survival rates compared to those with different subtypes. However, patients classified as C4 and C6 subtypes have the lowest survival rates [27]. With a p-value less than 0.001, this study found a strong association between HDAC expression and each of the six immunological subtypes (Fig. 8A). Further analysis reveals that high expression of HDAC4/5/6/9/11 was closely correlated with C5 subtype, suggesting that these members may serve as a cancer-suppressing effect. In contrast, HDAC7 appears to function as a pro-oncogene in the C6 subtype compared to other subtypes. Similarly, high expression of HDAC1/2/3 was more relevant in C1, C2, and C6 infiltrates compared to other subtypes, illustrating that these members have a tumor-promoting role, and predicting poorer survival in these patient categories. Several studies have shown that the presence of lymphocytes infiltrating the TME is a significant indicator in predicting both tumor sentinel node status and survival duration [28]. The association between HDACs and CD8 + T cell infiltration levels determined by seven distinct algorithms has been studied using TIMER. The level of expression of HDACs proved to be positively correlated with CD8 + T cell recruitment in multiple tumor types, as depicted in Fig. 8B. Estimate method was used to examine the correlation between HDAC expression levels and stromal and immune scores in tumor samples. Most notably, the results of immune Score and stromalScore showed almost positive correlation with the expression of both HDAC7 and HDAC9 (Fig. 8C and D). Specifically, HDAC1 exhibited a strong positive correlation with the immuneScore of UVM and LGG, and HDAC4 demonstrated a similar correlation with DLBC. In the case of TGCT, the immune score had a significant unfavorable tandem with both HDAC2 and HDAC11. Regarding StromalScore, in GBM and SARC, most tumors were negatively correlated except for HDAC7, which showed a remarkably positive associated with StromalScore. Moreover, HDAC6 and HDAC10 were consistently negatively related to stromal core in almost the tumors. The association between the expression of genes and tumor stemness (RNAss and DNAss) was determined using Spearman correlation analysis (Fig. 8E and F). For RNAss, we illustrated that a vigorously positive association exist between HDAC1, HDAC2, HDAC3 and most tumors. Meanwhile, an appreciably negative correlation was seen between HDAC1 and KIRC, LGG, and PCPG. Particularly, a positive correlation was observed in HDAC1, HDAC2, HDAC3, HDAC7 and HDAC10 in the THYM. We also identified a significant negative interrelation between THYM and HDAC1/2/3/7/10 in DNAss.

Fig. 8

Outcomes of correlation analysis of HDACs with immune subtypes, stemness index, and microenvironmental score. (A) Transcript levels of HDACs expression in multiple cancer types in C1-C6 immune subtypes. (B) Applying seven different algorithms, we examined the probable link between HDACs expression and the CD8+ T-cells infiltration. (C, D) A bubble plot illustrating the association between immunology value and stemness index with the proliferation of HDACs mRNA. (E, F) The association between stemness indices and HDAC gene expression across different cancer types is displayed by a bubble graph

Association analyses of HDACs with TMB, MSI and immune checkpoints

Given the widespread acceptance of immuno-surveillance in the prognosis of different tumor forms [29]. We studied the association between HDACs expression and 47 immune-mediated checkpoints that are frequently observed and are strongly linked to the manifestation of immunotherapy (Additional file 8: Figure S6), for instance, CTLA-4, PDL-1, and PD-1 (Fig. 9A). For instance, in LIHC, the expression of HDAC1-5 and HDAC7-9 was positively correlated with the mRNA levels of PD-L1, PD-1 and CTLA-4. HDAC1/3/8 was also identified to be positively associated with this gene expression in UVM. These results demonstrated that overexpression of HDACs may mediate immune evasion in some tumors. TMB and MSI are widely recognized as significant indicators of tumor development and advancement [24]. Additionally, investigations were conducted to explore the associations between HDAC expression, TMB, and MSI. A negative correlation was observed between HDAC1-4, HDAC6, HDAC7, and HDAC10 and TMB in THYM, whereas there was a positive association between TMB and HDAC11 in THYM, ESCA, and KIRP. It was observed that HDAC1-4, HDAC6, HDAC7, and HDAC10 are adversely associated to the TMB in THYM, while HDAC11 is positively correlated with the TMB in THYM, ESCA, and KIRP. HDAC4 was negatively associated with TMB in BLCA, BRCA, DLBC, etc., and positively correlated in ACC only (Fig. 9B, Additional file 9: Figure S7). Furthermore, HDAC1 showed an adverse correlation with MSI in CESC and a significant association with MSI in COAD, ESCA, LAML, etc. Moreover, there is merely an association between HDAC9 activity and MSI in MESO, and an unfavorable association with MSI in DLBC, ESCA, HNSC, etc. Notably, no malignancies were shown to be adversely correlated with MSI, whereas HDAC8 expression had a favorable relationship with MSI in ACC, KIRP, HNSC, etc. (Fig. 9C, Additional file 10: Figure S8).

Fig. 9

The relationship between HDAC expression and immune checkpoint genes, such as CTLA4, PD-1, and PD-L1 (A), as well as TMB (B) and MSI (C) in 33 tumor types

The role of HDACs expression in immunotherapy of independent cohorts

A correlation analysis was conducted between HDAC expression and immunotherapy effectiveness in renal cell carcinoma, bladder cancer, and melanoma. IMvigor210, GSE111636, GSE176307, GSE78220, GSE67501, as well as our own TRUCE-01 dataset were used for this analysis. The study shows that there was a significant correlation (P < 0.09) between the expression of HDAC1/2/3 and the therapeutic objective response to anti-PD-L1. Interestingly, the GSE111636 cohort also demonstrated comparable results for HDAC3, with a statistically significant level of correlation (P = 0.017). We found that the expression levels of HDAC6 and HDAC10 had a poor correlation (P = 0.03) with the response rate to immunotherapy in the GSE111636 and TRUCE-01 datasets.

Other members should be regarded as of borderline statistical significance (0.1 ≤ P < 0.3). In summary, the findings suggest that assessing the expression of HDACs may serve as a valuable predictor for determining the responsiveness of anti-PD-1/PD-L1 immunotherapy.

Furthermore, we explored the possible mechanisms behind HDAC manifestations and immunotherapy response. In the present single-arm phase 2 clinical trial (TRUCE-01, NCT04730219), a notable increase in HDAC4 or HDAC9 expression was observed in BC patients who received tislelizumab combined with nab-paclitaxel therapy, regardless of their response to treatment. In non-responsive cases, however, there was a decrease in expression of HDAC7 after treatment (Fig. 10B–G). Therefore, these findings suggest a strong association between the HDAC gene expression and observed immunotherapy responses. This aspect will be further investigated in our future research endeavors.

Fig. 10

The association between HDACs expression and clinical immunotherapy effect from several immunotherapy datasets. (A) Expression of HDACs in different cohorts based on their response to immunotherapy. (B-G) Changes of HDAC1/4/6/7/8/9 expression level before and after immunotherapy based on our Truce01 dataset

Drug sensitivity analysis

In order to examine whether HDAC expression may be associated with anti-cancer drug sensitivity in different types of cancer, we obtained relevant data from the CellMiner™ database. The results indicated that sensitivity to various chemotherapeutic drugs correlates significantly with HDAC expression (Additional file 11: Table S2, P < 0.05). HDAC7 expression was found to be inversely correlated and the sensitivity of Selumetinib, Cobimetinib, Trametinib and PD-98059, in contrast to a positive correlation between HDAC7 expression and the sensitivity of Everolimus, Rapamycin, Temsirolimus. The expression of HDAC4, PX-316, and chelerythrine drug sensitivity was found to be favorably correlated. Chelerythrine and Acrichine exhibited a favorable correlation with HDAC1 activity. Moreover, pharmaceutical tolerance to Ifosfamide, Oxiliplatin, Imexon, Lomustine, BN-2629, and Eribulin mesylate was inversely correlated with HDAC11 expression. Beside this, correlations results between HDAC expression and the sensitivity of commonly used chemotherapeutic agents against tumors are presented herein (Fig. 11B,  P < 0.05). Paclitaxel, Vinblastine, Doxorubicin, or Docetaxel IC50 have an adverse association with HDAC7/11 activity. However, there was a beneficial association between gemcitabine's IC50 and HDAC3.

Fig. 11

Pan-cancer drug sensitivity and HDACs. (A) We identified the top 25 drugs that were significantly correlated with HDACs expression (P < 0.01). (B) The association between the expression of HDACs and commonly employed chemotherapeutic drugs in cancer treatment is noteworthy (P < 0.05)

GSEA analysis of HDACs expression in bladder cancer

To assess the biological function and signaling pathways, we employed GSEA to identify enriched KEGG pathways within the groups with high and low levels of each HDAC gene family (Additional file 12: Table S3). High- or low-expression groups of different HDACs in bladder cancer showed significantly differential enrichment of tumor- and immune-related pathways including “CYTOKINE CYTOKINE RECEPTOR INTERACTION”, “OLFACTORY TRANSDUCTION”, “JAK STAT SIGNALING PATHWAY”, “CALCIUM SIGNALING PATHWAY”, “TGF BETA SIGNALING PATHWAY”, and “WNT SIGNALING PATHWAY”, and so on (Additional file 13: Figure S9).

Immune subgroups and clinical traits correspond to HDAC expression in BC

We compared the levels of HDACs across the five immune subtypes in BLCA. The findings highlighted significant differences in HDAC expression levels across the subtypes (C1–C4 and C6), with the exception of HDAC4/5/6 and HDAC8 (Additional file 14: Figure S10A). Furthermore, we explored the association between HDACs expression and bladder cancer clinical traits, such as stage of tumor node metastasis (TNM), tumor grade, and categorization of tumors into non-papillary or papillary subtypes. Several HDACs showed differential expression in bladder tumor tissues based on tumor grade, TNM stage, and tumor subtype, especially HDAC1, HDAC4, HDAC10, and HDAC11 (Additional file 14: Figure S10B-G).

Correlation between HDACs and stemness index, TMB, MSI, or TME in BC

As shown in Fig. 12, most HDAC family gene expression levels in BLCA patients were positively or negatively related with the RNAss/DNAss/immuneSore/stromalScore, except for HDAC7. Furthermore, the TIMER database was utilized to investigate the relationship between the expression of HDACs and the infiltration of six major immune cell types in BLCA. Findings demonstrated a substantial correlation between HDAC expression and the variety of lymphocytes penetrating BLCA (Fig. 13). The above results illustrated that HDACs might impact the survival outcomes of individuals with cancer by influencing immune cells infiltration in the TME.

Fig. 12

Correlations between the expression of HDACs and RNAss, DNAss, StromaScore, ImmuneScore, ESTIMATEScore, TMB, and MSI in BLCA

Fig. 13

Correlations between the expression of HDACs and immune infiltration by six frequently focused immune cells in BLCA based on TIMER database

Correlations between HDAC family CNVs, mutations, and immune infiltration in BLCA

To investigate the factors contributing to alterations in HDACs expression and the immune microenvironment, we examined the association between CNVs or mutations in HDACs and the immune microenvironment using data from the TIMER database. Figure 14A shows the relationships between changes in the copy number of HDACs and six types of immune infiltration in BLCA. In particular, copy number high-amplification or arm-level gain of HDAC3 and HDAC8 was associated with substantially lower levels of immune infiltrates (CD4 + T cells, dendritic cells, neutrophils, and CD8 + T cells) in BLCA. These findings may indicate a potential mechanism by which gene alterations in the HDAC family gene may predict responses to immunotherapy. Subsequently, all BLCA samples were classified according to the mutation status of HDACs, and difference analysis was employed to investigate their potential relationship with the level of immune cell infiltration (Fig. 14B). Mutation groups of HDAC8 or HDAC5 were enriched in NK cells and CD8 + T cells; however, mutation groups of HDAC6 or HDAC11 exhibited lower levels in CD8 + T cells compared to non-mutation groups.

Fig. 14

The correlations between variations in the copy numbers of HDACs (A) or mutations (B), and the presence of immune cell infiltrations in BLCA

Validation of HDACs protein expression in BC from Human Protein Atlas (HPA)

The HPA database was used to compare HDAC protein levels between BC and normal bladder tissues. Compared with expression in adjacent tissues, cancer tissues exhibited an increase in the expression of HDAC1/3 and a decrease in HDAC4/5/9. These findings align with our previous analysis of transcriptome differences for these genes (Fig. 15A-I).

Fig. 15

Verification of protein-level expression of HDACs. (A-I) Immunohistochemical staining for HDAC1-6 and HDAC8-10 protein expression based on HPA database