Dynamics of chromatin marks and the role of JMJD3 during pancreatic endocrine cell fate commitment

Induction of mature pancreatic β cells in vitro is an effective way to treat diabetes mellitus. Because the developmental mechanisms of pancreatic development, especially pancreatic endocrine development, are not fully understood, artificially generated β cells remain functionally immature (Johnson, 2016; McKnight et al., 2010; Pagliuca and Melton, 2013). Therefore, defining the regulatory mechanism that controls endocrine cell fate commitment is extremely important in the study of pancreatic development.

During pancreatic development, all endocrine lineages arise from multipotent pancreatic progenitor cells. The developmental potential of these progenitors is gradually constrained to the establishment of exocrine and islet endocrine lineages (Shih et al., 2013). The basic helix-loop-helix transcription factor (TF) neurogenin 3 (Neurog3 or Ngn3) plays a central role in endocrine lineage specification, and lineage-tracing experiments have shown that all endocrine lineages, and even some exocrine cells, are derived from Ngn3-expressing progenitors (Gu et al., 2002; Schonhoff et al., 2004). Levels of Ngn3 are also important for driving endocrine lineage specification, as demonstrated by the fact that genetically ablating Ngn3 expression or modifying the Ngn3 dose can affect the induction of other key TF genes involved in endocrine cell specification, including Neurod1, Isl1, Pax4 and Pax6 (Gradwohl et al., 2000), or can cause Ngn3-induced cells to follow exocrine cell fate (Beucher et al., 2012; Wang et al., 2010). These findings indicate that early-stage Ngn3-expressing cells retain their plasticity to differentiate into multiple pancreatic lineages.

During mouse embryogenesis, Ngn3 is expressed at a relatively low level in the dorsal endoderm domain beginning at embryonic day (E) 8.5, and this expression is maintained until E11.0. The second wave of Ngn3⁺ cell specification in both the dorsal and ventral pancreas begins at ∼E12.0, peaks at E15.5, and then abruptly decreases (Gradwohl et al., 2000; Schwitzgebel et al., 2000). However, Ngn3 is transiently expressed in individual cells and undergoes a transition from low expression (Ngn3-low) to high expression (Ngn3-high) during the second wave of endocrine differentiation. In addition, Ngn3-low cells are cell-cycle-active progenitors, whereas most Ngn3-high cells are cell-cycle-silent endocrine-committed cells (Bankaitis et al., 2015; Bechard et al., 2016). A recent EdU pulse-chase study demonstrated that pancreatic endocrine progenitor cells take ∼12 h to pass through the Ngn3⁺ state (Bankaitis et al., 2015).

Chromatin modifications play important roles in regulating cell fate transitions. During early developmental stages, both activating (H3K4me3) and repressive (H3K27me3) epigenetic marks are frequently located in the promoter regions of developmentally poised regulatory genes (Bernstein et al., 2006; Mikkelsen et al., 2007). Upon differentiation, these regions, which are known as bivalent domains, must be resolved to facilitate either rapid transcriptional activation or stable gene silencing (Voigt et al., 2013). The histone modifications H3K4me1 and H3K27ac are both chromatin marks of enhancers and show dramatic dynamics during cell fate differentiation. However, it is unclear how these chromatin modifications coordinate and promote the quick cell fate transition from Ngn3-low to Ngn3-high.

In this study, we performed transcriptomic analysis of sorted Ngn3-low and Ngn3-high cells from E13.5 Ngn3-GFP mouse embryos and discovered dramatic changes in gene expression between these two populations. Based on low-cell-number ChIP-seq analysis, we identified that the removal of a repressive H3K27me3 histone mark was linked to the activation of key TFs for the specification of Ngn3-high cells. The activation of one TF, NeuroD1, is sequentially accompanied by the deposition of histone modifications associated with active and primed enhancers. Moreover, we demonstrate that the histone H3K27 demethylase JMJD3 regulates the efficiency of the transition from Ngn3-low to Ngn3-high.

To profile the global changes in gene expression associated with dorsal pancreatic endocrine lineage specification in mice, we performed RNA-seq analyses on E13.5 Ngn3-low and Ngn3-high cells from Ngn3-GFP pancreata (Lee et al., 2002), and on E17.5 nascent β cells from Ins1-RFP pancreata (Piccand et al., 2014). All of these cells were purified via fluorescence-activated cell sorting (FACS) (Fig. S1A,B and Table S1). Although passing through the Ngn3⁺ state of pancreatic endocrine progenitor cells only takes ∼12 h (Bankaitis et al., 2015), hierarchical clustering revealed a striking difference in the gene expression profiles between Ngn3-low and Ngn3-high cell populations. The gene expression profile of E13.5 Ngn3-high cells more closely resembled those of differentiated E17.5 β cells than those of E13.5 Ngn3-low cells (Fig. S1C). These findings indicate that a significant cell fate transition towards pancreatic endocrine development occurred between Ngn3-low and Ngn3-high stages.

A comparative transcriptomic analysis identified 5197 genes that were differentially expressed between Ngn3-low and Ngn3-high cell populations (Fig. S1D and Table S2), among which 154 genes encode TFs (Fig. S1E). The expression levels of Sox9, Onecut1 and Hes1 (Fig. S1E) were enriched in Ngn3-low cells, which is consistent with the progenitor property of Ngn3-low cells (Shih et al., 2013), whereas the expression levels of many key TFs for endocrine development, such as Ngn3, Neurod1, Insm1, Rfx6, Isl1 and Pax6 were elevated in Ngn3-high cells (Ahlgren et al., 1997; Gierl et al., 2006; Mellitzer et al., 2006; Naya et al., 1997; Sander et al., 1997; Smith et al., 2010; St-Onge et al., 1997) (Fig. S1E). Cell cycle regulatory and DNA replication genes were enriched in Ngn3-low cells but displayed a drastic change during the transition from Ngn3-low to Ngn3-high (Fig. S1F,G and Table S2). KEGG pathway (Kanehisa and Goto, 2000) enrichment analysis revealed that Ngn3-high cells were enriched for pathways related to pancreatic functions (Fig. S1F and Table S2). Taken together, these transcriptomic analyses indicate that during the transition from Ngn3-low to Ngn3-high, endocrine progenitor cells lose their progenitor properties and begin specification.

To obtain a genomic view of the dynamics of chromatin bivalency during the endocrine fate transition, we mapped the modifications of H3K4me3 and H3K27me3 in E13.5 Ngn3-low and Ngn3-high cells using low-cell-number ChIP-seq (Brind'Amour et al., 2015) (Fig. S2A,B and Table S1). In Ngn3-low cells, 2131 genes displayed bivalency at their promoter regions [transcription start sites (TSSs) ±2 kb] (Fig. 1A). During the transition from Ngn3-low to Ngn3-high, bivalency was resolved in 590 of these genes, whereas 226 genes gained bivalent marks (Fig. 1B, top). Of the genes showing resolved bivalency in Ngn3-high cells, 80% (471 genes) lost only the H3K27me3 mark (Fig. 1B, top). However, only 49 genes became bivalent via the sole gain of H3K27me3 (Fig. 1B, top). To determine whether endocrine cells maintain this unbalanced loss and gain of bivalency-associated H3K27me3 during later development, we compared the dynamics of bivalent marks between E13.5 Ngn3-high and E17.5 β cells (Fig. 1A). Curiously, the numbers of genes gaining H3K27me3, H3K4me3 or both marks were similar to the numbers of genes losing those marks (Fig. 1B, bottom). This finding indicates that the removal of bivalency-associated H3K27me3 is a specific characteristic of the transition from Ngn3-low to Ngn3-high.

The composite plot shows that the removal of bivalency-associated H3K27me3 did not change the overall corresponding H3K4me3 level (Fig. 1C). However, according to the RNA-seq analyses of Ngn3-low and Ngn3-high cell populations, the genes associated with the loss of bivalency-associated H3K27me3 showed a net increase in their expression levels (Fig. 1D). Notably, this group of genes included key TFs for endocrine development, such as Neurod1, Rfx6 and Insm1 (Gierl et al., 2006; Mellitzer et al., 2006; Naya et al., 1997; Smith et al., 2010) (Fig. 1E,F and Table S3), but did not include Ngn3, which indicates that the regulation of Ngn3 is independent of bivalency resolution (Fig. 1G). These studies indicate that the resolution of bivalent domains, mainly by the removal of H3K27me3, is related to pancreatic endocrine cell specification.

We also examined the changes of H3K27me3 levels at promoters (TSSs ±2 kb) marked by H3K27me3 alone. To our surprise, during the transition from Ngn3-low to Ngn3-high, 124 promoter loci showed decreases (>2-fold changes) in H3K27me3 levels, whereas only six regions displayed increases in H3K27me3 levels (Fig. 2A,B). However, during the transition from Ngn3-high cells to β cells, we observed that the number of genes adding H3K27me3 (137 genes) was higher than the number removing H3K27me3 (36 genes) (Fig. 2B), which indicates a substantial increase in the number of promoters marked by H3K27me3 alone during β cell differentiation. Furthermore, the decrease in the H3K27me3 level at promoters during the transition from Ngn3-low to Ngn3-high was associated with an increase in H3K4me3 level (Fig. 2C) and a net increase in the expression level of corresponding genes (Fig. 2D). This class of genes included key TFs for endocrine lineage differentiation, such as Pax6, Arx, Fev and Neurod2 (Sander et al., 1997; St-Onge et al., 1997; Collombat et al., 2003; Gasa et al., 2008; Ohta et al., 2011) (Fig. 2E,F and Table S3). Taken together, these results indicate that the removal of bivalency-independent H3K27me3 is also associated with the regulation of endocrine cell specification.

Enhancers, which function as integrated TF-binding platforms, play central roles in controlling cell type- and cell state-specific gene expression patterns (Calo and Wysocka, 2013; Long et al., 2016). Active enhancers are characterized by the presence of both H3K4me1 and H3K27ac modifications, and are generally associated with active genes. However, primed enhancers are marked solely by H3K4me1 prior to the activation of gene transcription (Creyghton et al., 2010). To understand how enhancer-associated chromatin modifications change during the pancreatic endocrine cell fate transition, we mapped H3K4me1 and H3K27ac profiles in Ngn3-low and Ngn3-high cells from E13.5 dorsal pancreata (Fig. S2A and Table S1). We identified 9366 active and 30,600 primed putative enhancers in Ngn3-low cells, and 13,953 active and 54,514 primed putative enhancers in Ngn3-high cells (Fig. 3A). Genes associated with active enhancers had higher net expression levels than genes associated with primed enhancers in both Ngn3-low and Ngn3-high cell populations (Fig. 3B). Of the active enhancers in Ngn3-high cells, 20.7% were primed and 49.0% were unmarked in the Ngn3-low population (Fig. 3C). These 9718 newly formed active enhancers were associated with the upregulation of gene expression levels during the transition from Ngn3-low to Ngn3-high (Fig. 3D,E). Furthermore, MOTIF enrichment analysis revealed that the newly formed active enhancers in Ngn3-high cells were enriched for binding motifs, such as those for NEUROD1, FOXA2 and FOXO1 (Fig. 3F). KEGG pathway enrichment analysis showed that the genes associated with this class of enhancer were related to the regulation of pancreatic development and several cell signaling pathways (Fig. 3G and Table S4).

Of the 54,514 primed enhancers identified in Ngn3-high cells, 57.9% were unmarked in the Ngn3-low cells (Fig. 3C). The selection of new enhancer loci during cell fate transitions may be guided by ‘pioneer’ factors (Calo and Wysocka, 2013; Iwafuchi-Doi and Zaret, 2014), which are capable of directly binding to even compacted chromatin, thus making DNA-binding sites accessible to chromatin modifiers and other TFs. In Ngn3-high cells, the newly formed primed enhancers predominantly enriched the binding motif of NEUROD1 (Fig. 3H), which has been demonstrated to be a pioneer factor responsible for the establishment of a fraction of new enhancer sites during neuronal cell differentiation (Pataskar et al., 2016). Curiously, the newly primed enhancers in Ngn3-high cells were associated with genes involved in processes specific to the islet lineage, such as pancreatic secretion and metabolic pathways (Fig. 3I and Table S4). To investigate whether these newly formed primed enhancers in Ngn3-high cells contribute to the proportion of active enhancers in β cells, we performed H3K4me1 and H3K27ac ChIP-seq analyses on E17.5 β cells, and detected 5031 active enhancers newly formed from primed enhancers during the transition from Ngn3-high cells to β cells (Fig. S3A). Of this class of active enhancers, 63.5% originated from the newly primed enhancers in Ngn3-high cells (Fig. S3A). The genes associated with this enhancer state change showed an increase in net expression level during the transition from Ngn3-high cells to β cells (Fig. S3B,C).

To verify whether the occupancy of NeuroD1 was accompanied by the establishment of primed and active enhancers during the transition from Ngn3-low to Ngn3-high status, we performed NeuroD1 ChIP-seq in cells from E14.5 whole pancreas (Fig. S4A). We identified 1609 peaks, representing the NeuroD1-binding sites in endocrine lineages expressing Neurod1. Notably, 65.4% of the 1609 peaks contained the NEUROD1 motif (Fig. 4A), which indicated that the peaks identified in this study were reliable. After normalization for the small size of gene regulatory regions within the entire genome, it was clear that NeuroD1 binding-site peaks were preferentially enriched in enhancer and promoter regions (Fig. 4B,C). Further comprehensive analysis revealed that NeuroD1 bound to 569 Ngn3-high active putative enhancers, among which 79.4% were newly formed enhancers, and associated genes showed an increased net expression level during the transition from Ngn3-low to Ngn3-high (Fig. 4D-G). We also found that 274 Ngn3-high primed putative enhancers were occupied by NeuroD1, and 77.4% of those were newly formed primed enhancers (Fig. 4D). Gene ontology (GO) analysis revealed that the genes associated with these newly formed active and primed enhancers were involved in developmental processes and pancreatic endocrine functions (Fig. 4H and Table S5). Therefore, these data indicate that the pancreatic endocrine-specific TF NeuroD1 may function as a pioneer factor accompanied by the establishment of a subset of active or primed enhancers in Ngn3-high cells. However, NeuroD1 may not be the only pioneer factor during the transition, given that many other newly formed enhancers contain other motifs, such as FOXA2 and FOXO1, which suggests that the generation of Ngn3-high cells probably requires multiple pioneer factors (Fig. 3F,H). In addition to the enhancer loci, we identified NeuroD1-binding peaks at 67 promoter regions (Fig. S4B). The related genes, including Hes6 and Neurod1 (which encode TFs) showed increased net expression levels during the transition from Ngn3-low to Ngn3-high (Fig. 4C and Fig. S4C).

Taken together, these data characterize the dramatic dynamics of primed and active enhancers, which are accompanied by the transition from Ngn3-low to Ngn3-high. Our study also indicates that the activation of Neurod1 after the resolution of H3K27me3 from its bivalent promoter is accompanied by the establishment of active and primed enhancers in Ngn3-high cells. Given that the removal of H3K27me3 was associated with the activation of other key TFs related to endocrine development, in addition to the activation of NeuroD1 (Fig. 1E,F), we hypothesized that a disturbance in the removal of H3K27me3 would perturb the cell fate transition from Ngn3-low to Ngn3-high.

Only two histone demethylases, JMJD3 and UTX, have been proved to be responsible for the removal of H3K27me3 mark (Agger et al., 2007; Lan et al., 2007). Studies have indicated that JMJD3 regulates promoter-specific H3K27me3 levels in human ES cells (Akiyama et al., 2016). Hence, we examined whether the precocious loss of Jmjd3 expression would affect cell fate commitment from Ngn3-low to Ngn3-high. We used a Jmjd3 conditional allele (Fig. S5A) and a Pdx1-Cre transgene to delete Jmjd3 from pancreatic progenitor cells on an Ngn3-GFP background (Hingorani et al., 2003). The expression level of Jmjd3 was reduced in sorted Ngn3-low and Ngn3-high cells from E13.5 Pdx1-Cre;Jmjd3^fl/fl;Ngn3-GFP dorsal pancreata (Fig. S5B). Flow cytometric analysis also showed that in E13.5 Pdx1-Cre;Jmjd3^fl/fl pancreata, the number of Ngn3-low cells increased 1.8-fold, whereas the number of Ngn3-high cells decreased 1.5-fold, with the ratio of Ngn3-high to Ngn3-low cells decreasing 2.6-fold compared with the ratio in Pdx1-Cre;Jmjd3^+/+ pancreata (Fig. 5A,B). Later, at E14.5, similar analysis also revealed that the ratio of Ngn3-high to Ngn3-low cells in knock-out (KO) pancreata decreased 1.7-fold compared with the ratio in wild-type pancreata (Fig. S5C). To determine whether cell death occurred in Jmjd3 KO cells, we performed TUNEL detection of E14.5 FACS-sorted Ngn3-low and Ngn3-high cells from wild-type or KO pancreata. We observed few apoptotic cells in Jmjd3 wild-type and KO Ngn3-low and Ngn3-high cells (Fig. S5D). Therefore, the decreased ratio of Ngn3-high to Ngn3-low cells in Jmjd3 KO pancreata was not induced by cell death. These findings indicate that the cell fate transition from Ngn3-low to Ngn3-high is constrained in Jmjd3-ablated mice.

We then examined whether the deletion of Jmjd3 by Pdx1-Cre in pancreatic progenitor cells affected the properties of individual Ngn3-low and Ngn3-high cells. We performed single-cell transcriptomic analyses of 49 Ngn3-low and 59 Ngn3-high cells purified from wild-type pancreata and of 69 Ngn3-low and 44 Ngn3-high cells from Pdx1-Cre;Jmjd3^fl/fl pancreata at E13.5 (Fig. S6A-D). Principal component analysis (PCA) (Pearson, 1901; Hotelling, 1933) revealed a developmental pathway from Ngn3-low to Ngn3-high cells (Fig. 5C). However, the sources of Ngn3-low and Ngn3-high cells, i.e. wild-type or Pdx1-Cre;Jmjd3^fl/fl pancreata, could not be distinguished based on the PCA plot (Fig. 5C). Hierarchical clustering analyses on Ngn3-low and Ngn3-high populations showed that the single cells from wild-type and Jmjd3 KO pancreata were intermingled in the heat maps (Fig. S5F). Furthermore, a comparative transcriptomic analysis barely identified differentially expressed genes between wild-type and Jmjd3-deleted Ngn3-low or Ngn3-high cells (Fig. 5D). This result indicates that Jmjd3 KO decreased the rate of transition from Ngn3-low to Ngn3-high but did not affect the cell properties before or after cell fate transition.

Considering the dramatic change in the ratio of Ngn3-low to Ngn3-high cells, we next assessed the associated developmental consequence by quantifying the number of nascent β cells from E17.5 Pdx1-Cre;Jmjd3^fl/fl pancreata on an Ins1-RFP transgenic background. Using a different genetic background can eliminate the effects of the Ngn3-GFP knock-in mouse line, which has been shown to have an increased capacity for exocrine differentiation (Wang et al., 2010). We observed a 62.3% decrease in the proportion of β cells in Pdx1-Cre;Jmjd3^fl/fl dorsal pancreata compared with that in wild-type mice (Fig. 6A,B). In addition, pancreatic islet areas were reduced and the shape of islet was deficient in 6-week-old Pdx1-Cre;Jmjd3^fl/fl male mice (Fig. 6C). In Jmjd3 KO islets, the level of H3K27me3 modifications were higher than those in wild-type islets (Fig. S5G). Moreover, Pdx1-Cre;Jmjd3^fl/fl animals were smaller than wild-type mice (Fig. 6D and Fig. S5H) and had diabetes, demonstrating defects in glucose metabolism and serum insulin secretion after overnight fasting and glucose stimulation (Fig. 6E,F). Our data demonstrate that Jmjd3 promotes the cell fate transition from Ngn3-low to Ngn3-high cells, whereas the deletion of Jmjd3 by Pdx1-Cre at an earlier stage impairs islet development.

In this study, we compared the transcriptomic profiles of cell populations at the key stages of pancreatic endocrine cell development. To our surprise, the dramatic changes in gene expression occurred during the transition from Ngn3-low to Ngn3-high, which was completed within ∼12 h (Bankaitis et al., 2015). During this process, the cells release the suppression of key TFs, mainly by removing the H3K27me3 mark. Some of these TFs might in turn help establish the primed and active enhancers, which facilitate the fate commitment to Ngn3-high cells (Fig. 7). Therefore, our studies suggest the existence of sequential chromatin modification events during lineage specification.

Cell fate decisions are usually linked to changes in cell cycle progression, and loss of pluripotency is associated with an extension of the G1 phase of the cell cycle (Dalton, 2015). However, whereas pancreatic endocrine progenitor cells lose multipotency and become specified in the islet lineage, cell cycle regulatory genes are broadly repressed, culminating in G1 arrest in Ngn3-high cells. During the G1 cell phase, the key TF genes for endocrine differentiation, such as Neurod1, Pax4, Insm1 and Rfx6, start to be expressed. Considering the putative role of NeuroD1 in establishing enhancer-related chromatin marks, we presume that a significant number of dynamic changes in H3K4me1 and H3K27ac occur during the G1 phase. However, based on current findings, we cannot speculate whether the resolution of H3K27me3 and H3K4me3 is related to the cell cycle phase. Therefore, this system is suitable for investigating the mechanisms related to the regulation of the cell cycle and chromatin modifications during the loss of multipotency and lineage specification of progenitor cells.

The histone H3K27me3 balance is regulated by the histone methyltransferases EZH2 and EZH1, and their histone demethylase counterparts JMJD3 and UTX (Voigt et al., 2013). Our previous work demonstrated that EZH2 restricts the specification of pancreatic progenitors to an Ngn3-low population (Xu et al., 2014), but the current study, using the same Pdx1-Cre system, revealed that Jmjd3 promotes the transition from Ngn3-low to Ngn3-high. Therefore, rather than playing opposing roles during the same step of cell lineage commitment, Ezh2 and Jmjd3 function at different stages of endocrine cell specification.

In our study, we performed single-cell RNA transcriptomic analyses in Ngn3-low and Ngn3-high cells from Jmjd3 wild-type or KO embryos. Surprisingly, we found no obvious differences between wild-type and KO cells in terms of the transcriptome and developmental pathways (Fig. 5C,D). This result indicated that knocking out Jmjd3 blocked the transition from Ngn3-low to Ngn3-high. We discussed that in Ngn3-low cells, the histone modification H3K27me3 produces a chromatin barrier to pancreatic endocrine specification by repressing the expression of key TFs. During the transition from Ngn3-low to Ngn3-high, Jmjd3 reduced this chromatin barrier, thereby facilitating the cell fate transition. In Jmjd3 KO pancreas, the proportion of cells that could overcome this chromatin barrier was significantly decreased; however, once these cells overcame the barrier, they developed normally, which was demonstrated at the single-cell level. In addition, in E17.5 Jmjd3 KO pancreata, we observed differentiated β cells (Fig. 6A,B). This finding is consistent with our earlier studies of foregut endoderm differentiation, which indicated that the histone modifications, the histone modifier itself, or both play a modulatory role rather than an absolute governance role in cell fate induction (Xu et al., 2011). This finding also suggests that transiently regulating the activities of chromatin modifiers using small molecules at specific times in vitro might be a safe way to regulate cell fate transition.

The development of endocrine cells in vivo is controlled in a coordinated manner by cell signaling pathways, cell cycle regulators, key TFs and histone modifiers. However, how the cells integrate all of these different levels of regulation into the determination of cell fate at a precise developmental stage needs further investigation. Understanding the mechanisms underlying endocrine cell specification will undoubtedly facilitate the in vitro generation of functional β cells.

Pdx1-Cre (Hingorani et al., 2003), Ngn3-GFP (Lee et al., 2002), Ins1-RFP (Piccand et al., 2014) and Jmjd3^fl/fl mice were used. Jmjd3^fl/fl mice were generated at the Nanjing Biomedical Research Institute of Nanjing University. All protocols were approved by the Institutional Animal Care and Use Committees of Peking University. All mice were housed and maintained under specific pathogen-free conditions on a 12 h/12 h day/night cycle at 23±2°C and received an autoclaved standard diet and water ad libitum.

Dissected E13.5 and E14.5 dorsal pancreata were dissociated into single-cell suspensions by digestion with 0.25% trypsin for 5 min at 37°C. The digestion was terminated by adding 0.4 volumes of fetal bovine serum. E17.5 pancreata were digested with 0.5 mg/ml of collagenase P (Roche, 11213873001) for 2 min at 37°C, followed by trypsin treatment. Cells were sorted using a BD FACSAria II cell sorter.

Total RNA was isolated from ∼4-6×10⁴ sorted cells using an RNeasy Micro Kit (Qiagen, 74004). mRNA was purified from ∼200 ng of total RNA using an NEBNext Poly(A) mRNA Magnetic Isolation Module (NEB, E7490). Libraries were prepared using an NEBNext Ultra Directional RNA Library Prep Kit for Illumina (NEB, E7420), according to the manufacturer's instructions. The quality of the cDNA libraries was assessed using an Advanced Analytical Fragment Analyzer. RNA-seq analyses were performed on two independent biological replicates. Correlations between the replicates are shown in Fig. S1B.

Single cells were manually picked up from FACS-purified bulk cells by using a mouth pipette. cDNA synthesis was performed following the Smart-seq2 method described by Picelli et al. (2014). Briefly, 0.1 µl of ERCC spike-in RNA (1:1,000,000, Life Technologies, 4456740) was added to 4 µl of lysis buffer. The libraries were prepared with 2.5 ng of cDNA using a TruePrep DNA Library Prep Kit (Vazyme, TD502) (Qiu et al., 2017).

Low-cell-number ChIP-seq was performed following the method of Brind'Amour et al. (2015). We used ∼1.5×10⁴ cells per ChIP with primary antibodies (1 ng each) recognizing H3K4me3 (Abcam, ab8580, GK9535-2), H3K27me3 (Millipore, 07-449, 2318778), H3K4me1 (Millipore, 07-436, 2343581) and H3K27ac (Active motif, 39133, 01613007). NeuroD1 ChIP-seq was performed following the method of Johnson et al. (2007). We used ∼8×10⁶ cells per ChIP with an antibody against NeuroD1 (Santa Cruz Biotechnology, sc-1084, A0616), which has now been discontinued at Santa Cruz Biotechnology. The library was prepared using the Ultra DNA Library Prep Kit for Illumina (NEB, E7370). Correlations between the replicates are shown in Fig. S2A.

FACS-sorted E14.5 Ngn3-low and Ngn3-high cells were attached on poly-prep slides (Sigma, P0425-72EA), fixed with 4% paraformaldehyde and stained with a TUNEL apoptosis detection kit (Alexa Fluor 647) (Yesen, 40308), following the manufacturer's instructions. Pancreata from 6-week-old male mice were fixed with 4% paraformaldehyde and embedded in paraffin wax. Thereafter, 5 µm paraffin sections were stained with guinea pig anti-insulin (1:500, Abcam, ab7842, GR119778-1) and mouse anti-amylase (1:500, Santa Cruz Biotechnology, sc-46657, B0911) or goat anti-Pdx1 (1:2000, Abcam, ab47383, GR16180-1) and rabbit anti-H3K27me3 (1:800, Millipore, 07-449, 2318778) primary antibodies. The staining of Pdx1 and H3K27me3 was with antigen recovery using Vector Antigen Unmasking Solutions (Vector Labs, H-3300) at high pressure for 8 min. The secondary antibodies used were Alexa Fluor 488 donkey anti-mouse (Invitrogen, A21202), Alexa Fluor 594 goat anti-guinea pig (Invitrogen, A11076), Alexa Fluor 488 donkey anti-goat (Invitrogen, A11055) and Alexa Fluor 594 donkey anti-rabbit (Invitrogen, A21207). Fluorescent images were obtained using a Zeiss Axioimager M2 or a Zeiss Axio Scan Z1.

Total RNA was isolated from FACS-sorted cells using the RNeasy Micro Kit (Qiagen) and reverse-transcribed using HiScript II Q RT SuperMix (Vazyme, R223). The RT-qPCR primers are: Gapdh forward, 5′-ATGGTGAAGGTCGGTGTGAAC-3′; Gapdh reverse, 5′-GCCTTGACTGTGCCGTTGAAT-3′; Jmjd3 forward, 5′-CCTCACCGCCTATCAGTA-3′; Jmjd3 reverse, 5′-GCCATTCTCACTTGTAACG-3′.

Six-week-old male Jmjd3 KO mice and their littermates were fasted for 16 h and then intraperitoneally injected with glucose at 2 mg/g body weight. Blood glucose levels were monitored before injection and at 15, 30, 60 and 120 min after injection. For insulin secretion assays, blood was drawn both before and 30 min after glucose injection, and the serum insulin levels were assayed using a Rat/Mouse Insulin ELISA Kit (Millipore, EZRMI-13K), following the manufacturer's instructions. Jmjd3-deleted mice with or without the Ngn3-GFP background showed similar phenotypes; thus, we pooled their results for statistical analyses.

Adult islet cells were lysed in RIPA buffer with protease inhibitor cocktail (Sigma, P8340) and subjected to standard immunoblotting analysis. Goat anti-NeuroD1 (1:100, Santa Cruz Biotechnology, sc-1084, Lot# A0616) and donkey anti-goat IgG-HRP (1:500, Santa Cruz Biotechnology, sc-2020) antibodies were used.

Bulk-cell RNA-seq libraries were sequenced as 101 bp paired-end or 51 bp single-end clusters using an Illumina HiSeq 2000 System to a depth of 2-4×10⁷ reads. The sequencing reads were aligned to the Mus musculus genome (mm10) using TopHat (v2.1.0) (Kim et al., 2013) with the parameters ‘-o out_dir -G gtf --transcriptome-index trans_index bowtie2_index input_fastq1 input_fastq2’. Reads aligned to genes were counted using HTSeq (v0.6.0) (Anders et al., 2015) with the parameters ‘htseq-count -f bam -r pos -s reverse -a 30 input_bam gtf’. Gene expression levels were quantified as reads per kilobase of the longest transcript per million mapped reads (RPKM). The normalization of RPKM values was performed by size factor (Anders and Huber, 2010). Ward's hierarchical clustering was performed with Spearman's correlation distance metric (Fig. S1C). Genes shown in Fig. S1C were identified based on the following criteria: coefficient of variation >0.3 and expressed in at least two samples (normalized RPKM >5). Differential expression analysis was performed using the R package DESeq2 (v1.12.4) (Love et al., 2014) with adjusted P-value (padj) <0.05 and fold change >2. Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000) pathway and Gene Ontology (GO) (Ashburner et al., 2000) enrichment analyses were performed using the R package GOstats (v2.38.1) (Falcon and Gentleman, 2007).

ChIP-seq libraries were sequenced as 51 bp single-end reads on an Illumina HiSeq 2500 system to a depth of 1-4×10⁷ reads. The sequencing reads were aligned to the Mus musculus genome (mm10) using bowtie2 (v2.2.5) (Langmead and Salzberg, 2012) with default parameters, and only reads with a mapping quality ≥30 were retained for subsequent analyses. PCR duplicate reads were removed using samtools (v1.3.1) (Li et al., 2009) with the commands ‘samtools sort in_bam|samtools rmdup -s - out_bam’. Reads from biological replicates were pooled. TP10M (tags per 10 million) values were calculated and stored in the bedGraph format using Homer (Heinz et al., 2010) with the parameters ‘makeUCSCfile out_dir –o out.bdg –name sample_name -color track_color –fragLength 150 –avg -fsize 1e20’. To calculate signal intensity, we subtracted the input TP10M from the ChIP TP10M. Regions around the transcription start sites (TSSs) ±2 kb were considered as promoters. Promoter regions that overlapped with ENCODE blacklisted genomic regions (Encode Project Consortium, 2012) were excluded from subsequent analyses. The average signal intensity of H3K27me3 and H3K4me3 at the promoters was determined using the ScoreMatrixBin function of the R genomation package (v1.4.2) (Akalin et al., 2015). For genes with multiple transcript isoforms, the start sites of the isoform with the highest average H3K4me3 signal intensity was selected as the TSSs of these genes. Promoters for which the average H3K4me3 (or H3K27me3) signal intensity was greater than 1 (or 0.7) were considered to be marked by H3K4me3 (or H3K27me3), and promoters marked by both H3K4me3 and H3K27me3 were considered to be bivalent (Fig. 1A). A comparison of the averaged H3K27me3 signal intensities at the promoters is shown in Fig. 2A, in which genes with an H3K27me3 signal intensity >0.7 and a fold change >2 were identified.

Highly enriched regions (peaks) of H3K4me1 were identified using MACS (v1.4.2) (Zhang et al., 2008) with default settings. We merged H3K4me1 peaks that did not overlap with ENCODE blacklisted genomic regions and regions ±3 kb in length around known TSSs as potential enhancers. The average signal intensities of H3K4me1 and H3K27ac at the potential enhancers were calculated using the ScoreMatrixBin function of the R genomation package (v1.4.2) (Akalin et al., 2015). Potential enhancers with an H3K4me1 signal intensity >1 and an H3K27ac signal intensity ≤1 were defined as primed enhancers, whereas potential enhancers with an H3K4me1 signal intensity >1 and an H3K27ac signal intensity >1 were defined as active enhancers. The closest gene to a given enhancer was regarded as its putative target. Motif enrichment analysis was performed using the Homer (Heinz et al., 2010) ‘findMotifsGenome.pl’ script with default parameters.

NeuroD1-binding sites were identified using MACS (v1.4.2) (Zhang et al., 2008) with parameters ‘-f BAM -g mm -n sample_name -t ip_bam -c input_bam -p 1e-5’. The regions of enhancers (center ±2 kb) overlapping with NeuroD1 peaks were considered as NeuroD1-bound enhancers. Promoter regions (−800 to +200 bp flanking the TSSs) overlapping with NeuroD1 peaks were considered as NeuroD1-targeted promoters. The percentage of NEUROD1 motifs occurring at NeuroD1 peaks was calculated using Homer (Heinz et al., 2010) with commands ‘findMotifsGenome.pl Neruod1_peaks.bed mm10 motif_dir -size given’.

Single-cell RNA-seq libraries were sequenced as 51 bp single-end reads on an Illumina HiSeq 2500 system to a depth of 0.5-2×10⁶ reads. The sequencing reads were aligned to the Mus musculus genome (mm10) using TopHat (v2.1.0) (Kim et al., 2013) with the parameters ‘-o out_dir -G gtf --transcriptome-index trans_index bowtie2_index input_fastq’. Reads aligned to genes were counted using HTSeq (v0.6.0) (Anders et al., 2015) with the parameters ‘htseq-count -f bam -r pos -s no -a 30 input_bam gtf’. Gene expression levels were quantified as RPKM, which were correlated with the number of mRNA molecules (Fig. S6A). Saturation analysis show that 0.5 million mapped reads were sufficient to detect most genes in the single-cell libraries (Fig. S6B). Therefore, the cells having more than 0.5 million mapped reads were retained for further analyses (Fig. S6C and Table S1). Cells expressing the leukocyte marker gene Spi1 (PU.1) (Hromas et al., 1993) were also excluded from subsequent analyses. In general, 6000-9000 genes were detected in each cell (Fig. S6D). We accounted for technical noise and identified highly variable genes based on the ERCC spike-in as described previously (Brennecke et al., 2013). PCA was performed using the ‘PCA’ function of the FactoMineR package (v1.34) (Le et al., 2008) with the log2(RPKM+0.1) values of highly variable genes. The graphic displays for the PCA results were generated using the R package ggplot2 (v2.2.1) (Hadley, 2009). Ward's hierarchical clustering was performed with Spearman's correlation distance metric on highly variable genes (Fig. S5F).

The standard error of the mean (s.e.m.) and P-values of unpaired two tails t-tests were calculated for analyses of FACS, blood glucose levels, insulin secretion levels, imaging, RT-qPCR and body weights. The numbers of independent experiments/biological replicates are indicated in the figures. P<0.05 was considered significant. Gene expression data are shown in box plots; unpaired or paired two-tailed Wilcoxon rank-sum tests were performed, and the P-values are presented above the box plots. The chi-square test was performed for Fig. 4B, and the P-value is provided in the figure legend.

We thank Dr Gérard Gradwohl for the Ins1-RFP transgenic mouse; Dr Fuchou Tang for single-cell RNA-seq advice; Drs Erfei Bi, Ken Zaret, Fuchou Tang and Wei Xie for comments on the manuscript; Ms Si Si for technical assistance; members of the Xu laboratory for advice and comments; Mr Guopeng Wang from the Core Facility of the School of Life Sciences, Peking University, for fluorescent image acquisition with the Zeiss Axio Scan Z1; and the Peking-Tsinghua Center for Life Sciences Computing Platform.