An Unbiased Approach to Identify Cellular Reprogramming Inducible Enhancers

1. Introduction

Cellular reprogramming of somatic cells towards induced Pluripotent Stem Cells (iPSCs) via the over-expression of Oct4, Sox2, Klf4 and c-Myc (OSKM) provides an excellent model system to study the mechanisms by which a small set of transcription factors (TFs) could reactivate the pluripotency gene network to drive cell fate specification [1,2]. Reprogramming is a stochastic process, achieved only in a few cells, which manage to alter their chromatin structure to efficiently and appropriately regulate the necessary gene modules [3,4,5]. During the early stages of reprogramming, somatic cells need to undergo a series of morphological and functional changes, including among others an increase in the proliferation rate, a transition towards glycolytic metabolic pathways and a decrease in size [4,6]. In addition, mesenchymal cells undergo the Mesenchymal to Epithelial Transition process (MET) between day 2 to day 4 of reprogramming [7,8,9]. We have previously identified a gene regulatory network reconstructed from nine transcriptional regulators (TRs), the 9TR-GRN, which is required for the establishment of pluripotency [10]. The 9TR-GRN consists of Taf1c, Tead4, Tfap4, Rcan1, Cbfa2t3, Gli2, Irf6, Ovol1 and the Nanog TFs, all of which are direct OSKM targets. The 9TR-GRN is assembled in a stochastic and stepwise manner between days 1 and 6 of reprogramming only in a small number of cells expressing all factors, and is required for the gradual establishment of pluripotency [10]. These and other studies have suggested that OSKM-induced cellular reprogramming is highly plastic, as it is accompanied by the transient expression of genes from unrelated lineages, regardless of the starting cell type or species [11].

Previous landmark studies have provided insights for the OSKM-driven epigenetic and transcriptional changes occurring in the cells undergoing reprogramming [12,13,14,15,16,17,18,19,20,21]. Oct4, Sox2 and Klf4 function as pioneer transcription factors by binding to MEF-closed chromatin during the first days of reprogramming with or without c-Myc [14,18,19,22,23]. OSKM trigger chromatin opening in these sites leading to the stepwise induction of pluripotency genes either alone or in combination with other transcriptional regulators [12,13,14,15,17,18,19,20,21,24,25]. In parallel, OSKM participate in the deactivation of somatic and mesenchymal gene networks regulated by the presence of macroH2A-containing nucleosomes [26]. In addition, OSKM bind accessible and active somatic enhancers causing repression of their associated genes through various complex mechanisms including interactions with the somatic factors, whereas their concomitant binding to pluripotency enhancers activates the expression of the linked genes. OSKM-driven somatic gene repression involves chromatin deacetylation, displacement and/or down-regulation of the somatic TFs and up-regulation of negative transcriptional regulators [17,18,19,24].

Given that cellular reprogramming is a purely stochastic process, predictions regarding which cells will be reprogrammed and the timing for completion of the process are currently impossible. However, various gene markers have been proposed for the identification of pure and efficiently reprogrammable cell populations in vivo [27]. Among them, the best markers for mouse cells involve the early down-regulation of the mesenchymal surface antigen Thy1, the activated expression of alkaline phosphatase (Alpl), the upregulation of the embryonic antigen SSEA1, and the late activation of the Oct4 and/or the Nanog regulatory regions and the reconstruction of the 9TR-GRN [10,28,29,30,31,32,33]. Despite the existing methods to study cell populations entering the efficient reprogramming path(s), for identifying the actual mechanisms responsible for the gradual acquisition of pluripotency remain still elusive. Given the molecular complexity of the reprogramming process and the development of various reprogramming systems (different TF combinations, chemical molecules, genetically modified cells, etc) [27,34,35,36,37] it will be beneficial to establish new efficient ways to isolate cells during their transition to pluripotency in real-time, in order to illuminate the shady processes controlling reprogramming.

In this study, we have performed a large-scale analysis of integrated ChIP-seq, ATAC-seq and RNA-seq data and revealed the spatiotemporal highly dynamic and complex pattern of OSKM binding to the mouse genome during the course of cellular reprogramming. The combinatorial integration of multiple datasets ensures a more reliable and complete picture of the genomic OSKM binding activity, free of putative system-dependent discrepancies. We found that OSKM show distinct temporal patterns of binding at different groups of ESC enhancers. These highly dynamic and complex DNA binding events are associated with the function and expression of the neighboring genes by affecting the local chromatin dynamics. We found that genes essential for reprogramming (MET/EMT and pluripotency genes) are regulated by the coordinated activity of multiple enhancers, which are sequentially bound by OSKM to warrant strict transcriptional control. Based on these findings we developed an unbiased approach to identify Reprogramming Inducible Enhancers (RIEs) aiming to generate enhancer-traps for isolating in real-time cells undergoing reprogramming. We used a construct bearing the RIE taken from the Upp1 gene fused to GFP as a marker to isolate cells undergoing reprogramming in real-time. We demonstrated that these cells have acquired unique developmental capacities since they are reprogrammed with high efficiency due to their distinct molecular signatures, such as the high-level expression of 9TR GRN components, which are co-expressed in pre-iPSCs cells. In conclusion, our experiments have led to the development of an unbiased method to identify and isolate reprogrammable cells in real-time by exploiting the functional dynamics of OSKM that controls the transcriptional potential of RIEs, which can be used as efficient reprogramming biomarkers.

2. Results

2.1. Creation of an Integrated ChIP-Seq Dataset to Monitor Global Binding of OSKM During Cellular Reprogramming

To investigate the complex molecular mechanisms involved in TF-induced cellular reprogramming, we examined the binding of the OSKM factors to the mouse genome by integrating publicly available time-course ChIP-seq datasets [10,16,17,18,38] (Figure 1A, Table S1, Table S2) followed by computational analysis. We anticipated that a large-scale analysis of such integrated data acquired from a number of different experiments using the OSKM MEFs reprogramming system is an important parameter for obtaining a representative picture of the OSKM binding activity. Furthermore, the OSKM ChIP-seq time-course dataset (day 0, day 1, day 3, day 6 of reprogramming and in ESCs) was integrated with ATAC-seq and RNA-seq experiments in order to monitor both the requirements and the functional consequences of OSKM genomic binding. The ESCs RNA-seq data used in our analysis were obtained from our previously published dataset [26] (Table S1).

2.2. An Unprecedented Highly Dynamic OSKM Binding During Cellular Reprogramming

Our large scale analysis revealed that OSK tend to co-bind in a gradually increasing mode from day 1 to day 6 of reprogramming (Figure 1B, Oct4, Sox2 and Klf4 rows), while Myc preferentially co-binds with Klf4 only (Figure 1B, Myc rows) [14,17]. Interestingly, we discovered that the highly preferential Oct4-Sox2 co-binding in ESCs represents a stable state established late during reprogramming, because during the long reprogramming process, Sox2 shares more binding sites with Klf4 than with Oct4 (Figure 1B, compare Sox2 row with Oct4 and Klf4 columns at the different time points). Furthermore, comparison of the common Oct4 (“Oct4 to Oct4”), Sox2 (“Sox2 to Sox2”), Klf4 (‘Klf4 to Klf4”) and Myc (“Myc to Myc”) binding sites between days 1 and 3 and between days 3 and 6, revealed that ~50% of their total binding sites on day 1 are lost on day 3 and similarly, their binding sites were lost from day 3 to day 6 (Figure 1C, “D1→D3” and “D3→D6” columns). Importantly, with the exception of Klf4, OSM, especially Oct4 and Sox2, are highly mobile preserving only 10–30% of their day 6 sites in ESCs (Figure 1C, 1D, “D6→ES” column). These observations strongly suggest that the highly dynamic OSKM binding generates transient chromatin landscapes bearing unique functional potentials via the sequential occupation and abandonment of sites, which altogether, are not similar to the landscape found at the terminal and stable ESC state (Figure 1C,D). Overall, our data suggest that on the road to achieve pluripotency, OSKM are in a constant mobility phase onto chromatin rewiring various enhancer landscapes, thus providing a reasonable molecular explanation for the generation of cellular plasticity in reprogramming [39].

2.3. The Combinatorial Binding of OSKM to Early and Late Elements Prefigures the Induction of Pluripotency-Related Genes During Reprogramming

To investigate how the continuously changing OSKM DNA binding landscape during reprogramming culminates in the establishment of the unique stable DNA binding pattern in ESCs, we grouped the total sites occupied by OSKM from the beginning of reprogramming until the acquisition of pluripotency (ESCs) according to their temporal DNA binding pattern. We identified four major classes of sites: the early bound ESC sites (ESCs sites occupied from day 1), late-bound ESCs sites (ECS sites occupied after day 3 of reprogramming), transient sites (sites that are occupied transiently during reprogramming), and the MEF sites (sites bound by KM in MEFs and OSKM during reprogramming, but not in ESCs) (Figure 1E). Collectively, during reprogramming the majority of the total sites (>200,000 sites) bound by OSKM are occupied by these factors transiently only. This finding further underscores the dynamic and the stochastic nature of reprogramming.

More specifically, we found that of the total 83,674 OSKM binding sites in ESCs (Figure 2A, top panel), 31,188 (37.3%) sites are bound by OSKM from day 1 (early-bound sites) (Figure 2A, left panel), while the remaining 52,486 sites are sequentially occupied at days 3, 6, and later (late-bound sites) (Figure 2A, right panel). Furthermore, we found that the early-bound sites are preferentially located at ±0-5kb from the TSSs of the closest genes (more proximal regions), whereas the late-bound sites are primarily located between 5 and 500kb from the TSSs (distal regions) (Figure 2B,E). OSKM early binding events coincide with the rapid down-regulation and up-regulation of neighboring genes mainly involved in the somatic state and developmental and housekeeping processes, respectively, whereas the late-bound sites correlate to a more delayed activation (after day 3) of genes involved in developmental and differentiation processes (compare Figure 2C-D with Figure 2F-G).

Intriguingly, we found that the vast majority of the OSKM-bound genes in ESCs (9,992 genes) bear both early- and late-bound OSKM sites (Figure S1A, Table S3), thus suggesting that these genes could be regulated in a combinatorial manner by OSKM via widely separated cis-regulatory elements occupied at different times during reprogramming. Some of these genes are induced during reprogramming, such as the critical pluripotent regulators Pou5f1 (Oct4), Sox2 and Nanog (OSN) and various other reprogramming TFs and cofactors, like Bhlhe40, Ehf, Elf3, Esrrb, Etv5, Lin28a/b, Nr5a2, Sall1/4, Tbx3, Tfcp2l1, Utf1 and Zic3 (Table S3, Figure S1C,F). Importantly, 6 out of the 9 TRs (i.e. Rcan1, Tead4, Tfap4, Gli2, Ovol1 and Nanog) previously shown to reconstruct the 9TR gene regulatory network required for cellular reprogramming [10], are also regulated by the combinatorial action of early- and late- OSKM binding at proximal and distal regulatory elements, respectively. Of note, this group of genes includes also TFs which are down-regulated during reprogramming, like Snai1, Twist1/2 and Zeb1/2, all of which are involved in maintaining the mesenchyme phenotype [26,40], and thus inhibit reprogramming. Interestingly, although the OSKM ESC early-bound genes display a relatively stable expression pattern during the first 6-8 days of reprogramming, their expression is dramatically decreased towards ESCs, even though OSKM are still bound at their respective DNA elements (Figure S1B), whereas the late-bound OSKM ESC genes are induced at the latest phases of reprogramming (compare Figure S1C and Figure S1D) including among others genes involved in extracellular matrix (ECM) organization (Figure S1G and Table S3).

As predicted from our analysis, a significant fraction of the early-bound OSKM sites (37.4%) lie in open chromatin in MEFs (Figure 2A, left panel heatmap and Figure 2H, top panel) bearing also motifs for somatic TFs (e.g. Fosl, Jun and Tead families) (Figure S2A), whereas the OSKM late-bound elements are in a closed chromatin configuration in MEFs (Figure 2A, right panel heatmap and Figure 2H, lower panel), but contain various motifs for pluripotent TFs, such as Nanog and are linked to genes involved exclusively in cell fate decision processes (Figure 2G, Figure S2B). The vast majority of the late-bound sites (90%) remain closed throughout reprogramming (low ATAC-seq signal, Figure 2A right panel heatmap, and Figure 2H lower panel). We speculate that since these genes are bound by the master pluripotency factors Oct4, Sox2 and Nanog are in a “poised” state in ESCs, that is, ready to be expressed during differentiation upon receiving the right developmental signals [41].

The pioneering activity of OSKM in initiating reprogramming begins with their binding to the early-bound sites, a significant fraction of which is constitutively open in MEFs and remain open throughout reprogramming (23,282 sites), and a fraction of these lie within ESC-specific super-enhancers [17,42] (Figure 2I, Table S4). We also note that stable binding of OSKM in these open sites is also marked by the co-presence of somatic TF motifs, suggesting that these regions could play a binary role in both the mesenchymal and the pluripotent state (Figure 2I, Venn diagram, Figure S2C). We assume that OSKM binding to these sites eventually replaces pre-bound somatic TFs (Figure 2M, see a list of genes in Table S3). Intriguingly, a small percentage of early-bound ESCs sites, closes and re-opens at various stages of reprogramming (Figure 2J-L), thus causing dynamic changes in the gene expression programs of the associated genes (Figure 2N-P). This finding further underscores the continuously changing chromatin landscape controlled by OSKM.

2.4. Identification of Reprogramming-Inducible Enhancers in the Mouse Genome

So far, we have provided a comprehensive compilation of putative genomic regulatory regions bound by OSKM during different stages of cellular reprogramming by performing genome-wide occupation analytical studies integrated with gene expression and open chromatin signatures (Figure 1E, Figure 2). Previous studies have identified a large number of biochemical and functional interactions between O/S/K/M, thus implying that they collaboratively initiate cellular reprogramming [17,22]. This unprecedented and carefully-coordinated OSKM dynamic binding along with the binding of other TFs, including the 9TR network, lead to a sequential occupation of markedly different genomic regions between early and later stages of reprogramming, thus driving cells to abandon their mesenchymal phenotype and induce their transition towards pluripotency. For example, only the small percentage of cells in which the 9TR network is successfully reconstructed and activated by the direct action of OSKM will be reprogrammed into iPSCs [10] (Figure S3A-I). For an initial assessment of the existence and the biological relevance of putative Reprogramming Inducible Enhancers (RIEs), we searched the genome for such enhancers with in vivo reprogramming activity. Furthermore, the identification of RIEs will allow the isolation and study of the precious rare cell populations undergoing successful transition from the somatic to the pluripotent phenotype in real time. Our previous analysis identified early-bound ESC genomic sites (Figure 2A, left panel) implicated in the up-regulation of genes related to the transition to pluripotency (Figure 2C-D, Figure S1C,F). Importantly, 19,226 of these sites are bound by OSKM only after the initiation of reprogramming (Figure 3A, right panel de novo sites), and therefore we assumed that they could be contained within putative RIEs.

Based on the above, we carefully inspected these specific genomic sites to assess the value of these elements in functioning as true RIEs. More specifically, we searched the de novo acquired OSKM early-bound sites focusing on regions residing up to ±2.5kb from the TSS of up-regulated genes during reprogramming (Figure 3B). We did not take into account the biological function of the genes close to these elements, or any data regarding genomic regions previously characterized as regulatory elements. Thus, our RIEs identification approach was based on an unbiased logic. The resulting 702 sites were further filtered to remove loci bearing 10 or less Oct4, Sox2, Klf4 or Myc peaks in total. Based on this distribution, we defined 66 co-bound regions of adjacent OSKM peaks that are no more than 800bp apart on average, which we assumed that could function as putative RIEs (Table S5). Of these, three test enhancers were selected as representative examples for further functional validation analyses: a 600bp-region located ~700 bp upstream of the Lefty1 gene (“Lefty1₇₀₀” element, Figure 3C), a 300bp region located ~1800 bp upstream of the Pou5f1 gene (“Pou5f1₁₈₀₀” element, Figure 3D) and a 300bp region residing ~800 bp upstream of the Upp1 gene (“Upp1₈₀₀” element, Figure 3E). Lefty1 encodes for a protein belonging to the TGF-β family of ligands which plays a role in the determination of the left-right symmetry in the developing embryo and it is a stemness marker [43,44]. Pou5f1 encodes for the Oct4 transcription factor protein. Oct4 has a vital role in regulation of pluripotency by forming the core of the pluripotency transcription factor network along with Sox2 and Nanog [45]. Finally, Upp1 encodes for a uridine phosphorylase, an enzyme catalyzing the reversible phosphorolysis of uridine, often found up-regulated in rapidly dividing malignant cells [46,47]. All three elements are stably-bound by OSK from day 1 throughout reprogramming to ESCs, with the exception of Upp1₈₀₀, where Klf4 binding occurs from day 3 and afterwards. Importantly, all three elements are bound by Nanog in ESCs (Figure 3C-E). Furthermore, OSK early binding (day 1) correlates with an increase of the corresponding gene expression within the first 24 hours (Figure 3F-H). Taken together, these observations support the idea that the above genomic elements could indeed function as RIEs. In addition, ATAC-seq analysis revealed that these sites are occupied by nucleosomes in MEFs, but become accessible early in reprogramming (Figure S4A-C). More specifically, OSK binding at Lefty1₇₀₀ correlates with an immediate opening of the local chromatin, while the Pou5f1₁₈₀₀ and Upp1₈₀₀ become accessible after day 3. Accordingly, although at the beginning of reprogramming all three putative RIEs are marked by the repressive histone modification H3K27me3, during reprogramming there is a gradual replacement of the repressive mark by the H3K27ac modification, a transcriptionally active promoter/enhancer marker, a finding consistent with their function as RIEs (Figure S4D-F) [17]. Interestingly, all three RIEs coincide with ESC-related Multiple Transcription-factor binding Loci (MTLs) - regions short in size, bound by multiple TFs, which are believed to act as sites of enhanceosome assembly in ESCs (Figure 3C-E, bottom). The elements have been also annotated as ES-specific Super-enhancer regions (Figure 3C-E, bottom) [42,48]. Taken together, the Lefty1₇₀₀, Pou5f1₁₈₀₀ and Upp1₈₀₀ identified by our unbiased approach (Figure 3B), share many structural characteristics with functional enhancer elements and they are part of previously known regulatory elements in ESCs and therefore, they could indeed function as true RIEs.

2.5. The Upp1₈₀₀ Element Functions as a Reprogramming Inducible Enhancer Marking Cells Undergoing Reprogramming to Pluripotency

To investigate whether the three putative regulatory elements act as true RIEs, we generated reporter constructs in which a minimal promoter driving the expression of the GFP gene was fused to the test enhancers immediately upstream, that is, in a position similar to their natural genomic location (Figure 4A, left panel). The resulting constructs were cloned into a vector suitable for the generation of lenti-viral particles. Next, the produced lenti-viruses bearing the putative RIE reporters were used to transduce MEFs with the OKSM transgene incorporated into the Col1a1 locus and the rtTA*M2 transactivator under the ROSA 26 promoter [49,50]. Interestingly, we found that initiation of cellular reprogramming with the addition of DOX is accompanied by the activation of the putative RIEs. We observed a robust GFP expression driven by all three putative RIEs during reprogramming (Figure 4B, Figure S5A-B). More specifically, during the first 3 days of reprogramming, a gradually increasing number of GFP-expressing cells (GFP(+) cells) appeared within the population of cells undergoing reprogramming. Morphological examination revealed that the majority of the GFP(+) cells had abandoned their mesenchymal phenotype as they appeared smaller, round, with an increased number of intercellular contacts (especially since day 3) (arrows in Figure 4B and Figure S5A-B). This is a profound alteration characterizing cells that undergo MET [8]. Importantly, many of the GFP(+) cells cluster within the early iPSC colonies, especially those cells expressing the Lefty1₇₀₀ and Upp1₈₀₀ reporters (see day 6 in Figure 4B and Figure S5A). In contrast, expression of the Pou5f1₁₈₀₀ reporter appears to be less restricted to the early iPSC colonies. These experiments strongly suggest that the Lefty1₇₀₀, Pou5f1₁₈₀₀ and Upp1₈₀₀ function as ESC RIEs and they mark cells undergoing phenotypical transitions during reprogramming. Significantly, the Lefty1₇₀₀ and Upp1₈₀₀ reporters, not only specifically mark activation of GFP gene expression after initiation of reprogramming, but they also mark the cells that will form the early iPSCs colonies (Figure 4B, Figure S5A). On the other hand, the Pou5f1₁₈₀₀ RIE suffices to warrantee transcriptional activation in the context of reprogramming only, without marking specifically the early iPSCs formations (Figure S5B).

Next, we asked whether these putative RIEs suffice to regulate the expression of their adjacent genes during reprogramming. Thus, we compared the expression pattern of the endogenous genes (Lefty1, Pou5f1 and Upp1) with the expression of the exogenous GFP controlled by the respective regulatory elements (Figure 4C, Figure S5C-D). Remarkably, we found that both the Lefty1 and Upp1 GFP reporters display expression patterns that are similar to their endogenous counterparts, thus implying that these elements are the main regulatory enhancers controlling the expression of these genes. Of note, the small temporal difference (1-2 days) observed between the GFP and the endogenous gene expression pattern could be explained by the fact that as the exogenous RIE-GFP construct can be incorporated in various random open chromatin genomic sites, it may lack an appropriate epigenetic landscape, which normally could delay the endogenous gene expression. Since the Lefty1₇₀₀ element simulates the expression pattern of the endogenous gene during the first 6 days only, followed by a gradual decrease in its activity, it implies that this element is not sufficient to support Lefty1 gene expression throughout the entire course of reprogramming, but it suffices only for the initial stages of the process (Figure S5C). On the other hand, the Pou5f1₁₈₀₀ element did not mimic the expression program of the endogenous gene (Figure S5D), thus suggesting that Pou5f1 expression is controlled by additional regulatory elements during reprogramming.

2.6. The Cells “Marked” by the Reprogramming Inducible Enhancer Upp1₈₀₀ Element Achieve Earlier and More Robust Induction of the 9TR Network Leading to Efficient Reprogramming

We described above the integration of multiple features that led us to conclude that the Upp1₈₀₀ and Lefty1₇₀₀ enhancers function as autonomous RIEs activated early during reprogramming marking the cells undergoing MET in early iPSC colonies (Figure 4B, Figure S5A). Therefore, we used the RIE Upp1₈₀₀ GFP reporter as a tool to isolate the cells undergoing reprogramming by monitoring GFP expression in order to explore their transcriptional and phenotypical properties (Figure 4A, right panel). We chose the Upp1₈₀₀ GFP reporter because it accurately recapitulates the expression of the endogenous Upp1 gene throughout reprogramming (Figure 4C). FACS-isolated Upp1₈₀₀-GFP(+) cells and GFP(-) cells on day 4 of reprogramming were re-plated in parallel and their reprogramming efficiency was calculated. We found that the GFP(+) cells produced two times more AP-positive iPSC colonies as compared to the GFP(-) cells (Figure 4D). This finding is in agreement with the reporter assays, where we observed high level expression of Upp1₈₀₀ driven GFP expression in early iPSC colonies (Figure 4B, Figure S5E).

To investigate the molecular signature of the Upp1₈₀₀-GFP(+) cells that enables them for efficient reprogramming, as opposed to the GFP(-) cells, we performed RNA-seq experiments side by side with the GFP(-) cells used as a control to “capture” the underlying molecular pathways supporting highly efficient reprogramming. Comparison of the RNA-seq data between GFP(+) and GFP(-) cells revealed that these two subpopulations are characterized by distinct transcriptiome profiles, since we identified a significant number of up- and down-regulated Differentially Expressed Genes (DEGs) (Figure 4E-F). More specifically, the GFP(+) cells, where the Upp1₈₀₀ RIE is activated, bear 718 up-regulated and 394 down-regulated DEGs as compared to the GFP(-) cells (Figure 4F, Table S6). Functional enrichment analysis (GSEA) revealed that the up-regulated genes are related to intercellular adhesion and various signaling processes, whereas down-regulated genes are involved in chemotaxis and response to chemokines, processes indicative of the complex cell-cell interactions and the phenotypic alterations occurring during the early days of reprogramming (Figure 4G).

To further investigate the reconstruction of molecular networks marking the reprogrammable cells, we used the MCODE algorithm [51] to cluster groups of DEGs involved in the various biological processes described above (Figure S6). This detailed analysis verified that indeed in the GFP(+) cells, clusters of down-regulated genes are formed which are involved in chemotaxis, in responses to extracellular signals (Figure S6, Clusters 1-2), in cell migration and EMT (Figure S6, Clusters 1-2, 4-5, 9). For example, Bmp4, an EMT-promoting gene functioning as a reprogramming inhibitor [52,53] is down-regulated in GFP(+) cells, thus relieving the cells from a reprogramming barrier (Figure S6, Cluster 5). On the contrary, genes involved in the establishment of cell adhesion and epithelial functions are up-regulated, including among others Dcn and various basement membrane components like laminin and collagen IV genes (Figure S6, Clusters 5, 9). These genes function by facilitating MET, a requirement for reprogramming. Furthermore, genes involved in embryonic development and morphogenesis are also up-regulated in Upp1₈₀₀-GFP(+) cells (Figure S6, Clusters 5, 7-8). Of these, an important example includes the telomerase reverse transcriptase gene (Tert) (Figure S6, Cluster 5). This finding further validates our experimental approach since regulation of telomere length by the induced Tert is critical for the achievement of the high pluripotency potential of ES and iPS cells [54,55,56,57]. In this context, we also identified that Tpp1, a gene regulating Tert binding to telomeres during MEF reprogramming [58], is also up-regulated in the Upp1₈₀₀-GFP(+) cells (Table S6). As expected, we showed that a number of Hox TFs, previously implicated in high efficiency reprogramming [59], are also up-regulated in the RIE-GFP(+) cells (Figure S6, Cluster 7). Taken together these data further validate our previous results indicating that the RIEs could mark specific cells in a population, which undergo global morphological and molecular alterations in the beginning of reprogramming. More specifically, we have been able to “tag” cells that undergo MET in which specific sets of genes involved in the transition have been activated.

We have previously shown that the small percentage of cells completing reprogramming reconstruct a dynamic gene regulatory network known as 9TR-GRN [10]. We found that the TFs Gli2 and Irf6, which are direct targets of OSKM and components of the 9TR-GRN are up-regulated in the Upp1₈₀₀ GFP(+) cells (Figure 4H, Figure S5F, S6I). Next, we sorted Upp1₈₀₀-GFP(+) cells and GFP(-) cells on day 4 of reprogramming followed by replating to evaluate their molecular characteristics at later time points. First, we examined the EMT/mesenchymal markers Snail1 and Vimentin (Vim) along with the 9TR-GRN members Irf6 and Nanog, 12 days after cell sorting and re-plating (Figure 4I, Figure S5G). We found that the GFP(+) cells bear decreased levels of the EMT master regulator Snail1, despite the fact that Snail1 expression was nearly equal between the two cell populations on day 4 before re-plating (Figure S5G, left panel). This observation suggests that the GFP(+) cells, which over-express Gli2 and Irf6 early after OSKM induction, abandon the mesenchymal phenotype more efficiently. Interestingly, the classic EMT marker Vim is down-regulated in the GFP(+) cells as early as from day 4 (when the sorting is performed) (Figure S5G, right panel), in agreement with our previous data indicating that the Upp1-GFP(+) cells have epithelial characteristics (Figure 4B, 4G, Figure S6E,I). Intriguingly, Irf6 is expressed at 4 times higher levels in Upp1₈₀₀-GFP(+) than in GFP(-) cells on day 4, but later these expression levels are similar (Figure 4I). This implies that the RIE Upp1₈₀₀ marks cells acquiring a temporal “advantage” over the rest of the population in reconstructing the 9TR-GRN network. As a result, these cells could acquire earlier and more efficiently various pluripotent molecular components. Consistent with this is our observation that 12 days after sorting, the GFP(+) cells express more Nanog transcripts than the GFP(-) cells (Figure 4I). Taken together, the data described above strongly suggest that the RIE Upp1₈₀₀ is an efficient enhancer-trap for the early identification and isolation of cells entering the reprogramming trajectory.

3. Discussion

Cellular reprogramming is a stepwise stochastic process characterized by the emergence of multiple intermediate cell states towards MET from which a small percentage of cells continues to be reprogrammed to achieve pluripotency [3,4,5,10]. Reactivation of the pluripotency gene regulatory network by OSKM requires the rewiring of many gene regulatory elements to permanently inactivate somatic enhancers, to transiently activate enhancers of genes required for initiation of the process and finally to activate the pluripotency enhancers. In this study, we have highlighted the complexity of OSKM action by creating a composite OSKM ChIP-seq dataset [10,16,17,18,38] used to map in detail the association of the reprogramming factors with the genome to reveal how they orchestrate the gradual acquisition of pluripotency. Contrary to many previous studies, we focused on the temporal and permanent OSKM binding patterns to different genomic elements aiming to examine how the OSKM ESC binding pattern is achieved and established during reprogramming. Within this context, we identified novel putative regulatory elements bound by OSKM in a complex and dynamic fashion, some of which are capable of functioning as RIEs. We have exploited the transcriptional regulatory property of these elements and isolated cells undergoing reprogramming in real-time. Their characterization not only independently validated our previously identified 9TR-GRN as a critical player of reprogramming, but importantly, it also revealed the simultaneous reconstruction of additional reprogramming-specific gene networks controlling the global morphological and molecular alterations required for achieving pluripotency.

We discovered that of the 83,674 genomic sites occupied by OSKM in ESCs, 11,962 sites are pre-bound in MEFs by the endogenously expressed KM, and then they are immediately co-occupied, at least some of them, by OS from day 1 (Early ESC sites). These results suggest that pre-bound KM recruit OS via cooperative interactions. This is an important observation because the role of Myc is usually overlooked in reprogramming and it is often omitted from the transcription-factor “cocktails” [60,61]. However, even in the absence of Myc co-expression with OSK, it is still present in MEFs and its inhibition has dramatic effects in the reprogramming efficiency [21,60].

In contrast to the stably bound OSKM throughout reprogramming, a large number of genomic sites are occupied by OSKM (individually or as a group) in a temporal manner. Apparently, during reprogramming, OSKM become highly mobile by constantly binding to new sites and/or abandoning others inducing alterations in chromatin structure and the expression of the nearby genes. Thus, the dynamic binding of OSKM is directly related to the re-organization of chromatin and the gradual establishment of the pluripotent epigenetic landscape. It is conceivable to assume that at least in part, the inherent stochasticity of cellular reprogramming is due to the low percentage of cells that will be able to establish the correct pattern of dynamic OSKM binding suitable for orchestrating the correct transcriptional programs towards pluripotency.

An interesting implication of the dynamic OSKM binding patterns described herein stems from the fact that Klf4 appears to regulate the formation of “hubs” composed of pluripotency genes and ESC Super-enhancers [38]. Thus, we propose that the temporal pattern of Klf4 binding during reprogramming shown here may regulate the instruction of different hubs consisting of different enhancers associated transiently with the same or different genes during reprogramming. Based on the above, we propose that single genes may “use” multiple OSKM-bound enhancer elements, loaded at different time-points, to ensure the execution of an accurate reprogramming expression pattern.

Our extensive characterization of the various classes of OSKM-bound elements and their association with pluripotency genes prompted us to test whether some of these elements can be used as enhancer traps to isolate and study the small percentage of cells that are stochastically “chosen” to be reprogrammed. We developed and applied an unbiased method based on the OSKM global binding pattern for the identification of putative RIEs among the early-bound regions. We identified 66 putative RIEs, including the Lefty1₇₀₀, Pou5f1₁₈₀₀ and Upp1₈₀₀ elements. We showed that all three elements function as true reprogramming inducible enhancers. Interestingly, these elements map within ESC-specific Super-enhancers [42] and within loci that have been previously implicated in regulation of gene expression [62,63,64], thus validating our unbiased approach used to identify RIEs.

We found that, at least one of these RIEs, the Upp1₈₀₀ element when fused to a minimal promoter-GFP cassette can accurately reproduce the expression of the endogenous gene specifically in the cells undergoing reprogramming. Therefore, we used the Upp1₈₀₀-GFP reporter to isolate GFP(+) and GFP(-) cells during the first days of reprogramming. Importantly, the GFP(+) isolated cells have twice as much the capacity to complete their reprogramming process as compared to the control GFP(-) cells, since we showed that they acquire more efficiently the epithelial phenotype. This observation is also supported by our finding that genes involved in MET are significantly enriched in the Upp1₈₀₀-GFP(+) cells (e.g. Dcn, Col4a1/2, Lama2, Pdgfrb, Serpine1/2, Slit3, Tagln, Thbs2 etc.). Importantly, we found that the GFP(+) cells express at high levels the 9TR-GRN members Gli2 and Irf6, as early as on day 2 of reprogramming, leading to a more robust activation of the 9TR network. As Nanog induction is the final component required for reconstructing the 9TR-GRN network its efficient up-regulation paves the route to pluripotency (Figure 4J). In summary, by using the OSKM complex DNA binding patterns during reprogramming together with gene expression data as guides, we were able to accurately identify functional enhancers induced in reprogramming (RIEs). We provided evidence that RIEs can be used as effective biomarkers for the identification and isolation of rare cell populations that are being appropriately “prepared” to enter reprogramming successfully very early in the process. This may have important implications in the development of methods to achieve high reprogramming efficiencies for medical applications.

4. Materials and Methods

4.1. Experimental Protocols

4.1.1. MEFs Isolation from Mouse Embryos

For our cellular reprogramming experiments we used mouse embryonic fibroblasts (MEFs) from cross-bred homogygous transgenic mice originating from the B6;129S4-Col1a1^{tm1(tetO-Pou5f1,-Klf4,-Sox2,-Myc)Hoch}/J (RRID:IMSR_JAX:011001) and B6.Cg-Gt(ROSA)26Sor^{tm1(rtTA*M2)Jae}/J (RRID:IMSR_JAX:006965) strains (“4F-MEFs”) [49,50]. MEFs were isolated as described before [10]. In brief, mouse embryos E13.5-14.5 were harvested from euthanized pregnant female mice. After removing the head, the heart and the liver, the emrbyonic tissue was chopped up using a razor blade and turned into single-cell suspension by incubation in 0.25% Trypsin solution in PBS (Gibco™, #15090046) and multiple passes through 18G and 21G syringes. MEFs were cultured at 37^οC and 5% CO2 concentration, in high-glucose DMEM (Sigma-Aldrich, D6429) supplemented with heat-inactivated FBS at 10% final concentration (Pan Biotech, P30-1985), 1x GlutaMAX™ Supplement (Gibco™, #35050-061), 1x MEM Non-Essential Amino Acids Solution (Gibco™, #11140-050) and Penicillin-Streptomycin mix diluted at 1/100 (Gibco™, #15140-122).

4.1.2. Cellular Reprogramming Protocol

An appropriate number of “4F-MEFs” (depending on the protocol) was placed in plates (day 0) and cellular reprogramming (OKSM over-expression) was induced with the addition of doxycycline (DOX) at 2ug/ml (Sigma-Aldrich, D9891) in MEF medium as described above, but supplemented with 15% FBS. The medium was refreshed every two days. Upon the emergence of the first early iPSC colonies (around day 8) the medium was replaced with ESC medium: KnockOut™ DMEM (Gibco, #10829018), 20% Pansera FBS designed for ESCs (Pan Biotech, P30-2602), 1x GlutaMAX™ Supplement (Gibco™, #35050-061), 1x MEM Non-Essential Amino Acids Solution (Gibco™, #11140-050), Penicillin-Streptomycin mix diluted at 1/100 (Gibco™, #15140-122) and 20ng/mL mLIF (Santa Cruz Biotechnology, sc-4378).

4.1.3. Nanog Chromatin Immunoprecipitation Followed by High-Throughput Sequencing

ChIP was performed as described before [10]. Cells were cross-linked with 1% formaldehyde at room temperature for 10 min, and then quenched with 0.125M glycine for 5 min. Cells were washed with PBS and harvested with scraper. Lysis was performed in Lysis buffer [50mM Hepes (pH 7.9), 140mM NaCl, 1mM EDTA, 10% glycerol, 0.5% NP-40, 0.25% Triton X-100], washed in Wash buffer [10mM Tris-HCL (pH 8.1), 200mM NaCl, 1mM EDTA (pH 8.0), 0.5mM EGTA (pH 8.0)] and then resuspended in Sonication buffer [0.1% SDS, 1mM EDTA, 10mM Tris (pH 8.1)]. Chromatin was sheared to 200-500bp fragments in the Covaris S2 sonicator using TC12×12mm tubes (Tube AFA Fiber and Cap, Covaris) for 12 min (200 cycles per burst). After the sonication, we added Triton X-100 and NaCl to reach a final concentration of 1% and 150mM, respectively. The chromatin was then centrifuged and filtered via a 0.2um filter syringe. 10ug chromatin was incubated with 3,6ug rabbit IgG (crude serum) and 10ug anti-Nanog (Cell Signaling Technology Cat# 8785, RRID:AB_11220438) at 4^oC overnight. The next day pre-equilibrated protein G DynabeadsTM (Thermo Fisher Scientific, 10004D) in IP buffer [0.1% SDS, 1mM EDTA, 10mM Tris (pH 8.1), 1% Triton X-100, 150mM NaCl] were mixed with the chromatin and antibody solutions. The mix was left to incubate for 1 hour at room temperature. Next, the beads were washed with Low salt [0.1% SDS, 1% Triton X-100, 2mM EDTA, 20mM Hepes-KOH (pH 7.9), 150mM NaCl], High salt [0.1% SDS, 1% Triton X-100, 2mM EDTA, 20mM Hepes-KOH (pH 7.9), 500mM NaCl] and LiCl buffer [100mM Tris-HCl (pH 7.5), 0.5M LiCl, 1% NP-40, 1% Sodium deoxycholate]. The immunoprecipitated chromatin was digested with Proteinase K (Roche Life Sciences, 03115828001) at 50^οC for 15 min, followed by incubation with RNase A at 65°C overnight and an extra step of Proteinase K for 1 hour at 50^οC. The released DNA was cleaned-up using NucleoMag clean-up and selection kit (Macherey-Nagel, 744970) and the elution was performed in TE/10 buffer. Eluted DNA was then used for downstream analyses.

The Nanog ChIP-seq library construction and sequencing were carried out at the Greek Genome Center (GGC) of BRFAA. Libraries were generated with NEBNext® Ultra II™ DNA Library prep Kit for Illumina® (E7645L, NEB), as indicated by the manufacturer’s guidelines using 10ng of DNA (chromatin that was not incubated with the antibody was used as input DNA). The quality of the generated libraries was validated with the Agilent Bioanalyzer DNA 1000 chip. Subsequently, the libraries were quantitated with the Qubit™ High Sensitivity (HS) spectrophotometric method and pooled in equimolar amounts for sequencing. Approximately, 25 million, 101bp-long, single-end reads were generated with Illumina NovaSeq 6000 sequencer. The data of the Nanog ChIP-seq experiment are deposited in the Gene Expression Omnibus data repository under the GEO accession ID GSE274131.

4.1.4. ATAC-Seq

ATAC-seq libraries from different time-points were prepared starting from 50,000 cells, following previously published protocols [65,66,67,68]. Briefly, cell pellets were lysed in Lysis Buffer [10mM Tris-HCl, pH 7.5, 10mM NaCl, 3mM MgCl2, 0.1% NP-40, 0.1% Tween-20, 0.01% Digitonin (Promega, G9441)] for 3 min on ice, washed with Wash Buffer [10mM Tris-HCl, pH 7.5, 10mM NaCl, 3mM MgCl₂, 0.1% Tween-20] and centrifuged at 500g for 10 min at 4^oC. Each pellet (nuclei) was resuspended in 50ul Transposition reaction mix [1x Tagment DNA Buffer (Illumina, 15027866), 16.5ul 1x PBS, 0.1% Tween-20, 0.01% Digitonin (Promega, G9441), 2.5ul Tn5 Transposase (Tagment DNA Enzyme 1, Illumina, 15027865)] and incubated at 37^oC for 30 min on a thermomixer at 1000rpm. DNA purification was performed using the Qiagen MinElute Reaction Cleanup Kit (Qiagen, 28204). Library purification was performed with the AMPure XP beads (Beckman Coulter, A63880). The quality of libraries was assessed using an Agilent 2100 Bioanalyzer and libraries were quantified with the Qubit™ High Sensitivity (HS) spectrophotometric method (Invitrogen). Libraries were sequenced paired-end using a NovaSeq 6000 sequencer at the Greek Genome Center at BRFAA. The data of the ATAC-seq experiment are deposited in the Gene Expression Omnibus data repository under the GEO accession ID GSE274130.

4.1.5. Construction of Enhancer Reporters and Production of Lenti-Viral Particles

Cloning of regulatory elements and construction of enhancer reporters

We isolated the three identified enhancer elements (Lefty1₇₀₀, Pou5f1₁₈₀₀, Upp1₈₀₀) from the mouse genome via PCR using specific primer pairs (Table S7). The primers were designed to enclose the identified OSKM peaks of the respective RIE element in the smallest amplicon size possible. For the construction of our enhancer reporters we utilized the plasmidic vector LeGO-G2 [69]. LeGO-G2 was a gift from Boris Fehse (Addgene plasmid #25917; http://n2t.net/addgene:25917; RRID: Addgene_25917). The strong regulatory element SFFV upstream of the EGFP gene was removed and replaced each time by one of the three isolated enhancer elements and a minimal promoter sequence. As minimal promoter was used either the endogenous proximal promoter sequences of each of the three genes (Table S7), or the SCP1 element [70] which was isolated from the pSTARR-seq_human vector [71]. pSTARR-seq_human was a gift from Alexander Stark (Addgene plasmid #71509; http://n2t.net/addgene:71509; RRID: Addgene_71509). The two types of promoters (SCP1 and endogenous) did not affect the pattern of EGFP expression (data not shown) and thus, were used interchangeably in our reporter assays. LeGO-G2 was used as a positive control in our experiments since it showed stable GFP expression over-time.

4.1.6. Generation of Lenti-Viral Particles

To introduce our constructs into the MEF genome we produced lenti-viral particles via transfection in HEK293T cells (Human Embryonic Kidney 293T; RRID:CVCL_0063). The HEK293T cells were transfected with a combination of the plasmids pMD2.G, psPAX2 and each one of our constructs (separately) using the classic calcium phosphate method. The day after the transfection the supernatant was replaced with fresh medium and 3 days after the medium change, the lenti-viral particles were harvested. HEK293T were cultured at 37^οC and 5% CO2 concentration, in high-glucose DMEM (Sigma-Aldrich, D6429) supplemented with heat-inactivated FBS at 10% final concentration (Pan Biotech, P30-1985), 1x GlutaMAX™ Supplement (Gibco™, #35050-061), 1x MEM Non-Essential Amino Acids Solution (Gibco™, #11140-050) and Penicillin-Streptomycin mix diluted at 1/100 (Gibco™, #15140-122). Plasmids pMD2.G (Addgene plasmid #12259; http://n2t.net/addgene:12259; RRID:Addgene_12259) and psPAX2 (Addgene plasmid #12260; http://n2t.net/addgene:12260; RRID:Addgene_12260) were a gift from Didier Trono.

4.1.7. Reporter Assays

• MEFs transduction and cellular reprogramming

For the reporter assays we transduced ~0.8-1x10⁶ early-passage “4F-MEFs” with the above described lenti-viral particles in the presence of polybrene at a final concentration of 6ug/ml (Sigma-Aldrich, H9268) overnight, followed by change of medium. After 2-3 days, an appropriate number of transduced-MEFs (depending on the assay) was re-plated and cellular reprogramming was induced with the addition of DOX (day 0), as described above.

2.

Microscopy

~600.000 “4F-MEFs” transduced with the enhancer reporters were plated in 10cm plates and cellular reprogramming was induced after 24 hours by the addition of DOX, as described above. Cells were observed regularly under the microscope for the monitoring of the reprogramming process and the detection of GFP(+) cells. A Leica DM IRE2 inverted microscope was used (phase contrast and fluorescence microscopy) with an ORCA-Flash 4.0 LT digital camera (Hamamatsu, C11440-42U, SN 001432). The images were captured with the HCImage Live software (Hamamatsu) and processed using FIJI/ImageJ [72]. Brightness adjustments were made uniformly to all the images of a single experiment in order to correctly depict the increase/decrease of GFP intensity between the different time-points of each experiment. The pixel intensities of the images in Figure S5E were transformed using the LOG function to be able to depict both time-points with very low and very high GFP intensity.

3.

Reporter assay for the comparison of exogenous-GFP and endogenous gene expression pattern (Figure 4C, Figure S5C-D)

~70.000 “4F-MEFs” transduced with the enhancer reporters were plated in 6-well plates and cellular reprogramming was induced the next day by the addition of DOX, as already described. For each enhancer-reporter two biological replicates were prepared. Cells were collected using 1x Trypsin solution (Gibco™, 15090046) at different time-points for RNA extraction.

4.

Sorting of Upp1₈₀₀-GFP(+) and GFP(-) cells for transcriptome analysis and reprogramming efficiency calculation

1.5-3x10⁶ transduced “4F-MEFs” were plated in 10cm plates and reprogramming was initiated after 24 hours (day 0). At least two biological replicates were prepared for each experiment. At day 2 (for day 2 Upp1₈₀₀ RNA-seq) or day 4 (for efficiency calculation and transcriptome analysis in Figure 4I and Figure S5G) cells were harvested using 1x Trypsin solution (Gibco™, 15090046), washed with PBS and stained with 0.3uM DAPI for 5-10 min (Roche, 236276). Next, cells were washed again with PBS and resuspended in cold Sorting buffer [1x PBS, 5% FBS, 1mM EDTA] to a final concentration of 5-10x10⁶ cells/ml. Cell aggregates were removed prior to sorting with 50um strainer cap filters (BD Biosciences, 340629). GFP(+) and GFP(-) sorted cells aimed for day 2/day 4 transcriptome analyses (RNA-seq and analysis in Figure 4I and S5G, respectively) were used directly for RNA extraction (100-200.000 cells per sample). For reprogramming efficiency calculation alive GFP(+) and GFP(-) cells were counted using Trypan blue and plated in equal numbers on 12-well plates (70-80.000 cells per well). After letting cells to overcome stress and expand for almost 10-12 days, they formed healthy iPSC colonies. To compensate for inequalities in cell duplication rates caused by the stress of the sorting process, before estimating the reprogramming efficiency we re-plated equal numbers of GFP(+) and GFP(-) cells in 6-well plates (100-200.000 cells per well) covered with gelatin (Sigma-Aldrich, G1393) and feeders (C57Bl/6 MEFs treated with 10ug/ml Mitomycin C for 2.5-3 hours). iPSC colonies were left to grow and were cultured 1 week without DOX, prior to AP staining. GFP(+) and GFP(-) sorted cells aimed to be used for transcriptome analysis at a later time-point (“12 days after sorting” in Figure 4I and S5G) were re-plated on 12-well plates (100.000 cells per well) and used for RNA-extraction 12 days later.

5.

Alkaline Phosphatase (AP) staining for reprogramming efficiency calculation

Reprogramming cell cultures were subjected to Alkaline Phosphatase (AP) staining for the identification of iPSC colonies as described before [10]. In brief, cells growing on plates were washed with PBS and fixed with 4% formaldehyde at 4^oC for 10 min. The cells were washed with NTMT buffer [100mM Tris-HCl, 100mM NaCl, 50mM MgCl₂ and 0.1% Tween-20, pH 9.5] and then stained for AP activity, using NBT/BCIP substrate solution (Roche, 11681451001) diluted 1/100 in NTMT buffer. Reprogramming efficiency was calculated by dividing the number of AP(+) iPSC colonies with the total number of cells which were seeded on the feeder-covered plates and then, by multiplying by 100. Reprogramming efficiency evaluation was performed in biological triplicates.

6.

RNA isolation and real-time PCR

Cells were lysed with NucleoZOL (740404, Macherey-Nagel) and their RNA was extracted following the manufacturer’s instructions. Reverse transcription followed to produce cDNA with the PrimeScript™ RT Reagent Kit (Perfect Real Time) (TaKaRa, RR037A) using oligodT primers and random 6mers. The cDNA quantification was performed in the CFX96 Touch™ Real-Time PCR Detection System (Bio-Rad) using the 2x KAPA SYBR® FAST Universal mix (KAPA, KK4618) and the appropriate primer pair (Table S7). Calculation of expression values was performed with the 2^(-ΔCt) method and Gapdh served as the reference gene. GFP expression values of the enhancer-reporters in Figure 4C and Figure S5C-D were normalized to the GFP values of LeGO-G2 at the respective time-points, to compensate for putative fluctuations in GFP expression caused by the intense epigenetic alterations occuring in reprogramming at the random chromatin integration sites of the reporters.

7.

RNA sequencing

• Reprogramming time-course

“4F-MEFs” (day 0) and cells undergoing reprogramming were isolated at different time-points (day 1, day 2, day 3, day 4, day 6, day 8) and their RNA was extracted using NucleoZOL (740404, Macherey-Nagel), as by the manufacturer’s instructions. Two biological replicates were prepared for each sample. The RNAs were DNAseI-treated and cleaned-up using the classic phenol-chloroform procedure. RNA-seq experiments were carried out in the Greek Genome Center (GGC) of the Biomedical Research Foundation of the Academy of Athens (BRFAA). RNA-seq libraries were prepared using the NEBNext® Ultra II™ Directional RNA Library prep Kit for Illumina® (E7760, NEB) with 1ug of total RNA input. Library QC was performed with the Agilent Bioanalyzer DNA 1000 kit and quantitation with the Qubit™ High Sensitivity (HS) spectrophotometric method (Invitrogen). Sequencing was performed in the Illumina NovaSeq 6000 sequencer and we obtained on average 22x10⁶ single-end, reverse-stranded, 101bp-long reads per sample. The data of the RNA-seq experiment are deposited in the Gene Expression Omnibus data repository under the GEO accession ID GSE274132.

• Upp1₈₀₀-GFP sorted cells

GFP(+) and GFP(-) cells collected after sorting (as described above) were subjected to RNA extraction using the NucleoSpin RNA kit (740955, Macherey-Nagel). The DNA was removed by a DNAseI on-column step, as instructed by the kit’s manual. The library preparation and sequencing was performed as described above. Two biological replicates were sequenced per sample. The data of the RNA-seq experiment are deposited in the Gene Expression Omnibus data repository under the GEO accession ID GSE279976.

4.2. Bioinformatics Analyses

4.2.1. ChIP-Seq Data Analysis

8.

Selection and primary analysis of published O/S/K/M and Nanog ChIP-seq datasets

For our O/S/K/M ChIP-seq analysis we utilized multiple datasets derived either from our laboratory (current publication and [10]), or by others [16,17,18,38], deposited under the GEO IDs GSE114581, GSE274131 GSE90893, GSE67520, GSE101905 and GSE113429 (Table S1). We selected datasets which use the 4-factor reprogramming system (OKSM) and MEFs as the starting population. Prior to merging with other datasests, we re-analyzed our O/S/K/M ChIP-seq data under the GEO ID GSE114581 [10], wishing to take advantage of the latest adjustments to the existing algorithms. In more detail, Fastq files were aligned against the mm10 mouse genome with Bowtie2 aligner [73] and the output files were filtered for low-quality reads, duplicates, blacklist regions and mitochondrial genome entries with custom-made pipeline from BEDtools [74] and SAMtools [75]. Peak calling was performed with MACS2 callpeak command [76]. Bigwig files were constructed with the deepTools bamCompare tool [77], using the “subtract“ option to remove Input signal from the respective IP signal, and choosing RPKM for normalization. The same pipeline was performed for our Nanog ChIP-seq experiment analysis, but first, the quality of the Fastq files was examined with the FastQC software [78]. Finally, for the production of bigwig files for Figure 3 and Figure S3, raw data of the O/S/K/M published datasets which were used in this study (GSE90893, GSE67520, GSE101905, GSE113429) were downloaded from the SRA database and analyzed as described above.

9.

Merging of datasets

Our O/S/K/M ChIP-seq dataset (GSE114581) was re-analyzed as described above, while the other ChIP-seq peaks were used as published (GSE90893, GSE67520, GSE101905, GSE113429) (see Table S1 for more details). The peaks (bed files) were combined per time-point for each individual transcription factor (using cat and sort functions) and overlaps between peaks from different datasets were combined using the BEDtools mergeBed tool (option “-d 300”) [74]. For the day 1 time-point we combined the 18 hours and day 1 data and for the day 6 time-point we combined the day 5 and day 6 data. We also combined data from both high- (SSEA1+) and low- (Thy1+) efficiency cells -when available, since we are interested in studying the bulk population of cells undergoing reprogramming. For the analysis of OSKM binding as a group (Figure 1E and Figure 2), we merged the Oct4, Sox2, Klf4 and Myc peaks per time-point, as we did above (cat, sort, mergeBed -d 300). The peaks of the combined dataset are presented in Table S2.

For the visualization of the combined dataset in the Integrative Genomics Viewer (IGV) [79] (Figure 3 and S3) we used the normalized bigwig files from each dataset, as produced by the analysis described above. Then, the different bigwig files were merged and averaged for each factor per time-point, using the WiggleTools mean tool [80].

10.

ATAC-seq analysis

Bioinformatics analyses were performed using The Galaxy Suite [81]. The quality of the sequencing reads was evaluated using the FastQC algorithm [78]. Reads were trimmed to 50bp using the Trimmomatic tool [82]. Sequencing reads were mapped to the mm10 version of the mouse genome using the Bowtie2 algorithm with the “very sensitive end to end” option [73]. The Samtools Fixmate command was used [75] prior to duplicate removal with the MarkDuplicates command from Picard tools using the option “do not write duplicates to the output file”. Reads from all samples were downsampled to the same sequencing depth using the Downsample BAM option from Picard tools. Peaks were called using the MACS2 algorithm with the options “no model”, “extension-200”, “shift-100” and a q-value cut-off of 0.05 [76]. Bigwig files were constructed using the bamCoverage command from the Deeptools suite [83].

4.2.2. Investigation of OSKM Binding Sites

11.

Calculation of common binding sites between Oct4, Sox2, Klf4 and Myc

For the analysis presented in Figure 1B-D we utilized our combined dataset after downsampling the Oct4 day 1, day 3 and day 6 datasets. This was necessary in order to avoid a bias over Oct4 overlaps, due to the higher number of Oct4 peaks in our combined dataset. We downsampled randomly the number of Oct4 peaks to match the average number of the Sox2 and Klf4 peaks of the respective time-point, using a custom-made pipeline (sed command). For the calculation of the common sites between each pair of the OSKM factors we used BEDTools intersect (-wa -u options) [74].

12.

Identification of ESC, transient and MEF sites

For this analysis (Figure 1E and Figure 2) we used the merged OSKM peaks and not the individual Oct4, Sox2, Klf4 and Myc sites (see ChIP-seq analysis above). All sites where OSKM were present at the ESC stage were considered as “ESC OSKM binding sites”. For the identification of the “early-bound” and “late-bound” categories of ESC sites presented in Figure 1E and 2A we performed intersection of all the ESC sites with the day 1 OSKM binding sites using the BEDTools intersect tool (-wa -u and -v options, respectively) [74]. For the early-bound ESC subcategories in Figure 2I-L, we performed successive intersections as described above (-wa -u or -v options), with the day 3 and day 6 OSKM binding sites, accordingly. The same was done for the identification of MEF “pre-bound” and “de novo sites” presented in Figure 3A.

For the OSKM “Transient sites” (Figure 1E) we followed a similar logic of successive intersections, as described above, in order to combine all the different categories of transiently-bound sites (sites where OSKM bind at least at one intermediate time-point, but neither in MEFs, nor in ESCs). In more detail, for the transient sites where OSKM bind only at one time-point (i.e. day 1, day 3, or day 6) we removed (-v option) from the respective bed files (i.e. the total day 1, day 3, or day 6 OSKM peaks) the overlaps with the ESC OSKM sites and subsequently, the overlaps with the other two remaining intermediate time-points (i.e. we removed the day 3-day 6, day 1-day 6 and day 1-day 3 OSKM peaks, respectively). Finally, we removed the common peaks of the resulting files with the KM MEF binding sites. The remaining sites were considered the day 1-, day 3- and day 6-only transient OSKM sites (data not shown). The transient sites where OSKM bind at two intermediate time-points (i.e. day 1-day 3, day 1-day 6, and day 3-day 6 sites, data not shown) were calculated by pairwise intersections of the respective time-point OSKM peaks. In each case, for every overlap we kept the original sites of both intersected files and merged them (cat, sort and mergeBed with default options). Then, we removed the overlaps with the remaining intermediate time-point (i.e. day 6, day 3 or day 1, respectively), with ESCs and subsequently with MEFs. Finally, for the transient sites where OSKM bind at all intermediate time-points (day 1-day 3-day 6 sites, data not shown), we intersected in a stepwise manner the three respective OSKM files and for every overlap we kept and merged the original sites of all intersected files. From the resulting file we kept only the sites where OSKM were not also present in ESCs and MEFs. The final file containing all “Transient sites” (shown in Figure 1E) was produced by merging all the above categories of sites (cat, sort and mergeBed with default options).

For the “MEF sites, lost at ESCs” (Figure 1E) we combined the different categories of sites where (OS)KM bind at MEFs but not in ESCs. In more detail, the sites bound only at MEFs (data not shown), were calculated by keeping the KM MEF peaks which did not overlap with ESCs and neither with any of the other time-points in reprogramming. The rest of the MEF sites, which are lost progressively in reprogramming (data not shown) were calculated by merging all the transient site files before the removal of the overlaps to the MEF peaks and intersecting them with the MEFs KM peaks. The final file containing all “MEF sites, lost at ESCs” (shown in Figure 1E) was produced by merging the two above categories of sites (cat, sort and mergeBed with default options).

Using multiIntersect tool, we ensured that the identified (sub)categories of the ESC, transient and MEF OSKM-bound sites did not overlap with each other. Thus, all OSKM sites presented in this study belong to a single (sub)category of sites.

13.

Overlaps of OSKM sites with ESC Super-enhancers

For the calculation of overlaps between the ESC Super-enhancers [42] and the various types of OSKM ESC sites presented in Table S4 we used the BEDTools intersect tool, as above (-wa -u and -c options) [74].

14.

Epigenetic characterization of OSKM binding sites

The percentage of statistically significant ATAC-open and ATAC-closed chromatin regions (heatmaps and Venn diagrams in Figure 2A,I-L) was calculated using the BEDTools intersect tool for each category of sites and for all ATAC-seq examined time-points. The tornado plots were constructed with the plotHeatmap tool of deepTools (v3.5.2) [77]. The matrices that were used as input for the plotHeatmap were calculated with computeMatrix reference-point mode and the options “--referencePoint center” “--averageTypeBins mean”, “--missingDataAsZero”, “-bs 5”.

15.

Gene assignment to OSKM binding sites

Assignment of the various OSKM binding sites to genes was performed with GREAT on-line tool (v4.0.4) [84,85] using the whole genome as background and the “Single nearest gene” option for 1000kb maximum distance. For the calculation of common gene sets between the different OSKM ESC site subcategories (Figure S1, Table S3) we used the Venny tool (v2.1) [86].

16.

Motif analysis

Our de novo motif analysis was performed with HOMER [87]. We used the options “-size given” and “-preparse”. The motifs presented in this paper have a final enrichment p-value<1e-10, appear at least at 5% of the target sequences and have a score ≥0.9.

4.3. Method for the Identification of Putative Reprogramming Inducible Enhancers

To identify Reprogramming Inducible Enhancers (RIEs) in the mouse genome we used the de novo acquired early-bound ESC sites (Figure 3A, right panel). Using GREAT (v4.0.4) [84,85] and the “Single nearest gene” option for 2.5kb maximum distance, we found that 2865 genes are associated with the above regions. Among them, we chose only the sites close to the 610 genes which are up-regulated between MEFs and ESCs, with Log₂FC>2 and p-adjusted<0.05 (702 sites). Then, we intersected the resulting sites with all the individual Oct4, Sox2, Klf4 and Myc peaks in reprogramming using BEDTools (-c option) [74] to identify and remove regions with a low number of binding events (=<10 peaks). The remaining 66 sites were considered as putative RIE elemements (Table S5).

4.4. RNA-Seq Analysis

The quality of the sequenced samples was evaluated by the FastQC tool [78]. Then, the FASTQ files were aligned to the UCSC mm10 genome version with the HISAT2 tool (v2.1.0) using the options “-5 4 -3 4” to trim the 5’ and 3’ ends of the reads and “--rna-strandness R” [88]. Annotation of the SAM files to produce the raw counts was performed with the HTSeq htseq-count tool (v0.6.1p1 and 2.0.3) [89] and the options “-s reverse” “-m intersection-nonempty”, using the UCSC mm10refGene.gtf genome annotation file. Normalized counts were calculated using DESeq2 (v1.30.1, estimateSizeFactors) [90]. To calculate the normalized counts for the reprogramming time-course we used also the raw counts of the sample “ESCs wild type” as they are uploaded in the GEO database (GSE215883, Table S1). Genes with zero normalized counts at all time-points were excluded when calculating the median expression value of multiple genes. DEGs were considered the genes with Log₂FC>0.58 or <-0.58 and p-adjusted<0.05.

4.5. Functional Enrichment of Gene Sets

For the functional enrichment analyses presented in Figure 2D, 2G and 4G we used Webgestalt [91], while the analyses in Figure S1 were performed with EnrichR [92,93,94]. Finally, for the networks analysis in Figure S6 we utilized the STRING Enrichment tool of the stringApp [95] in Cytoscape [96]. We utilized the libraries GO Biological Process 2023 and GO Molecular Function 2023 [97,98], KEGG Pathways [99,100,101], Reactome Pathways [102], MSigDB Hallmark 2020 [103,104] and PanglaoDB Augmented 2021 [105].

4.6. Network Analysis

For the network analysis of the day 2 Upp1800-GFP(+)/(-) RNA-seq dataset presented in Figure S6 we inserted all the DEGs (Log₂FC>0.58 or <-0.58 and p-adjusted<0.05) in Cytoscape STRING protein query [96] to create a “full STRING network” with confidence cut-off 0.4. Then, the network was clustered with the MCODE application [51] using the default parameters and without the “haircut” option. The Log₂FC value of each DEG node was assigned to a blue-red color gradient (blue for down-regulated genes, red for up-regulated genes). The color scale has the same minimum and maximum values for all networks appearing in Figure S6.

4.7. Other Graphical Representations

Ggplot2 package [106] was used in R (v4.0.5) for the construction of the gene expression line-plots, the functional enrichment and motif analyses dot-plots, the PCA and volcano plots. All heatmaps were designed using the pheatmap package in R. Statistical analyses and plots in Figure 4 and S5 were performed using GraphPad Prism version 8.0.2 for Windows, GraphPad Software, Boston.