PLZF targets developmental enhancers for activation during osteogenic differentiation of human mesenchymal stem cells

The PLZF transcription factor is essential for osteogenic differentiation of hMSCs; however, its regulation and molecular function during this process is not fully understood. Here, we revealed that the ZBTB16 locus encoding PLZF, is repressed by Polycomb (PcG) and H3K27me3 in naive hMSCs. At the pre-osteoblast stage of differentiation, the locus lost PcG binding and H3K27me3, gained JMJD3 recruitment, and H3K27ac resulting in high expression of PLZF. Subsequently, PLZF was recruited to osteogenic enhancers, influencing H3K27 acetylation and expression of nearby genes important for osteogenic function. Furthermore, we identified a latent enhancer within the ZBTB16/PLZF locus itself that became active, gained PLZF, p300 and Mediator binding and looped to the promoter of the nicotinamide N-methyltransferase (NNMT) gene. The increased expression of NNMT correlated with a decline in SAM levels, which is dependent on PLZF and is required for osteogenic differentiation.


Introduction
Human mesenchymal stem cells (hMSCs) possess self-renewal and multi-lineage differentiation potential toward osteogenic, chondrogenic and adipogenic specification (Pittenger et al., 1999;Prockop, 1997). Therefore, hMSCs represent a promising resource for regenerative medicine. However, for treatment efficiency and safety, it is instrumental to obtain a thorough understanding of basic molecular mechanisms that control the orchestrated activation of lineage-specific genes, determining cell fate of hMSCs and factors that ensure maintenance of the terminal differentiated state(s).
The adipogenic differentiation of MSCs have been extensively studied and have revealed that it consists of a cascade of transcriptional events that are driven by a series of transcription factors (Rosen and Spiegelman, 2014). However, the knowledge regarding osteogenic differentiation of MSCs is mainly limited to the activity of the key transcription factors RUNX2 and SP7 (osterix), and markers representing the different stages of differentiation such as proteins involved in matrix formation and mineralization (Karsenty and Wagner, 2002;Neve et al., 2011). Although it is well established that Polycomb (PcG)-and Trithorax (Trx)-group protein complexes are intimately linked to cell specification through gene repression and activation, respectively (Piunti and Shilatifard, 2016;Schuettengruber et al., 2017), very limited information regarding these complexes are available in naive hMSCs and cells undergoing osteogenic lineage specification. Moreover, there is a considerable lack in understanding of the chromatin changes associated with enhancer activation that leads to lineage-specific gene transcription during differentiation of mesenchymal stem cells.
Enhancers are cis-acting DNA regulatory elements that function as transcription factor (TF) binding platforms, increasing the transcriptional output of target genes through chromatin topological changes involving enhancer-promoter looping (Calo and Wysocka, 2013;Plank and Dean, 2014). Binding of lineage-specific TFs to enhancers is key to cell specification during stem cell differentiation and organism development (Spitz and Furlong, 2012). However, the enhancer landscape that becomes active and regulates lineage-specific gene expression during osteogenic commitment has not been fully explored.
The ZBTB16 gene locus encodes a BTB/POZ domain and zinc finger containing TF known as PLZF (Li et al., 1997). Targeted deletion of Zbtb16 in mice disrupts limb and axial skeleton patterning (Barna et al., 2000;Fischer et al., 2008), and although PLZF has been shown to be involved in osteogenic differentiation (Djouad et al., 2014;Ikeda et al., 2005), the underlying mechanisms of its function is only partly understood. PLZF has been described to have a dual function in transcription 1) as a gene repressor through interaction with HDAC1, mSin3a, SMRT and NCOR (David et al., 1998;Hong et al., 1997;Melnick et al., 2002;Wong and Privalsky, 1998), and PcG proteins (Barna et al., 2002;Boukarabila et al., 2009), 2) as a gene activator due to its positive impact on transcription (Doulatov et al., 2009;Hobbs et al., 2010;Labbaye et al., 2002;Xu et al., 2009). These studies have been performed in other tissues and cell types than MSCs, such as hematopoietic and germline cells. In these cell types, PLZF is already expressed at the stem cells stage and found to be required for the maintenance of the stem cell pool. However, in MSCs, PLZF is expressed only in differentiating MSCs and not in naive cells (stem cells), and its molecular function is so far unknown in these cells.
Through the use of genome-wide ChIP-sequencing and expression analysis, we now present evidence for a novel function of PLZF at developmental enhancers directing osteogenic differentiation of hMSCs. Interestingly, we find that the ZBTB16 (PLZF) locus is bound and repressed by PcG proteins and extensively H3K27 tri-methylated (H3K27me3) in naive hMSCs. Upon osteoblast commitment of progenitor cells (pre-osteoblast stage), H3K27me3 demethylation takes place through JMJD3 (also known as KDM6B) recruitment to the ZBTB16 locus, followed by extensive H3K27 acetylation and PLZF expression. Interestingly, upon expression, PLZF binds to a subset of enhancers, which gain H3K27 acetylation and correlates with induced expression of nearby genes important for osteogenic differentiation and function. Intriguingly, an initially PcG-repressed enhancer within the ZBTB16 locus becomes active in pre-osteoblasts, gaining PLZF, p300 and Mediator binding. By use of Chromosome Conformation Capture combined with high-throughput sequencing (4C-seq) (van de Werken et al., 2012b), we show that this ZBTB16 intragenic enhancer 'EnP' physically contacts the promoter of the nicotinamide n-methyltransferase gene (NNMT) and regulate its expression during osteogenic differentiation.
Taken together, our data reveal that the activation of the ZBTB16 (PLZF) gene locus is an important event during osteogenic differentiation. PLZF regulates osteogenic differentiation of hMSCs through binding at gene enhancers and thereby it positively affects transcription of lineage-specific genes. The enhancer, EnP, localized within the ZBTB16 locus is an example of a latent developmental enhancer that becomes active upon differentiation and regulates NNMT expression in a PLZFdependent manner through Mediator and RNA-PolII recruitment and enhancer-promoter looping.

Results
The ZBTB16 locus is activated during early osteogenic differentiation To understand the transcriptional events that control the osteogenic differentiation of hMSCs, we performed RNA-seq analysis for global expression levels and ChIP-seq analysis for histone marks known to be associated with transcriptionally repressed (H3K27me3) and active (H3K4me3, H3K27ac) gene loci. We looked for the changes taking place at an early time point of differentiation (10 days) as compared to naive hMSCs (day 0). Figure 1A shows the timeline of osteogenic differentiation of hMSCs (Kulterer et al., 2007;Neve et al., 2011;Qi et al., 2003). The time point chosen for our analyses was based on prior microarray analyses (Supplementary file 1, Figure 1-figure supplement 1B), and an in situ calcium staining assay and western blots, which revealed day 10 as a timepoint where progenitor cells have committed to the osteoblast lineage, but cells are still far from the mineralization stage that takes place between day 15 and day 21 (Figure 1-figure supplement 1A-C). We will refer to the cells at this stage as immature osteoblasts (day 10) ( Figure 1A). Two days after adding the osteogenic differentiation medium, the hMSCs are still proliferating (data not shown), which is supported by the expression of cyclin E and the presence of hyperphosphorylated pRB2 (p130) as revealed by western blotting (Figure 1-figure supplement 1A). Cells at this stage can be referred to as osteogenic committed progenitor cells (pre-osteoblasts) ( Figure 1A) (Neve et al., 2011). While cells are still proliferating to some extent at day 5 (Cyclin E expression), they were completely arrested 10 days after induction of osteogenic differentiation (immature osteoblasts).
By analyzing the genome wide data revealing the most pronounced collective changes in histone marks and gene transcription, we observed a correlation between the loss of H3K27me3, a gain in H3K4me3 and increased gene expression at several genomic loci (Figure 1-figure supplement 1D-H, Supplementary file 2). The GO-terms for these groups of genes were associated to bone formation such as 'extracellular matrix', 'calcium' and 'osteogenesis' among others (Figure 1-figure supplement 1I). Examples of induced genes include 1) BMP4, BMP6, IGF2, IGFBP2; gene products known to promote osteogenic differentiation (Hamidouche et al., 2010;Lavery et al., 2008), 2) Leptin, COMP, PPL; genes encoding proteins involved in matrix formation (Ishida et al., 2013;Turner et al., 2013), 3) ZBTB16, HIF3a, examples of transcription factors involved in bone development (Ikeda et al., 2005;Zhu et al., 2014). Examples of ChIP-seq tracks and aligned RNA-seq data are shown in Figure 1C and Figure 1-figure supplement 1J. In contrast to differentiation of other cell types such as neuronal lineage, where a substantial number of gene loci showed dynamic changes in Polycomb repression (Prezioso and Orlando, 2011), H3K27me3 levels were surprisingly stable in differentiating hMSCs, and only few genes have significant changes in H3K27me3 levels at their TSS ( Figure 1B, Supplementary file 3).
The most striking loss of H3K27me3 was found at the ZBTB16 locus including intragenic exons and 3'UTRs ( Figure 1B and C). Interestingly, the loss of H3K27me3 was accompanied by strong gain of H3K27ac across the whole ZBTB16 locus and gain of H3K4me3 at the TSS, which was furthermore confirmed by ChIP-QPCR ( Figure 1C and D).

PLZF is expressed in osteoblast committed progenitor cells
The ZBTB16 gene encodes the transcription factor PLZF, which has been shown to be critical for limb and axial skeleton patterning (Barna et al., 2000;Fischer et al., 2008). Intriguingly, our gene expression and histone modification data revealed that the ZBTB16 gene locus was highly transcribed in response to osteogenic commitment of hMSCs, which is consistent with the loss of H3K27me3 in immature osteoblasts (day 10 of differentiation) (Figure 2A,B, and Figure 1-figure supplement 2C). Observing earlier timepoints, we found that the expression of ZBTB16 was strongly induced (! 4000 fold analyzed by QPCR) already at day 2 of differentiation which corresponds to pre-osteoblasts. We therefore analyzed the ZBTB16 locus for the presence of PcG proteins and H3K27me3 levels at day 2 of osteogenic differentiation. Interestingly, we observed a loss of H3K27me3 and binding of the PRC2 subunit SUZ12 already in pre-osteoblasts ( Figure 1D and E and data not shown). Furthermore, recruitment of the H3K27me3 demethylase JMJD3 was observed at the locus correlating to the loss of H3K27me3 ( Figure 1E). These data revealed that derepression and activation of the ZBTB16 locus and high expression of the PLZF (Figure 2A and     Scatter diagram of base mean (X-axis) and log2 fold differences in H3K27me3 levels (Y-axis) in immature osteoblasts relative to naive hMSC. Values were quantified from two replicates within ± 5 kbp of each TSS, and genes with a significant gain or loss were colored green or orange, respectively. ChIP-seq signal for H3K27me3 (biological replicates) in naive hMSCs and immature-osteoblasts (day 10 of osteogenic Figure 1 continued on next page event in osteogenic differentiation of hMSCs corresponding to the transition of stem cells to osteoblast committed progenitor cells (pre-osteoblasts). Several markers that are expressed during osteogenic differentiation of hMSCs have been described. To understand how the induction of PLZF expression was timed relative to these, we investigated our RNA-seq data for known markers of osteogenesis. As summarized in Figure 1-figure supplement 2A, we observed that RUNX2, often referred to as the master regulator of bone formation (Liu and Lee, 2013;Long, 2011) was expressed already in proliferating naive hMSCs and that the expression was quite stable throughout the time course. Importantly, SOX9 an essential factor for chondrocyte specification and known to antagonize RUNX2 transcriptional activity (Loebel et al., 2015;Zhou et al., 2006), was highly expressed in naive hMSCs, but was strongly reduced at day 2 of osteogenic differentiation. This down regulation of SOX9 likely reflects the transition to osteoblast committed progenitor cells. Furthermore, FOXO1, previously shown to be involved in osteogenic differentiation (Teixeira et al., 2010) was induced at day 2. Importantly, transcription factors of other lineage fates such as MYOD1 (myocytic commitment) and SOX5/8 (chondrocytic commitment) (Chimal-Monroy et al., 2003) were repressed and marked by H3K27me3 in naive hMSCs and still at day 10 of osteogenic differentiation (Figure 1-figure supplement 2A). In addition, factors known to mark the transition from the pre-osteoblast stage to immature osteoblasts such as ALPL (alkaline phosphatase) and SPP1 (osteopontin) (Aubin et al., 1995;Huang et al., 2007) were strongly induced at day 10 of osteogenic differentiation (Figure 1-figure supplement 2A,B), while SP7 (osterix) and BGLAP (osteocalcin), markers of the mineralization stage (Aubin et al., 1995;Harada and Rodan, 2003) were still repressed and marked by H3K27me3 in immature osteoblasts (Figure 1-figure supplement 2A,B,C and data not shown). In conclusion, our analyses position PLZF expression as a very early event during osteogenic differentiation and likely implicated in the transition from naive hMSCs to osteoblast committed progenitor cells (preosteoblasts).
PLZF is recruited to distal regulatory regions and affects expression of osteogenic-specific genes during differentiation of hMSCs To understand the function of PLZF in osteogenic differentiation of hMSCs, we performed ChIP-seq against PLZF in naive proliferating (day 0) and immature osteoblasts (day 10 of differentiation). We identified 2,282 PLZF peaks that were observed only in immature osteoblasts ( Figure 2C, Supplementary file 4). The analysis revealed that the majority of PLZF-binding sites localized at intergenic (45.5%) and intronic regions (34.7%) combined representing 1,830 out of 2,282 PLZFbinding sites ( Figure 2D). Very few PLZF-binding sites were observed at gene promoters (± 1 kb of TSS: 3.0%) and upstream (À10 kb from TSS: 8.5%) or downstream (+ 10 kb from TSE: 4.8%) regions ( Figure 2D). Discriminative de novo motif analysis of the DNA sequences within the PLZF peaks, compared to sequences from negative control regions, revealed a significant occurrence of a TACAGC motif (E = 1.1 Â 10 À81 ), which is highly similar to the TAC(T/A)GTA PLZF motif identified Figure 1 continued differentiation) were quantified at 8,822 TSS (±5 kbp), and analyzed by Deseq2. Green dots represent gene loci with gain of H3K27me3, while orange dots represent gene loci with loss of H3K27me3 (see details in Materials and methods and Supplementary file 3). As indicated, the largest loss of H3K27me3 was observed at the ZBTB16/PLZF locus in immature osteoblasts. (C) ChIP-seq tracks for histone marks and RNA-seq tracks for transcripts are shown for naive hMSCs (day 0) and immature osteoblasts for the ZBTB16/PLZF gene locus. (D) Upper part is a schematic presentation of the ZBTB16 locus (an upstream region; exons 1-7; coding exons 2-7) showing the position of primers used for ChIP-QPCR (marked A-F). Lower panel represents ChIP-QPCR validation of H3K27me3 loss, and H3K27ac gain. The data shown represents three biological replicates shown as average of triplicate values from QPCR ± SD. (E) Loss of Polycomb (SUZ12) binding and increased JMJD3 enrichment after 2 days of osteogenic differentiation analyzed by ChIP-QPCR. The data shown represents three independent experiments and are averages of triplicate samples from QPCR ± SD. General IgG was used as control in ChIP-QPCR. The online version of this article includes the following source data and figure supplement(s) for figure 1: Source data 1. Validation of PLZF antibody in ChIP by PLZF kncokdown. Figure supplement 1. After induction of differentiation, the hMSCs cells undergo a hyperproliferative phase within the first 2 days followed by growth arrest and differentiation.   (Ivins et al., 2003). Gene ontology analyses of closest genes to PLZF bound regions enriched for terms highly relevant for osteogenic differentiation such as 'cytoskeleton' and 'extracellular matrix' among others ( Figure 2E). The specificity of the PLZF antibody in ChIP was confirmed by PLZF knockdown in hMSCs followed by induction of differentiation and ChIP-QPCR at selected loci (

PLZF binds to active chromatin regions in immature osteoblasts
To study the relationship of PLZF binding at genomic sites and histone modifications, we performed a cluster analysis relating PLZF-binding sites to ChIP-seq data for H3K27me3, H3K4me3, H3K4me1 and H3K27ac histone marks in immature osteoblasts. It has previously been suggested that PLZF functions as a repressor of gene transcription and has a link to PcG recruitment. Therefore, we were surprised to observe a very limited co-occurrence of PLZF with PRC1 (RNF2) and PRC2 (SUZ12) as well as the H3K27me3 repressive mark catalyzed by PRC2 ( Figure 2F and Figure 2-figure supplement 2A). In contrast, there was a clear co-occurrence of PLZF binding, with H3K27ac and H3K4me1, histone marks characterizing active gene enhancers (H3K27ac, H3K4me1) (Creyghton et al., 2010) and promoters (H3K27ac). However, the number of gene promoters bound by PLZF was very low ( Figure 2D), indicating that the majority of PLZF-bound regions represent active gene enhancers.

PLZF is recruited to gene enhancers in immature osteoblasts
Activating transcription factors often recruit histone acetyl transferases (HATs) leading to H3K27 acetylation as part of the gene activation mechanism (Spitz and Furlong, 2012). To investigate if PLZF during osteogenic differentiation would lead to increased H3K27ac at its binding sites, we next analyzed the levels of H3K27ac before and after PLZF recruitment genome wide. As evident in the ratio-metric heat map ( Figure 2G), approximately 18% of the 2,282 genomic sites binding PLZF, gained H3K27ac (! 2 fold, Supplementary file 5) upon osteogenic differentiation (red-colored (n = 2,282) that overlapped with different genomic regions; promoter (À1 kb to TSS, n = 70); upstream (À10 kb to TSS excluding promoters, n = 195); downstream (+10 kb from the end of the gene, n = 110); exon (n = 75); intron (n = 792); and intergenic (> ± 10 kb from TSS or end of the gene, n = 1,040). (E) GO term analyses of genes nearby to PLZF peaks (n = 2,282) observed by DAVID (The Database for Annotation, Visualization and Integrated Discovery) (Huang da, Bailey et al., 2009). The key words enriched are displayed in the plot. The Y-axis represents the Benjamini Hochberg corrected -LogP value. Grey parts of the bars are above a 0.05 cut-off. (F) Heat-maps at PLZF peaks (n = 2,282) clustered according to local densities of the H3K27me3, H3K4me3, H3K27ac and H3K4me1 histone marks in immature-osteoblasts. ChIP-seq signals from each antibody were normalized to library sizes and visualized at PLZF peaks ± 5 kb. Only PLZF peaks specific for day 10 of osteogenic differentiation were included, and peaks with PLZF signal at day 0 were excluded (for threshold and genomic peak positions please see Supplementary file 3, and Supplementary methods). (G) Heat maps showing the changes in H3K27ac and H3K4me1 at PLZF peaks (n = 2,282), and expression of the nearest gene in naive hMSCs (day 0) and immature osteoblasts. The distribution of normalized H3K27ac or H3K4me1 ChIP seq signal in naïve hMSCs or immature osteoblasts is shown in the heat maps, centered at PLZF peaks identified in immature osteoblasts (± 5 kb). The ratiometric heatmap depicts the changes in H3K27ac and H3K4me1 levels in immature-osteoblasts compared to naive hMSCs (log 2 fc day10/naïve) and all heatmaps were sorted according to the H3K27ac ratio with the highest day10/naive H3K27ac ratio at the bottom. Gain is colored red (log 2 fold ! 1, corresponds to 403 PLZF peaks), no change is colored green (log 2 fold>-1 and < 1, n = 1,844 PLZF peaks) and loss is colored blue (log 2 fold 1, n = 35 PLZF peaks). The RNA-seq based heat-map shows the log 2 fc in expression (day10/naive) of nearest transcript mapped to the center of individual PLZF peaks. Transcripts with induced expression are colored red, while blue coloring reflects decreased expression between naïve hMSCs and immature osteoblasts. Examples of genes with increased expression correlating with increased H3K27ac and PLZF binding at nearest genomic region (enhancer) are indicated on the right side.    population), whereas 81% of PLZF-bound regions retained their H3K27ac pattern (green-colored population). A small subset of PLZF-bound regions (approximately 1%) showed reduction in H3K27ac ((! 2 fold, blue-colored population). In contrast, the H3K4me1 mark was relatively constant at the PLZF-bound regions for the majority of sites (n = 2,123; 93%) when comparing immature osteoblasts to naïve hMSCs ( Figure 2G). Interestingly, when comparing the PLZF bound enhancers that gained H3K27ac to those that retained H3K27ac levels and related it to the expression of nearest transcript, we observed a strong correlation between gain of H3K27ac and increased expression of proximal genes. Gene Ontology analyses revealed that many of these genes were related to 'bone development', 'osteoclast differentiation' (rightmost panel in Figure 2G and examples of ChIP-seq tracks in Figure 2-figure supplement 2B,C and Supplementary file 5). This suggested that a subset of PLZF-bound enhancers become active for the regulation of nearest gene that relates to an increase in H3K27ac. Furthermore, some of the PLZF-binding sites with unchanged levels of H3K27ac, showed increased expression of nearest gene upon induction of osteogenic differentiation. This might be related to other chromatin features, for example, other histone modifications and/or cooperative TF-binding effects at individual enhancers ( Figure 2G) (Spitz and Furlong, 2012).

PLZF localizes to eRNA expressing active enhancers
Active enhancers produce so-called enhancer RNAs (eRNAs), which can be detected by global 5' end RNA sequencing (CAGE). In the FANTOM5 project, this technique was used to create an atlas of 43,011 CAGE-predicted active enhancer regions across numerous cell types (432 primary cell types and 241 cell lines) and tissues (135 tissues) (Andersson et al., 2014). When comparing the PLZF-binding sites in immature osteoblasts with the CAGE-derived enhancer candidates, we observed that 603 out of 2,282 PLZF-binding sites (26%) localized within ± 1 kb of a predicted FAN-TOM5 enhancer, corresponding to 740 enhancers (Supplementary file 6). Interestingly, PLZF-bound enhancers in our hMSCs ChIP-seq data were highly enriched for enhancers that were significantly expressed in FANTOM CAGE libraries, of naive hMSCs (p < 5*10 À8 , Fisher's exact test) and osteoblasts (p < 4*10 À39 , Fisher's exact test) ( Figure 2H) (it should be noted that the time point 10 days of osteogenic differentiation, immature osteoblasts, does not exist in the FANTOM5 data). Interestingly, almost equal numbers of PLZF-bound enhancers identified in our ChIP-seq analysis in immature osteoblasts were CAGE positive (FANTOM5 data) in naive hMSCs or in osteoblasts (13.2%: 98 out of 740% and 13.6%: 101 out of 740, respectively) and approximately 42% (42 out of 98 or 101) of these enhancers are in common between naive hMSCs and osteoblasts. This suggests that 1) in naive hMSCs, where PLZF is absent, other TFs define a subpopulation of active enhancers that later on during osteogenic differentiation is bound by PLZF and 2) that PLZF likely continues to be bound to a significant number of such enhancers in terminal differentiated osteoblasts. Moreover, the data suggest that there exist a number of lineage-specific enhancers that become active and gain PLZF binding in the process of differentiation (

PLZF affects histone H3K27ac at bound enhancers and an osteogenic gene expression signature
To further investigate the relationship between PLZF recruitment to enhancers and changes in H3K27ac genome wide, we performed siRNA-mediated PLZF knockdown in naïve hMSCs followed by induction of osteogenic differentiation, and H3K27ac ChIP-seq. We tested three different siRNA oligos targeting PLZF (Figure 3-figure supplement 1). All three oligos showed significant knockdown of PLZF mRNA compared to non-targeting control siRNA. We chose oligo #57 for further experiments. Since it is not possible to maintain a stable knockdown of PLZF for 10 days of differentiation using siRNA, we analyzed the changes in H3K27ac at day 2 of differentiation (pre-osteoblasts), where PLZF knockdown was still pronounced ( Figure 3A). Clustering of the data revealed a highly reproducible increase in H3K27ac in pre-osteoblasts, which was reduced when PLZF expression was down-regulated by siRNA ( Figure 3B).
Examples of regions showing a gain of H3K27ac in a PLZF-dependent manner are shown as ChIPseq tracks in Figure 3C. In addition to regions that gained H3K27ac, there was furthermore a fraction of regions that lost H3K27ac upon osteogenic differentiation (blue in left panel). However, only To further gain insight to the transcriptional changes taking place in the absence of PLZF during osteogenic differentiation, we performed a whole transcriptome analysis by RNA-seq after siRNAmediated PLZF knockdown in hMSCs. There were 1,897 genes differentially regulated upon osteogenic induction in control siRNA-transfected hMSCs compared to 1,184 genes differentially regulated in PLZF knockdown cells. Only 44% (950/1,897) of these genes overlapped between control siRNA and PLZF knockdown cells (Figure 3-figure supplement 1B). As shown in the Figure 3D and E, the PLZF knockdown significantly reduced the expression of genes that were normally induced in osteoblast committed progenitor cells ( Figure 3D and E, Supplementary file 8). The heatmaps in Figure 3D, representing the RPKM-values of differentially regulated genes from RNAseq of PLZF knockdown or control siRNA (in biological triplicates) in hMSCs. The vertical order of genes is similar to the clustered heatmap shown in Figure 1-figure supplement 1B, (differentially regulated genes from microarray analyses). The RNA-seq data clearly indicate; i) a high degree of reproducibility in biological replicates, ii) pronounced impact of PLZF knockdown on gene regulation at early stage (2 days) of osteogenic differentiation. The expression of upregulated genes is mainly affected by the PLZF knockdown. Furthermore, in Figure 3E, the plot shows that PLZF knockdown significantly affected the osteogenic transcriptional program (p < 0.0001) when comparing the log2 fold differences between PLZF knockdown or negative control siRNA in hMSCs, harvested at day 0 and day 2. The Gene Ontology (GO) analyses of differentially regulated genes in control siRNA transfected cells enriched the terms such as 'extracellular matrix organization, regulation of cell differentiation and ossification'. Whereas genes regulated in PLZF knockdown hMSCs, enriched terms including 'Interferon signaling and chondrocyte differentiation' (Figure 3-figure supplement 1C). Moreover, we looked for the influence of PLZF knockdown on expression of genes proximal to the PLZF-bound enhancers that gain H3K27ac. Interestingly, we observed that lack of PLZF in many cases correlated with the reduced expression of nearby genes to the regions that gain H3K27ac in PLZF-dependent manner (examples are shown in Figure 3C and F).
The PLZF-bound developmental enhancer EnP becomes active upon osteogenic differentiation of hMSC Interestingly, among the PLZF-bound genomic regions, we identified the ZBTB16 locus itself ( Figure 4A). Multiple PLZF peaks appeared across the intronic region covering around 50 kb, and flanking exons 3 and 4. We also observed a remarkable gain of H3K27ac at these regions overlapping with PLZF peaks in immature osteoblasts. In contrast, H3K27ac was absent in naive hMSCs and   the region was heavily marked by H3K27me3. Intriguingly, two out of three observed peaks overlapped with CAGE defined enhancers. We refer to the enhancer with the most significant PLZF peak as 'EnP' (Enhancer within PLZF locus bound by PLZF) and decided to further characterize it as an example of a PLZF-bound osteogenic enhancer. To validate, that the PLZF-bound intragenic EnP element correspond to a functional enhancer, we examined the presence of histone marks and gene regulatory factors known to bind enhancers by ChIP QPCR. First, we analyzed a number of relevant histone modifications and confirmed that the EnP element gained H3K4me1 and H3K27ac in preosteoblasts as compared to naive hMSCs, but contained low levels of H3K4me3 and H3K36me3 compared to other regions within the ZBTB16 locus ( Figure 4B). The gain in H3K4me3 and H3K27ac was observed at the TSS, and H3K36me3 was more enriched at exon two within the ZBTB16 locus in pre-osteoblasts correlating with transcription of the gene ( Figure 4B). We obtained similar results when investigating another PLZF peak within the ZBTB16 locus (Figure 4-figure supplement 1A). It is well established that binding of the CBP/p300 histone acetyl transferases (HATs) marks enhancers and catalyzes H3K27 acetylation (Visel et al., 2009). Moreover, the binding of lineagespecific TFs together with the Mediator co-activator complex at enhancers facilitates RNA polymerase II (RNA PolII) recruitment to the promoter of target genes (Carlsten et al., 2013). We therefore, performed ChIP for p300 and the Mediator co-activator complex (MED1, MED12) in naive hMSCs and in pre-osteoblasts. The data revealed enrichment of all three factors at EnP ( Figure 4C). Taken together, our data strongly suggest that the EnP element that lies within the ZBTB16 locus represents a latent (no H3K4me1 and no H327ac) PcG repressed enhancer in naive hMSCs which gains the characteristics of an active enhancer in pre-osteoblasts.
To further validate EnP as a differentiation-induced enhancer, we analyzed the enhancer activity of EnP in a GFP reporter assay ( Figure 4D). As a control element of similar size, we cloned another upstream distal region from the ZBTB16 locus (CtE) that did not show PLZF binding. There was a complete absence of GFP fluorescence in hMSCs infected with the CtE-GFP element or empty vector control ( Figure 4E and F). Interestingly, EnP showed activity in a differentiation-specific manner, being completely inactive in naive hMSCs ( Figure 4E and F). These results demonstrated that the Figure 4 continued three biological replicates, ****p < 0.0001, **p = 0.0025 for H3K27ac, not significant for other regions; **p = 0.0092 for H3K4me1, not significant for other regions; **p = 0.0071 for H3K4me3 at TSS, not significant for other regions; ****p < 0.0001, *p = 0.027 for H3K36me3, not significant for other regions; calculated by two-way ANOVA with Sidak's multiple comparison tests. The regions (promoter, TSS, EnP) amplified using primers A, B and E shown in Figure 1E and primers binding at Exon 2 of the ZBTB16 locus. (C) ChIP-QPCR for the Mediator components MED1 and MED12 as well as p300 in naive or in pre-osteoblasts (2 days of osteogenic differentiation) at the EnP element. Data shows mean ± SD from triplicates of QPCR from three biological replicates, **p < 0.001 calculated by multiple t-tests with FDR 1% and two-stage step-up method of Benjamini, Krieger and Yekutieli. (D) Schematic presentation of the lentiviral enhancer GFP-reporter system used (pSINMIN). EnP (region corresponding to the largest PLZF peak in Figure 4A) or a control region (CtE; region within ZBTB16 locus distal to EnP but without a PLZF peak) was cloned upstream of the minimal TK promoter driving the expression of GFP. (E) Histogram represents GFP expression measured by flow cytometer in hMSCs cells transduced with lentivirus encoding EnP-GFP followed by induction of osteogenic differentiation (2 days). Data is a representative of five biological replicates. The percentage of GFP positive cells from gated live cells are indicated in the upper right corner. (F) Box plot represents the GFP expression obtained by flow cytometer measurements in hMSCs transduced with empty vector (EV), control element (CtE) or enhancer (EnP) cloned in the pSINMIN GFP reporter, before and after induction of osteogenic differentiation for two days. The data shows median calculated from five biological replicates. Median is shown by horizontal line. ****p < 0.0001 calculated by 2-way ANOVA with Sidak's multiple comparison tests. (G) Integration of the GFP coding sequence was analyzed by QPCR using primers for GFP on genomic DNA isolated from each group of samples. Data shown are mean ± SD from three biological replicates. Untransduced hMSCs were used as a negative control and shown as right most bars in the plot. The online version of this article includes the following source data and figure supplement(s) for figure 4: Source data 1. ChIP for Med1, Med12 and P300 in hMSCs. Source data 2. H3K4me3 ChIP in hMSCs. Source data 3. H3K4me1 ChIP in hMScs. Source data 4. H3K27ac ChIP in hMSCs. Source data 5. H3K36me3 ChIP in hMSCs. Source data 6. Med12 ChIP in hMSCs. Source data 7. Med1 ChIP in hMSCs. Source data 8. GFP reproter integration analyses. EnP element likely represents a latent repressed enhancer in naive hMSCs, which upon induction of osteogenic differentiation acquires the properties of an active enhancer and can drive the expression of a nearby gene. To ensure that the integration of the lentiviral reporter constructs was comparable for EnP-GFP and CtE-GFP, we performed a QPCR for the GFP coding part on genomic DNA (gDNA) isolated from transduced cells. The data confirmed that the integration efficiency of the two reporter constructs were highly comparable and therefore the increase in GFP fluorescence in the EnP-GFP infected hMSCs was due to the functional activity of the EnP enhancer in response to induction of osteogenic differentiation ( Figure 4G).
The PLZF-bound intragenic enhancer (EnP) loops to the promoter of the neighboring nicotinamide N-methyl transferase gene (NNMT) gene It is widely accepted that enhancers interact with gene promoters through chromatin looping and thereby regulate the transcriptional activity of associated genes (de Laat and Duboule, 2013; Levine et al., 2014). Therefore, we decided to use 4C-seq in order to obtain an unbiased genomewide screen for DNA contacts made by the EnP element upon induction of osteogenic differentiation. The experiment was performed using naïve and pre-osteoblasts in biological replicates. Interestingly, the 4C-seq data revealed several contacts between the EnP element and other genomic sites observed within a 1 Mb window. We found that the EnP element looped to the promoter of the NNMT gene, which is located~100 kb downstream of the EnP element ( Figure 5A). Furthermore, this contact was specific for pre-osteoblasts, since no looping was detected in naïve hMSCs ( Figure 5A). This is in agreement with the ZBTB16 locus being repressed by PcG and H3K27me3 in naive hMSCs and the inaccessibility of the EnP enhancer. To test the reliability of our 4C-seq data, we repeated the experiment under similar conditions, but using the NNMT promoter as viewpoint for the analysis. The results confirmed the contact of the NNMT promoter with the EnP element within the ZBTB16 locus and that it was specific for hMSCs undergoing osteogenic differentiation ( Figure 5B). Interestingly, a rather extended contact area around EnP (> 10 kb) was observed when the NNMT promoter was used as viewpoint. Additionally, a NNMT upstream region and the promoter of ZBTB16 were used as a negative 4C-seq controls. The NNMT upstream region showed limited, but much less than NNMT-promoter enriched contact frequencies with the EnP element, whereas the ZBTB16 promoter showed no enrichment of contact frequencies interaction with EnP. Thus this may represent a Polycomb repressed compacted region, since the whole locus was covered by H3K27me3 ( Figure 5-figure supplement 1A and B). In conclusion, these data support the existence of what we believe could be a cluster of enhancers localized within the ZBTB16 gene locus, where we observed a strong and broad H3K27ac enrichment in ChIP-seq and multiple PLZF peaks (intron 2-4 in Figure 4A). This suggests, that beside EnP several other intragenic elements within the ZBTB16 locus have the potential to act as gene enhancers.

The ZBTB16 intragenic enhancer EnP regulates the NNMT gene promoter
In order to reveal whether the contact of the EnP enhancer with the NNMT promoter had an impact on NNMT expression, we measured mRNA and protein levels. Indeed, NNMT expression was increased upon induction of osteogenic differentiation of hMSCs ( Figure 6A and B). Combined with the 4C-seq analyses these data suggested a positive influence of the distal EnP enhancer located within the ZBTB16 locus on NNMT expression.
The NNMT gene encodes the nicotinamide N-methyl transferase, a metabolic enzyme that regulates SAM (S-adenosyl-methionine) homeostasis. NNMT functions in a biochemical reaction where it transfers the methyl group from SAM to nicotinamide and converts SAM to SAH (S-adenosyl-homocysteine) ( Figure 6C) (Aksoy et al., 1995). To investigate whether the increase in NNMT expression upon osteogenic differentiation affected the overall SAM levels in the cells, we measured SAM using a fluorescence-based assay. Interestingly, we observed a decrease in SAM levels upon differentiation, suggesting a correlation between increase in NNMT expression and a decrease in SAM levels (t test, p = 0.01, Figure 6D). However, the drop in SAM levels was only minor, although significant (p = 0.01 by t-test), suggesting a delicate balance on SAM homeostasis during osteogenic differentiation. SAM being a critical co-factor for methylation of DNA through the action of DNA methyl transferases (DNMTs) and histones through the activity of histone methyl transferases (HMTs)  obviously needs to be tightly regulated. However, our data suggest that fine-tuning of global methylation during osteogenic differentiation could partly result from small changes in SAM levels via NNMT expression regulated through PLZF and the intragenic enhancer EnP. Overall, these results provide evidence for a dynamic activation and function of a lineage specific developmental enhancer EnP present within the ZBTB16 locus, which regulates the expression of NNMT and thereby SAM homeostasis upon osteogenic differentiation of hMSCs.
To further characterize the effect of NNMT on osteogenic differentiation of hMSCs, we used siRNA to knock-down NNMT expression. We tested three different siRNA oligos to target NNMT and found all three oligos showed significant knockdown efficiency (Figure 6-figure supplement 1). As shown in Figure 6 reduction in NNMT levels (panel E) inhibited the increase in ZBTB16 and ALPL expression (panels F and G, respectively) normally observed during osteogenic differentiation suggesting a block in differentiation.

PLZF is required for the recruitment of Mediator and RNA PolII at the EnP enhancer
Next, we wanted to explore the function of PLZF at the EnP enhancer and its involvement in NNMT expression. Therefore, we first investigated if PLZF knockdown had any influence on NNMT expression. Interestingly, the knockdown of PLZF negatively influenced NNMT gene induction ( Figure 7A, B). This was also reflected in SAM levels, since the reduction in SAM levels previously observed during osteogenic differentiation ( Figure 6D) was abrogated in the absence of PLZF ( Figure 7C). We next tested if the looping between EnP and the NNMT promoter was affected in the absence of PLZF. We performed 4C-seq experiments in the PLZF knockdown hMSCs using the EnP enhancer element as the viewpoint in naïve and in pre-osteoblasts. To our surprise, the EnP element looped to the NNMT promoter with the same frequency irrespective of PLZF knockdown (Figure 7-figure  supplement 1A and B). This suggested that PLZF does not affect the chromatin looping between EnP and the NNMT promoter.
Earlier, we have shown an enrichment of Mediator at EnP upon osteogenic differentiation ( Figure 4C). Therefore, we wondered if PLZF played a role in facilitating the binding of Mediator and RNAPolII at EnP. To address this question, we next performed ChIP for MED12 and RNAPolII in the presence or absence of PLZF and induction of osteogenic differentiation. Interestingly, the recruitment of MED12 and RNAPolII at the EnP enhancer was dependent on PLZF expression upon induction of osteogenic differentiation ( Figure 7D).
Taken together, these results suggest that PLZF influences the expression of the NNMT gene via recruitment of Mediator and RNAPolII at the EnP enhancer, without regulating enhancer-promoter looping per se.

Discussion
In this study, we have shown that the induction of the ZBTB16 gene locus encoding the transcription factor PLZF is an important event during osteogenic differentiation of hMSCs. This is in agreement with the genetic knockout of the Zbtb16 in mice, showing severe defects in limb and axial skeleton Figure 5 continued enhancer in naive hMSCs and in pre-osteoblasts (2 days of osteogenic differentiation). The plot represents the overlay between naïve hMSCs and pre-osteoblasts. Data are displayed as reads per million (RPM). The arrow indicates the view point for depicted track and asterisks indicate contact points revealed by 4C-Seq. The plots reveal that EnP physically contacts the NNMT promoter and vice versa (shown in B). Gray rectangles indicate areas with higher contact frequency in pre-osteoblasts over naive hMSCs (X 2 -test; False Discovery Rate < 0.01). (B) The NNMT promoter as a view point in 4C-Seq revealed a high frequency of contact at a 30 kb region in intron 2 and 3 of the ZBTB16 locus after 2 days of osteogenic differentiation. The data presented here is representative of two biological replicates. Gray rectangles indicate the EnP-enhancer area with higher contact frequency in preosteoblasts over naive hMSCs (X 2 -test; False Discovery Rate < 0.01). The online version of this article includes the following figure supplement(s) for figure 5:  NNMT expression during osteogenic differentiation as analyzed by RT-QPCR and western blot. RT-QPCR data represents mean ± SD from three biological replicates. *p = 0.0223, **p = 0.0011, ****p < 0.0001, significance was calculated by two-way ANOVA using Bonferroni's multiple comparisons test. (C) Schematic presentation of the biochemical reaction whereby NNMT transfers a methyl group from SAM (S-adenosyl methionine) to nicotinamide, producing Figure 6 continued on next page patterning (Barna et al., 2000). Unexpectedly, however, we found that PLZF, rather than functioning as a gene repressor by binding at promoters as described in other cellular models (Liu et al., 2016), mainly localized to gene enhancers in immature osteoblasts. Increased H3K27ac at PLZF-binding sites often correlated with higher expression of nearby transcripts in a PLZF-dependent manner, many of which represent protein encoding genes known to be important in osteogenic differentiation, and thereby supporting a more general role of PLZF at osteogenic lineage-specific enhancers (Supplementary file 4).
We also revealed the presence of an intragenic enhancer (EnP) within this locus that is repressed by PcG and H3K27me3 in naive hMSCs but becomes highly active at an early stage of osteogenic differentiation of hMSCs. Our data also demonstrate the differentiation-specific looping of EnP to its target promoter of the NNMT gene. This is now one of the few known examples where dynamics of enhancer promoter interaction has been described during mesenchymal stem cell differentiation (Dixon et al., 2015).

PLZF marks the differentiation onset of naive hMSCs
In line with previous studies (Djouad et al., 2014), we observed that PLZF expression was not specific for the osteogenic lineage since induction of chondrogenic-and adipogenic differentiation of hMSCs also led to induction of ZBTB16 expression (data not shown). These observations suggest that PLZF functions at an early state of hMSCs differentiation and promotes the transition from the naive stem cell stage to a more committed state for all three lineages. Therefore, to obtain specific gene expression patterns that characterize the three different lineages and their functions, other cell-type-specific TFs must operate together with PLZF at enhancers in a context dependent manner. Recently, the genome-wide binding profile for Runx2 was characterized in mouse pre-osteoblastic cell lines and suggested a possible link of this TF to enhancer function (Meyer et al., 2014;Wu et al., 2014). Future studies using human MSCs should be designed to reveal if RUNX2 and FOXO1 bind enhancers and cooperate with PLZF or work independently in gene regulation during osteogenic commitment. In our de novo motif analyses of PLZF-binding sites, we identified the TACAGT motif, that is identical to previously published data (Ivins et al., 2003). The motif is, furthermore, a recognition site for the transcription factor OSR2, implicated in proliferation of osteoblasts and in osteogenic differentiation of dental follicle cells (Kawai et al., 2007;Park et al., 2015). Whether OSR2 and PLZF compete or cooperate for binding to enhancers and if OSR2 functions together with PLZF to regulate osteogenic gene expression is not known and would require further studies. We observed that MED12 binding at the osteogenic enhancer EnP was dependent on PLZF, although no direct interaction has been observed between PLZF and MED12 (data not shown), suggesting involvement of other cofactor(s). It should be considered that, the influence of PLZF on H3K27ac at distal regulatory elements might facilitate MED12 binding at these regions which is reverted in the absence of PLZF. 1MNA (one methyl nicotinamide) and SAH (S-adenosyl homocysteine). (D) SAM levels measured by a fluorescence based assay (Mediomics) in hMSCs before and after osteogenic differentiation for 2 days. The data shown are averages of three biological replicates (p value = 0.01, two tailed t-test). (E) NNMT knockdown using siRNA in naive and osteogenic differentiated hMSCs, assessed by RT-QPCR. Two different siRNAs were used for knockdown of NNMT expression. Experiments shown are mean ± SD from three biological replicates. (F) Expression of ZBTB16 after NNMT knockdown analyzed by RT-QPCR. Experiments shown are mean ± SD from two biological replicates. (G) Expression of ALPL after NNMT knockdown analyzed by RT-QPCR. Experiments shown are mean ± SD from two biological replicates. ****p < 0.0001, statistical significance was calculated by two-way ANOVA using Tukey's multiple comparisons test. The online version of this article includes the following source data and figure supplement(s) for figure 6: Source data 1. NNMT expression by RT-QPCR. Source data 2. ALPL expression by RT-QPCR. Source data 3. PLZF expression by RT-QPCR.  A latent developmental enhancer EnP within the ZBTB16/PLZF locus EnP that lies within the ZBTB16 locus represents a latent PcG-repressed developmental enhancer in naive hMSCs. Upon induction of osteogenic differentiation, H3K27me3 was removed from EnP possibly through JMJD3 recruitment and there was a subsequent gain of chromatin marks characterizing active enhancers (H3K4me1 and H3K27ac) and gene regulatory factors such as Mediator, p300 and PLZF. Intriguingly, a reporter assay in which EnP was integrated in the genome of hMSCs cells showed that EnP functions only upon induction of differentiation. This suggests that the EnP enhancer element in the context of a GFP reporter also require lineage-specific factor(s) such as PLZF for its activity.
The region of the ZBTB16 locus that we found to be in contact with the promoter of the NNMT gene covered an area around 50 kb and included introns 3-5 which showed a strong enrichment for H3K27ac and H3K4me1 upon osteogenic differentiation and gained PLZF binding at several positions. These characteristics suggest that this ZBTB16 intragenic region could represent a cluster of several active enhancers (Whyte et al., 2013) that likely binds other TFs besides PLZF. Since increased expression of PLZF was also observed during induction of chondrogenic and adipogenic differentiation (data not shown), we assume that EnP might be functional in these lineages as well. However, whether all potential intragenic enhancers in the ZBTB16 locus would be regulated similarly between the three different cell lineages and contact similar or different target promoters would need further investigation.

Regulation of SAM levels through NNMT expression and activity
We found that EnP loops to the NNMT promoter and correlated with induced expression of the NNMT gene during osteogenic differentiation, and that NNMT expression is required for the differentiation of hMSCs. The enzyme NNMT regulates nicotinamide (vitamin B3) levels through methylation, using SAM as a methyl donor ( Figure 6C). NNMT is overexpressed in many cancers (Sartini et al., 2013;Ulanovskaya et al., 2013) and can result in an altered epigenetic state. This has been proposed to happen due to draining of methyl groups from the methionine cycle leaving very stable 1MNA and SAH byproducts, thereby changing the methylome of cancer cells (Shlomi and Rabinowitz, 2013;Ulanovskaya et al., 2013). Based on our data, it is tempting to speculate that the ZBTB16 intragenic enhancer EnP and PLZF, helps to fine tune global methylation patterns through transcriptional regulation of the NNMT gene during osteogenic commitment of progenitor cells. These important aspects require further detailed investigations.
In conclusion, our results identify the transcription factor PLZF as a novel component that marks the enhancer landscape in hMSCs during osteogenic differentiation, acting as a positive regulator of enhancer function. Importantly, we find that the derepression of the ZBTB16 locus, besides causing increased PLZF expression, furthermore, exposes an intragenic developmental enhancer EnP that Figure 7 continued ****p < 0.0001). The impact of PLZF knockdown on ALPL expression is shown in the lower part of the panel B. The data shown are mean ± SD from three biological replicates, ****p < 0.0001, significance was calculated by two-way ANOVA using Tukey's multiple comparisons test. (C) Fluorescence based SAM assay revealed that PLZF knockdown prevented the decline in SAM levels as observed in control siRNA cells (tested as significant p = 0.01 two tailed t-test) during early osteogenic differentiation (day2). The data shown here is an average of three biological replicates. (D) ChIP followed by QPCR for PLZF, RNA PolII, and MED12 at the EnP element before and after PLZF knockdown in hMSCs. Osteogenic differentiation was induced one day after siRNA transfection and cells were fixed for ChIP two days later. The data represents two biological replicates and show averages of triplicate values from QPCR ± SD. ****p < 0.0001, significance was calculated by two-way ANOVA using Tukey's multiple comparisons test. (E) Model to show that the ZBTB16 locus is repressed by Polycomb protein complexes (PcG) and marked by H3K27me3 in naive hMSCs. Upon induction of osteogenic differentiation, the ZBTB16 locus gets derepressed by losing PcG binding, and H3K27me3 through JMJD3 recruitment, gain H3K27ac which eventually results in high expression of PLZF. The intragenic enhancer 'EnP' gets exposed, gains H3K27ac and H3K4me1 histone marks and binds PLZF as well as P300 and the Mediator complex. Subsequently, the EnP element loops to the promoter of the NNMT gene (100 kb downstream) and induce its expression and as a consequence regulates SAM homeostasis during osteogenic differentiation of hMSCs. The online version of this article includes the following source data and figure supplement(s) for figure 7: Source data 1. ChIP for enhancer binding proteins in hMSCs, in the absence of PLZF. Figure supplement 1. 4C-sequencing map revealing that the PLZF knockdown (using siRNA) did not affect the contact frequency between the EnP element and NNMT promoter observed during early osteogenic differentiation (day 2).

ChIP sequencing
ChIP DNA from three parallel ChIPs were pooled or individual ChIPs from three biological replicates were used and 10 ng each was used for making ChIP-seq libraries. The libraries were prepared using 'ChIP seq DNA sample preparation kit' from Illumina following manufacturer's instructions. Individual samples were runsequenced in a single lanes on HiSeq2000 (Illumina). H3K27ac samples were or multiplexed (using NEB kit) and run on a Illumina Genome Analyzer IIx, HiSeq2000 (Illumina), or NextSeq500 as single end sequencing.

4C-Sequencing
Templates for 4C were prepared as described previously (van de Werken et al., 2012a;van de Werken et al., 2012b). The enzymes used were four base pair cutters DpnII and Csp6I. Osteogenic induced (2 days; preosteoblasts) or naive hMSCs (10 Â 10 7 cells per group) were harvested and fixed in 2% (v/v) formaldehyde in PBS, 10% (v/v) FBS, for 10 min with rotation at room temperature. The procedure was continued as described earlier (van de Werken et al., 2012a).
The primer sequences are given in the Table with all primers. The final PCR amplification was performed using 3 mg of 4C DNA, and PCR conditions were 98˚C (2 min), 98˚C (15 s), 58˚C (1 min), 72˚C (3 min), repeated for 30 cycles and final amplification for 7 min at 72˚C using Phusion High Fidelity DNA Polymerase (Thermo Scientific).
S-adenosyl methionine (SAM) assay SAM was measured using the Bridge-It S-Adenosyl Methionine (SAM) Fluorescence Assay Kit from Mediomics according to the manufacturer's instructions (Mediomics, LLC, St. Louis, Missouri, U.S. A.). Osteogenic differentiation of hMSCs was induced and cells were harvested after 2 days together with naïve cells and processed for SAM assay. The procedures was followed as described in the assay kit manual.

RT-PCR, microarray and RNA sequencing
Total RNA representing the indicated time points were purified from hMSCs using RNeasy Plus Mini kit (QIAGEN 74106), and cDNA was generated by RT-PCR using the TaqMan Reverse Transcription Reagents (Applied Biosystems#N808-0234). Quantifications were done using the Fast SYBR Green Master Mix (Applied Biosystems#4385612) and an ABI Step One Plus instrument. GAPDH was used for normalization. The sequences of the primers used are listed in the Supplementary file 9 with primers. The RNA samples at days 0, 2, 5, and 10 of osteogenic differentiation (in three biological replicates) were sent for microarray hybridization using affymetrix gene chips (HT_HG-U133_Plus_PM).
For RNA sequencing, samples were processed using a True Seq RNA Sample preparation kit (Illumina # 15025062) following the manufacturer's instructions. Samples were multiplexed and sequenced in HiSeq2000.

Enhancer activity reporter assay
The pSINMIN lentiviral 'enhancer reporter vector' was kindly provided by Johanna Wysoca's lab (Rada-Iglesias et al., 2011). The enhancer fragment (EnP) (corresponding to the biggest PLZF peak within the ZBTB16 locus) or control element region without a PLZF peak within the ZBTB16 locus (CtE) were PCR amplified from genomic DNA (gDNA) isolated from hMSCs using primers described in the primer Table and cloned into the PCR8 TOPO TA vector (Invitrogen). The relevant fragments were cut out from the TOPO vector using EcoRI sites and subsequently cloned into the pSINMIN vector using EcoRI sites. Positive clones were confirmed by sequencing.
Lenti-virus were produced by transfecting 20 mg reporter plasmid together with 15 mg PAX8 and 6 mg of VSV in 293FT cells seeded at a density of 5 Â 10 6 cells per 15 cm dish (Nunc cell culture) using the CaPO 4 method. 48 hr following transfection and media change the supernatant containing viral particles were collected and filtered through 0.45 mm filters. Polybrene was added to the supernatant 8 mg/ml and hMSCs were infected with virus overnight. Cells were trypsinised and reseeded 48 hr later and Puromycin (0.5 mg/ml) was added for selection. Two days later, osteogenic differentiation was induced and cells were harvested after 2 days for flowcytometric analyses or gDNA isolation. Genomic DNA was isolated using the DNeasy Blood and tissue kit (QIAGEN: 69504). Quantitative PCR was performed using GFP primers described in the primer Table. QPCR on the HOXD locus was used for normalization. Untransduced hMSCs were used as a control.

Alizarin red staining
Naive or osteogenic induced hMSCs were washed two times with PBS followed by fixation in 96% ethanol for 15 min at room temperature. Alizarin red solution 1% (w/v, dissolved in milliQ water) was added to the fixed cells and incubated for 60 min at room temperature with gentle rotation. Finally, cells were carefully (not to lose precipitates) washed three times with water and dried. Pictures were taken with a Leica camera (DFC295) fitted to a microscope at 20X magnification.

Bioinformatic analyses ChIP sequencing
Reads were trimmed for low-quality nucleotides and mapped to the hg19 human genome sequence using Bowtie v0.12.7 (Langmead et al., 2009) using the parameters '-S -m 1'. When nothing else is mentioned, subsequent processing, signal quantitation, peak-finding, distance calculation, normalization, clustering and visualization was done using EaSeq and default settings (http://easeq.net Lerdrup et al., 2016). Values in 2D-histograms in tracks or heatmaps were counted and normalized to fragments per million reads per kbp. The number of fragments was derived from the count by dividing it with (1 + DNA fragment size / bin size). Aligned datasets were filtered for multiple reads at the same position likely to occur from PCR amplification. The reads were extended to a DNA fragment size of 300 bp. Genes subject to differential H3K27me3 levels were identified as follows: Read counts from two biological H3K27me3 ChIP-seq replicates from day 0 (naïve hMSCs) and day 10 (immature osteoblast) samples were quantified within 5five kbp windows surrounding 30,716 nonredundant TSS in Refseq Hg19 gene annotation data (O'Leary et al., 2016), which were downloaded from the refflat table at UCSC (Kent et al., 2002). TSSes were filtered for low abundance by requiring a read count >=15 from each sample and >=100 for all four samples collectively. The 8,222 TSSes meeting these criteria were analysed for fold change and significance using DESeq2 (Love et al., 2014 ). For 2D-histograms of ChIP-seq data, Tthe number of reads overlapping with an area of +/-1 kbp of transcript start to end for H3K27me3 and ±1 kbp of TSS for H3K27ac and H3K4me3 were quantified. Each data set was quantile normalized in order to compensate for variations in ChIP efficiency between the samples. Regions with low H3K27me3 levels (<0.5) or low H3K4me3 (<5) in both samples were filtered and two fold changes in each histone mark was calculated (day10/day0 (naïve)). Regions with gain or loss in each histone mark were calculated by selecting the population with twofold change or higher in either direction. Values in 2D-histograms in tracks or heatmaps were counted and normalized to fragments per million reads per kbp. The number of fragments was derived from the count by dividing it with (1 + DNA fragment size/bin size).

Peak-finding
Specific PLZF peaks were identified using the IgG library data set as negative control. The procedure resembles that of MACS with some modifications (Zhang et al., 2008). Each dataset was divided into 100 bp windows, and the reads within each window scanned genome-wide. A normalization coefficient (NCIS) serving to normalize the background levels of the two datasets was analyzed in accordance with Liang and Keles, (Liang and Keleş, 2012). Window size was manually set to 100 bp. Global thresholds were calculated based on a Poisson Distribution using the genome-wide average number of reads in the windows and a p-value of 1E-05. Adaptive thresholds were modeled the following way: The average number of negative control reads in areas corresponding to 10x, 50x and 250x window size (100 bp) was calculated for each position. This number was used as lambdas for Poisson distributions and thresholds that matched a p-value of 1E-05 were calculated. The most conservative threshold was chosen from the three local thresholds of the control and the global threshold of the sample. Thresholds from the negative control were scaled according with the NCIS normalization factor. The position and statistics of windows passing the most conservative threshold, and having a NCIS-normalized log2-fold Sample/Control-ratio above 2, as well as less than a 3:1 difference between the signal on the plus and minus strands were extracted into a separate list, and the entire procedure was done four times where the windows were shifted 25 bp each time. Windows within 100 bp of each other and overlapping windows were merged. For each region in the resulting the borders were refined by sliding a window of 100 bp from one window-size upstream to downstream of the temporary border. The exact position where the number of samples reads within the window fell below the threshold was defined as the new border of that region. Shoulders were excluded at values below m +2 SD. After border refinement and peak-merging, peaks were positively selected for an FDR value of 1E-05 Benjamini or better and minimum a NCIS normalized log2 fold difference of 2. To focus on those PLZF peaks from day 10 that were unique for differentiated cells, those peaks that had more than 0.76 fragments/kbp/M in the day 0 (naïve hMSCs) sample (quantified ±500 bp of the peak center) were excluded.

Motif analysis
De novo motif analysis was done using the MEME suite at http://meme-suite.org (Bailey et al., 2009) and the MEME-ChIP tool set (Machanick and Bailey, 2011) to search for new motifs in discriminative mode. The central enrichment of output from MEME and DREME (Bailey, 2011) was Central enrichment was visualized using the integrated Centrimo tool (Bailey and Machanick, 2012). Negative control loci were generated using EaSeq (Lerdrup et al., 2016), and the tool 'Controls' to make a set of control loci that on population level matched the PLZF peaks in distance and orientation to the closest TSS. Fasta files of sequences from the two sets of regions were generated using the 'Get sequences' tool in EaSeq.

RNA sequencing
Data received from RNA sequencing was analyzed by using Genomatix software (Expression Analysis for RNASeq Data, GGA) (Genomatix Software GmbH, Germany) as per their guidelines. Normalized gene expression was calculated using following formula (by the software): NE = c * #reads / (#reads * length) where NE is the normalized expression or enrichment value, #reads: the reads (sum of base pairs) of falling into either the transcript or the cluster region, #reads: all mapped reads (in base pairs),length: the transcript or cluster length in base pairsand c a normalization constant set to 10 7 .
ChIP-seq, RNA-seq, and microarray integration Significantly regulated genes from the microarray analysis were imported as a Regionset in EaSeq ( (Lerdrup et al., 2016) and translated into genomic coordinates using the 'Find coordinates in an already loaded geneset file' option based on gene names and a Refseq (O'Leary et al., 2016) hg19 annotation file downloaded and imported as Geneset on October 30 2018 from the UCSC table browser (Karolchik et al., 2004). Of 1,286 imported genes the gene names corresponding coordinates were not found in 187 cases, which were depicted as white in the heatmaps. The list of upregulated genes at d2 in Supplementary file 1 was imported similarly, resulting in 441 coordinates with non-zero read counts in the control. Heatmaps were made using the 'ParMap' plot type, and the color scale based on http://www.ColorBrewer.org (Harrower and Brewer, 2003). Beeswarm plotting and Mann-Whitney U-testing was done using R (R Development Core Team, 2014) and the beeswarm package (Eklund, 2016). Similarly, the list of transcript coordinates and mRNA quantities derived from the RNA-seq analysis was imported into EaSeq and used for quantitation of ChIP-seq signal at TSS of the gene encoding each transcript.

Colocalization of PLZF peaks and gene expression
To visualize expression changes of genes near PLZF peaks, For each peak the the nearest list of transcripts was identified from the imported RNA-seq analysis using the 'Colocalize'-tool was searched for the one with the closest transcript and the fold change in normalized expression of this transcript was used for visualization.

4C-Sequencing
The 4C-seq mapping and normalization was carried out using the mapping and normalization strategy as previously described (van de Werken et al., 2012b). In brief, 4C-seq reads were demultiplexed and the primers sequences without the first restriction enzyme recognition site, were removed from the reads. Subsequently, the trimmed reads were mapped to a database with in silico digested genomic fragment-ends of the human reference genome build hg19. The view point, the non-cut fragment-end and the self-ligation fragment-end were discarded. After linear interpolation of the center of the fragment-ends, the distribution of the blind-fragment reads and the non-blind fragment reads were quantile normalized using R's limma package (Smyth, 2005) . Furthermore, the 4C-seq contact frequencies within the locus of interest were normalized to its library size. The 4Cseq contact frequencies were used to calculate the trimmed mean (10%) with a running window size of 21 fragment-ends. Since 4C-seq contact frequencies have a PCR amplification bias, we carriedout a non-parametric approach by ranking each fragment-end and applying a 2 -test on the summed rank position for both samples on each running window. Subsequently, we determined the sample with the highest contact frequencies per window and corrected for multiple hypothesis testing using the Benjamini-Hochberg method (Benjamini and Hochberg, 1995). The R statistical package version 3.1.1 was used for the statistical calculations and for generating the 4C-seq plots (R Development Core Team, 2014). The R Gviz package was used for plotting the annotation data (Hahne et al., 2018).