Differential Expression of Zinc-Dependent HDAC Subtypes and their Involvement in Unique Pathways Associated with Carcinogenesis.

OBJECTIVE: The present study aims to identify the effect of ZnHDACs expression on the survival of the patients. Further, reveal the unique and common genes associated with each ZnHDACs and their associated pathways.
METHODS: The patient data was obtained from the Cancer Genome Atlas Program (TCGA) database and was analyzed using cBioportal and Gene Expression Profiling Interactive Analysis 2(GEPIA2) online tools. Protein-protein interactions and functional interactomic analysis were done using STRING, DAVID, and KEGG pathway databases.
RESULTS: HDAC1, 2, 8, 11 were over-expressed and, HDAC4, 5, 6, 7, and 10 were down-regulated in all the cancer types, but there are few exceptional expression patterns such as HDAC7 and HDAC10 overexpression in HNSC, HDAC3 down-regulation in LUAD, and PRAD. The unique genes interacting with each ZnHDACs provided a better understanding of ZnHDAC's putative role in carcinogenesis. The present study reported that JARID2, stem cell regulation gene uniquely interacts with HDAC1, BPTF-CHRAC-BAZIA axis, enzymes for chromatin modeling selectively interacting with only HDAC2, HDAC3 in H2A acetylation via DMAP1 and YEATS4. HDAC6 associated unique genes regulate protein stability, HDAC7 in subnuclear localization and splicing, HDAC8 in telomere maintenance, HDAC9 in chromosomal rearrangements, and HDAC11 in maintaining histone core and folding.
CONCLUSION: The unique genes and pathways associated with a particular ZnHDACs could provide a wide window for interrogating these genes for obtaining putative drug targets.

Introduction

According to previous data and predicted statistics for cancer incidence and mortality, the most common cancer types worldwide are prostate, lung and Bronchus, Breast and head and neck cancers (Siegel et al., 2019). These cancers are studied extensively from a genomics, transcriptomics, and proteomics perspective. The data from thousands of patient tumors and blood samples have been deposited into databases like The Cancer Genome Atlas Program (TCGA). TCGA is a cancer-specific database rich in genetic, proteomic, histologic, and clinical data (Sun et al., 2018). Exploring these datasets may provide promising avenues for cancer research, several tools for systematic analysis of these big datasets such as cBioportal (Cerami et al., 2012), Genomic data portal (Grossman et al., 2016), Xena browser (Goldman et al., 2018), MEXPRESS (Koch et al. 2015), and Gene Expression Profiling Interactive Analysis 2 (GEPIA2) (Tang et al., 2019) are available. GEPIA2, a recently launched server that is helpful in data mining by filling the gap between big data and its integrated information for cancer research. The target identification for designing new drugs for chemotherapy and immune therapy is the growing field for cancer research. HDAC is one of the promising therapeutic targets. HDAC inhibitors are potential chemotherapeutic drugs, a few pan HDACi’s are already approved, for example, Vorinostat, Panobinostat, Chidamide (Kavanaugh et al., 2010; Bandolik et al., 2019). HDAC isoforms have four classes with eighteen subtypes. Class I, II, and IV HDACs are zinc-dependent whereas, Class III, also called sirtuins, are NAD+ dependent and have seven subclasses.(Lane and Chabner, 2009). HDAC 1, 2, 3, and 8 are four isoforms belonging to class I further, class II is sub-divided into class IIa (HDAC 4, 5, 7, 9) and Class IIb (HDAC6, 10). Class IV has only one member, HDAC11 (Ho et al., 2020). The ZnHDAC subtypes are already known for their significant role in cell cycle regulation, apoptosis, and other carcinogenic processes (Ropero and Esteller, 2007). Therefore, the ZnHDACs have shown great potential in anti-cancer drug discovery. Despite various merits in using HDAC inhibitors, researchers worldwide are encountering a challenge in the global use of pan HDACi for cancer therapeutics because these cause serious side effects. Literature reports several studies with superiority of subtype-selective HDACi (SSI) over Pan HDACi (PHI) such as HDAC1 SSI (Etinostat) over PHI (Panobinostat) in ovarian cancer (Bandolik et al., 2019). To develop ZnHDAC SSI’s various efforts are made using different approaches such as high throughput screening and subtype-specific filters designing computationally (Bertrand, 2010; Ukey et al., 2021), chemical synthesis of inhibitor molecules and their bioassays to validate activities (Ganai, 2018), studies at the molecular level such as mutagenesis (Dowling et al., 2008), and identification of active site properties by various experiments (Dowling et al., 2008; Ganai, 2018; Bandolik et al., 2019). Some structure-based studies report that despite sequence similarity in residues that make up the catalytic pocket of the binding sites, HDAC subtypes exhibit their unique structural features (Hai and Christianson, 2016).

HDAC subtype selectivity and specificity of inhibitors remained the concerned problem for decades. Big-data analyses and interactomics are emerging fields in drug target identification, prognostic and diagnostic biomarker identification, and a step towards personalized medicine.

The present study attempted to analyze the differential expression patterns of each ZnHDACs, their correlation to overall and disease-free survival in various cancer types. Further, protein-protein interaction, interatomic and functional enrichment analyses were done to obtain the genes and pathways associated with these ZnHDACs. The present study aims to identify the effect of ZnHDACs expression on the survival of the patients. Further, reveal the unique and common genes associated with each ZnHDACs and their associated pathways.

Materials and Methods

ZnHDACs expression analysis

In the present study, we used GEPIA2 and cBioportal online analysis tools to investigate the mRNA expression distribution of different ZnHDACs in the most common cancer types from the data uploaded on the public databases. GEPIA2 and cBioportal works on the UCSC Xena project (http://xena.ucsc.edu) and TCGA. TCGA database contains genome-wide expression profiling data from various cancer types for almost 9,736 tumors and 8,587 control samples. The data for prostate adenoma (PRAD), lung squamous cell carcinoma (LUSC), lung adenoma (LUAD), breast invasive carcinoma (BRCA), colon adenocarcinoma (COAD), rectum adenocarcinoma (READ), and head and neck squamous cell carcinoma (HNSC) used in the present study. Both male and female patients with ages ranging from 16 to 85 years were included. The number of samples for each cancer is reported in Table 1. We also used the Metabolic-gene Rapid Visualizer (MERAV), comprising human gene expression data from normal tissues, primary tumors, and cancer cell lines, to further validate the expression distribution of various ZnHDACs.

Table 1

The Number of Tumor and Normal Samples Available in Each Datasets

TCGA Datasets	Tumor	Normal
PRAD	492	52
LUSC	486	50
LUAD	483	59
BRCA	1085	112
COAD	275	41
READ	92	10
HNSC	519	44

Comprehensive survival analysis of patients against ZnHDACs expression in various cancer

Overall survival and disease-free survival analysis were done using the GEPIA2 functions under the head of the “Expression Analysis” survival analyses panel. GAPDH normalized ZnHDACs expression was plotted against survival for each cancer dataset. Overall and disease-free survival are plotted individually with median group cutoff and 95% confidence interval. The hazards ratio (HR) was calculated based on the Cox PH model.

Gene set enrichment analysis

The interacting genes/ proteins for each ZnHDACs were listed using the STRING database (Szklarczyk et al., 2019). STRING generates the network of interacting genes/proteins depending upon the evidence and the strength of the data support. The STRING network settings were set to a minimum required interaction score of 0.7, which is considered high confidence. The genes of up to five generations were selected and downloaded as TSV files which contain all the interacting nodes with their interaction scores. The visualization of these interacting genes was done using Cytoscape (Shannon et al. 2003). The genes obtained interacting with each HDAC were subjected to Database for Annotation, Visualization, and Integrated Discovery (DAVID ) (Jiao et al., 2012) for gene set enrichment analyses.Overall workflow has been shown in Figure1.

Figure1

The Overall Workflow of the Study

Results

Expression profile of ZnHDACs in various cancer types

The mRNA expression analysis of ZnHDACs between various cancers has shown differential expression patterns shown in Figure 2. Class I HDACs HDCA1, HDAC2, and HDAC8, are overexpressed in all cancers but, the extent of upregulation is different such as, in BRCA and LUAD, the transcripts per million (TPM) for HDAC1 are around 80 or less whereas, COED and READ have TPM as high as 150. HDAC2 also shows a similar kind of trend. HDAC3 has similar expressions in cases compared to controls in most cancers, but there is observable upregulation in HNSC and downregulation in PRAD. Class II HDACs, HDAC5, HDAC6, HDAC7 and, HDAC10 mostly were downregulated, except in HNSC where HDAC7 and HDAC10 show upregulation and HDAC6 has a similar expression as controls. HDAC4 and HDAC9 have a negligible overall expression in all cancer types, still in COAD and READ a downregulation observed compared to controls. Class IV, HDAC11 have similar expression patterns as controls in each cancer type except for, BRCA where it shows a significant upregulation.

Figure 2

Expression of ZnHDACs in Different Cancer Types between Patients Having Cancer and Normal Controls

Comparative analysis of ZnHDAC expression with the survival of patients in various cancer types

mRNA expressions of ZnHDACs are compared with patient survival to see the effect of aberrant expression patterns clinically on overall (OA) and disease-free survival (DF). The analysis was done by constructing Kaplin Meyer plots for each ZnHDACs in all cancer types. The parameters, such as log-rank P values (LR) and hazard ratio (HR) reported in Table 2. LR less than 0.05 is considered significant. Although class I ZnHDACs upregulate in various cancers, only HDAC2 overexpression adversely affects breast cancer patient survival with an HR of 1.7. HDAC 3 shows slight upregulation in BRCA, COAD, HNSC, and down-regulation READ, LUAD, LUSC, and PRAD. But, HDAC3 downregulation is responsible for short DF survival of PRAD significantly by an HR of 2. In Class II ZnHDACs, HDAC4 expression was low compared to HDAC5, HDAC6, and HDAC7 still, HDAC4 downregulation is affecting DF survival adversely for BRCA patients with HR of 1.4. There is HDAC4 and HDAC6 upregulation in HNSC patients, and interestingly it increases survival rate with HR of 0.73 and 0.77. HDAC5, HDAC6, HDAC7, and HDAC10 are downregulated in all cancers but, only affect the DF survival of PRAD with HR of 1.6, 1.5, 2.2, and 1.9 respectively. In LUSC, HDAC11 upregulation reduces the overall survival significantly with an HR of 1.5.

Table 2

The Overall and Disease Free Survival Data for Each HDAC in Various Cancer Types

		BRCA		COAD		READ		HNSC		LUAD		LUSC		PRAD
		OA	DF	OA	DF	OA	DF	OA	DF	OA	DF	OA	DF	OA	DF
	HR	0.99	0.82	1	0.93	1.7	1.1	0.97	0.51	1.3	0.98	1.1	1.1	0.94	0.95
HDAC1	LR	0.93	0.29	0.86	0.78	0.29	0.86	0.84	0.89	0.09	0.89	0.41	0.48	0.92	0.79
	HR	1.7	1.4	0.74	1.1	0.43	0.75	0.39	1	1.3	1.2	0.91	1.1	1.2	1.1
HDAC2	LR	0.00085*	0.08	0.22	0.59	0.09	0.54	1.1	0.96	0.55	0.2	0.48	0.75	0.78	0.76
	HR	1.2	0.83	1.1	1.5	0.51	0.79	1.1	1.3	1	1	1.1	1.4	2.4	2
HDAC3	LR	0.35	0.32	0.61	0.1	0.17	0.6	0.55	0.087	0.85	1	0.63	0.07	0.2	0.00094*
	HR	1.2	1.4	0.7	0.66	0.85	1.4	0.73	0.69	1.3	1.2	1.2	1.2	0.68	1.3
HDAC4	LR	0.17	0.05*	0.15	0.093	0.74	0.49	0.021*	0.027*	0.11	0.29	0.29	0.3	0.55	0.17
	HR	1	1.3	1.5	1.1	0.51	1.3	0.99	0.81	0.9	0.87	1.1	1.3	0.91	1.6
HDAC5	LR	0.99	0.12	0.1	0.59	0.17	0.54	0.93	0.2	0.49	0.37	0.68	0.11	0.89	0.039*
	HR	1.1	0.76	1.2	1.3	0.83	1.6	0.77	1.1	1.2	0.87	1.1	1.3	4.4	1.5
HDAC6	LR	0.75	0.1	0.55	0.28	0.72	0.35	0.05*	0.46	0.34	0.35	0.61	0.19	0.047	0.039*
	HR	1.11	1.1	1.2	1.1	1	1.8	0.78	1	1.1	1.2	1.1	1.3	2.6	2.2
HDAC7	LR	0.55	0.73	0.4	0.82	0.92	0.23	0.06	0.83	0.46	0.18	0.39	0.18	0.15	0.00022*
	HR	1.1	1.1	1	1.2	0.65	1	0.92	1.2	1.3	1.2	0.94	1.1	1.4	1.3
HDAC8	LR	0.59	0.5	0.94	0.37	0.37	0.94	0.55	0.22	0.11	0.18	0.65	0.55	0.58	0.2
	HR	0.86	0.94	1	1.5	1	1.5	1.2	1.2	0.94	0.8	1.1	0.89	0.34	0.92
HDAC9	LR	0.35	0.72	0.97	0.076	0.93	0.42	0.26	0.39	0.66	0.15	0.6	0.53	0.16	0.68
	HR	0.97	1.4	1.3	1.5	0.69	2.6	0.84	0.84	0.96	0.88	1	1.1	2.3	1.9
HDAC10	LR	0.87	0.09	0.31	0.11	0.45	0.064	0.19	0.3	0.78	0.4	0.88	0.67	0.22	0.0045*
	HR	0.92	1	0.77	0.68	0.46	0.5	0.84	1	0.87	1	1.5	1.3	1.7	1.3
HDAC11	LR	0.59	0.99	0.29	0.12	0.14	0.16	0.21	0.92	0.35	0.96	0.0035*	0.2	0.41	0.16

HR, Hazard ratio; LR, logrank p value; OA, overall survival; DF disease free survival; *logrank p value less than 0.05 are considered significant

Functional enrichment analyses of ZnHDACs

The STRING networks for five generations gave us numerous genes interacting with these ZnHDACs. These genes were subjected to DAVID functional annotation tools to understand the biological meanings behind this list of genes. The significantly enriched Gene Ontology (GO) terms and pathways in all ZnHDACs and their neighboring genes are further analyzed.

The enrichment terms in common for all ZnHDACs and their interacting genes are listed in Supplementary Table 1. Common functional enrichment associated with the binding is either sequence-specific, transcription regulatory region, double-strand binding or chromatin, transcription factor binding, and histone deacetylase binding. ZnHDACs have a common role in catalytic activity, transcription regulatory activity, protein binding, NAD-dependent histone deacetylase activity, protein N terminus binding, and nuclear receptor binding. The common biological processes are chromatin organization, histone modification such as deacetylation, regulation of RNA metabolic, and other macromolecular metabolic processes. These ZnHDACs are common in the negative regulation of gene silencing, miRNA processing, and gene expression. The enriched common pathways are the notch signaling pathway, thyroid signaling pathways, transcriptional misregulation in cancer, and alcoholism. ZnHDACs are also commonly enriched in Epstein Barr virus infection, viral carcinogenesis, myeloid leukemia, Human papillomavirus infection, cell cycle, and miRNAs in cancer. The –log FDR values for these common enrichment terms components for each ZnHDACs and its interacting genes are different, shown in Figure 3. That means ZnHDACs and associated genes are common in functions but, their involvement is more significant in some HDACs and less in others. For example, histone deacetylase complex and nuclear chromatin components are highly enriched in HDAC1, HDAC2, HDAC9, and least on HDAC4 and HDAC6. Similarly, chromatin organization, Histone deacetylation, Histone modification are more significant in HDAC3, HDAC7, and HDAC9 than other ZnHDACs.

Figure 3

Plot of -log FDR Values of each ZnHDACs in Common Enrichment Terms

Further, we obtained the unique and significant enrichment terms for each ZnHDACs these enrichment terms are reported in Supplementary Table 2. These unique interaction terms show the uniqueness of each ZnHDAC in terms of functions and gene interaction. The genes that interact with a ZnHDAC directly or in a cascade and not with others are unique. We obtained a few uniquely interacting genes for each ZnHDACs shown in Figure 4. Some of these are JARID2 for HDAC1, BPTF, CHRAC1, and BAZIA for HDAC2, YEATS4, and DMAP1 for HDAC3, RCOR1, KDM1A, TAF10, and ATXN7 for HDAC7, etc.

Figure 4.

Unique Genes for Each ZnHDAC Subtypes

The unique enrichment terms significant for class I ZnHDACs are stem cell proliferation, negative regulation of transcription for HDAC1. Chromatin remodeling, and positive regulation of transcription for HDAC2. H2A acetylation, and transcriptional corepressor activity for HDAC3. Regulation and maintenance of telomere, and chromatin assembly and disassembly for HDAC8. The unique enrichment terms significant for class II ZnHDACs are viral carcinogenesis, positive regulation of apoptosis, and negative regulation of cell proliferation for HDAC4. Negative regulation of transcription from RNA polymerase II promoter, and negative regulation of receptor biosynthetic process for HDAC5. Regulation of protein stability complex and binding of ubiquitin complex and heat shock protein for HDAC6. STAGA complex and nucleoplasm translocation for HDAC7.Histone methyltransferase and B cell differentiation for HDAC9. The unique enrichment terms significant for class III ZnHDAC, HDAC11 in alcoholism and extracellular exosome.

Discussion

The differential expression of ZnHDACs in different cancer types and within the same cancers supports using SSI’s over PHI (Li and Seto, 2016; Bandolik et al., 2019). Many studies support the differential expression of these ZnHDACs in various cancers (Lucio-Eterovic et al., 2008; Müller et al., 2013). Even the ZnHDAC’s belonging to the same class shows differential expressions in the present study. Not only class-specific but subtype-specific inhibitors are the need of cancer therapeutics. The survival analysis with ZnHDAC expression revealed that a single HDAC might affect the overall and disease-free survival adversely or favorably. HDAC 1 has high expression in almost all cancers and has a null effect on survival in the present study whereas, an earlier report shows significance with lung adenoma (Minamiya et al., 2011). The survival data from the present study reported in Table 2 shows that only HDAC2 overexpression is adversely affecting the survival of the patients having breast cancer. In support, few studies report that increased expression of HDAC2 is associated with aggressive progression and poor prognosis of breast cancer (Müller et al., 2013; Zhao et al., 2016 p. 2). Therefore, HDAC 2 is considered a drug target in breast cancer, and oleuropein is a therapeutic lead for HDAC2 inhibition (Bayat et al. 2019). The results from our study and an earlier study have an agreement upon combined overexpression of class I HDACs is in colorectal carcinomas(Li and Seto, 2016). Mostly, Class II HDACs are down-regulated in all cancers except HNSC. A recent study has shown HDAC6 role in regulating stem cells and differentiated cancer cells (Sharif et al., 2019). HDAC10 downregulation was associated with KRAS-associated lung adenocarcinomas as its downregulation is significant in LUSC and LUAD in the present study (Li et al., 2020). HDAC11 upregulation in LUSC shows adverse effects on the survival of the patients and causes stem cell renewal in LUAD (Wang et al., 2009; Bora-Singhal et al., 2020). Further, the unique and significant enrichment terms and genes for each ZnHDACs segregate them depending on their function despite their similar structure. The gene unique for HDAC1 is JARID2, a transcriptional regulator for methyltransferases, and is responsible for nucleation and propagation of heterochromatin (Ge et al., 2019). JARID2 is important in stem cell and embryonic development (Cho et al. 2018). The JARID-HDAC1 axis represses the pluripotency of stem cells (Wiśnik et al., 2017). BAZIA, BPTF, CHRAC1 are crucial nucleosome modeling enzymes unique to HDAC2 (Migliori et al., 2012). HDAC3 and its co-repressor DMAP1 are involved in transcriptional repression by co-repressor mediated path that is elucidated and reported in the HDAC3 crystal structure (Watson et al., 2012). HDAC4 and HDAC5 are downregulated, have common enrichment terms, these are negative regulation of proliferation, receptor biosynthesis process, DNA binding regulation. HDAC4 and HDAC5 associated genes such as MYOCD, SRF, KAT5 are coactivators of proliferation (Wei et al., 2019 p. 5; Tong et al., 2020). Therefore, in some studies, HDAC4 and 5 are upregulated in cancers (Wang et al., 2014; Li et al., 2016 p. 5). Some genes such as FOXO1, XRCC6, KAT5, and SAP30 are negatively regulating cell cycle, and proliferation (Unoki et al., 2009; Zhu et al., 2016; Jiang et al., 2018) therefore, HDAC4 and HDAC5 expression is varying in different cancers. HDAC6 is uniquely associated with various heat shock proteins, poly ubiquitin-binding, and aggresome formation. HDAC6 directly affects protein folding and stability, earlier studies reported that HDAC6 regulates oxidative stress by interacting with HSP’s and stress granules (Aldana-Masangkay and Sakamoto, 2011 p. 6). In this study, HDAC7 was associated uniquely with the genes for STAGA complex formation, involved in pre-mRNA splicing and nucleoplasm translocation. HDAC is involved in splicing and subnuclear translocation but, HDAC7 involvement individually is not reported yet (Koutelou et al., 2010). HDAC8 unique enrichment terms are telomere maintenance and chromatin assembly and disassembly it’s in agreement with earlier reports (Chakrabarti et al., 2015). HDAC9 and 11 in methyltransferase activity and overall histone modification. HDAC10 was not having any unique enrichment term or unique associated genes in this study. ZnHDACs genes in unique interactions should be participating differently in various cell cycle, and proliferation-related pathways. Overlooking the present expression analysis and involvement of these ZnHDACs in fundamental processes, a careful selection of SSIs should be done in individual patients to avoid serious side effects. The correlation between survival and expression should also be considered for SSI selection. Expression patterns and functional analysis can also be used for other targets facing subtype selectivity issues.

In conclusion, the HDAC subtype-specific data-mining and interactome analyses have provided significant differences in their role in carcinogenesis even if these subtypes share similar structures. The problem of subtype selectivity issue of HDAC inhibitors can be taken towards the functional aspects of these subtypes. The unique interactions by these subtypes may be explored more as some of the interactions we obtained are proven to be true experimentally. This may help to progress a step towards personalized medicine for cancer therapeutics.

Author Contribution Statement

Conceptualization, work plan and preliminary draft writing by SU, data procurement by AR, draft editing by CC and PP supervision by PS.