Changes in chromatin accessibility landscape and histone H3 core acetylation during valproic acid-induced differentiation of embryonic stem cells

Directed differentiation of mouse embryonic stem cells (mESCs) or induced pluripotent stem cells (iPSCs) provides powerful models to dissect the molecular mechanisms leading to the formation of specific cell lineages. Treatment with histone deacetylase inhibitors can significantly enhance the efficiency of directed differentiation. However, the mechanisms are not well understood. Here, we use CUT&RUN in combination with ATAC-seq to determine changes in both histone modifications and genome-wide chromatin accessibility following valproic acid (VPA) exposure. VPA induced a significant increase in global histone H3 acetylation (H3K56ac), a core histone modification affecting nucleosome stability, as well as enrichment at loci associated with cytoskeletal organization and cellular morphogenesis. In addition, VPA altered the levels of linker histone H1 subtypes and the total histone H1/nucleosome ratio indicative of initial differentiation events. Notably, ATAC-seq analysis revealed changes in chromatin accessibility of genes involved in regulation of CDK serine/threonine kinase activity and DNA duplex unwinding. Importantly, changes in chromatin accessibility were evident at several key genomic loci, such as the pluripotency factor Lefty, cardiac muscle troponin Tnnt2, and the homeodomain factor Hopx, which play critical roles in cardiomyocyte differentiation. Massive parallel transcription factor (TF) footprinting also indicates an increased occupancy of TFs involved in differentiation toward mesoderm and endoderm lineages and a loss of footprints of POU5F1/SOX2 pluripotency factors following VPA treatment. Our results provide the first genome-wide analysis of the chromatin landscape following VPA-induced differentiation in mESCs and provide new mechanistic insight into the intricate molecular processes that govern departure from pluripotency and early lineage commitment. Supplementary Information The online version contains supplementary material available at 10.1186/s13072-021-00432-5.

subject to developmental regulation. For instance, in mouse embryonic stem cells, chromatin is considered 'open, ' hyperdynamic and transcriptionally hyperactive compared to differentiated cells [5][6][7]. Chromatin remodeling alters nucleosome structure and function to accommodate various cellular processes, such as cell cycle progression, growth, and differentiation, and facilitates the response to environmental or pharmacological stimuli. Core histones as well as linker histone H1 are thereby subject to multiple histone posttranslational modifications (PTMs) that confer changes to the DNA binding capacity as well as interactions with other protein complexes, either consolidating or diminishing their inherent transcriptional repressive attributes. Thus, PTMs are essential to the biochemical mechanisms that remodel chromatin to accommodate fluctuating physiological requirements [8,9].
Histone acetylation is of utmost importance for the regulation of chromatin accessibility with far-reaching implications for transcriptional activity, DNA repair, and replication [10][11][12] and requires a dynamic balance between the addition of acetyl groups to lysine residues through the activity of histone acetyltransferases (HATs) and their removal catalyzed by histone deacetylases (HDACs) [13][14][15]. Disruptions to these homeostatic mechanisms are commonly associated with pathophysiological changes, particularly in several types of human cancers that show overexpression of HDAC enzymes [16][17][18]. Selective HDAC inhibitors (HDACi), such as valproic acid (VPA), an FDA approved class I deacetylase inhibitor, are therefore important pharmacological compounds for novel epigenetic cancer therapies [19].
Regulating the levels of global histone acetylation is also a powerful strategy to improve the efficiency of human and mouse embryonic stem cell reprogramming. Small molecule HDACi, such as valproic acid [20,21], are of critical importance in facilitating the efficiency of reprogramming of somatic cells into induced pluripotent stem cells by more than 100-fold and to promote transcription factor free reprogramming [20,[22][23][24][25], and disruption of histone deacetylation interferes with maintenance of pluripotency in mESCs [26]. In addition, differentiation of human adipose-derived stem cells into cardiomyocyte-like cells that express bona fide cardiac marker proteins can be enhanced by VPA supplementation [27]. Collectively, these effects have been suggested to rely on changes induced by histone hyperacetylation [28]. However, the precise mechanisms of action and molecular pathways involved are not fully understood and little is known regarding the effects of VPA treatment on chromatin structure and the regulation of specific genome-wide changes in the chromatin landscape during induced differentiation of embryonic stem cells.
We show here that mESCs treated with VPA undergo striking changes in both nuclear accumulation as well as genome-wide occupancy of histone H3 acetylated at lysine 56 (H3K56ac), an important lateral surface PTM that increases the rate of local DNA unwrapping at the DNA entry/exit site of the nucleosome [29][30][31]. Highresolution confocal microscopy and high-content image analyses demonstrate a global increase in H3K56ac nuclear staining per cell as well as an increase in the number of cells with elevated H3K56ac nuclear accumulation. In addition, chromatin profiling with CUT&RUN revealed an increase in locus-specific H3K56ac genome occupancy. Changes in H3K56 acetylation were associated with differential expression levels of linker histone H1 subtypes as well as critical changes in DNA accessibility to nicking enzymes. ATAC-seq analysis demonstrated genome-wide changes to the chromatin landscape following VPA-induced differentiation. Importantly, our results indicate that changes in H3K56ac promoter binding occur concomitantly with changes in chromatin accessibility at key genomic loci involved in pluripotency as well as cardiomyocyte differentiation. Our studies provide novel insight into the specific molecular changes to the chromatin environment brought about by VPA treatment and inform the critical intersections of remodeling processes, such as changes in histone acetylation, linker histone expression, and locus-specific chromatin accessibility during early lineage commitment.

Valproic acid induces global changes in histone H3K56 acetylation
Treatment with VPA (2 mM) for 48 h significantly enhances the efficiency of reprogramming into cardiomyocyte cell lineages [28] as well as neuronal precursors [32]. This process has been linked to changes in histone acetylation [28]. However, the molecular mechanisms involved are not known. To determine the upstream epigenetic changes associated with the onset of directed differentiation of ES cells, we used high-content image analysis [33] to compare the patterns of global histone H3K56 acetylation in control cells with the patterns observed after VPA exposure. Confocal high-content image analysis following immunochemical detection of H3K56ac (green) revealed basal levels of nuclear staining in control cells (Fig. 1A). However, treatment with VPA (2 mM) induced a striking upregulation of H3K56ac nuclear localization within 48 h in karyotypically stable, female PGK12.1 mESCs compared to untreated control cultures (Fig. 1A). This was also confirmed by Western blotting (Additional file 1: Fig. S1) and line scan analysis using high-resolution confocal microscopy of control and VPA-treated mESCs, demonstrating elevated H3K56ac fluorescence intensity (green) within euchromatic and to some extent heterochromatic chromocenter domains (shaded blue) (Fig. 1B), while major changes were not detectable for several histone N-terminal tail modifications, including H4K5ac and H4K8ac (Additional file 1: Fig. S1).
To quantify the degree of H3K56 hyperacetylation in response to VPA treatment in large cell populations, we used threshold maps to segment individual cell nuclei in control and VPA-treated cultures and measured fluorescence intensities using an automated system for high-content confocal image analysis (ImageXpress, Molecular Devices, LLC). This analysis demonstrated that VPA treatment induces a significant increase in H3K56ac fluorescence intensity, while the majority of control mESCs show only baseline levels of acetylation (Fig. 1C). Nuclear circumference measurements revealed no significant differences between cultures (Fig. 1C), suggesting that VPA treatment does not induce detectable changes in nuclear size. These results indicate that VPA treatment causes global hyperacetylation of an atypical histone posttranslational modification that plays a critical role in the regulation of nucleosome stability and higher-order chromatin structure [34]. Notably, hyperacetylation of H3K56 has been previously observed in ES cells carrying HDAC1 knockout alleles [35], indicating a possible molecular mechanism by which VPA may regulate H3K56 acetylation levels. However, whether VPAinduced changes in H3K56ac are associated with altered levels of chromatin compaction is not known.

Profile of linker histone H1 subtypes during induced differentiation of ES cells
The 'open' chromatin status of mouse embryonic stem cells is thought to confer flexibility toward dynamic transitions in nuclear architecture in response to developmental or environmental cues [7,36]. In turn, higher-order chromatin compaction is essential for the epigenetic regulation of pluripotent stem cell differentiation and critically depends on appropriate nuclear levels of individual linker histone H1 subtypes [37,38]. Linker histones bind nucleosome core particles at the site of DNA entry/exit to confer supra-nucleosomal chromatin structure. Due to the proximity of hyperacetylated lysine 56 within the globular domain of histone H3 to the site of DNA entry/exit, histone H1 expression patterns may provide critical insight into changes to higher-order chromatin structure and elucidate potential changes in the expression patterns of specific histone H1 subtypes in response to VPA.
We first conducted reverse-phase HPLC experiments of isolated histones to determine the protein levels of individual H1 subtypes. Notably, H1 0 was nearly undetectable in the chromatograms of control mESCs. However, this subtype showed a significantly larger peak in VPA-treated cells ( Fig. 2A). To analyze the H1d/H1e fraction collected from HPLC eluates in more detail, this fraction was subjected to mass spectrometry and revealed a striking increase in the subtype H1e peak following VPA exposure ( Fig. 2A, inset). Quantitative analyses revealed that the total histone H1/nucleosome ratio rose by 15% from 0.58 in control to 0.67 following valproic acid exposure (Fig. 2B). This increase in total H1 was not attributable to one particular subtype, but rather the result of composite differential expression patterns. For instance, subtypes H1a, H1b, H1e, and H1 0 all showed significantly increased ratios (P < 0.05) in VPAtreated cells, while the H1d/nucleosome ratio was significantly reduced and H1c showed no differences between groups (Fig. 2B). These results demonstrate that treatment with histone deacetylase inhibitors, such as VPA, has the potential to affect total linker H1 levels and alters subtype-specific expression patterns in mESCs that may impact chromatin structure and function.

Changes in the genome-wide chromatin landscape during induced differentiation of ES cells
Hyperacetylation of histone tail residues reduces the electrostatic affinity between histone proteins and DNA, and thereby promotes a more 'open' chromatin structure that is permissive to active gene transcription [8,15]. Due to the proximity of the lysine 56 residue of the histone H3 globular domain to the DNA entry/exit point of the dyad at the lateral surface of the nucleosome, its acetylation has the capacity to directly influence DNA-histone interactions [39][40][41]. Global hyperacetylation of H3K56ac may, thus, contribute to changes in mesoscale chromatin (See figure on next page.) Fig. 1 VPA treatment induces global histone H3K56 hyperacetylation. A Treatment of mESCs with 2 mM VPA for 48 h increases the levels of H3K56 acetylation (green) compared to untreated control cultures. Scale bar = 20 μm. B High-resolution confocal imaging of individual control and VPA-treated nuclei stained with H3K56ac (green). DNA is counterstained using DAPI (blue). Fluorescence intensity across nuclei is visualized by line scan (red arrows) graphs. DAPI bright chromocenters are highlighted (shaded blue). C High-content confocal analysis using threshold maps for the quantification of H3K56ac fluorescence intensity and nuclear circumference reveals a higher proportion of nuclei with bright H3K56ac staining (yellow) in VPA-treated cultures compared to controls. Nuclear circumference measurements indicate no significant differences (Cntl: n = 98; VPA: n = 108). Scale bar = 20 μm. Data are expressed as the mean ± SD of three independent replicates. Unpaired t-test, two tailed resulted in P < 0.001. organization or affect the genome-wide chromatin landscape in VPA-treated mESCs.
To determine whether altered expression of Histone H1 subtypes and changes in H3K56 acetylation are associated with genome-wide or locus-specific changes in chromatin accessibility, we next conducted assays for transposase-accessible chromatin using sequencing (ATAC-seq) in control and VPA-treated mESCs in two biological replicates (Additional file 1: Fig. S2). Proportional ATAC-seq peak distributions were essentially indistinguishable relative to genomic features in VPAtreated and control mESCs (Additional file 1: Fig. S3) and heatmaps of tag distributions across transcription start sites (TSSs) (Fig. 3A), merged regions, or genebodies (Additional file 1: Fig. S3) exhibited no major differences between treatment groups. Notably, differential ATACseq peaks (P < 0.001) were detected in 5537 genomic loci between control and VPA-treated cells (Fig. 3B), with 3923 loci showing reduced (red) and 1614 loci exhibiting increased (green) accessibility. Comparative analysis of The inset shows the relative signal intensity of H1d and H1e mass spectral peaks in the H1d/H1e fraction collected from HPLC eluates. Prominent changes are indicated in red arrows. B Quantification of H1/nucleosome ratios of individual H1 subtypes as well as of total H1 in control versus VPA-treated mESCs as determined by RP-HPLC and Mass Spectrometry. Data are expressed as the mean ± SD over three independent replicates. Unpaired t-test, two tailed; *P < 0.05, **P < 0.01, ***P < 0.005, ****P < 0.001 peak numbers demonstrated significantly (P < 0.05) fewer peaks associated with specific genomic features, such as exons, 3'UTRs, proximal and distal gene loci, promoters, and genebodies following VPA treatment for 48 h (Additional file 1: Fig. S3). 5'UTRs, CpG islands, as well as distal intergenic regions exhibited no significant differences. These results illustrate that, rather than inducing a unidirectional increase in chromatin accessibility, VPA treatment leads to a complex chromatin response involving both opening and closing of genomic loci following VPA exposure. VPA treatment is, moreover, associated with changes in chromatin accessibility not only at specific protein coding regions but also at lncRNA loci (Additional file 1: Fig. S3) and various classes of transposable elements (TE), such as long terminal repeats (LTRs) and short interspersed nuclear elements (SINEs) (Fig. 3C).
Importantly, the top 10 genomic loci with VPAinduced gain of chromatin accessibility encode for factors involved in chromatin remodeling, cardiomyocyte, and neuronal cell differentiation (Fig. 3D). For example, Troponin T3 (Tnnt3) as well as the Coiled-Coil Domain Containing Protein 80 (CCdc80), Alpha-3-Catenin (Ctnna3), and the Polycomb Group Ring Finger Protein (Pcgf5) play critical roles in cardiac/myogenic cell lineages and are crucial for striated muscle contraction [42][43][44][45][46]. Additionally, Zinc Finger Matrin-Type 4 (Zmat4) and the Paternally Expressed Gene 10 (Peg10) are implicated in mammalian neurogenesis [47,48], suggesting that VPA treatment induces changes in chromatin accessibility that favor onset of cell differentiation. This notion is supported by the fact that VPA exposure led to a net loss of ATAC-seq peaks and hence decreased accessibility at 3923 loci ( Conversely, a striking loss of chromatin accessibility was notable at transcription start sites of the bona fide pluripotency marker Pou5f1 (Oct4) locus that was also associated with a significant reduction in protein expression as detected by immunochemistry and Western blotting (Fig. 4A). Other pluripotency genes, such as Nanog (Fig. 4A) and cMyc, also showed significantly reduced accessibility at their respective TSSs, and the Sox2 locus demonstrated loss of ATAC-seq peaks at a putative enhancer upstream of the TSS (data not shown), further suggesting initiation of differentiation and departure from pluripotency in response to VPA exposure. Moreover, we observed changes in chromatin accessibility at putative enhancer sites upstream of the Tnnt2 locus, encoding cardiac dominant Troponin T2 (Fig. 4A), as well as in several regions, including the TSS of the homeodomain-only protein homeobox, Hopx (Fig. 4A), a crucial regulator of cardiac cell proliferation and differentiation [49] that functions through interactions with HDAC2 downstream of Nkx2.5 to mediate repression of cardiac-specific genes [50]. These results demonstrate that VPA treatment leads to altered chromatin accessibility at specific genomic loci associated with loss of pluripotency and induction of cell differentiation and lineage commitment, favoring cardiac and neuronal pathways.
Digital genomic footprinting (DGF) has recently emerged as valuable approach to overcome limitations of ChIP-based methods [51] and allows for the massive parallel interrogation of ATAC-seq datasets to predict differential binding of transcription factors. Accordingly, bound TFs, similar to nucleosomes, hinder DNA cleavage resulting in characteristic footprints-regions with diminished signal strength within a region of higher signals [51]. DGF analysis of our ATAC-seq datasets revealed a striking switch in TF occupancy from a chromatin landscape dominated by pluripotency-associated factors, such as POU domain/Sox-related factors, like OCT4 and SOX2, to a chromatin environment that is controlled by TF binding characteristic for differentiating cells (Fig. 4B). A complete list of TFs with enrichment in control and VPA-treated cells, respectively, is shown in Additional file 1: Tables S1 and S2. For instance, Jun/Fos and Smad2/3 transcription factor occupancy increased strikingly following VPA exposure (Fig. 4B), suggesting initiation of differentiation programs into mesoderm and/or endoderm lineages [52]. Our DGF findings, thus, further complement the notion that VPA treatment induces complex changes in both chromatin accessibility at coding genes as well as TF occupancy that support the departure form pluripotency and entry into lineage commitment.

Changes in H3K56 acetylation genome occupancy during VPA-induced differentiation of embryonic stem cells
Acetylation of H3K56 is linked to the core transcriptional network in human and mouse embryonic stem cells and is enriched within nucleosomes at promoter sites of actively transcribed genes regulated by at least one of the pluripotency factors, SOX2, NANOG, or OCT4 [53,54]. Importantly, in human ESCs, H3K56 is also acetylated at poised promoters bound by RNA Pol-II in a manner that allows switching from pluripotency genes to promoters of tissue-specific developmental genes [53]. However, the effects of VPA-induced hyperacetylation on genome-wide or locus-specific occupancy of H3K56ac in mouse embryonic stem cells is not known.
To undertake a genome-wide analysis of the effects of HDAC inhibition on this critical epigenetic modification, we conducted Cleavage Under Targets and Release Using Nuclease (CUT&RUN) assays. CUT&RUN is an efficient and robust strategy for quantitative mapping of protein-DNA interactions and local chromatin environment analysis [55]. First, we used an antibody against H3K27me3 to validate our experimental procedure on mESCs. CUT&RUN resulted in specific enrichment of reads in bona fide H3K27me loci, such as the Hoxd cluster and the Wnt5a locus (Additional file 1: Fig. S4) [56][57][58]. We used Peak calling by Sparse Enrichment Analysis for CUT&RUN (SEACR) to identify peaks and enriched regions and validated our results with previously published data [59]. Next, we set out to conduct CUT&RUN assays using an antibody against H3K56ac in control and VPA-treated mESCs. Notably, while peak widths were indistinguishable between control and VPAtreated samples, H3K56ac CUT&RUN revealed higher mean peak numbers (135,194) in VPA-treated samples compared to controls (130,907; Additional file 1: Fig. S4). Chromatin immunoprecipitation experiments have previously shown an association of H3K56ac with the Thy1 locus in mouse embryonic fibroblasts [60], a locus that was also found enriched for H3K56ac in control mESCs in our CUT&RUN assays (Additional file 1: Fig. S4). Representative signal alignment heatmaps visualizing the midpoints of signal blocks returned form SEACR analysis in control and VPA-treated cells show similar alignment patterns (Fig. 5A). In contrast, analysis of annotated genomic regions with H3K56ac enrichment revealed a 22% increase in peaks within genes in VPA-treated cells compared to control mESCs. A large majority of peaks was found within gene bodies, although promoter regions made up to 6.5% and 5.7% in control and VPA-treated samples, respectively (Fig. 5B). Collectively, these findings demonstrate that H3K56ac hyperacetylation observed by immunocytochemistry is mirrored by increased genome-wide mapping of this epigenetic mark following HDAC inhibition. However, data also reveal that H3K56ac peaks, albeit globally increased in VPA-treated cells, show gain or loss of enrichment in a locus-dependent manner. Functional annotation shows an over-representation of GO Terms associated with cell morphogenesis, cytoskeleton organization, and neuronal differentiation among H3K56ac-enriched loci in VPAtreated samples (Fig. 5C, top 15 processes out of 659 shown), suggesting that locus-specific gain of H3K56ac may be linked to local chromatin remodeling processes in response to VPA-induced directed differentiation. For example, enrichment of H3K56ac was observed within the promoter region of Pax6 (Fig. 5D), suggesting potentially increased nucleosome dynamics and facilitated/ poised expression of this important developmental regulator [61]. A striking reduction in H3K56ac across the entire gene body was evident across the genomic region encoding the histone chaperone ASF1A involved in deposition of histone H3 variants (Fig. 5E). Intriguingly, ablation of ASF1A expression by short hairpin RNA interference has been shown to decrease the levels of pluripotency-related markers and increases differentiation-associated factors in mESCs [54]. In addition, ASF1A is also essential for maintenance of pluripotency and reprogramming of human stem cells [62]. It is thus plausible that the loss of H3K56ac at the Asf1a locus observed in our study may be important to decrease the expression of pluripotency markers during induced differentiation.
To investigate whether chromatin accessibility patterns in VPA-treated mESCs correlate with changes in H3K56ac, we performed integrative analyses. Interestingly, we found overlapping gain in chromatin accessibility and H3K56ac enrichment within the TSS of  11 coding genes ( Fig. 6A and Additional file 1: Fig.  S5), while the opposite pattern, i.e., loss of accessibility associated with loss of H3K56ac enrichment was detectable at the TSS of 73 genes, revealing a gene regulatory network that is subject to correlative changes in accessibility and H3K56 acetylation. Enriched GO Terms involving loci gaining accessibility and H3K56ac encompassed critical myogenic processes, such as skeletal muscle contraction, sarcomere, and cytoskeletal organization as well as regulation of chromatin integrity (Fig. 6B). Notably, the fast-twitching muscle troponin T3 (Tnnt3) encoding locus, a major genomic target undergoing VPA-induced chromatin opening at the TSS, also showed a correlative enrichment of H3K56ac as demonstrated by SEACR analysis of CUT&RUN profiles (Fig. 6C). Similar correlative changes in chromatin accessibility and H3K56ac occupancy were also observed at the Lefty locus (Fig. 7A), that encodes two members of the transforming growth factor (TGF)-β superfamily, Lefty 1 and Lefty 2. Lefty ligands are known to play crucial roles in the maintenance of pluripotency in mouse and human ESCs [63][64][65]. Interestingly, conserved enhancer features have been described in mouse and human synthenic regions of the Lefty locus [64,65], that are characterized by binding of multiple pluripotency factors, such as Nanog, OCT4, SOX2, KLF4, and KLF2 in ESCs [65]. Our analyses demonstrate a striking loss of chromatin accessibility at these enhancer regions and the TSS within the Lefty locus that are also accompanied by a reduction in H3K56ac as well as OCT4 occupancy following VPA treatment as indicated by CUT&RUN and ChIP-qPCR analyses (Fig. 7A, C). These findings indicate changes in the chromatin landscape at the Lefty locus that may contribute to departure from pluripotency upon VPA-induced directed differentiation as depicted in the model in Fig. 7B. Collectively, our genome-wide analyses demonstrate that VPA treatment-induced directed differentiation is associated with significant changes in locus-specific chromatin accessibility patterns that encompass genes encoding pluripotency factors as well as genes involved in cell cycle control, and myogenic and neurogenic differentiation, thus providing critical molecular insight into the mechanisms by which HDAC inhibitors induce directed differentiation of embryonic stem cells.

Discussion
Modulating the levels of global histone acetylation is a powerful strategy to improve the efficiency of human and mouse embryonic stem cell reprogramming. Both small molecule histone deacetylase inhibitors, such as valproic acid [20,21], as well as biophysical signals transmitted to the nucleus by topographical and mechanical properties of functionalized biomaterials enhance the efficiency of ESC reprogramming into cardiomyocyte cell lineages [28] and neuronal precursors [32]. mESCs exhibit profound changes in histone modifications and chromatin structure during epigenetic reprogramming [66]. However, little is known regarding the mechanisms of reprogramming into distinct cell lineages and the effects on the chromatin landscape.
Here, we provide a genome-wide map of epigenetic signatures and changes in chromatin accessibility during induced differentiation of mESCs. Our results indicate that valproic acid induces dramatic changes in both global H3K56 acetylation in the nucleus as well as H3K56ac genome occupancy at promoter regions and gene bodies of factors associated with cytoskeleton organization, cell junction assembly, and cellular morphogenesis. ATACseq analysis revealed significant changes in genome-wide chromatin accessibility with the top 10 changes detected in genes involved in cardiomyocyte differentiation, chromatin remodeling, and neuronal differentiation. Consistent with the role of histone acetylation in chromatin structure, we detected 1614 genes that gained chromatin accessibility following VPA exposure. However, we also observed a loss of chromatin accessibility in 3923 loci. Notably, ontology analysis revealed changes in genes involved in CDK serine/threonine kinase activity and DNA duplex unwinding. Moreover, we found evidence for a significant increase in the levels of several linker histone subtypes, such as H1 0 and H1e, and the total H1/nucleosome ratio consistent with chromatin remodeling events.
Human and mouse ESCs are characterized by the presence of several atypical PTMs, such as H3K56ac [53,67]. In contrast with the majority of histone acetylation marks that occur at histone tails, H3K56ac affects the core (globular) domain of histone H3 and has the capability for altering nucleosome organization and chromatin structure [34]. A previous study used mass spectrometry and ChIP-on-Chip to demonstrate that H3K56ac is expressed in human ESCs, where it binds regulatory sequences at pluripotency genes [53]. However, the patterns of genome-wide occupancy of this PTM in response to chemically induced differentiation in mouse ESCs remained to be determined. Consistent with previous studies [53,54], we found evidence for the binding of H3K56ac at promoter regions of Pou5f1 (Oct4) and Nanog. In addition, we detected binding of H3K56ac at enhancer and promoter sites of the Lefty locus in mouse ES cells, also encoding pluripotency factors. However, our study did not reveal binding at regulatory regions of Sox2.
Notably, although OCT4 exhibits a direct interaction with H3K56ac, Sox2 does not co-immunoprecipitate with H3K56ac in mouse ESCs [54]. Our results suggest that the genome occupancy of this histone mark at pluripotency factors has been mostly conserved in mouse ESCs. Importantly, VPA induced elevated occupancy at genebodies and promoter regions of factors involved in cytoskeleton organization, cell junction assembly, and cell morphogenesis, suggesting that changes in histone acetylation of multiple cytoskeleton-related genes are an important component of the cellular reprogramming mechanisms activated by both chemically induced reprogramming of mESCs (this study) as well as biophysical properties of functionalized biomaterials [68]. Notably, retinoic acid-induced differentiation of human ESCs also caused a redistribution of H3K56ac occupancy from pluripotency genes toward factors involved in organ morphogenesis, signal transduction, and tissue remodeling [53]. Thus, changes in global H3K56ac and its genome occupancy at cytoskeleton-related genes constitute an early response to VPA-induced differentiation. Enrichment of H3K56ac within the promoter region of the neuroectoderm cell fate determinant Pax6 is in agreement with previous studies demonstrating that, aside from an activating role in neural gene expression, Pax6 is also involved in the repression of pluripotency genes [61]. Loss of H3K56 binding in response to VPA treatment, in contrast, supports the notion of a complex regulatory, and context-dependent response in mESC. For example, the histone chaperone ASF1A is essential for the maintenance of pluripotency, and changes in Asf1 expression in mouse embryos alter their developmental potential [69]. However, whether H3K56ac-induced changes to ASF1A expression would result in differences in H3 variant incorporation or affect genome stability remains to be determined.
Linker histone H1 is a major chromatin architectural protein that plays crucial roles in global chromatin organization and nucleosome stability [70]. The potential of H3K56ac to affect nucleosome structure led us to determine the effect on linker histone subtypes. VPA treatment altered the expression patterns of histone H1. For instance, we observed significant differences in both the total H1 to nucleosome ratio as well as that of several individual H1 subtypes. Of these differentially expressed subtypes, all but H1c and H1d exhibited significantly higher ratios, suggesting that VPA treatment leads to shifts in H1 subtype expression and may impact nucleosome structure. We previously reported that differentiation of mESCs into embryoid bodies for up to 10 days is associated with increasing total H1/nucleosome ratios from approximately 0.46 in mESCs at d 0 to 0.55 (day 7) and 0.62 (day 10) [37]. Reverse-phase HPLC analyses in our current study revealed ratios of 0.58 in control undifferentiated mESCs and 0.67 in VPA-treated cultures after 48 h and show significant increases in response to a 48 h VPA exposure with total H1/nucleosome ratios superseding the levels previously observed following 10 days of embryoid body (EB) differentiation [37]. However, in contrast to EB differentiation, H1c protein levels remained unchanged and H1d levels declined in mESCs upon VPA treatment. Based on these findings we conclude that a brief supplementation of culture media with VPA is sufficient to induce changes in H1/nucleosome ratios at levels at least comparable to those achieved following a 10-day standard differentiation protocol, with distinct changes in profiles of individual H1 subtypes. For instance, previous reports have shown a progressive increase of H1c, H1d, H1e, and H1 0 during ESC differentiation [37], whereas we find a rapid increase of H1a, H1b, H1e, and H1 0 upon acute, highdose VPA treatment, suggesting that EB differentiation and VPA-mediated directed differentiation exert distinct effects on the expression of specific H1 subtypes. Notably, changes in H1 subtype expression levels were not associated with altered H3K56ac occupancy at their respective loci, suggesting that H3K56ac is not directly involved in the regulation of H1 subtype expression in response to VPA. This notion is supported by reports indicating that other HDAC inhibitors, such as Trichostatin A (TSA) and sodium butyrate, have the capacity to induce expression of the differentiation-specific histone H1 0 in mouse and human tumor-derived cell lines [71][72][73][74][75]. Recent fluorescence spectroscopy studies have also identified direct interactions between VPA, histone H1, and chromatin, and demonstrated that VPA can bind to H1 with even higher affinity than to chromatin [76]. Using high-performance polarization microscopy and Fourier-transform infrared (FTIR) microspectroscopy, studies have shown that such interactions of VPA extend to involve histone H1 and H3 conformations and suggested that VPA induces changes to the suprastructure of helical DNA aggregates during crystallization [77]. Thus, VPA treatment-induced changes in histone H1 expression presented here provide support for the concept of a wide-spectrum pharmacological potential of VPA in exerting activity not only indirectly by inhibiting class I HDAC enzymes but also directly on chromatin structure through effects on major chromatin architectural proteins, such as histone H1 and H3, during induced differentiation. Our results also provide mechanistic insight into the pathways regulating chromatin remodeling and epigenetic reprogramming during VPA-induced differentiation of ES cells on a locus-specific scale. Consistent with a role in changing the chromatin landscape, we detected changes in chromatin accessibility within genes involved in regulation of CDK serine/threonine kinase activity and DNA duplex unwinding. For example, an increase in chromatin accessibility was evident in the Annexin-1 locus. However, we observed a loss of accessibility at regulatory regions within the Bloom helicase (Blm), the suppressor of variegation RNA helicase (Supv3L1), and the cyclin-dependent kinase (Mnat-1) genes. Notably, we also detected a loss of both chromatin accessibility and H3K56ac near the transcriptional start site of Wdr5, a subunit of the MLL complex that is required for ES self-renewal and maintenance of active chromatin at pluripotency genes [78]. Departure from a state of pluripotency at the beginning of differentiation involves highly complex and dynamic chromatin reorganization events to allow for proper cell lineage commitment [79][80][81][82]. Our analysis revealed loss of chromatin accessibility at TSS in several key pluripotency factors, such as Oct4 (Pou5f1), Nanog, and Sox2 after exposure to VPA, consistent with induced differentiation of ES cells. Interestingly, we also detected significantly lower chromatin accessibility at known enhancer sites [83] within the Klf4 locus (Additional file 1: Fig. S5), further supporting a mechanism by which VPA treatment induces changes in chromatin accessibility that promote a phenotype consistent with the onset of differentiation. Similarly, the striking changes in chromatin accessibility and H3K56ac enrichment at the pluripotency locus Lefty [25,64,65] provide mechanistic insight into the early differentiation events by demonstrating how acute, high-dose VPA treatment induces specific alterations at key loci involved in loss of pluripotency and initiation to differentiation programs. Transcription factor footprinting and protein expression analyses support this concept by demonstrating a striking shift form a chromatin landscape dominated by pluripotency factor binding in undifferentiated control cultures to a pattern reflective of TF occupancy in differentiating cells following VPA treatment. Intriguingly, both permissive epigenomes as well as nucleosome flexibility have recently been associated with changes in the binding capacity of pioneering pluripotency factors to DNA motifs in nucleosomal DNA during reprogramming [84] and may conceivably also be affected by H3K56-induced changes in nucleosome stability and chromatin accessibility during directed differentiation as shown for OCT4 (this study). Whether changes in H3K56ac occupancy at pluripotency loci coincide with similar changes in classical bivalent histone marks, such as H3K4me3 and H3K27me3, remains a critical question for future analyses.
In stark contrast to pluripotency-related genes, chromatin accessibility increased in response to VPA treatment at specific loci involved in differentiation programs, particularly myogenic, neuronal, and chromatin remodeling pathways. For example, isoforms for troponin T (TnT), a key regulator of actin thin filament function that is essential for contraction of striated muscle, demonstrated prominent differences in chromatin accessibility. Notably, while accessibility was altered at putative regulatory sequences upstream of cardiac muscle-specific Tnnt2, fast skeletal muscle-specific Tnnt3 gained accessibility at the TSS suggesting a transition toward (poised) activation. Importantly, in response to VPA, we also observed a significant increase in H3K56ac enrichment the TSS of Tnnt3 that overlapped with the observed ATAC-seq peak gain in the same genomic location, suggesting a coordinate transition toward a transcriptionally permissive chromatin environment. In addition, highly significant increases in chromatin accessibility at loci, such as alpha-T-catenin (Ctnn3) and polycomb group ring finger 5 (Pcgf5), further suggest that VPA exposure leads to changes in chromatin structure at specific key regulatory factors of cardiac and neuronal differentiation programs [43,44]. Intriguingly, loss of chromatin accessibility at putative downstream enhancer regions within (See figure on next page.) Fig. 7 Chromatin landscape at the Lefty locus in mESCs and during induced differentiation. A IGV browser alignment of H3K56ac CUT&RUN enrichment and ATAC-seq patterns at the Lefty locus. Coherent broad domains of H3K56ac enrichment were identified by SEACR analysis and are denoted as blue bars. Differential enrichment and correlative accessibility change at known enhancers (shaded light green) and TSS (shaded gray). B Putative model of key transitions in chromatin accessibility, nucleosome stability and transcription factor binding* at enhancers* # and TSS sites within the Lefty locus in response to HDAC inhibitor treatment (adapted from * [65] and # [64]). C H3K56ac and OCT4 (POU5F1) ChIP-qPCR data following immunoprecipitation and amplification of regulatory sequences of the Lefty1 locus (P < 0.05). IgG was used as negative control the Hopx locus suggests that HOPX expression levels may be subjected to HDACi effects in mESCs. HOPX is an unusual homeodomain-only homeobox protein involved in regulating the activity of the transcription factor GATA4 and is critical for the modulation of embryonic cardiomyocyte proliferation and differentiation [85].
Notably, Hopx-mediated regulation of GATA4 depends on HDAC2 activity and confers repression of cell cycle gene expression and cardiac progenitor cell differentiation [49]. It remains to be tested whether VPA-induced changes in chromatin accessibility, in concert with inhibition of the catalytic activity and/or proteasomal degradation of HDAC2 [86,87], lead to hyperacetylation of GATA4 and regulation of cardiomyocyte proliferation. To this extent, the susceptibility of the Hopx locus to VPAinduced changes in chromatin accessibility observed in this study may have functional implications for current differentiation protocols. For example, Hopx is critical for the regulation of late stages of cardiomyocyte maturation [88]. Thus, abnormal Hopx regulation may be associated with the high proportion of immature cardiomyocyte cell types generated through directed differentiation of iPSCs [88]. Notably, the striking potency of VPA treatment in directed differentiation is illustrated by reports demonstrating that differentiation of stem cells of mesenchymal origin into cardiac lineage-type cells improves with VPA supplementation and is accompanied by the expression of many cardiac markers, such as GATA4, troponins, as well as several cell surface markers [27]. Analysis of chromatin accessibility provides mechanistic insight into the factors involved in directed differentiation of ES cells and is also essential to identify key transcription factors and chromatin changes involved in driving cell fate (epigenetic state and cell potency). Our massive parallel genomic footprinting analysis elucidates a prominent switch from a chromatin landscape dominated by pluripotency TFs in undifferentiated control mESCs, while VPA-induced directed differentiation triggered, for instance, changes in TF occupancy known to be associated with mesoderm/endoderm lineage commitment via TGF-β signaling [52]. Notably, these findings are in agreement with the patterns of chromatin accessibility changes detected at genes associated with pluripotency (loss of accessibility) and cardiac development, neuronal differentiation, and chromatin remodeling (gain of accessibility), and shed further light onto the pathways and mechanisms responsible for induced differentiation toward the cardiomyocyte cell lineage. Our study establishes the importance for generating maps of chromatin structure dynamics during differentiation events to improve our mechanistic understanding and, ultimately, aiding in advancing lineage-specific differentiation protocols.

Conclusion
Our results provide novel information on the early regulatory networks controlling chromatin organization in response to histone deacetylation-induced directed differentiation. Focusing specifically on shifts in chromatin accessibility and markers of nucleosome stability, we demonstrate how these regulatory processes act on a locus-specific scale. Given the integral roles of nucleosome dynamics in a wide variety of cellular processes, additional efforts are needed to enhance our understanding of the critical factors involved in key cellular transitions. The specific mechanisms triggering the expression of cardiac or neuronal markers during directed lineage differentiation remain to be established. However, our results indicate that redistribution of H3K56ac and changes in chromatin accessibility and TF occupancy during VPA-induced differentiation are important components of the epigenetic reprogramming mechanism during lineage commitment. Dissecting the chromatin pathways that directly affect transcription and maintenance of lineage commitment is essential to harness the potential of induced differentiation protocols to direct and maintain specific lineage differentiation.

Image acquisition, processing, and quantitative measurements
Laser scanning confocal microscopy was conducted using a Nikon Eclipse Ti-U/D-Eclipse C1 laser scanning confocal microscope equipped with 40× objective lens following excitation of AlexaFluor488 fluorochromes with a 488 Coherent Sapphire laser. Image acquisition was conducted using EZC1 3.91 software (Nikon) with a step size of 1 μm and a Z-stack range of 15 μm. Imaging data were subsequently analyzed by maximum intensity projections using NIH Elements 4.0 software (Nikon) and fluorescence intensity quantification was performed using the multidimensional imaging capabilities of the Automated Measurements Module of NIS Elements 4.0 software (Advanced Research, Nikon Instruments). Confocal Z-stacks of control and VPA-treated samples were imaged using identical laser power and gain parameters.

HPLC and Mass-Spec
Histone H1 expression was analyzed in control and VPAtreated mESCs by reverse-phase HPLC and Mass-Spec as described previously [37]. Briefly, cells were harvested and snap-frozen before nuclei extraction using 0.5% Nonidet P-40 in RSB (10 mM NaCl, 3 mM MgCl2, 10 mM TrisHCl, pH 7.5, protease inhibitors) and in a Dounce homogenizer at 4ºC. Released nuclei were then pelleted and resuspended in RSB buffer before extraction of chromatin and histone proteins as described previously [90]. 50-100 mg of total histone preparations was injected into a C18 reverse-phase HPLC column (Vydac) on an A¨ KTA UPC10 system (GE Healthcare). The effluent from the column was monitored at 214 nm (A214), and the peaks areas were recorded and determined with AKTA UNICORN 5.11 software. Relative amounts of total H1s were determined by ratio of the total A214 of all H1 peaks to half of the A214 of H2B peak. The A214 values of the H1 and H2B peaks were adjusted to account for the differences in the number of peptide bonds in each H1 subtype and H2B. Fractions corresponding to the H1d/H1e peak from HPLC analysis were collected and subjected to mass spectrometry analysis on a Qstar XL MS/MS system (Applied Biosystems) with electrospray ionization (ESI) as the ionization method. Analyst QS software (Applied Biosystems) was used for data acquisition and analysis.

ATAC-seq
To map chromatin accessibility genome wide, DNA was probed with hyperactive Tn5 transposase and sequencing adapters were inserted into accessible regions of chromatin as described previously [91]. Sequencing reads were used to infer regions of increased accessibility. 100,000 control and VPA-treated cells per replicate were used for ATAC-seq sequencing with the Ilumina NextSeq 500. Paired-end 42 bp sequencing reads (PE42) were mapped to the genome using the BWA (Burrows-Wheeler aligner) algorithm with default settings. Only reads that passed Illumina's purity filter, aligned with no more than two mismatches, and mapped uniquely to the genome were used in subsequent analyses and duplicate reads were removed. Genomic regions with high levels of transposition/tagging events were determined using the MACS2 peak calling algorithm [92]. Differentially enriched regions were identified using DESeq2 in GUAVA [93] as described previously [94], and regions with a P-value < 0.001 and a fold change greater than four upward or downward are considered differentially enriched peaks. Peaks within 5000 bp upstream and 3000 bp downstream of the TSS were considered for the analysis. The peaks identified were associated and annotated with the nearest genes in GUAVA's peak annotation algorithm, and these genes were then used to determine over-represented Kegg pathways and GO Terms [93].
Differential digital TF footprinting was conducted using the TOBIAS framework tools [51]. Replicate ATAC-seq files were merged into joined bam files of reads prior to analysis and the TF motifs file was prepared using all Mus musculus motifs from the JASPAR 2020 database [95]. In the Tobias pipeline, the combined bam files were subjected to bias correction, footprint score calculation, and finally the estimation of differentially bound motifs. TF motifs with −log10(p-value) above the 95% percentile and/or differential binding scores smaller than the 5% percentile and larger than the 95% percentile were highlighted as differentially bound motifs. The volcano plot was generated from the BINDetect function and the two aggregated ATAC-seq signals plots for two representative motifs were generated using the PlotAggregate function. TE analysis of the ATAC-Seq data was carried out by overlapping the differential peak region coordinates from GUAVA analysis with the mm10 RepeatMasker TE coordinates file from the UCSC table browser [96]. Any differential peak overlapping with a TE region was recorded as a TE-overlapped peak, with the possibility of overlap with multiple TEs. Finally, the proportion of peaks overlapping each type of TE was calculated relative to total number of gained-opened or gained-closed peaks.

CUT&RUN
Cleavage Under Targets and Release Using Nuclease (CUT&RUN) was conducted using H3K27me3-and H3K56ac-targeted cleavage by micrococcal nuclease in situ to release specific histone--DNA complexes into the supernatant for paired-end DNA sequencing. Approximately 2.5 × 10 5 control and VPA-treated mESCs were harvested in two experimental replicates a magnetic bead (BioMag Plus, cat #86057)-based protocol as described before [55] using chip-seq validated antibodies against H3K56ac (EMDMillipore, cat# 07-677) or H3K27me3 (Cell Signaling Technologies, cat# C36B11). Resulting DNA fragments were used for the construction of sequencing libraries using SWIFT ACCEL-NGS 2S Plus library preparation kit protocol. The Bowtie2 [97] method was first used to align the paired-end FASTQ sequencing reads to the Bowtie2 reference sequences of Mus musculus genome assembly GRCm38 (mm10). Browser Extensible Data (bed) files with read pairs on the same chromosome with fragment length less than 1000 bp were retained. The bed files were used to compute the genome coverage information on each position in the genome. Bedgraph files containing the genome coverage information were loaded into the SEACR [59] analysis pipeline to identify the sites of DNA binding with high signal to noise ratio on the genome-wide scale using the 'stringent' setting and a threshold of 0.05. SEACR output files were normalized by removing any peaks with signal value lower than half the average peak value in this file. Subsequently, output files of SEACR analyses were analyzed to identify high-signal regions consistent in either both replicate experiments within a minimum distance of 300 bp from high-signal blocks. In addition, the SEACR region overlap between replicates of the same treatment group was set to a minimum of 300 bp to be considered a unique differential peak. The unique differential peaks were then mapped to the mm10 genome annotation file from Ensemble [98] to identify peaks within a region 10,000 bp up-or downstream of genebodies. To elucidate concomitant changes in chromatin accessibility and H3K56ac occupancy, gained-open and gained-close differential peaks in the ATAC-Seq analysis were compared with unique differential peaks for control and treatment detected by CUT&RUN SEACR analysis. Only overlapping peaks within TSS of annotated genes were recorded and are classified as either I) gaining both accessibility and H3K56ac enrichment, II) losing both accessibility and H3K56ac enrichment, or III) losing chromatin accessibility while simultaneously gaining H3K56ac enrichment.

Chromatin immunoprecipitation (ChIP) and ChIP-qPCR
ChIP on control and VPA-treated mESCs was performed according to manufacturer's instructions using the ChIP-IT Express Enzymatic kit (Active Motif ) on 10 μg of chromatin per IP. Briefly, cells were first crosslinked with 1% formaldehyde for 10 min at room temperature before the reaction was stopped using glycine for 5 min. The cells were rinsed with PBS, resuspended in cell scraping buffer supplemented with PMSF. Cells were dislodged using rubber policemen before centrifugation and resuspension in ice-cold lysis buffer supplemented with PMSF and protease inhibitor cocktail. A dounce homogenizer was used to lyze cells. Solubilized chromatin was obtained by enzymatic digest at 37 °C for 8 min before ChIP using unconjugated IgG control antibodies (Invitrogen) as negative control, as well as ChIP validated anti-H3K56ac (EMDMillipore, cat# 07-677) and anti-OCT4 (Active Motif, cat# 39811) antibodies. Protein G magnetic beads were used to bind immune complexes before sequential washes and elution. Crosslinks were reversed at 95 °C for 15 min followed by Proteinase K treatment, and DNA clean-up with UltraPure Phenol: Chloform/ Isoamylalcohol. ChIP-qPCR was performed using the ChIP-IT qPCR analysis kit (Active Motif ) according to the manufacturer's instructions. Relative enrichment was calculated by dividing the normalized level of ChIP DNA to that of input DNA at the corresponding locus. ChIP-qPCR results are reported as mean ± s.d. for two experimental replicates. Primers used for qPCR to quantify the ChIP enrichment at the Lefty1 locus were previously described [100].

Statistical analyses
All data presented were collected from three independent biological replicates (except for ATAC-seq and CUT&RUN analyses, which were conducted in two independent experimental replicates) and statistically analyzed using GraphPad Prism software. Comparison of all pairs was conducted using parametric and non-parametric tests (Mann-Whitney or t-test) according to the sample distribution. Data are presented as means, and variation between individual replicates is indicated as the standard deviation (SD). Differences were considered significant when P < 0.05 and are indicated by *, P < 0.05; **, P < 0.01; ***, P < 0.005; and ****, P < 0.001.

Supplementary Information
The online version contains supplementary material available at https:// doi. org/ 10. 1186/ s13072-021-00432-5.  Figure S4. CUT&RUN validation and genome-wide distribution patterns of H3K56ac. (A) To test the CUT&RUN approach in our hands, control mESCs were used for profiling H3K27me3 genomic distributions according to [55,101]. IGV browser views of 200 kb of the Hoxd and 35 kb of the Wnt5a loci are shown. H3K27me3 enrichment is in agreement with previously published reports [56][57][58][59]. Coherent broad domains of enrichment were identified by SEACR analysis and are denoted as blue bars. (B) SEACR analysis comparison of peak numbers, % peaks reproduced and width of peaks of H3K56ac CUT&RUN assays on control (Cntl) and VPA-treated (2 mM, 48 h) mESCs in two biological replicates (R1, R2). Box pots represent the median value as well as upper and lower quartiles with whiskers representing the minimum and maximum values observed. IGV browser view of 35 kb of the Thy1 locus encoding a cell surface antigen is shown and is in accordance with previously reports [60]. Figure S5. Integrative ATAC-seq/ CUT&RUN analysis. (A) List of 11 genes found to show correlative gain in chromatin accessibility and gain in H3K56ac at the TSS. (B) Venn diagram illustrating total overlap between ATAC-seq peak changes and altered H3K56ac enrichment by SEACR analysis. (C) IGV browser view of the Klf4 locus, including E2/E3 enhancer sites. ATAC-seq chromatin accessibility tracks are aligned with H3K56ac CUT&RUN enrichment and differential patterns are indicated in shaded box. Coherent broad domains of H3K56ac enrichment were identified by SEACR analysis and are denoted as blue bars. Table S1. Transcription factor footprinting. List of TFs with high scores in control samples. Table S2. Transcription factor footprinting. List of TFs with high scores in VPA-treated samples.