Aggresome–Autophagy Associated Gene HDAC6 Is a Potential Biomarker in Pan-Cancer, Especially in Colon Adenocarcinoma

Background: Histone deacetylase 6 (HDAC6) regulates cytoplasmic signaling networks through the deacetylation of various cytoplasmic substrates. Recent studies have identified the role of HDAC6 in tumor development and immune metabolism, but its specific function remains unclear.

Methods: The current study determined the role of HDAC6 in tumor metabolism and tumor immunity through a multi-database pan-cancer analysis. The Cancer Genome Atlas (TCGA), Genotype-Tissue Expression (GTEx), and Cancer Cell Line Encyclopedia (CCLE) datasets were used to determine the expression levels, prognosis, tumor progression, immune checkpoints, and immune metabolism of HDAC6 in 33 tumors. Pathways, immune checkpoints, immune neoantigens, immune microenvironment, tumor mutational burden (TMB), microsatellite instability (MSI), DNA mismatch repair (MMR), and the value of methyltransferases. The R package was used for quantitative analysis and panoramic description.

Results: In the present study, we determined that HDAC6 is differentially expressed in pan carcinomas, and by survival, we found that HDAC6 was generally associated with the prognosis of pancreatic adenocarcinoma, Thymoma, and uveal melanoma, where low expression of HDAC6 had a significantly worse prognosis. Secondly, through this experiment, we confirmed that HDAC6 expression level was associated with tumor immune infiltration and tumor microenvironment, especially in PAAD. Finally, HDAC6 was associated with immune neoantigen and immune checkpoint gene expression profiles in all cancers in addition to TMB and MSI in pan-cancers.

Conclusion: HDAC6 is differentially expressed in pan-cancers and plays an essential role in tumor metabolism and immunity. HDAC6 holds promise as a tumor potential prognostic marker, especially in colon cancer.

Introduction

Living standards, dietary habits, and living environments of people is changing (1), and the incidence of various diseases is exploding (2). Globally, cancer has become a severe public health hazard; cancer incidence and mortality rates are increasing rapidly every year (3). Malignancy is one of the leading causes of death worldwide, and treatment success is low in developed and developing countries (4). Although breakthroughs in tumor diagnosis and tumor treatment have been made in recent years, conventional therapies, including surgery, chemotherapy, and radiation therapy, remain the first-line treatment for most cancer patients (5). However, patients with many cancers still do not achieve the desired level of 5-year survival with treatment (6). Currently, the application of tumor biomarkers has attracted much attention from scholars, among which pan-cancer analysis has emerged as a potential option to explore new tumor biomarkers (7).

Overexpression of histone deacetylases (HDACs) in cancer cells is an important cause of acetylation imbalance. A total of 18 HDACs have been identified in Ijsselsteijn mans, which can be classified into four categories according to their class namely, I (HDAC1, 2, 3, and 8), II (HDAC4, 5, 6, 7, 9, and 10), III (SIRT1, 2, 3, 4, 5, and 6) and IV (HDAC11), of which classes II and III play an essential role in the life course. (HDAC6) is the most widely studied class II histone deacetylase isoform (8), and studies have confirmed that HDAC6 plays an essential role in the liver (9), kidney (10), brain (11), and other areas of the brain (11); in pancreatic (12) tumors the expression of HDAC6 is upregulated. SIRTs, as class III histone deacetylases, play essential roles and broad cellular functions in aging. Recently, it has been shown that HDAC6 and SIRT2 act as deacetylases that regulate the acetylation status of KRAS in cancer cells. SIRT2 also promotes the differentiation and proliferation of intestinal epithelial cells by regulating Wnt-β-catenin signaling. Although HDAC6 plays a role in deacetylating histones, recent reports have found that HDAC6 can be involved in tumorigenesis and development through multiple pathways.

Cellular autophagy has emerged as a hot topic in recent years. It has been confirmed that autophagy is associated with tumor metabolic patterns and the formation of tumor heterogeneity. According to its properties, it can be divided into two aspects. On the one hand, autophagy antagonizes the inflammatory response and can inhibit the degree of infiltration of chronic inflammation, thus improving inflammatory cancer transformation. On the other hand, autophagy can act as an essential pathway of cellular energy metabolism and material recycling by degrading damaged mitochondria and assisting tumors in escaping from Reactive Oxygen Species (ROS) during aerobic glycolysis, thus ensuring the sustainability of the Warburg effect. The aggresome–autophagy pathway is a specific type of induced autophagy. In other words, misfolded and aggregated proteins are selectively recognized for reverse transportation to the center of the microtubule tissue, where they form aggresomes in the pericentriolar region and are subsequently cleared by autophagy. An increasing number of studies have identified the aggresome–autophagy pathway as a critical cellular defense system to prevent the accumulation of misfolded and aggregated proteins. Recent studies have found that HDAC6 is closely associated with aggresome–autophagy and that cells lacking HDAC6 cannot form an aggresome and are highly sensitive to the accumulation of misfolded proteins. It may be a potential mechanism for the involvement of HDAC6 in autophagy.

In addition, HDAC6 activity can also affect the gene expression of several critical immune molecules. These include programmed death receptor-1 (PD-1) and programmed death receptor ligand-1 (PD-L1), tumor-associated antigens, and these factors are central targets for cancer immunotherapy (13).

Unfortunately, the mechanism of HDAC6 in tumor metabolism and tumor immunity is not yet precise. In the present study, we used multiple datasets, such as Cancer Cell Lineage Encyclopedia (CCLE) and The Cancer Genome Atlas (TCGA) for pan-cancer analysis to reveal the precise mechanism of HDAC6. With the rise of high-throughput sequencing, histology technologies are gradually coming into the limelight. New perspectives for cancer research have been provided through histological technologies.

In this study, we analyzed HDAC6 expression, prognosis, TMB, and MSI in 33 tumors by pan-cancer analysis. In addition, we also examined the correlation of HDAC6 in an immune microenvironment, immune-related antigens, and checkpoint genes. We further confirmed that HDAC6 expression affected the expression of DNA repair genes and methyltransferases in pan-cancer. By gene set enrichment analysis, we found that HDAC6 regulated signaling pathways related to cancer apoptosis and tumor immunity.

Methods and Materials

Sample Source

Studies were based on The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) dataset (14), Genotype-Tissue Expression (GTEx, https://gtexportal.org/) dataset (15), Cancer Cell Line Encyclopedia (CCLE, https://portals.broadinstitute.org/) dataset (16) and the timer (https://cistrome.shinyapps.io/timer/) dataset (17). The TCGA (Table 1) and GTEx datasets were used to obtain gene expression information of tumor and normal samples and clinical information data and to analyze HDAC6 expression in 27 tumors. HDAC6 expression in 31 cancers in the GTEx dataset was also analyzed. To examine differential gene expression in cancers on a larger scale, RNA sequencing datasets for each cell line in the CCLE dataset were downloaded, while cancer immune infiltrating cell score data were downloaded through the timer dataset for tumor immune infiltration analysis. In this study, only open access data were used, which precluded the requirement for ethics committee approval.

TABLE 1

Table 1 TCGA Study Abbreviations in manuscript.

Analysis of HDAC6 Expression Levels and Prognosis in Pan-Cancer

In this study, the R package (edgeR) was used to analyze the differential levels of HDAC6 expression in the dataset. For the expression of HDAC6 in different tumor cells and other normal tissues, we used Kruskal–Wallis test line analysis, and the R package ggplot presented the figures. For the predictive analysis of HDAC6 in pan-cancer, a one-way Cox regression test was used to analyze the correlation between HDAC6 and patient survival. The Kaplan–Meier (K–M) test was also used to compare the survival of patients. The images were plotted using forest plot visualization software.

Gene Set Enrichment Analysis

Gene set enrichment analysis (GSEA) is a powerful analysis method for interpreting gene expression data and analyzing statistically significant and consistent differences between different groups with different biological states (18). The signaling pathway of HDAC6 was analyzed by gene set enrichment analysis. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed with the R package clusterProfiler. Another dataset used for GSEA analysis is the Molecular Signature Database (MsigDB) (19), in which the Hallmark gene set was used. Implementation criteria were | NES | >1, p-value <0.05, FDR <0.25.

Immune Checkpoint Genes and Immune Neoantigen Analysis

Immune checkpoints are a set of molecules expressed on immune cells that regulate the degree of immune activation by maintaining their normal immune function in vivo (20). We analyzed the relationship between HDAC6 expression levels and the expression of 47 common immune checkpoint genes. Correlations between HDAC6 expression levels and expression of various immune checkpoint genes in multiple cancers were assessed by the R package (limma, reshape2, RColorBrewer). The number of neoantigens in each tumor sample was counted and detected using a scanner.

Correlation Analysis of HDAC6 in Immune Infiltration and Tumor Microenvironment

To evaluate the performance of HDAC6 gene in immune infiltration, we performed correlation analysis between HDAC6 and six immune cell infiltrations with the help of purity-adjusted Spearman. In addition, we performed an ESTIMATE algorithm by standardized expression matrix to assess the tumor microenvironment-related scores of patients, including estimation score, stromal score, and immune score. p-value <0.05 was considered statistically significant.

Tumor Mutational Burden and Microsatellite Instability Analysis

The Perl language and R package (fmsb) were used to calculate the total TMB and analyze its correlation with HDAC6 expression levels in patients with pan-cancer. MSI has novel microsatellite alleles in tumors compared to normal tissue and can be used for screening, diagnosis, and prognosis of pancreatic cancer (21). The correlation between MSI and HDAC6 expression levels in pan-cancer was analyzed using the R language with the fmsb package.

DNA Mismatch Repair Gene Mutation and DNA Methylation Analysis

The correlation between MMR gene expression levels and HDAC6 expression levels was analyzed using the Spearman correlation method. In addition, DNA methylation is an essential factor affecting gene expression. In this study, we evaluated the expression levels of DNMT1, DNMT2, DNMT3A, and DNMT3B and assessed the correlation between these four methyltransferases and HDAC6 expression using the Spearman correlation method.

Results

HDAC6 Expression Was Significantly Upregulated in Pan-Cancer

In this study, to observe the expression level of HDAC6 in pan-cancer, the expression of HDAC6 in various cancer databases was first analyzed. The expression levels of HDAC6 in the GTEx (Figure 1A), CCLE (Figure 1B), and TCGA (Figure 1C) datasets were shown by a histogram. Considering the small number of normal sample data in TCGA, the expression of HDAC6 in each cancer in the TCGA and GTEx datasets was further integrated. By integration, 27 tumor modules were retained, and the final analysis showed that HDAC6 expression was downregulated in the combined 27 tumor modules (Figure 1D).

FIGURE 1

Figure 1 Expression levels of HDAC6 in pan-cancer. (A) Analysis of HDAC6 expression levels in pan-cancer based on the GTEx dataset. (B) Analysis of HDAC67nbsp;expression levels in various cancer cells based on the CCLE dataset. (C) Analysis of HDAC6 expression levels in pan-cancer based on TCGA database. (D) Integration of GTEx and TCGA databases to obtain 27 tumor modules and to analyze HDAC6 expression levels in tumor modules. * indicates P < 0.05 compared with control, ** indicates P < 0.01 compared with control, *** indicates P < 0.001 compared with control.

Prognostic Analysis of HDAC6 in Pan-Cancer

In the above study, the expression levels of HDAC6 in various tumor tissues were determined. To understand the relationship between HDAC6 and tumor prognosis, one-way Cox regression analysis was used to analyze HDAC6 expression levels and patient prognosis. Based on the TCGA database data, HDAC6 expressions were grouped into high and low expression groups according to the median value of HDAC6 expression in each tumor. It was further observed that HDAC6 was only expressed in 33 tumors in PAAD (HR = 0.960, 95CI%: 0.990–1.000, P = 0.032), THYM (HR = 0.870, 95CI%: 0.800–0.960, P = 0.005), and UVM (HR = 0.930, 95CI%: 0.870–0.990, P = 0.018) with of prognostic significance (Figure 2). In addition, the relationship between HDAC6 and PAAD, THYM, and UVM was further observed by plotting K–M survival curves (Figure 3). It suggested that there was an association between HDAC6 and tumor prognosis, especially with THYM. Considering the possible existence of non-tumor-related deaths during follow-up, we analyzed the relationship between HDAC6 expression levels and prognostic DSS (disease-specific survival) in 33 tumors of TCGA. The results of the analysis showed that HDAC6 was associated with KIRP (HR = 0.940, 95CI%: 0.900–0.990, P = 0.014), THCA (HR = 1.16, 95CI%: 1.010–1.330, P = 0.038), THYM (HR = 0.820, 95CI%: 0.700–0.960, P = 0.015), and UVM (HR = 0.870, 95CI%: 0.800–0.960, P = 0.005), indicating that HDAC6 and DSS had prognostic significance (Figure 4). Also, the K–M survival curve indicated that HDAC6 was a prognostic indicator of tumor DSS (Figure 5).

FIGURE 2

Figure 2 Analysis of the relationship between HDAC6 and OS in 33 tumors using one-way Cox regression, presented using forest plots.

FIGURE 3

Figure 3 Relationship between HDAC6 and PAAD, THYM, and UVM OS. (A) The relationship between HDAC6 high and low expression levels and PAAD OS was analyzed using K–M survival. (B) Relationship between HDAC6 high and low expression levels and THYM OS using K–M survival analysis. (C) The relationship between HDAC6 high and low expression levels and UVM OS using K–M survival analysis.

FIGURE 4

Figure 4 Analysis of the relationship between HDAC6 and DSS in 33 tumors using one-way Cox regression, presented using forest plots.

FIGURE 5

Figure 5 Relationship between HDAC6 and KIRP, THCA, THYM, UVM DSS. (A) The relationship between HDAC6 high and low expression levels and KIRP DSS was analyzed by K–M survival. (B) The relationship between HDAC6 high and low expression levels and THCA DSS was analyzed using K–M survival. (C) Relationship between HDAC6 high and low expression levels and THYM DSS using K–M survival analysis. (D) Relationship between HDAC6 high and low expression levels and UVM DSS using K–M survival analysis.

Tumor Metabolic and Immune Signaling Pathways Involved in HDAC6

To further understand the relationship between HDAC6 involvement in cancer metabolism and tumor immunity, a GSEA analysis was performed to analyze the signaling enrichment of KEGG and markers in both groups according to HDAC6 gene expression, which was divided into high and low groups. The results were ranked according to NES scores for the top ten most enriched signaling pathways or biological processes (Tables S1, S2). In addition, we mapped the top three most abundant signaling pathways in both databases. Among them, we found that the Notch signaling pathway, Wnt/β-catenin signaling pathway, and Hedgehog signaling pathway were the most enriched (Figures 6A, B).

FIGURE 6

Figure 6 Signaling enrichment of HDAC6 in KEGG and markers. (A) GSEA analysis of the top three correlations of HDAC6 with signaling pathways in the KEGG database. (B) GSEA analysis of HDAC6 with the top three results of signaling pathway correlations in the marker dataset. (C) String constructs HDAC6 protein co-expression network.

In addition, we constructed an HDAC6 PPI (protein–protein interaction) network and found that HDAC6 was related to eight proteins (Figure 6C).

Correlation Between HDAC6 Expression Levels and Immune Checkpoints and Immune Neoantigens in Pan-Cancer

Immune checkpoints are a series of molecules expressed in immune cells that regulate the degree of immune activation and play an essential role in regulating autoimmunity (22). In contrast, tumor neoantigens are nascent antigens encoded by mutated genes in tumor cells. The synthesis of neoantigen vaccines is facilitated by exploiting the immune activity of tumor neoantigens (23). This time, to explore the relationship between HDAC6 and immune regulation, it was preferred to analyze the correlation between HDAC6 and immune checkpoints and immune neoantigens. Our results found that HDAC6 expression levels in various tumors correlated with more than 40 checkpoints, with a positive correlation mainly with HNSC, THCA (Figure 7A). It suggested that HDAC6 had a role in modulating these immune checkpoints in some tumors and may have the ability to strengthen immunity. In addition, in the present study, we calculated the number of neoantigens for each tumor type. The results showed that HDAC6 was only correlated in KIPR (R = 0.173, P = 0.027), and CESE (R = −0.163, P = 0.025) (Figure 7B).

FIGURE 7

Figure 7 Correlation analysis of HDAC6 expression levels with immune checkpoints and immune neoantigens. (A) Correlation analysis of HDAC6 expression levels with immune checkpoint gene expression levels in various tumors. (B) Correlation analysis of HDAC6 expression levels in various tumors with the number of neoantigens in tumors. * indicates correlation difference of P < 0.05, ** indicates correlation difference of P < 0.01, *** indicates correlation difference of P < 0.001.

Correlation Between HDAC6 and Tumor Immune Infiltration and the Tumor Microenvironment in Pan-Cancer

Differences in the degree of infiltration of different immune cells are highly correlated with tumor progression and prognosis (24). In contrast, the tumor microenvironment consists of various cells, external matrix, and associated factors that significantly influence tumor diagnosis, survival outcome, and clinical treatment sensitivity (25). In the present study, we analyzed immune infiltration and tumor immune microenvironment of HDAC6 in different cancers. A total of 3 most relevant tumors (HNSC, KIRC, PAAD) were identified after analyzing 33 tumors in the TCGA database. Among these tumors, HDAC6 was positively correlated with B cells (R = 0.286, P < 0.001), CD4+ T cells (R = 0.424, P < 0.001), neutrophils (R = 0.293, P < 0.001), macrophages (R = 0.232, P < 0.001), and dendritic cells (R = 0.201, P < 0.001) in HNSC. In contrast, HDAC6 was positively correlated only with CD4+ T cells in KIRC (R = 0.231, P < 0.001). In addition, HDAC6 was positively correlated with CD4+ T cells (R = 0.211, P = 0.005), CD8+ T cells (R = 0.317, P < 0.001), neutrophils (R = 0.283, P < 0.001), macrophages (R = 0.414, P < 0.001) and dendritic cells (R = 0.351, P < 0.001) in PAAD (Figure 8A). In addition to analyzing the relationship between HDAC6 and tumor immune infiltration, the present study also analyzed the correlation between HDAC6 and tumor microenvironment in pan-cancer. Individual tumor samples were analyzed by R package (ESTIMATE), and the relationship between HDAC6 expression levels and immune score, stromal score, and the immune score of ESTIMATE was analyzed separately. We showed the top three tumors in which LGG (R = −277, p < 0.001), SARC (R = −0.359, p < 0.001), BGM (R = −0.412, p < 0.001) immune scores were negatively correlated with HDAC6 expression levels in 33 tumors. UCEC (R = −0.176, P < 0.001), THCA (R = −0.13, P = 0.003), BGM (R = −0.412, P < 0.001) stromal scores were negatively correlated with HDAC6 expression levels. In addition, the immune scores of SARC (R = −0.359, P < 0.001), UCEC (R = −0.176, P < 0.001), BGM (R = −0.412, P < 0.001) ESTIMATE were found to be negatively correlated with HDAC6 expression levels (Figure 8B). The above results suggested that HDAC6 was negatively correlated with both tumor immune scores.

FIGURE 8

Figure 8 correlation between HDAC6 and tumor immune infiltration and tumor microenvironment. (A) Correlation analysis between HDAC6 and immune cell infiltration in HNSC, KIRC, PAAD. (B) Relationship between HDAC6 expression levels and immune scores, stromal scores, and immune scores of ESTIMATE.

TMB and MSI Analyses of HDAC6 Expression Levels in Pan-Cancer

TMB has a high value in predicting the efficacy and prognosis of immune checkpoint therapy to assess the total number of substitutions and insertion/deletion mutations per megabase in the exon coding regions of genes in tumor samples (26). MSI is caused by defects in mismatch repair (MMR) genes and is strongly associated with tumorigenesis (21). In the present study, we analyzed the relationship between HDAC6 in TMB and MSI. Analysis by Spearman test revealed that HDAC6 expression levels were positively correlated with BLCA, COAD, and OV, and negatively correlated with BRCA, HNSC, PRAD, and THCA, respectively (Figure 9A). In addition, Spearman test analysis found HDAC6 expression levels in tumors given in the MSI. HDAC6 expression levels were positively correlated with BLCA, CESC, COAD, LGG, LUAD, LUSC, and SARC, respectively, and negatively correlated with DLBC by Spearman test analysis (Figure 9B).

FIGURE 9

Figure 9 TMB and MSI analysis of HDAC6 expression level vs. (A) Spearman test to analyze the relationship between HDAC6 and TMB in pan-cancer. (B) Spearman test to analyze the relationship between HDAC6 and MSI in pan-cancer.

Relationship Between HDAC6 Expression Level and MMR Gene and Methyltransferase in Pan Carcinoma

DNA methylation has a role in altering the structure and stability of DNA, and the appearance of MMR leads to mutations in somatic cells that are more deleterious (27). DNA methyltransferases are involved in the occurrence of methylation and play an essential regulatory role in the methylation process (28). In the present study, the relationship between HDAC6 and MMR genes (MLH1, MSH2, MSH6, PMS2, EPCAM) in various tumors was evaluated according to the TCGA database,. The analysis showed that HDAC6 expression levels were positively correlated with MMR genes in all 33 tumors except READ and UCS (Figure 10A). In addition, further investigation revealed that HDAC6 expression levels were significantly and positively correlated with both methyltransferase expression levels (Figure 10B). These results suggest that HDAC6 is involved in tumorigenesis development and has epigenetic properties that regulate various tumors.

FIGURE 10

Figure 10 Correlation of HDAC6 expression level with MMR gene and methyltransferase. (A) Spearman test to analyze the correlation between HDAC6 expression level and MMR gene in various tumors. (B) Spearman test to analyze the correlation between HDAC6 expression level and methyltransferase in various tumors. * Indicates that EPCAM, PMS2, MSH6, MSH2, MLH1 and HDAC6 are correlated in tumors, P<0.05. **Indicating that EPCAM, PMS2, MSH6, MSH2, MLH1 and HDAC6 are correlated in tumors, P< 0.01. *** means EPCAM, PMS2, MSH6, MSH2, MLH1 and HDAC6 are correlated in tumors, P<0.001.

Discussion

The number of cancer cases worldwide has exceeded 90 million, and cancer is currently a more prominent public health problem (3). A pan-cancer analysis mines genomic similarities and differences in tumors through major databases to provide better insights about cancer diagnosis and treatment to guide cancer-related research (29).

HDAC6 is a unique member of the HDAC family located on chromosome Xp11.23 and is the giant protein molecule in the HDAC family (30). The primary role of HDAC6 is to acetylate and deacetylate histones, which are modified by histone acetyltransferase (HAT) and histone deacetylase (HDAC). HAT has a role in promoting chromosome depolymerization and activating transcription, while HDAC has a role in blocking DNA and inhibiting transcription. Loss-of-nest apoptosis is a specific type of programmed cell death that results from the detachment of cells from the extracellular matrix and surrounding basement membrane. Recent studies have found that HDAC6 promotes oncogenic transformation and tumor formation by promoting anchorage-independent proliferation in transduced cells. In addition to histone regulation, HDAC6 also affects tumor cell motility by regulating non-histone substrates. α-microtubulin was the first non-histone substrate identified for HDAC6, and the reversible deacetylation of α-microtubulin by HDAC6 affects microtubulin. The reversible deacetylation of α-microtubulin by HDAC6 affects microtubule stability and function. It has been found that HDAC6 is highly expressed in melanoma, and by knocking down HDAC6 acetylated α-microtubulin increases, acetylated microtubules accumulate, causing CYLD to get around the nucleus, ultimately leading to a reduction in the interaction between CYLD and BCL3, preventing the transcriptional activity of the nuclear factor NF-κB, which in turn affects cell growth and metastasis.

Immunomodulation is of interest as a new idea in cancer therapy, especially in tumors lacking specific molecular targets. TME comprises pro and anti-cancer immune cells containing CD8-effective T cells, natural killer (NK) cells, macrophages, regulatory T cells, and myeloid-derived suppressor cells. Studies have shown a strong link between HDAC6 and the immune microenvironment. For example, inhibition of HDAC6 enhances antitumor immune signaling and reduces tumor load in ovarian cancer (31). HDAC6 has a regulatory antitumor immune response in breast cancer and plays an atypical role in disseminated and invasive breast cancer (32). Although anti-programmed death-1 (PD-1)/programmed death ligand-1 (PD-L1) drugs have achieved considerable clinical efficacy and low toxicity, they are not effective in all cancer types or do not achieve the desired effect in all patients (33). In a study, a novel HDAC6 inhibitor (MPT0G612), which induces apoptosis and inhibits IFN-γ-induced programmed death-ligand 1 in human colorectal cancer cells, was found to be a potential strategy for the combination of immune checkpoint inhibitors in the treatment of CRC (34). It suggests that HDAC6 holds promise as a new drug target for anti-cancer immunotherapy or in combination with known immune checkpoint inhibitors to enhance immune infiltration and response to cancer. Therefore, we analyzed HDAC6 expression in pan-cancer, prognosis, immune microenvironment, immune-associated antigens, and checkpoint genes for correlation analysis.

In this study, we first analyzed the expression and prognosis of HDAC6 in pan-cancer. After analyzing GTEx, CCLE, and TCGA databases, we found that HDAC6 was differentially expressed in the TCGA database for pan-cancers. However, considering the small number of control samples in the TCGA database, we integrated samples from the TCGA and GTEx databases. Here the results were reversed, and the expression of HDAC6 in the integrated database for pan-cancers was low, contrary to recent findings (31, 35). We speculate that this may be related to the recent addition of sample data to the TCGA database. Also in our prognostic analysis, we found that patients with low HDAC6 expression in PAAD, THYM, and UVM had significantly lower OS according to HDAC6 expression level which was further divided into high and low expression groups, and patients with low HDAC6 expression in KIRP, THYM, and UVM had significantly lower DSS. These results indicate that HDAC6 is expected to be a prognostic indicator in tumors as mentioned earlier, but further clinical validation is needed.

Cancer metabolism and tumor immune signaling pathways are important areas in basic tumor research (34, 36). To better observe the mechanism of HDAC6 in tumors, the analysis by GSEA revealed that HDAC6 is associated with the Notch signaling pathway (37), Wnt/β-catenin signaling pathway (38), and Hedgehog signaling pathway (39) were the most enriched. As the three most relevant pathways to tumor development, single-gene enrichment allowed identifying HDAC6 as involved in the feedback of these pathways and involved in tumorigenesis based on these pathways. In addition, by constructing a PPI protein co-expression network, eight genes were co-expressed with HDAC6, among which HDAC11 was of interest. HDAC11 belongs to class IV histone deacetylation, and earlier studies found that HDAC6 and HDAC11 act as a common transcriptional activator to regulate IL-10 expression, suggesting a possible regulatory role of HDAC6 in cellular immunity.

We found HDAC6 expression correlated with more than 40 checkpoints in the follow-up study, which were not lacking essential immune genes such as PD-1, TIGIT, TNFRSF9, and CTLA4. These results strongly suggest that HDAC6 may be a potential biomarker and play a crucial role in tumor immunity. The immune microenvironment is an essential component in tumor development, and tumor growth and metastasis can be altered by regulating the infiltration of immune cells and immune modifications and immune escape in the tumor microenvironment (40). In addition, the present study also found that HDAC6 was involved in the immune infiltration of tumors, and there was a positive correlation between the expression level of HDAC6 and immune score. They suggested that HDAC6 may be a potential immunosuppressant.

TMB and MSI scores have a high value in predicting the efficacy and prognosis of immune checkpoint therapy (41). With the continuous improvement of immune checkpoints in recent years, immunosuppressants are increasingly used in tumor immunotherapy and have played an influential role in clinical practice (42). In addition to this, PD-1/PD-L1 immunosuppressants have become a first- or second-line treatment option for some tumors (43). In our study, HDAC6 expression levels were shown to correlate with TMB and MSI in four cancer types, respectively. Among them, BLCA and COAD were positively correlated in TMB and MSI with overlapping relationships, and COAD was more significantly associated with TMB and MSI. In recent years, studies on COAD and HDAC6 have been increasing, and Yang et al. (12) found that HDAC6 in colon cancer cells regulates the deacetylation of K-RAS genes in the acetylated state and through the RAS/MAPK signaling pathway. KRAS-activating mutations are present in more than 40% of patients with colon cancer. The recent finding that KRAS mutations in LUAD may serve as a potential predictor for guiding anti-PD-1/PD-L1 immunotherapy needs to be further validated in colon cancer.

We concluded the study by analyzing the relationship between HDAC6 expression levels and the MMR gene and methyltransferase. Under normal circumstances, mismatches may occur during DNA replication, leading to genetic mutations, but regulation by the MMR genome in cells can identify and correct mutations (44). Moreover, mutations in the MMR genome can exacerbate the increased accumulation of genetic errors causing genomic instability or MSI (45). DNA methylation is an essential factor leading to altered tumor development, and an increasing number of studies show that hypermethylation of gene promoters is a common epigenetic feature of cancer (46). In the present study, we found that HDAC6 expression was strongly associated with mutation levels in five MMR genes in human pan-cancers, but there were some exceptions in tumors (READ and UCS). These results are consistent with our conclusion that high HDAC6 expression plays an essential regulatory role in tumor development by regulating MMR gene mutations through DNA methylation.

In the present study, we identified the potential value of HDAC6 in pan-cancer by primary analysis. However, limitations of the study still exist. First, as a primary analysis, we did not perform validation analysis on clinical samples or animal models, such as HDAC6 expression in multiple cancers. The results in the current study differed from recent studies. Second, although we determined that HDAC6 expression levels were associated with tumor immune cell infiltration and patient survival, we could not directly demonstrate that HDAC6 affects patient survival through immune infiltration. Finally, the current study used sequencing data from multiple databases and microarray data for analysis, which are subject to systematic bias. Therefore, we hope to collect clinical samples in future studies and conduct primary research to refine our findings.

Despite the flaws of the current study, we have to acknowledge the close association of HDAC6 in tumor immunity and cancer development. In the present study, we determined that HDAC6 was differentially expressed in pan-cancers and that abnormal expression was associated with tumor progression, especially in COAD. Abnormal expression of HDAC6 was associated with immune checkpoints, immune cell infiltration, tumor microenvironment, TMB, MSI, and DNA methylation. HDAC6 holds promise as a tumor potential prognostic marker, especially in colon cancer.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: https://portal.gdc.cancer.gov/, https://gtexportal.org/, https://portals.broadinstitute.org/, https://cistrome.shinyapps.io/timer/.

Author Contributions

ZZ: conception and design, writing the article, critical revision of the article. XZ: data collection, writing the article, critical revision of the article, analysis and interpretation. AH: conception and design, writing the article, critical revision of the article. All authors contributed to the article and approved the submitted version.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors acknowledge the sangerbox tool developed by Hangzhou Mugu science and technology Co., Ltd. The authors thank them for their contribution to this study.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2021.718589/full#supplementary-material

References

1. Budnik LT, Casteleyn L, Paschalidou AK, Kassomenos P. Pollution in Living and Working Environments, Climate Variability, and Their Impact on non-Communicable Disease Burden. Sci Total Environ (2019) 660:593–4. doi: 10.1016/j.scitotenv.2019.01.055

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Jemal A, Center MM, DeSantis C, Ward EM. Global Patterns of Cancer Incidence and Mortality Rates and Trends. Cancer Epidemiol Biomarkers Prev (2010) 19(8):1893–907. doi: 10.1158/1055-9965.EPI-10-0437

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global Cancer Statistics 2018: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin (2018) 68(6):394–424. doi: 10.3322/caac.21492

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Global Burden of Disease Cancer C, Fitzmaurice C, Abate D, Abbasi N, Abbastabar H, Abd-Allah F, et al. Global, Regional, and National Cancer Incidence, Mortality, Years of Life Lost, Years Lived With Disability, and Disability-Adjusted Life-Years for 29 Cancer Groups, 1990 to 2017: A Systematic Analysis for the Global Burden of Disease Study. JAMA Oncol (2019) 5(12):1749–68. doi: 10.1001/jamaoncol.2016.5688

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Seymour L, Bogaerts J, Perrone A, Ford R, Schwartz LH, Mandrekar S, et al. Irecist: Guidelines for Response Criteria for Use in Trials Testing Immunotherapeutics. Lancet Oncol (2017) 18(3):e143–52. doi: 10.1016/S1470-2045(17)30074-8

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Straker N. The Treatment of Cancer Patients Who Die. Psychodyn Psychiatry (2020) 48(1):1–25. doi: 10.1521/pdps.2020.48.1.1

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Gentles AJ, Newman AM, Liu CL, Bratman SV, Feng W, Kim D, et al. The Prognostic Landscape of Genes and Infiltrating Immune Cells Across Human Cancers. Nat Med (2015) 21(8):938–45. doi: 10.1038/nm.3909

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Miyake Y, Keusch JJ, Wang L, Saito M, Hess D, Wang X, et al. Structural Insights Into HDAC6 Tubulin Deacetylation and Its Selective Inhibition. Nat Chem Biol (2016) 12(9):748–54. doi: 10.1038/nchembio.2140

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Lv Z, Weng X, Du C, Zhang C, Xiao H, Cai X, et al. Downregulation of HDAC6 Promotes Angiogenesis in Hepatocellular Carcinoma Cells and Predicts Poor Prognosis in Liver Transplantation Patients. Mol Carcinog (2016) 55(5):1024–33. doi: 10.1002/mc.22345

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Ding Y, Kang A, Tang Q, Zhao Y. Inhibition of HDAC6 Expression Decreases Brain Injury Induced by APOE4 and Abeta Coaggregation in Rats. Mol Med Rep (2019) 20(4):3363–70. doi: 10.3892/mmr.2019.10571

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Zhang Z, Cao Y, Zhao W, Guo L, Liu W. HDAC6 Serves as a Biomarker for the Prognosis of Patients With Renal Cell Carcinoma. Cancer Biomark (2017) 19(2):169–75. doi: 10.3233/CBM-160298

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Yang Y, Sun Y, Wang H, Li H, Zhang M, Zhou L, et al. Microrna-221 Induces Autophagy Through Suppressing HDAC6 Expression and Promoting Apoptosis in Pancreatic Cancer. Oncol Lett (2018) 16(6):7295–301. doi: 10.3892/ol.2018.9513

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Knox T, Sahakian E, Banik D, Hadley M, Palmer E, Noonepalle S, et al. Selective HDAC6 Inhibitors Improve Anti-PD-1 Immune Checkpoint Blockade Therapy by Decreasing the Anti-Inflammatory Phenotype of Macrophages and Down-Regulation of Immunosuppressive Proteins in Tumor Cells. Sci Rep (2019) 9(1):6136. doi: 10.1038/s41598-019-42237-3

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Blum A, Wang P, Zenklusen JC. Snapshot: TCGA-Analyzed Tumors. Cell (2018) 173(2):530. doi: 10.1016/j.cell.2018.03.059

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Consortium GT. The Genotype-Tissue Expression (Gtex) Project. Nat Genet (2013) 45(6):580–5. doi: 10.1038/ng.2653

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Nusinow DP, Szpyt J, Ghandi M, Rose CM, McDonald ER 3rd, Kalocsay M, et al. Quantitative Proteomics of the Cancer Cell Line Encyclopedia. Cell (2020) 180(2):387–402.e16. doi: 10.1016/j.cell.2019.12.023

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, et al. Liberzon. Cancer Res (2017) 77(21):e108–10. doi: 10.1158/0008-5472.CAN-17-0307

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Powers RK, Goodspeed A, Pielke-Lombardo H, Tan AC, Costello JC. GSEA-Incontext: Identifying Novel and Common Patterns in Expression Experiments. Bioinformatics (2018) 34(13):i555–i64. doi: 10.1093/bioinformatics/bty271

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Liberzon A, Subramanian A, Pinchback R, Thorvaldsdottir H, Tamayo P, Mesirov JP. Molecular Signatures Database (Msigdb) 3.0. Bioinformatics (2011) 27(12):1739–40. doi: 10.1093/bioinformatics/btr260

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Danilova L, Ho WJ, Zhu Q, Vithayathil T, De Jesus-Acosta A, Azad NS, et al. Programmed Cell Death Ligand-1 (PD-L1) and CD8 Expression Profiling Identify an Immunologic Subtype of Pancreatic Ductal Adenocarcinomas With Favorable Survival. Cancer Immunol Res (2019) 7(6):886–95. doi: 10.1158/2326-6066.CIR-18-0822

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Schrock AB, Ouyang C, Sandhu J, Sokol E, Jin D, Ross JS, et al. Tumor Mutational Burden is Predictive of Response to Immune Checkpoint Inhibitors in MSI-High Metastatic Colorectal Cancer. Ann Oncol (2019) 30(7):1096–103. doi: 10.1093/annonc/mdz134

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Zhou B, Gao S. Pan-Cancer Analysis of FURIN as a Potential Prognostic and Immunological Biomarker. Front Mol Biosci (2021) 8:648402. doi: 10.3389/fmolb.2021.648402

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Cerezo-Wallis D, Soengas MS. Understanding Tumor-Antigen Presentation in the New Era of Cancer Immunotherapy. Curr Pharm Des (2016) 22(41):6234–50. doi: 10.2174/1381612822666160826111041

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Sokratous G, Polyzoidis S, Ashkan K. Immune Infiltration of Tumor Microenvironment Following Immunotherapy for Glioblastoma Multiforme. Hum Vaccin Immunother (2017) 13(11):2575–82. doi: 10.1080/21645515.2017.1303582

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Gajewski TF, Schreiber H, Fu YX. Innate and Adaptive Immune Cells in the Tumor Microenvironment. Nat Immunol (2013) 14(10):1014–22. doi: 10.1038/ni.2703

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Merino DM, McShane LM, Fabrizio D, Funari V, Chen SJ, White JR, et al. Establishing Guidelines to Harmonize Tumor Mutational Burden (TMB): In Silico Assessment of Variation in TMB Quantification Across Diagnostic Platforms: Phase I of the Friends of Cancer Research TMB Harmonization Project. J Immunother Cancer (2020) 8(1):e000147. doi: 10.1136/jitc-2019-000147

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Baretti M, Le DT. DNA Mismatch Repair in Cancer. Pharmacol Ther (2018) 189:45–62. doi: 10.1016/j.pharmthera.2018.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Bennett MR, Shepherd SA, Cronin VA, Micklefield J. Recent Advances in Methyltransferase Biocatalysis. Curr Opin Chem Biol (2017) 37:97–106. doi: 10.1016/j.cbpa.2017.01.020

PubMed Abstract | CrossRef Full Text | Google Scholar

29. He D, Zhang X, Tu J. Diagnostic Significance and Carcinogenic Mechanism of Pan-Cancer Gene POU5F1 in Liver Hepatocellular Carcinoma. Cancer Med (2020) 9(23):8782–800. doi: 10.1002/cam4.3486

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Verdel A, Khochbin S. Identification of a New Family of Higher Eukaryotic Histone Deacetylases. Coordinate Expression of Differentiation-Dependent Chromatin Modifiers. J Biol Chem (1999) 274(4):2440–5. doi: 10.1074/jbc.274.4.2440

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Moufarrij S, Srivastava A, Gomez S, Hadley M, Palmer E, Austin PT, et al. Combining DNMT and HDAC6 Inhibitors Increases Anti-Tumor Immune Signaling and Decreases Tumor Burden in Ovarian Cancer. Sci Rep (2020) 10(1):3470. doi: 10.1038/s41598-020-60409-4

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Banik D, Noonepalle S, Hadley M, Palmer E, Gracia-Hernandez M, Zevallos-Delgado C, et al. HDAC6 Plays a Noncanonical Role in the Regulation of Antitumor Immune Responses, Dissemination, and Invasiveness of Breast Cancer. Cancer Res (2020) 80(17):3649–62. doi: 10.1158/0008-5472.CAN-19-3738

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Chen G, Huang AC, Zhang W, Zhang G, Wu M, Xu W, et al. Exosomal PD-L1 Contributes to Immunosuppression and is Associated With Anti-PD-1 Response. Nature (2018) 560(7718):382–6. doi: 10.1038/s41586-018-0392-8

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Chen L, Zhang YH, Lu G, Huang T, Cai YD. Analysis of Cancer-Related Lncrnas Using Gene Ontology and KEGG Pathways. Artif Intell Med (2017) 76:27–36. doi: 10.1016/j.artmed.2017.02.001

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Sharif T, Martell E, Dai C, Ghassemi-Rad MS, Hanes MR, Murphy PJ, et al. HDAC6 Differentially Regulates Autophagy in Stem-Like Versus Differentiated Cancer Cells. Autophagy (2019) 15(4):686–706. doi: 10.1080/15548627.2018.1548547

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Kim I, Choi S, Kim S. BRCA-Pathway: A Structural Integration and Visualization System of TCGA Breast Cancer Data on KEGG Pathways. BMC Bioinf (2018) 19(Suppl 1):42. doi: 10.1186/s12859-018-2016-6

CrossRef Full Text | Google Scholar

37. Deskin B, Lasky J, Zhuang Y, Shan B. Requirement of HDAC6 for Activation of Notch1 by TGF-Beta1. Sci Rep (2016) 6:31086. doi: 10.1038/srep31086

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Yin Z, Xu W, Xu H, Zheng J, Gu Y. Overexpression of HDAC6 Suppresses Tumor Cell Proliferation and Metastasis by Inhibition of the Canonical Wnt/Beta-Catenin Signaling Pathway in Hepatocellular Carcinoma. Oncol Lett (2018) 16(6):7082–90. doi: 10.3892/ol.2018.9504

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Steg AD, Burke MR, Amm HM, Katre AA, Dobbin ZC, Jeong DH, et al. Proteasome Inhibition Reverses Hedgehog Inhibitor and Taxane Resistance in Ovarian Cancer. Oncotarget (2014) 5(16):7065–80. doi: 10.18632/oncotarget.2295

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Lei X, Lei Y, Li JK, Du WX, Li RG, Yang J, et al. Immune Cells Within the Tumor Microenvironment: Biological Functions and Roles in Cancer Immunotherapy. Cancer Lett (2020) 470:126–33. doi: 10.1016/j.canlet.2019.11.009

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Picard E, Verschoor CP, Ma GW, Pawelec G. Relationships Between Immune Landscapes, Genetic Subtypes and Responses to Immunotherapy in Colorectal Cancer. Front Immunol (2020) 11:369. doi: 10.3389/fimmu.2020.00369

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Jasiak NM, Park JM. Immunosuppression in Solid-Organ Transplantation: Essentials and Practical Tips. Crit Care Nurs Q (2016) 39(3):227–40. doi: 10.1097/CNQ.0000000000000117

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Wang Z, Wu X. Study and Analysis of Antitumor Resistance Mechanism of PD1/PD-L1 Immune Checkpoint Blocker. Cancer Med (2020) 9(21):8086–121. doi: 10.1002/cam4.3410

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Ijsselsteijn R, Jansen JG, de Wind N. DNA Mismatch Repair-Dependent DNA Damage Responses and Cancer. DNA Repair (Amst) (2020) 93:102923. doi: 10.1016/j.dnarep.2020.102923

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Caja F, Vodickova L, Kral J, Vymetalkova V, Naccarati A, Vodicka P. DNA Mismatch Repair Gene Variants in Sporadic Solid Cancers. Int J Mol Sci (2020) 21(15):5561. doi: 10.3390/ijms21155561

CrossRef Full Text | Google Scholar

46. Gupta D, Heinen CD. The Mismatch Repair-Dependent DNA Damage Response: Mechanisms and Implications. DNA Repair (Amst) (2019) 78:60–9. doi: 10.1016/j.dnarep.2019.03.009

PubMed Abstract | CrossRef Full Text | Google Scholar