PLZF targets developmental enhancers for activation during osteogenic differentiation of human mesenchymal stem cells
Abstract
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 PLZF-dependent 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, HIF3α, 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 B) was an early 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 (pre-osteoblasts).
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 PLZF-binding 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 by Ivins et al. (Figure 2—figure supplement 1A and B) (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 (Figure 2—figure supplement 1C–E).
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 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 FANTOM5 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 (Figure 2—figure supplement 2B).
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 ChIP-seq 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 a small fraction of these regions showed PLZF-dependent changes that are in line with the previous reported function of PLZF as a transcriptional repressor (Supplementary file 7).
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 siRNA-mediated 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 RNA-seq 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 pre-osteoblasts 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 lineage-specific 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.
-
Figure 4—source data 1
- https://cdn.elifesciences.org/articles/40364/elife-40364-fig4-data1-v2.xlsx
-
Figure 4—source data 2
- https://cdn.elifesciences.org/articles/40364/elife-40364-fig4-data2-v2.xlsx
-
Figure 4—source data 3
- https://cdn.elifesciences.org/articles/40364/elife-40364-fig4-data3-v2.xlsx
-
Figure 4—source data 4
- https://cdn.elifesciences.org/articles/40364/elife-40364-fig4-data4-v2.xlsx
-
Figure 4—source data 5
- https://cdn.elifesciences.org/articles/40364/elife-40364-fig4-data5-v2.xlsx
-
Figure 4—source data 6
- https://cdn.elifesciences.org/articles/40364/elife-40364-fig4-data6-v2.xlsx
-
Figure 4—source data 7
- https://cdn.elifesciences.org/articles/40364/elife-40364-fig4-data7-v2.xlsx
-
Figure 4—source data 8
- https://cdn.elifesciences.org/articles/40364/elife-40364-fig4-data8-v2.xlsx
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 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 genome-wide 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.
-
Figure 6—source data 1
- https://cdn.elifesciences.org/articles/40364/elife-40364-fig6-data1-v2.xlsx
-
Figure 6—source data 2
- https://cdn.elifesciences.org/articles/40364/elife-40364-fig6-data2-v2.xlsx
-
Figure 6—source data 3
- https://cdn.elifesciences.org/articles/40364/elife-40364-fig6-data3-v2.xlsx
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 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.
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 contact the promoter of the NNMT gene through chromatin looping and affects SAM homeostasis in osteogenesis (see model Figure 7E).
Materials and methods
Human mesenchymal stem cells (hMSC) culture and differentiation
Request a detailed protocolBone-marrow-derived primary human mesenchymal stem cells (hMSC) were purchased from Lonza (Lonza, Cat. No. PT-2501). Cell purity and their ability to differentiate into osteogenic, chondrogenic and adipogenic lineages were tested by Lonza. Cells were positive for CD105, CD166, CD29, and CD44. Cells tested negative for CD14, CD34 and CD45. Cells were cultured according to manufacturer’s instructions using mesenchymal stem cell growth medium (MSCGM) (Lonza MSCGM: PT-3238) supplemented with one MSCGM SingleQuots (PT-4105). Differentiation was induced using StemPro osteogenesis (Gibco A10072-01), chondrogenesis (Gibco A10071-01) or adipogenesis (Gibco A10070-01) differentiation kit as per manufacturer’s instructions. Cells were harvested at the indicated time points for ChIP, RNA, 4C-seq and protein using appropriate buffers as described below.
PLZF knockdown using siRNA
View detailed protocolMission siRNA oligos were ordered from Sigma. Three different oligos were tested and sequences are provided in the Supplementary file 9.
The hMSCs were reverse transfected with siRNA oligos using Lipofectamin 2000 and Optimem media (Gibco). Briefly, 600 pmol of siRNA (in 500 μl Optimem was mixed with 10 μl Lipofectamin 2000 (Invitrogen) in 500 μl Optimem. After 15 min of incubation at room temperature, 1 × 106 cells were mixed in the transfection mix and incubated for another 10 min. Cells in the transfection mix were seeded in dishes (15 cm, Nunc cell culture) with normal HMSCB media. Osteogenic differentiation was induced 24 hr after transfection and cells processed for analyses at indicated time points.
Chromatin immunoprecipitation (ChIP) assay
View detailed protocolCells were processed for ChIP as previously described with slight modifications (Dietrich et al., 2012). Briefly, cells were fixed for 10 min in culture media containing 1% formaldehyde and were processed for ChIP as previously described with slight modifications (Dietrich et al., 2012). Briefly, formaldehyde fixed chromatin was sonicated in IP buffer (IP buffer = 2 volumes of SDS Buffer: 1 vol Triton Dilution Buffer; SDS Buffer- 50 mM Tris-HCl, (pH8.1), 100 mM NaCl, 5 mM EDTA, (pH 8.0), 0.2% (w/v) NaN3, 2% (w/v) SDS; Triton Dilution Buffer-: 100 mM Tris-HCl, (pH 8.6), 100 mM NaCl, 5 mM EDTA, (pH 8.0), 0.2% (w/v) NaN3,5.0% (v/v) Triton X-100) using a Branson Sonifier (4 cycles of 30 s each at 22% of max amplitude). Twenty μg DNA (sonicated chromatin) was used for each ChIP in IP buffer. Antibodies used: Rabbit general IgG (DAKO), rabbit monoclonal anti-H3K27me3 (Cell Signalling # 9733, Rabbit monoclonal,), anti-SUZ12 (Cell Signalling # 3737 Rabbit monoclonal), RING1B (produced in our lab), H3K4me3 (Cell Signalling # 9751, Rabbit monoclonal), H3K27ac (Abcam # ab4729), H3K4me1 (Abcam # 8895), (Cell Signalling # 5326S, Rabbit monoclonal used for ChIP-seq), H3K36me3 (Cell Signalling # 4909, Rabbit monoclonal), PLZF (Santa Cruz # sc-22839), P300 (Santa Cruz # sc-585), RNA POLII (Santa Cruz # sc 9001X), MED1/TRAP220 (BETHYL # A300-793A), MED12 (BETHYL # A300-774A), JMJD3 (produced in our lab)). Samples were first incubated with antibodies overnight at 4°C by end-over-end rotation and subsequently immune-complexes were enriched using incubation with protein A-Sepharose beads for 3 hr. After washing the samples were de-crosslinked overnight at 68°C (shaking) in ChIP Elution Buffer (20 mM Tris-HCl (pH7.5), 5 mM EDTA (pH 8.0), 50 mM NaCl, 1% (w/v) SDS, 50 μg/ml Proteinase K). Finally, enriched DNA was purified using a Qiagen PCR purification kit. The ChIPs were validated at Polycomb target genes by QPCR.
ChIP sequencing
Request a detailed protocolChIP 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
Request a detailed protocolTemplates 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 × 107 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 μg 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
Request a detailed protocolSAM 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
Request a detailed protocolTotal 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
Request a detailed protocolThe 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 μg reporter plasmid together with 15 μg PAX8 and 6 μg of VSV in 293FT cells seeded at a density of 5 × 106 cells per 15 cm dish (Nunc cell culture) using the CaPO4 method. 48 hr following transfection and media change the supernatant containing viral particles were collected and filtered through 0.45 μm filters. Polybrene was added to the supernatant 8 μg/ml and hMSCs were infected with virus overnight. Cells were trypsinised and reseeded 48 hr later and Puromycin (0.5 μg/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.
Western blotting
Request a detailed protocolWestern blotting was performed according to standard procedures. The cells were lysed in high salt lysis buffer (50 mM Tris-HCl (pH 7.2), 300 mM NaCl, 0.5% (w/v) Igepal CA-630, 1 μg/ml aprotinin, 1 μg/ml leupeptin, 1 mM PMSF, 1 mM EDTA (pH 8.0), 50 mM NaF, 20 mM β-glycerophosphate) followed by a brief sonication on a Branson Sonifier with a 2 mm probe (10% of max amplitude for 2 s on ice). Samples were left on ice for 30 min, then centrifuged at 22,000 g for 15 min. Ten percent Tris-glycine or 4–12% Bis-Tris SDS-PAGE gels (Invitrogen) were used to separate proteins for analyses. Seablue Plus two prestained marker (Invitrogen) was used as a molecular weight standard. Blotting was performed using the following antibodies: PLZF (Santa Cruz # sc-22839), NNMT (Abcam # ab58743), GAPDH (Abcam # ab8245), cyclinE (mab, HE12), pRB2 (Santa Cruz, sc317).
Alizarin red staining
Request a detailed protocolNaive 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
Request a detailed protocolReads 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 non-redundant 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
Request a detailed protocolSpecific 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 µ +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
Request a detailed protocolDe 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
Request a detailed protocolData 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 107.
ChIP-seq, RNA-seq, and microarray integration
Request a detailed protocolSignificantly 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
Request a detailed protocolTo 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
Request a detailed protocolThe 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 4C-seq 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 carried-out 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).
Data availability
Sequencing data have been made available via GEO under accession number GSE125168. All other data generated or analysed during this study are included in the manuscript and supporting, source files: Supplementary tables associated with Figure 1, Figure 2 and Figure 3 are uploaded with the submission. CAGE-derived enhancer candidates were taken from http://slidebase.binf.ku.dk/.
-
NCBI Gene Expression OmnibusID GSE125168. PLZF targets developmental enhancers for activation during osteogenic differentiation of human mesenchymal stem cells.
-
European Nucleotide ArchiveID DRA000991. A promoter-level mammalian expression atlas.
References
-
MEME SUITE: tools for motif discovery and searchingNucleic Acids Research 37:W202–W208.https://doi.org/10.1093/nar/gkp335
-
DREME: motif discovery in transcription factor ChIP-seq dataBioinformatics 27:1653–1659.https://doi.org/10.1093/bioinformatics/btr261
-
Inferring direct DNA binding from ChIP-seqNucleic Acids Research 40:e128.https://doi.org/10.1093/nar/gks433
-
Plzf regulates limb and axial skeletal patterningNature Genetics 25:166–172.https://doi.org/10.1038/76014
-
Controlling the false discovery rate: a practical and powerful approach to multiple testingJournal of the Royal Statistical Society: Series B 57:289–300.https://doi.org/10.1111/j.2517-6161.1995.tb02031.x
-
The PRC1 Polycomb group complex interacts with PLZF/RARA to mediate leukemic transformationGenes & Development 23:1195–1206.https://doi.org/10.1101/gad.512009
-
Modification of enhancer chromatin: what, how, and why?Molecular Cell 49:825–837.https://doi.org/10.1016/j.molcel.2013.01.038
-
The multitalented Mediator complexTrends in Biochemical Sciences 38:531–537.https://doi.org/10.1016/j.tibs.2013.08.007
-
PLZF is a regulator of homeostatic and cytokine-induced myeloid developmentGenes & Development 23:2076–2087.https://doi.org/10.1101/gad.1788109
-
ColorBrewer.org: an online tool for selecting colour schemes for mapsThe Cartographic Journal 40:27–37.https://doi.org/10.1179/000870403235002042
-
Inducible pRb2/p130 expression and growth-suppressive mechanisms: evidence of a pRb2/p130, p27Kip1, and cyclin E negative feedback regulatory loopCancer Research 60:2737–2744.
-
The promyelotic leukemia zinc finger promotes osteoblastic differentiation of human mesenchymal stem cells as an upstream regulator of CBFA1Journal of Biological Chemistry 280:8523–8530.https://doi.org/10.1074/jbc.M409442200
-
The UCSC Table Browser data retrieval toolNucleic Acids Research 32:493D–496.https://doi.org/10.1093/nar/gkh103
-
Zinc-finger transcription factor odd-skipped related 2 is one of the regulators in osteoblast proliferation and bone formationJournal of Bone and Mineral Research 22:1362–1372.https://doi.org/10.1359/jbmr.070602
-
Enrichr: a comprehensive gene set enrichment analysis web server 2016 updateNucleic Acids Research 44:W90–W97.https://doi.org/10.1093/nar/gkw377
-
An interactive environment for agile analysis and visualization of ChIP-sequencing dataNature Structural & Molecular Biology 23:349–357.https://doi.org/10.1038/nsmb.3180
-
Sequence-specific DNA binding and transcriptional regulation by the promyelocytic leukemia zinc finger proteinThe Journal of Biological Chemistry 272:22447–22455.https://doi.org/10.1074/jbc.272.36.22447
-
Normalization of ChIP-seq data with controlBMC Bioinformatics 13:199.https://doi.org/10.1186/1471-2105-13-199
-
Transcriptional regulatory cascades in Runx2-dependent bone developmentTissue Engineering Part B: Reviews 19:254–263.https://doi.org/10.1089/ten.teb.2012.0527
-
In vitro osteogenic potential of human mesenchymal stem cells is predicted by Runx2/Sox9 ratioTissue Engineering. Part A 21:115–123.https://doi.org/10.1089/ten.TEA.2014.0096
-
Building strong bones: molecular regulation of the osteoblast lineageNature Reviews Molecular Cell Biology 13:27–38.https://doi.org/10.1038/nrm3254
-
MEME-ChIP: motif analysis of large DNA datasetsBioinformatics 27:1696–1697.https://doi.org/10.1093/bioinformatics/btr189
-
Critical residues within the BTB domain of PLZF and Bcl-6 modulate interaction with corepressorsMolecular and Cellular Biology 22:1804–1818.https://doi.org/10.1128/MCB.22.6.1804-1818.2002
-
The RUNX2 cistrome in osteoblasts: characterization, down-regulation following differentiation, and relationship to gene expressionThe Journal of Biological Chemistry 289:16016–16031.https://doi.org/10.1074/jbc.M114.552216
-
Osteoblast physiology in normal and pathological conditionsCell and Tissue Research 343:289–302.https://doi.org/10.1007/s00441-010-1086-1
-
SoftwareVEENY. An interactive tool for comparing lists with Venn diagramsBioinfo GPCNB-CSIC.
-
SoftwareR: A language and environment for statistical computingFoundation for Statistical Computing, Vienna, Austria.
-
Nicotinamide N-methyltransferase in non-small cell lung cancer: promising results for targeted anti-cancer therapyCell Biochemistry and Biophysics 67:865–873.https://doi.org/10.1007/s12013-013-9574-z
-
Metabolism: Cancer mistunes methylationNature Chemical Biology 9:293–294.https://doi.org/10.1038/nchembio.1234
-
Booklimma: Linear Models for Microarray Data Linear Models for Microarray DataIn: Irizarry RA, Carey V. J, Huber W, Gentelman R, Dudoit S, editors. Statistics for Biology and Health, Bioinformatics and Computational Biology Solutions Using R and Bioconductor. Springer Link. pp. 397–420.
-
Transcription factors: from enhancer binding to developmental controlNature Reviews Genetics 13:613–626.https://doi.org/10.1038/nrg3207
-
Foxo1, a novel regulator of osteoblast differentiation and skeletogenesisJournal of Biological Chemistry 285:31055–31065.https://doi.org/10.1074/jbc.M109.079962
-
Peripheral leptin regulates bone formationJournal of Bone and Mineral Research 28:22–34.https://doi.org/10.1002/jbmr.1734
-
NNMT promotes epigenetic remodeling in cancer by creating a metabolic methylation sinkNature Chemical Biology 9:300–306.https://doi.org/10.1038/nchembio.1204
-
4C technology: protocols and data analysisMethods in Enzymology 513:89–112.https://doi.org/10.1016/B978-0-12-391938-0.00004-5
-
Components of the SMRT corepressor complex exhibit distinctive interactions with the POZ domain oncoproteins PLZF, PLZF-RARalpha, and BCL-6Journal of Biological Chemistry 273:27695–27702.https://doi.org/10.1074/jbc.273.42.27695
-
Model-based analysis of ChIP-Seq (MACS)Genome Biology 9:R137.https://doi.org/10.1186/gb-2008-9-9-r137
Article and author information
Author details
Funding
Forskerakademiet (09-073526)
- Shuchi Agrawal Singh
Novo Nordisk (NNF17CC0027852)
- Kristian Helin
Danish National Research Foundation (DNRF 82)
- Kristian Helin
- Klaus Hansen
Memorial Sloan-Kettering Cancer Center (Support Grant NIH P30 CA008748)
- Kristian Helin
Lundbeckfonden (R34-A3479)
- Klaus Hansen
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Joanna Wysocka for providing the lentiviral GFP reporter vector (pSINMIN) and Prof. Carsten Muller-Tidow for critical reading and discussion on the manuscript. Thanks to the members of Hansen group for support and fruitful discussions, and Kathryn Wattam for excellent technical assistance. HJGvdW was supported by a Zenith grant (93511036) from the Netherlands Genomics Initiative (NGI). SAS was supported by FSS (Danish Research Council). KH was supported by the Danish National Research Foundation (DNRF82) and through a center grant from the Novo Nordisk Foundation (NNF17CC0027852). KHA was supported by the Lundbeck Foundation and the Danish National Research Foundation.
Copyright
© 2019, Agrawal Singh et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 1,903
- views
-
- 326
- downloads
-
- 29
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Chromosomes and Gene Expression
Chromatin immunoprecipitation (ChIP-seq) is the most common approach to observe global binding of proteins to DNA in vivo. The occupancy of transcription factors (TFs) from ChIP-seq agrees well with an alternative method, chromatin endogenous cleavage (ChEC-seq2). However, ChIP-seq and ChEC-seq2 reveal strikingly different patterns of enrichment of yeast RNA polymerase II (RNAPII). We hypothesized that this reflects distinct populations of RNAPII, some of which are captured by ChIP-seq and some of which are captured by ChEC-seq2. RNAPII association with enhancers and promoters - predicted from biochemical studies - is detected well by ChEC-seq2 but not by ChIP-seq. Enhancer/promoter-bound RNAPII correlates with transcription levels and matches predicted occupancy based on published rates of enhancer recruitment, preinitiation assembly, initiation, elongation, and termination. The occupancy from ChEC-seq2 allowed us to develop a stochastic model for global kinetics of RNAPII transcription which captured both the ChEC-seq2 data and changes upon chemical-genetic perturbations to transcription. Finally, RNAPII ChEC-seq2 and kinetic modeling suggests that a mutation in the Gcn4 transcription factor that blocks interaction with the NPC destabilizes promoter-associated RNAPII without altering its recruitment to the enhancer.
-
- Chromosomes and Gene Expression
- Microbiology and Infectious Disease
Candida glabrata can thrive inside macrophages and tolerate high levels of azole antifungals. These innate abilities render infections by this human pathogen a clinical challenge. How C. glabrata reacts inside macrophages and what is the molecular basis of its drug tolerance are not well understood. Here, we mapped genome-wide RNA polymerase II (RNAPII) occupancy in C. glabrata to delineate its transcriptional responses during macrophage infection in high temporal resolution. RNAPII profiles revealed dynamic C. glabrata responses to macrophages with genes of specialized pathways activated chronologically at different times of infection. We identified an uncharacterized transcription factor (CgXbp1) important for the chronological macrophage response, survival in macrophages, and virulence. Genome-wide mapping of CgXbp1 direct targets further revealed its multi-faceted functions, regulating not only virulence-related genes but also genes associated with drug resistance. Finally, we showed that CgXbp1 indeed also affects fluconazole resistance. Overall, this work presents a powerful approach for examining host-pathogen interaction and uncovers a novel transcription factor important for C. glabrata’s survival in macrophages and drug tolerance.