Systems approach to the pharmacological actions of HDAC inhibitors reveals EP300 activities and convergent mechanisms of regulation in diabetes

ABSTRACT

Given the skyrocketing costs to develop new drugs, repositioning of approved drugs, such as histone deacetylase (HDAC) inhibitors, may be a promising strategy to develop novel therapies. However, a gap exists in the understanding and advancement of these agents to meaningful translation for which new indications may emerge. To address this, we performed systems-level analyses of 33 independent HDAC inhibitor microarray studies. Based on network analysis, we identified enrichment for pathways implicated in metabolic syndrome and diabetes (insulin receptor signaling, lipid metabolism, immunity and trafficking). Integration with ENCODE ChIP-seq datasets identified suppression of EP300 target genes implicated in diabetes. Experimental validation indicates reversal of diabetes-associated EP300 target genes in primary vascular endothelial cells derived from a diabetic individual following inhibition of HDACs (by SAHA), EP300, or EP300 knockdown. Our computational systems biology approach provides an adaptable framework for the prediction of novel therapeutics for existing disease.

KEYWORDS:

• epigenomics
• HDAC inhibitors
• meta-analysis
• SAHA
• histone acetylation
• diabetes
• endothelial dysfunction
• network analysis
• pathway analysis
• computation biology

Background

Drug repurposing is a process by which existing, clinically approved drugs are repurposed for different conditions. It offers advantages over traditional drug discovery and development, including being time- and cost-effective. However, in order for drug repurposing to be efficient, we must first have an extensive understanding of the cellular and molecular effects of these drugs.

Potential candidates for drug repurposing are histone deacetylase inhibitors (HDACi), which represent a structurally diverse group of compounds that have emerged as important anticancer therapeutics. Although early work focused almost exclusively on antineoplastic properties, HDACi are now being investigated for a host of non-oncological applications that are driven, at least in part, by their potent anti-inflammatory effects.¹⁵ For example, pre-clinical models show that HDAC treatment attenuates cardiac hypertrophy,^36,53 while SAHA pre-treatment attenuates cardiac ischemia injury.⁵¹ In addition, HDAC inhibitors SAHA and givinostat were shown to prevent the development of type 1 diabetes (T1D) in the non-obese diabetic (NOD) mouse model by modulation of the immune system.⁸

The initial paradigm was that HDACi activate gene expression by hyperacetylating histones, resulting in a more open chromatin confirmation⁴⁹ while gene suppression was usually attributed to non-histone effects of HDAC inhibition.⁴ However, the utilization of genome-wide technologies has revealed more complex interactions with the epigenome. In particular, histone deacetylation at the promoters of suppressed genes, mediated by loss of histone acetyltransferase (HATs) binding, was identified as a key feature of pan-HDAC inhibition.³⁸

HDAC inhibitors include drugs from different chemical families and can target different HDAC classes. Pan-HDAC inhibitors, such as hydroxamic acids suberoylanilide hydroxamic acid (SAHA, Vorinostat, Zolinza™) and trichostatin A (TSA), inhibit the activity of Class I, II, and IV, but not Class III, HDAC enzymes.^13,30,31 However, compounds such as romidepsin, largazole, and valproic acid (VPA) selectively inhibit Class I HDAC enzymes and there is currently interest in the further development of class- or isoform-specific compounds.^{13,22,30,31,52} However, despite the diversity in structure of these drugs, they appear to have overlapping clinical applications. For example, both SAHA and romidepsin are approved for the treatment of CTCL, although they target different HDAC enzymes. However, it remains unclear if this diverse group of drugs regulates common gene expression pathways.

To investigate the effects of diverse HDAC inhibitors on the transcriptome, we performed a comprehensive meta-analysis on gene expression profiles in normal, transformed and malignant cell lines, which we then integrated with extensive network analysis (Fig. 1). While network analysis is currently utilized to develop a better understanding of diseases and their progression, the application of this approach to understanding cellular effects of drugs, known as ‘systems pharmacology’, is less widely used. The integration of systems pharmacology with network-based analysis of disease provides a time- and cost-effective solution to repurposing clinically approved drugs for human diseases.⁷ By utilizing this approach, we observe highly correlated gene activation and suppression signatures, independent of drug class or cell type, which were grouped into broad pathways, including cell cycle, the immune system, and RNA and protein processing and metabolism. Integrative epigenomic analyses combined with transcription factor binding sites (TFBS) derived from the Encyclopedia of DNA Elements (ENCODE) project (2004) show histone deacetylation is strongly associated with the suppression of EP300 target genes. The integration of multiple datasets identified EP300 target genes were associated with complex diseases that were previously implicated in, but not limited to, cancer, and uncovered novel HDAC inhibition networks associated with the removal of acetylation at gene targets functionally linked with diabetes and vascular disease. Thus, to understand the mechanisms of action, we investigated regulatory patterns of expression using primary human aortic endothelial cells (HAECs) of diabetic and non-diabetic origin. We show reversal of diabetes-associated EP300 target genes in diabetic primary cells by SAHA. Our results suggest the integrative analyses we have used to investigate pharmacological HDAC inhibition show complex epigenetic interactions with EP300 and novel therapeutic links relevant to diabetic vascular pathology.

Figure 1. Broad study overview. The overall design of this study has been summarized in 4 steps. We first performed an extensive meta-analysis of public microarray datasets that involved treatment of diverse human cell types with various pharmacological HDAC inhibitors (Step 1). Network analysis was used to broadly profile the pathways altered by HDAC inhibition (Step 2). Of note, we identify consistent changes to the immune system, cell cycle, and signaling pathways. Transcription factor analysis identified EP300 as a key regulator of genes involved in immunity, lipid metabolism, trafficking and insulin receptor signaling following HDAC inhibition (Step 3), pathways that are dysregulated in diabetic vascular disease. Next, we performed ex vivo validation in a diabetic disease model to explore the therapeutic potential of HDAC inhibitors (Step 4).

Results

Pharmacological HDAC inhibition segregates gene expression

To determine a representative expression signature, we individually re-analyzed 33 microarray experiments derived from human cells exposed to pharmacological HDAC inhibitors. The studies include isotype specific and pan-HDACi in cancer and primary human cell types. We observed robust expression changes (Fig. 2a) that include gene activation and suppression (Fig. 2b). A heatmap shows strong clustering based on gene expression, irrespective of cell type and drug class (Fig. 2c). Spearman's ranks also show strong correlations in gene expression between datasets (Fig. 2d), even between structurally unrelated HDAC inhibitors (Figure S1a-d). This initial assessment of the data suggests common gene expression signatures induced by all classes of HDAC inhibitors.

Figure 2. Gene expression changes induced by HDAC inhibition in different human cell types. Differential gene expression for independent human microarray HDACi datasets are summarized. (a) Line plot showing the total number of genes differentially expressed at FDR P values between 0 and 0.05. (b) Pie chart summary of the number of differentially expressed genes (FDR P value <0.05). (c) Heatmap of gene expression changes conferred by HDAC inhibition in human cell types. Each column represents an independent microarray experiment. Blue represents gene activation while red indicates suppression. The intensity of the color correlates with the strength of the P value for any given gene. Columns were normalized by Z-score normalization. (d) Spearman's rank correlation (ρ) of gene expression changes conferred by HDAC inhibitors in different human cell types. Positive correlations are shown in blue while negative correlations are shown in red. A histogram of values and a color key are provided for each heatmap.

Meta-analysis using Fisher's method for combining P values was used to identify gene expression in response to HDAC inhibition. We identify genes that are either activated or suppressed (Figure S1e) in diverse human cell types. For example, HDAC inhibition increased the expression of Class I (HDAC1, HDAC3), IIA (HDAC5), and Class III (SIRT4, SIRT2) HDACs, while suppressing HATs, including CREBBP and EP300.

In order to determine broad classes of genes consistently regulated by HDAC inhibition, network analysis was performed in combination with Gene Set Enrichment Analysis (GSEA) using REACTOME. We identified common pathways altered by HDAC inhibition and the suppression of genes involved in the cell cycle, immune system, and RNA and protein processing and metabolism (Fig. 3). We observed activation of genes related to signaling and receptor binding (including insulin receptor signaling), telomere packaging, membrane trafficking, and lipid metabolism (Fig. 3). Gene lists were generated from each cluster to define functional gene sets that consistently respond to HDAC inhibition.

Figure 3. Network analysis identifies common pathways altered by HDAC inhibition. Gene set enrichment results for REACTOME pathways altered by HDAC inhibition were visualized using Enrichment Map in Cytoscape (FDR P value <0.05; overlap coefficient cut-off 0.8, combined constant 0.5). Clusters were broadly classified based on biological function. The node colors indicate pathway activation (red) or suppression (blue). Green lines indicate the overlap between two nodes.

To validate the predictive ability of the meta-analysis, genes identified from the microarray studies were compared to independent RNA-seq and ChIP-seq HDAC inhibitor datasets, obtained from GEO. The datasets include exposure of both cancer and primary cell types to various HDAC inhibitors, including SAHA, TSA, and panbinostat. Using GSEA, we show that top HDAC inhibitor genes derived from the microarray meta-analysis were reproduced in the independent RNA-seq studies (Figure S2a). Next, we compared the gene expression signature from the meta-analysis to histone acetylation data from ChIP-seq. In primary human aortic endothelial cells (HAECs) exposed to SAHA or TSA, genes activated in the meta-analysis were strongly associated with increased histone acetylation, but suppressed genes were only weakly associated with histone deacetylation (Figure S2b). Taken together, these results suggest the expression of genes activated by HDAC inhibition are strongly associated with histone acetylation, whereas gene suppression is likely to be driven by both acetylation dependent and independent mechanisms.

To further explore the role of histone acetylation by HDAC inhibition, we integrated the gene expression (RNA-seq) and histone acetylation (ChIP-seq) data from primary HAECs exposed to SAHA and then intersected this data with the top genes activated or suppressed by HDAC inhibition as determined by the meta-analysis. We confirmed that genes activated in the meta-analysis were indeed activated and acetylated by SAHA in HAECs (Fig. 4a, gated box of 426 genes). We also show genes suppressed in the meta-analysis are generally suppressed and deacetylated in HAECs exposed to SAHA (Fig. 4b, gated box of 460 genes). Genes in each gated box were classified based on the gene sets defined from the REACTOME network clusters from Fig. 3 (Fig. 4c). Genes involved in trafficking, neuronal systems, metabolism, and signaling (including insulin receptor signaling) were associated with increased expression and acetylation. In contrast, genes involved in RNA processing, DNA repair, nucleotide metabolism, and the cell cycle were more likely to be suppressed and have reduced promoter acetylation. Genes involved in immunity (including infection and MAPK related pathways), although representing the most number of genes, were equally likely to be increased or decreased by SAHA.

Figure 4. Top genes predicted by meta-analysis are consistent with NGS gene expression and histone acetylation HDACi datasets. RNA-seq and ChIP-seq data from non-diabetic HAECs exposed to SAHA was obtained from GEO (GSE37378). A density plots shows the relationship between gene expression (RNA-seq) and histone acetylation at promoters (ChIP-seq) in response to SAHA in HAECs for the top 1000 (a) activated and (b) suppressed HDACi response genes as defined by the meta-analysis (legend: relative gene density). The gated red boxes highlight genes with (a) increased gene expression and acetylation or (b) decreased gene expression and acetylation by SAHA in HAECs. These genes were then classified based on (c) broad functional classes defined by the network analysis.

HDAC inhibition suppresses EP300 target genes

Next, to determine pathways regulated by HDAC inhibition for each microarray dataset, we utilized GSEA with the MSigDB gene set database.⁴⁶ We observe a strong association of regulatory pathways altered by iso- and pan-selective HDAC inhibitors in multiple human cell types. Pathway meta-analysis shows a strong correspondence for activated and suppressed gene sets. Since HDAC inhibition alters co-regulators of gene transcription,³⁸ we integrated transcription factor binding sites (TFBS) using ChIP-seq data sets originating from the ENCODE project.¹ We show candidate transcription factors as drivers of gene expression following HDAC inhibition (Figure S2c). Transcription factors were ranked and short-listed by meta-analysis using Fisher's method for combining P values. We identified EZH2, a PRC2 subunit and methyltransferase, as a regulatory component of activated target genes. In contrast, EP300 was strongly associated with gene suppression (Fig. 5a) and consistent with the notion that SAHA and TSA regulate histone deacetylation of EP300/CREBBP target genes.³⁸

Figure 5. EP300 target genes are associated with gene suppression following HDAC inhibition in human cell types. GSEA ENCODE-TF results associated with gene expression changes induced by HDAC inhibitors in multiple cell lines for (a) gene sets associated with epigenetics modifiers. GSEA was used to compare changes in expression induced by HDAC inhibition in human cell types to microarray gene expression signatures association with the knock down of (b) EP300 (GSE31873) and CREBBP (GSE50588). Blue represents gene sets associated with gene activation while red indicates gene sets associated with gene suppression. The intensity of the color correlates with the NES scores. Pathways with a FDR P value >0.05 are shown in white. (c) The relationship between gene expression (RNA-seq) and histone acetylation (ChIP-seq) for EP300 target genes (from the ENCODE project) are represented as density plots for HAECs exposed to SAHA (legend: relative gene density). (d) EP300 target genes with suppressed expression and reduced acetylation (red box) by SAHA in HAECs were classified based on biological function as defined by the network analysis.

Since EP300 regulates histone acetylation and gene expression, our observation of deacetylation at EP300 target genes was unexpected. Using gene sets derived from re-analyzed gene expression patterns derived from EP300 and CREBBP knockdown (KD) studies obtained from GEO, we used GSEA to show that genes suppressed following EP300 and CREBBP KD are strongly correlated with genes suppressed by pharmacological HDAC inhibition (Fig. 5b). We found suppressed EP300 target genes were strongly associated with histone deacetylation by the hydroxamic acids SAHA (Fig. 5c) and TSA (Figure S2d-e). Classifications of genes by function show biological processes that include RNA and protein processing, immunity and inflammation, the cell cycle and lipid metabolism, suggesting that EP300 is a critical driver of gene expression of key pathways altered by HDAC inhibition (Fig. 5d, Figure S2f).

Deacetylation of diabetes-associated EP300 target genes

Based on the pathways associated with HDAC inhibition and published pre-clinical and clinical trials, we identified a number of diseases that could potentially be treated by HDAC inhibition. These include cancer (cell cycle, telomeres, DNA repair), epilepsy, Huntington's and Alzheimer's disease (neuronal systems), inflammatory disorders and infection (immunity, signaling), cardiovascular disease (inflammation, developmental biology), and metabolic disorders such as metabolic syndrome and diabetes (insulin receptor signaling, lipid metabolism, immunity and trafficking). We questioned whether HDAC inhibition could regulate the expression of genes implicated in diabetic vascular disease, in an EP300-dependent manner. To do this, primary human aortic endothelial cells derived from non-diabetic (nHAEC) and type 1 diabetic (dHAEC) individuals were exposed to HDAC inhibition. We selected SAHA, an FDA approved drug, previously shown to regulate inflammatory gene expression via EP300-dependent histone deacetylation in endothelial cells.³⁸ Network analysis confirmed suppression of pathways related to the cell cycle, the immune system, and RNA and protein metabolism in diabetic HAECs treated with SAHA (Fig. 6).

Figure 6. Network analysis identifies pathways altered by SAHA in diabetes HAECs. Gene set enrichment results for REACTOME pathways altered by SAHA treatment in HAECs cells of diabetic origin were visualized using Enrichment Map in Cytoscape (FDR P value <0.05; overlap coefficient cut-off 0.8, combined constant 0.5). Clusters were broadly classified based on biological function. The node colors indicate pathway activation (red) or suppression (blue). Green lines indicate the overlap between two nodes.

Next, we examined the role of EP300 using the small molecule inhibitor C646 as well as somatic knockdown using EP300 siRNA. Gene expression profiles from non-diabetic and diabetic HAECs exposed to SAHA are highly correlated (r = 0.9) (Fig. 7a, Figure S3a). Gene expression response following C646 is also highly correlated between non-diabetic and diabetic HAECs (r = 0.89) (Fig 7a). In contrast, the correlation (r = 0.79) between non-diabetic and diabetic HAECs following EP300 KD by siRNA was not as strong (Fig. 7a).

Figure 7. Inhibition of HATs and HDACs mediates robust changes in gene expression. (a) Correlation plots with the linear model and Pearson's correlation for non-diabetic (nHAEC) and diabetic HAECs (dHAEC) exposed to SAHA, C646 (pharmacological EP300/CREBBP inhibitor) or EP300 siRNA. (b) Density plots show the relationship between gene expression (RNA-seq) and histone acetylation at promoters (ChIP-seq) for all genes and the top 1000 genes suppressed by SAHA, C656 or EP300 KD in diabetic HAECs (legend: relative gene density).

Furthermore, we identify genes suppressed by both HDAC and HAT inhibition (Figure S3b). To confirm deacetylation of EP300 target genes, we compared gene suppression patterns from SAHA, C646 and EP300KD experiments in diabetic HAECs with independent RNA-seq and histone acetylation ChIP-seq derived from non-diabetic HAECs stimulated with SAHA (Fig. 7b). These data show genes suppressed by SAHA and C646 in diabetic HAECs strongly correlate with histone deacetylation in non-diabetic HAECs. We also observe a similar association for EP300 KD, supporting the notion that EP300 target genes deacetylated by SAHA are implicated in diabetes.

We hypothesize that pharmacological HDAC inhibition may reverse gene expression patterns in diabetic HAECs by an EP300 dependent mechanism. To test this hypothesis, we inhibited HDAC activities using SAHA and examined genes activated in diabetes in dHAECs (Fig. 7c). We identified three distinct signatures for diabetes-associated genes. The first cluster shows diabetes genes suppressed by SAHA, C646, and EP300 KD. These results support a possible role for acetylation by EP300 in diabetic gene activation and the notion of functional reversal by SAHA. The second cluster consists of genes suppressed by SAHA and C646, but not EP300 KD, suggesting a role for HATs other than EP300, such as CREBBP. Finally, in clusters 3 and 4 we observe genes suppressed by SAHA that are regulated independently of EP300 (Fig. 8a). In summary, these results suggest that SAHA inhibits the inflammatory expression signature in diabetic HAECs by an EP300 dependent mechanism.

Figure 8. SAHA reverses gene expression in diabetic HAECs by EP300 dependent mechanism (a) A heatmap showing genes that are upregulated in diabetes and that are suppressed by SAHA in this dataset. Clustering identifies diabetic gene expression likely to be EP300-dependent or independent. Validation of gene expression was performed by RT-PCR. (b) Relative expression of genes associated with diabetes in diabetic HAECs relative to non-diabetic HAECs. (c) Relative expression of genes associated with diabetes in response to SAHA or C646 (relative to DMSO control) or (d) EP300 siRNA (relative to non-target siRNA). Gene expression was determined by normalizing against the housekeeping gene HPRT1, which was stably expressed. Error bars represent SEM. n = 3, (*) P value <0.05, (**) P value <0.005, (***) P value <0.0005, (****) P value <0.0001, unpaired t-test.

SAHA reverses the expression of EP300 target genes in diabetic HAECs

To validate the role for EP300 in diabetic HAECs identified from our integrative meta-analysis, we selected candidate genes that were activated in diabetes and suppressed by EP300KD in diabetic HAECs. The expression of VCAN, CD58, CD200, and SIGIRR genes were upregulated in diabetic HAECs (Fig. 8b). For example, VCAN gene expression was highly ranked with an almost 40-fold increase using qRT-PCR. This gene encodes versican, a large extracellular matrix proteoglycan involved in inflammation, and implicated in diabetes.^2,54 We observed increased H3K9/14ac at the VCAN gene promoter (Figure S4). The cell adhesion molecule CD58, which is involved in immunoregulation³⁴ and CD200, a cell surface molecule involved in endothelial T-cell interactions,²⁴ were also activated, as was the immunoregulatory protein SIGIRR.¹⁶ Exposure of dHAECs to SAHA and C646 reversed VCAN, CD58, CD200, and SIGIRR gene expression (Fig. 8c). Indeed, we observe greater than 50% reduction in the expression of CD200, VCAN, and SIGIRR by SAHA and C646. The expression of CD58 was suppressed more strongly by C646 compared to SAHA (Fig. 8c). Furthermore, we also observe strong suppression of the adhesion molecule ICAM2²⁹ and pro-inflammatory chemokine CCL2⁴² by both SAHA and C646. These findings show that both SAHA and C646 suppress the expression of pro-inflammatory genes in diabetic cells, suggesting a EP300-dependent mechanism.

To test this hypothesis, we examined whether knockdown of EP300 in dHAECs would suppress diabetes-related genes. We show suppression of the genes CD58, SIGIRR, CD200, ICAM2, and CCL2 is consistent with gene suppression induced by SAHA exposure (Fig. 8d). In particular, VCAN was more susceptible to suppression by EP300KD than by SAHA or C646 exposure. Collectively, these findings suggest that gene suppression mediated by SAHA is dependent on EP300 in diabetic endothelial cells.

Discussion

HDAC inhibitors are a diverse group of epigenetic modifying drugs with potent gene regulatory properties. Three compounds, SAHA, romidepsin, and panbinostat, are currently FDA approved.³⁵ Many other drugs are in phase II clinical trials for both hematological and solid tumors. While these compounds are most well known for their antineoplastic properties, there is less data available for their applicability in non-oncological diseases. For example, VPA has a long-standing history as an antiepileptic drug, while other agents that inhibit HDAC activity show pre-clinical applicability in diabetes and cardiovascular disease.^8,25,36

Current clinical use indicates some consistency in the therapeutic potential of diverse HDAC inhibitors. For example, SAHA and romidepsin, two structurally unrelated HDAC inhibitors, are both currently approved for the treatment of cutaneous T-cell lymphoma, while panbinostat and belinostat, both hydroxamic acids like SAHA, are in late stage clinical trials for the same malignancy. This suggests that HDAC inhibitors may regulate common pathways, despite diversity in drug class or chemical structure. To test this hypothesis, we integrated thirty-three microarray studies for gene expression changes conferred by a diverse group of pharmacological agents. Despite the diversity of drug and cell types, we observe comparable gene expression changes in response to HDAC inhibition. The strong overlap of gene expression changes conferred by pan- and Class I-specific HDAC inhibitors implicates Class I HDAC enzymes, rather than Class II, as potential drivers of gene expression. This is likely associated with the ubiquitous expression of Class I HDACs in human cell types. In contrast, Class II enzymes are cell type specific. For example, Class IIA HDAC expression is generally limited to smooth muscle, the brain and heart. Therefore, detection of expression signatures mediated by inhibition of this class of HDACs would not be possible using the publically available datasets access from GEO.

Consortium-based studies have shown the power of integrative approaches to multiomics. For example, the ENCODE project is the quintessential integrative analysis of transcriptomic and epigenomic datasets that have been highly informative of chromatin biology and gene function.¹ The Genotype-Tissue Expression (GTEx) project has characterized inter-individual transcriptomic variation across different tissues³² and has identified how genetic variation can affect gene expression.⁹ Collectively, such consortia are a powerful resource that can be leveraged by computational biologists. For example, one study utilizing data from the ENCODE project identified that changes in gene expression between individuals were associated with genetic variants that altered chromatin accessibility and transcription factor binding.¹¹

Currently, meta-analyses are widely used for combining clinical datasets and to gain deeper insights into biological questions. For example, a recent GWAS meta-analysis of more than 200,000 individuals identified novel genetic loci implicated in body fat distribution.⁴⁴ While there is no doubt that large datasets are informative for GWAS studies, there is still very little integration investigating epigenomic and regulatory mechanisms. The rapid generation of microarray and high throughput sequencing data available on public depositories in recent years provides a valuable resource to researchers, especially in the context of meta-analyses, pathway analysis, and systems biology.

Complementary to a systems biology understanding of disease, is the extension of network-based approaches to gain a deeper appreciation of the transcriptional and cellular effects of drugs. The integration of this ‘systems pharmacology’ approach with network-based analysis of disease provides a rational, cost-effective way and evidence based solution to repurposing clinically approved drugs for human diseases.⁷ This is a viable and cost-effective alternative to the costly and time-consuming development of new drugs. Furthermore, this approach could also provide a rational approach to drug combinatorial approaches, as well as the prediction of drug-drug interactions and off-target effects.²⁸ Furthermore, the rapidly declining cost of high throughput sequencing presents an opportunity for personalized medicine.³ For example, an in silico computational approach was utilized to predict that the anticonvulsant drug topiramate could treat inflammatory bowel disease.¹⁴ More recently, a similar approach identifies that dexamethasone may inhibit the activity of EGR, an oncogenic transcription factor.¹⁷

HDAC inhibitors, in particular, are excellent candidates for integrative meta-analysis given their extensive molecular characterization and the diverse compounds that exist in this class. Our meta-analysis approach integrated data from the ENCODE project with microarray expression profiles to show EP300 target genes suppression by HDAC inhibition, including genes implicated in diabetes.

Our findings of histone deacetylation suggest altered EP300 activity following HDAC inhibition. Currently, the mechanism by which HDAC inhibition alters EP300 activity remains unclear. However, we know that EP300 can auto-acetylate and acquire acetylation by other HATs,^6,45 and is deacetylated by HDAC enzymes.^6,48 While early studies suggest that auto-acetylation of EP300 enhances its activity and stimulates interactions with other proteins,^18,45 more recent evidence suggests that auto-acetylation of EP300 impairs its ability to acetylate K18 on histone H3.⁵⁵ Although the exact effects of EP300 acetylation status remain unclear, we propose that altered EP300 acetylation following the inhibition of HDAC enzyme activity may account for the histone deacetylation patterns observed in our study. Furthermore, EP300 auto-acetylation may also indirectly regulate gene expression via its altered interactions with other transcription factors.

As for the importance of chromatin modification, several regulatory determinants involved in histone acetylation have been described in the regulation of insulin resistance and diabetes.⁴³ For example, the activation of mouse insulin resistance genes is associated with both increased H3K27 acetylation and increased glucocorticoid receptor recruitment at gene enhancers in adipocytes.²¹ Furthermore, the activity of HDAC3 is elevated in peripheral blood mononuclear cells derived from individuals with type 2 diabetes when compared to non-diabetic controls.⁴¹ These results suggest changes to histone and non-histone substrates conferred by HDAC3 could be involved in the pathogenesis of diabetic complications. In addition, HDAC inhibition is associated with increased angiogenesis in endothelial cells.²⁰ In a previous study, we show the suppression of pro-inflammatory genes is mediated in part by EP300 activity.³⁸ In the current study, we extend these results showing pharmacological inhibition of HDAC activity is associated with histone deacetylation at EP300 target genes implicating the reversal of pro-inflammatory gene expression in diabetic endothelial cells. While the clinical implications for direct EP300 inhibition remains to be functionally confirmed, the current study adds to the growing body of literature of the potential therapeutic benefit of HDAC inhibitors in non-malignant disease.

While the datasets assessed in this study were derived from malignant cell types, our findings were informative of therapeutic potential in non-malignant disease. Recent studies have shown therapeutic benefit of HDAC inhibitors in diseases other than cancer. For example, SAHA and givinostat prevent the development of type 1 diabetes in a pre-clinical animal model.⁸ In that study, acetylation of the NFKB subunit p65 is linked to suppression of genes involved with inflammation. While that study focused on the prevention of type 1 diabetes in mice, our approach examined primary human vascular endothelial cells derived from diabetic and non-diabetic individuals. Endothelial cells form a critically important and highly functional barrier between the circulation and the blood vessel wall.³⁹ Endothelial cells are also among the first non-hematopoietic cells exposed to extracellular stimulation derived from blood, including medication, but also disease-causing stimuli such as elevated blood glucose. In particular, prolonged hyperglycemia can lead to endothelial dysfunction and is important in the progression of diabetic complications and the development of atherosclerosis.

Our studies identified pathways that are regulated by HDAC inhibition, with implications for novel drug repurposing for different diseases. For example, we show strong regulation of neuronal systems, implicating that HDACi may have therapeutic potential as treatments for neurological disorders. Indeed, recently, in an animal model of Alzheimer's disease involving severe amyloid deposition, treatment with SAHA improved spatial reference memory without altering plaque load.⁴ Further integrative analysis of genome-wide expression and histone acetylation data revealed that SAHA rescued the expression of genes associated with synaptic plasticity and inflammation. HDAC inhibitors are also protective in models of Huntington's disease.¹⁹ Collectively, these findings suggest that our approach to studying HDAC inhibition may be informative for identifying potential drug-diseases pairs.

While this study provides important insights into HDAC inhibition, the lack of diverse datasets limited the scope of our meta-analysis. For example, HDAC inhibitors are known to induce apoptosis at much higher rates in cancer and immortalized cell types when compared to primary cells. However, due to a lack of studies in primary cell types, comparison between HDAC inhibition in primary versus malignant cells was not feasible. We could not compare Class I and II HDAC inhibitors because many of the studies in GEO used pan- or Class I-specific HDAC inhibitors. RNA-seq and ChIP-seq datasets from these studies were also limited, and while we had access to H3K9/14ac profiles following HDAC inhibition, other histone acetylation marks such as H3K27ac were not assessed. Unlike H3K9/14ac, which is predominantly localized at gene promoters, H3K27ac regulates gene expression via acetylation of histones localized at both promoters and enhancers. In addition, further studies are required to investigate the acetylation status of EP300 post-HDAC inhibition. With the recent advances in ChIP sequencing technologies in non-malignant disease we believe more comprehensive meta-analysis will be achievable.

In summary, our integrative systems biology approach examining pharmacological HDAC inhibition identified key pathways regulated by SAHA with implications for novel therapeutic applications. As major consortia continue to grow, the integrative approach to meta-analysis and transcription factor analysis presenting in this paper will become increasingly relevant to further elucidating mechanistic pathways.²³ Furthermore, while cellular and molecular pathways, such as the cell cycle and immunity, are well defined by databases, such as REACTOME and KEGG, extensive disease-specific databases that integrate new and existing expression datasets are still lacking. As the field moves forward, establishing a better understanding of genomic and epigenomic mechanisms of drug action and disease etiology via systems biology approaches will be key for the repurposing of pharmaceutical therapies for existing diseases as well as rational drug design.

Experimental procedures

Screening for public HDAC inhibitor datasets

GEO (gene expression omnibus) was mined for independent HDAC inhibitor microarray or next generation sequencing using the ‘series’ filter. The search terms “HDAC inhibitors”, “HDAC inhibition” and “HDACi” were used for human datasets (searches performed between the 1st and 6th of May, 2014 and on the 21st of April 2015).

General inclusion criteria

Independent experiments in primary and transformed human cells were included. Tissue types include whole organs, organ components, or tumors. Pharmacological exposure to isotype specific and pan-HDAC inhibitors was included, with no restrictions on length or dose of exposure. Experiments involving exposure to a combination of HDACi were also included.

General exclusion criteria

HDAC gene knockout/down and co-exposure to a non-HDACi drug or siRNA/shRNA were not included.

Additional exclusion criteria for microarray

Gene expression datasets that match the following criteria were excluded from analysis: two-channel microarrays, experiments with no replicates (minimum n = 2), data that is not normalized, and arrays other than gene arrays (exon, miRNA, ChIP arrays). Datasets with less than 50 genes with a FDR P value of less than 0.05 were also excluded. Further details are described in the Supplementary Experimental Procedures.

Public microarray and NGS datasets included for meta-analysis

A total of 13 different isotype specific and pan-HDAC inhibitors are represented in this study (Table S1). The six Class I specific HDAC inhibitors are VPA, romidepsin, largazole, CHR-3996, ITF-A and ITF-B. The seven pan-HDAC inhibitors include TSA, SAHA, panbinostat, givinostat, HC-toxin, sodium butyrate (NaB), and ST7612AA1.

Differential gene expression of public microarray data

Thirty-three human microarray datasets were included for meta-analysis (Table S2). Normalized microarray data were obtained from GEO and analysis was performed using the limma⁴⁷ and GEO2R (default parameters). Differential gene expression used the Benjamini-Hochberg adjusted FDR P value <0.05⁵ with fold change greater than log2(1.2). Probes not annotated by gene name were removed from the analysis.

Differential gene expression and histone acetylation of public NGS data

NGS datasets analyzed are summarized in Table S3. For each study, data obtained in SRA format was converted to fastq using fastq-dump (sratoolkit.2.3.5-2-ubuntu64). RNA-seq: fastq files were aligned to the hg19 RefSeq transcriptome (download from UCSC) using the Burrows-Wheeler Aligner, (BWA).²⁶ A matrix was generated to summarize transcript reads and edgeR software was used to determine differential gene expression (Robinson et al.,⁴⁰ Significance is defined as FDR P value <0.05. ChIP-seq: fastq files were aligned to the hg19 (human) or mm9 (mouse) RefSeq genomes using BWA. Reads from each sample were counted at gene promoters (2.5 kb either side of the TSS) using bedtools³⁷ and edgeR was used to determine differential histone acetylation. Significance is defined as FDR P value <0.05.

Gene set enrichment analysis

GSEA was used for pathway and transcription factor analysis.⁴⁶ Transcripts from microarray and NGS expression data were assigned a score based on the negative log10 of the P value multiplied by the sign of the fold change (–logP*signFC). GSEA was run using classical scoring with 1000 permutations. Only gene sets between 30 and 2000 genes were included in the analysis, with the exception of GSEA on MSigDB pathways (version 4) (<1000 genes). For transcription factor analysis, gene sets of transcription factor target genes were generated from ENCODE transcription factor binding site (TFBS) ChIP-seq bed files. TFBS bed file represents a list of genome regions targeted by transcription factor binding. Target genes were defined as possessing a TFBS within 3 kb either side of the TSS.

Network analysis

GSEA of REACTOME pathways was performed on a single rank file generated from the gene expression results from the meta-analysis. REACTOME pathways altered by HDAC inhibition were visualized using Enrichment Map³³ in Cytoscape (FDR P value <0.05; overlap coefficient cut-off 0.8, combined constant 0.5). The REACTOME database was selected because it is a comprehensive curated database of human pathways used for pathway analysis.^10,50

Analysis

Integrative analyses for gene expression, network pathways and ENCODE (TFBS) transcription factor binding sites (using GSEA) were performed using Fisher's method for combining P values from the R package MADAM. To perform Fisher's test, a matrix of P values for genes derived from all datasets was generated. To test for gene activation, P values for suppressed genes were replaced by a P value of 1 (to indicate that they are not activated) and Fisher's method was performed on the matrix. To test for gene suppression, all P values for activated genes were replaced with 1 (to indicate that they are not suppressed) and Fisher's method was performed on the matrix. Because this test is performed twice, some genes may be identified as both activated and suppressed, and were excluded. Genes that are not represented in at least 70% of arrays (i.e., are not included in the microarray platform) were excluded from the test. Furthermore, genes not represented in pan- and iso-HDACi datasets were excluded from analysis (minimum requirement of 30% of datasets pan-HDACi and 30% iso-HDACi). Results are reported as the S-statistic and the FDR P value.

Primary human cells

Human aortic endothelial cells (HAEC) from type 1 diabetic (CC-2919, Lonza) and non-diabetic (CC-2535, Lonza) individuals, obtained at autopsy, were maintained in endothelial cell basal growth media-2 (EBM-2, Lonza) with additional EBM-2 growth factors and supplements and with 10% heat-inactivated fetal bovine serum (Gibco). HAECs used in this study were between 5–9 passages. Cells were exposed to DMSO, SAHA (2 μM) or C646 (10 μM) for 12 hours. EP300 knockdown was generated using ON-TARGETplus Human EP300 (2033) siRNA – SMARTpool (L-003486-00-0005, Dharmacon). Cells were cultured in Opti-MEM (Life Technologies) media with Lipofectamine RNAiMAX Reagent (13778-075, Life Technologies) with either non-target (ON-TARGETplus Non-targeting Pool, D-001810-10-05, Dharmacon) or 25 nM EP300 siRNA for four hours. Media was replaced with EBM-2 and RNA was harvested after 48 hours.

Quantitative Reverse Transcriptase-PCR (qRT-PCR) amplification was performed using a 7500 Fast Real-Time PCR System (Applied Biosystems). cDNA was synthesized from extracted RNA using a High-Capacity cDNA Reverse Transcription Kit using manufacturers protocol (Applied Biosystems). FAST SYBR Green Master Mix (Roche), cDNA template and 5 pmoles of forward and reverse primer were mixed to a final volume of 20 μl. Reactions were incubated for 10 minutes at 95°C, followed by 40 cycles of 95°C for 10 seconds and 60°C for 30 seconds. Gene expression changes were determined by normalizing against HPRT1.

RNA sequencing and analysis

Total RNA was extracted in TRIzol (Invitrogen) and purified using the Direct-zol RNA MiniPrep kit (cat. no: R2052, Zymo Research) according to the manufacturer's instructions. RNA was quantified on the MultiNA Bioanalyzer (Shimadzu). NEBNext Poly(A) mRNA Magnetic Isolation Module was used to enrich mRNA from 1 μg of total RNA. The NEBNext Ultra Directional RNA Library Prep Kit for Illumina was used to generate barcodes. Libraries were quantified on the MultiNA Bioanalyzer (Shimadzu) and pooled to equimolar ratios for sequencing. Cluster generation was performed at a concentration of 12 pM (TruSeq SR Cluster Kit v4-cBot-HS) and the flow cell was run on Illumina HiSeq2500 generating 100 nt reads. Fastq files were aligned to the hg19 human genome (ensembl GRCh37 release 75) using STAR.¹² The tool featureCounts²⁷ was used to generate a matrix of read counts per gene and the edgeR software was then used to determine differential gene expression.⁴⁰ Sex-linked genes were removed from the analysis. Significance is defined as FDR P value <0.05.

Conflict of interest

There is no conflict of interest to disclose.

Data access

Original RNA-seq data can be accessed from the NCBI Gene Expression Omnibus (GEO) with the accession number GSE77108.

Consent for publication

Not applicable.

Author contributions

H.R. and A.E. designed the study. H.R., J.O. and M.Z. carried out the experiments. H.R. performed the bioinformatics analysis with support from A.K. and M.Z. T.C.K. supported experimental design and manuscript writing. H.R. and A.E. wrote the manuscript.

Acknowledgments

We thank Dr. Ross Lazarus for bioinformatics support. We also acknowledge grant and fellowship support from the National Health and Medical Research Council (NHMRC). Supported in part by the Victorian Government's Operational Infrastructure Support Program.