Natural killer (NK) cells are highly specialized cytotoxic lymphocytes of the innate immune system that play a critical role in the control of microbial infections and several types of tumors (1). Unlike T or B cells, NK cells exert cytolytic functions without prior sensitization by utilizing germline-encoded, non-rearranged activating receptors (2, 3). They also modulate the function of other innate and adaptive immune cells by secreting several cytokines after activation (4, 5). Recently, key molecules involved in multiple signaling cascades downstream of major activating receptors were identified (6, 7). Moreover, the development and function of NK cells rely on the precise regulation of gene expression mediated by a complex matrix of transcription factors in response to the environmental cues (8, 9). Therefore, a better understanding of the molecular mechanisms underlying the transcriptional regulation of NK cells is important for developing new treatment strategies to combat infectious diseases and cancers.

Gene expression of eukaryotic cells is controlled by transcription factors (TFs) that bind to regulatory elements such as promoters and enhancers (10, 11). Moreover, epigenetic features, such as histone modifications and DNA methylation, are thought to play a crucial role in shaping chromatin structures and accessibilities of these regulatory elements (12). Enhancers are often located far from their target genes, and three-dimensional (3D) genome structure provides a mechanism to precisely regulate gene expression. The action range of enhancers is restricted by topologically associating domain (TAD) boundaries; intra-TAD chromatin looping facilitates the interaction of distal enhancers with their target promoters (13, 14).

Histone acetylation is a post-translational modification that can neutralize the positive charge of lysine residue, thereby reducing the interaction between the histone and negatively charged DNA, making DNA more accessible to the transcriptional machinery, and promoting gene expression (15). Histone acetylation is strictly regulated by the concerted activities of histone acetyltransferases (HATs) and histone deacetylases (HDACs) (16). Thus far, quite a few HDAC inhibitors have been reported, which have been shown to increase the acetylation of core histones, resulting in altered chromatin structure and gene expression (17).

Recent studies have demonstrated that the drug Vorinostat, an inhibitor of class I/II HDAC, could be used in a promising novel approach for treating leukemias and other malignancies (18). Moreover, Vorinostat has shown potential as an immunotherapeutic enhancer, improving the ability of NK cells to recognize tumor cells that have evaded immune surveillance (19, 20). However, it may also inhibit activating receptors’ expression and NK cell lytic activity against leukemic cells (21, 22). Therefore, further extensive studies are required to optimize the immunotherapeutic application of Vorinostat for the treatment of cancers and other diseases. Particular focus should be made on understanding how Vorinostat regulates the complex transcriptional network of NK cells in the 3D chromatin architecture. We analyzed the global transcriptional, epigenetic, and topological changes in NK cells in response to HDAC inhibitor and demonstrated the regulatory role that Vorinostat exerts in 3D enhancer network via acetylation of RUNX3 protein.

Genome-wide changes in gene expression and chromatin state of NK-92 cells in response to Vorinostat

Human NK cell line NK-92 cells display characteristics of activated NK cells in vivo and are currently in clinical trials as an off-the-shelf cellular immunotherapy for various cancers due to its prolonged activation, non-toxicity, low immunogenicity, and can be easily and rapidly expanded in culture (23). Exposing NK-92 cells to up to 1 μM Vorinostat for 24 hours did not significantly inhibit cell viability (Fig. 1A) or induce apoptosis (Fig. 1B). To understand how HDAC affects gene expression in NK cells, NK-92 cells were cultured for 24 hours in the absence/presence of 1 μM Vorinostat, and the global transcriptional, epigenetic, and topological changes occurring in NK-92 cells were analyzed (Fig. 1C). Upon Vorinostat treatment, differential RNA expression analysis revealed 1,829 upregulated and 1,991 downregulated genes (Fig. 1D). Gene ontology analysis indicated that the genes upregulated by Vorinostat were primarily enriched for immune-related pathways, while the downregulated genes were mostly associated with ribosome biogenesis and protein folding (Fig. 1E).

H3K27 acetylation (H3K27ac) is closely linked to the epigenetic control of gene expression and is a marker for active enhancers and promoters. To investigate the impact of Vorinostat treatment on H3K27ac levels, immunoblot assays were performed on total cell lysates, revealing higher levels of H3K27ac in the cells treated with Vorinostat (Fig. 1F). To further explore the effects of Vorinostat on enhancer activities, chromatin immunoprecipitation sequencing (ChIP-seq) was conducted to compare genome-wide H3K27ac occupancies (n = 2 per group) (Fig. 1G). Differential H3K27ac regions were identified: 7,061 were hyperacetylated and 785 were hypoacetylated by Vorinostat (Fig. 1G). We next performed assay for transposase-accessible chromatin with sequencing (ATAC-seq) to map the location and accessibility of regulatory elements genome-wide. Compared to the dynamic changes in H3K27ac levels, the degree of chromatin accessibility was marginally affected by Vorinostat in NK-92 cells (Fig. 1H). Representative genomic regions where Vorinostat induced changes in gene expression and chromatin state were shown in Fig. 1I.

Limited impact of Vorinostat on the 3D chromatin architecture in NK-92 cells

Mammalian genomes are packaged into intricate 3D structures consisting of hierarchical components such as compartments, TADs, and chromatin loops (24). Chromosomes separate into either active “A” or inactive “B” compartments, which can occur due to homotypic interactions between regions of the genome that share similar transcription and chromatin characteristics (25). TADs, the structural basis for organizing chromatin at a scale of hundreds of kilobases or below, are defined by preferential interactions between their chromatin (26, 27). We then investigated the extent to which Vorinostat affects higher-level genome organization of NK-92 cells by performing in situ Hi-C. Contact maps (Fig. 2A, Supplementary Table 1), the first eigenvector values from principal component analysis (Fig. 2B), and compartmentalization strength represented by saddle plots (Fig. 2C) indicated that Vorinostat had a negligible effect on segmentation of active and inactive chromosome domains into A and B compartments. Later, we defined TADs using insulation scores and found that the number of TADs (Fig. 2D) and the insulation capacity at TAD boundaries (Fig. 2E) were comparable between control and Vorinostat-treated cells. Moreover, genome-wide intra-TAD chromatin interactions (Fig. 2F) showed minimal changes upon Vorinostat treatment. We further explored intra-TAD chromatin looping by performing H3K27ac HiChIP in NK-92 cells. When statistically significant H3K27ac-based chromatin interactions were called using FitHiChIP, we identified a comparable number of high-confidence H3K27ac HiChIP loops (Fig. 2G, Supplementary Table 2). The major fraction (∼60%) of H3K27ac HiChIP loops represented chromatin interactions between promoter pairs, enhancer pairs, and promoter-enhancer pairs (Fig. 2G). Differential analysis for H3K27ac HiChIP loops demonstrated that the spatial proximity of nearly all H3K27ac HiChIP loops was preserved following Vorinostat treatment (Fig. 2H). Based on these topological data, it appears that Vorinostat has a limited impact on the three-dimensional chromatin architecture at the resolution of compartments, TADs, and enhancer-centric chromatin looping.

Vorinostat-induced RUNX3 acetylation reprograms enhancer landscape in NK-92 cells

To investigate the molecular mechanism underlying the effects of Vorinostat on gene expression, we analyzed TF binding motifs on ATAC-seq peaks that overlapped with differential H3K27ac peaks. Runt family motifs were enriched at hyperacetylated regions, while ERG and ETV family motifs were enriched at hypoacetylated regions in Vorinostat-treated cells (Fig. 3A). The Runt-related TF family comprises RUNX1, RUNX2, and RUNX3, which shares the evolutionarily conserved Runt domain. Tissue-specific and overlapping expression patterns indicate both exclusive and redundant roles for these three RUNX genes (28). Of these, RUNX3 is crucial in NK cell development and effector function, and considering its higher expression level compared to RUNX1 and RUNX2 in NK-92 cells (Fig. 3B) (28-31), we prioritized the analysis of RUNX3 and conducted ChIP-seq to examine their genome-wide binding pattern (Fig. 3C-F). The ChIP-seq analysis uncovered the presence of approximately 1,500 distinct genomic regions enriched with RUNX3 peaks, distributed across diverse genomic loci such as promoters, introns, and intergenic regions (Fig. 3D). While the levels of RUNX3 occupancy were similar in NK-92 cells cultured with or without Vorinostat (Fig. 3E, F), there was a notable increase in H3K27ac enrichment at the sites where RUNX3 was bound in Vorinostat-treated cells compared to control cells (Fig. 3G, H). Specifically, we identified a subset of RUNX3 binding sites (n = 604) with significantly elevated H3K27ac levels upon Vorinostat treatment, while only 36 RUNX3 peaks showed reduced H3K27ac levels (Fig. 3H). We next examined whether the increase in H3K27ac levels upon Vorinostat treatment was associated with posttranslational modifications of RUNX3. Previous studies have shown that acetylation plays a crucial role in enhancing RUNX3-dependent transactivation activity (32). The acetylation of RUNX3 by the histone acetyltransferase p300 is critical for the formation of a complex with bromodomain-containing protein 2 (BRD2), which in turn facilitates the recruitment of SWI/SNF and TFIID (33, 34). To test this hypothesis, we immunoprecipitated lysates from NK-92 cells using anti-acetyl-lysine antibodies and performed Western blotting with an anti-RUNX3 antibody. Our results indicate that Vorinostat treatment increases the acetylation level of RUNX3 protein, as shown in Fig. 3I.

RUNX3 acetylation enhances immune response-related gene expression via long-range enhancer-promoter chromatin interactions

Next, we sought to identify candidate genes whose expression could be regulated by Vorinostat-induced RUNX3 acetylation (i.e., activated RUNX3). To do so, we analyzed RUNX3 binding sites with enhanced H3K27ac levels and used H3K27ac HiChIP loops to link the activated enhancers with their target genes (Fig. 4A). Remarkably, only a few genes were directly bound by activated RUNX3 at their promoter regions, while most RUNX3 target genes were regulated by distal enhancers through long-range chromatin interactions (Fig. 4A). Analysis of the differential RNA expression of genes associated with activated RUNX3 found that 154 genes were deregulated by Vorinostat treatment (FDR < 0.05, FC > 1.3). Among these, 111 and 43 genes were upregulated and downregulated, respectively (Fig. 4B, Supplementary Table 3). Notably, gene ontology (GO) analysis revealed that the upregulated RUNX3 target genes (n = 111) were primarily enriched in immune response and cytokine production pathways (Fig. 4C). Specifically, Vorinostat treatment upregulated genes involved in transcription and signaling activity, such as NFKBIZ and MXI1, mediated by long-range enhancer-promoter chromatin interactions connecting their promoters with activated RUNX3-dependent distal enhancers (Fig. 4D, E).

The essential role of enhancers in the control of gene transcription during NK cell differentiation, activation, and function has been extensively studied in recent years (35, 36). However, there is limited research on the 3D chromatin architecture in NK cells and how it relates to enhancer activity (37). In this study, we systematically mapped the landscape of enhancer histone acetylation, chromatin accessibilities, higher-order chromatin structure, and enhancer-promoter regulatory interactions in human natural killer cells (Fig. 1C) to characterize the underlying mechanisms of the effect of Vorinostat on the gene expression profile. The ChIP-seq and ATAC-seq analysis revealed that HDAC inhibition by Vorinostat dramatically reprogrammed enhancer activities through changes in histone acetylation rather than chromatin accessibility in NK-92 cells (Fig. 1F-H). In contrast, in situ Hi-C analysis demonstrated that Vorinostat treatment did not induce changes in global genome organization at the compartment and TAD levels (Fig. 2A-F). Interestingly, Vorinostat-induced enhancer activation did not change spatial proximity between enhancer and promoter, given that most H3K27ac HiChIP loops remain stable after Vorinostat treatment (Fig. 2G, H). Taken together, gene expression changes following Vorinostat-mediated HDAC inhibition are controlled mainly at the level of enhancer histone acetylation rather than 3D chromatin architecture or chromatin accessibility.

Our analysis also revealed an overrepresentation of binding motifs for RUNX family transcription factors in Vorinostat-induced activated enhancer regions, with RUNX3 having the highest expression level in Vorinostat-induced NK-92 cells (Fig. 3A, B). Previous studies have demonstrated that expression of RUNX3 increases with NK maturation, and defects in its expression or activity can impede NK cell development (28, 38). Additionally, RUNX3-deficiency alters expression of genes involved in proliferation, survival, migration, and effector function of NK cells (39). The post-translational regulation of RUNX3’s transcriptional activity involves acetylation by p300 acetyltransferase, which protects it from degradation by the ubiquitin ligase Smurf and facilitates its interaction with BRD2 (32-34). The level of RUNX3 acetylation is influenced by environmental cues and is downregulated by HDAC activity. Therefore, RUNX3 may function as a sensor of HDAC inhibitor Vorinostat, and its acetylation potentially induces transcriptional reprogramming by elevating enhancer activities (Fig. 4F). However, it cannot be ruled out that other RUNX family proteins may exhibit compensatory effects on Vorinostat-induced enhancer reprogramming. To regulate its target genes, activated RUNX3 primarily depends on long-range chromatin interactions (Fig. 4A). The target genes upregulated by RUNX3 were significantly enriched in pathways related to immune response (Fig. 4B, C). Thus, RUNX3 acetylation induced by Vorinostat elevates the expression of immune response-related genes of NK cells through a spatial enhancer network.

It is important to note that mRNA expression levels are often used as a proxy for protein abundance. However, they do not always reflect the true abundance of the corresponding protein due to various post-transcriptional regulatory mechanisms that can affect mRNA stability, splicing, and translation (40). Given that Vorinostat downregulates the expression of genes involved in protein folding (Fig. 1G), further studies are warranted to investigate the regulatory effect of Vorinostat on the actual effector function of NK cells.

In conclusion, our study provides a comprehensive analysis of the molecular mechanisms underlying the effect of HDAC inhibitors on the transcriptional profile of NK cells. Our findings shed light on the intricate interplay between 3D chromatin structure, histone modifications, and gene expression in NK cells and offer insights into regulatory mechanisms that control NK cell function. This understanding may facilitate the development of novel therapeutic approaches to treat diseases resulting from dysregulated NK cell function.

Materials and methods are available in supplementary material.

This work was supported by the National Research Foundation of Korea (NRF) grants funded by the Korean government (MSIP) (2018M3A9D3079290, 2020R1A2C2013258, and 2022M3A9 B6017424 to H.-P. Kim).

The authors have no conflicting interests.