Comprehensive single cell mRNA profiling reveals a detailed roadmap for pancreatic endocrinogenesis

In rodents, pancreatic endocrine cells are generated at two distinct stages (transitions) during embryonic development. The first transition [embryonic day (E) 9.0-12.5] produces mainly glucagon-producing α-cells, whereas all other endocrine cell types (insulin-producing β-cells, somatostatin-producing δ-cells, pancreatic polypeptide-producing PP-cells and ghrelin-producing ε-cells) are mainly generated during the secondary transition (E12.5-15.5) and thereafter. At both stages, endocrine cells are derived from the endocrine progenitor-precursor (EP) pool (Bastidas-Ponce et al., 2017a). EPs are marked by the transient expression of the transcription factor (TF) neurogenin 3 (Neurog3, hereafter called Ngn3) and are derived from Sox9⁺ bipotent ductal/endocrine progenitors, which are located within the pancreatic epithelium (Gradwohl et al., 2000; Gu et al., 2002; Kopp et al., 2011; Seymour et al., 2007; Shih et al., 2012). Upon lineage priming and specification, EPs express low levels of Ngn3 and by receiving additional signals these mitotic and non-committed Ngn3^low cells (which we refer to as endocrine progenitors) increase the levels of Ngn3 to become determined Ngn3^high EPs (to which we refer to as endocrine precursors) (Bechard et al., 2016; Pan and Wright, 2011; Wang et al., 2010). Further differentiation of Ngn3^high precursors generates distinct endocrine subtypes of the aforementioned hormone-expressing lineages. Changes in the gene regulatory networks, extracellular signaling cues from the respective niche and epithelial states (non-polarized versus polarized and epithelialized) cooperatively orchestrate the specification and determination of EPs (Apelqvist, 1999; Arda et al., 2013; Bankaitis et al., 2015, 2018; Cortijo et al., 2012; Löf-Öhlin et al., 2017). Yet, during progressive lineage restriction from specified EP progenitors to determined EP precursors the molecular signature changes over time and the mechanisms underpinning lineage allocation towards a specific endocrine subtype fate have remained largely elusive. Although several recent studies (Byrnes et al., 2018; Krentz et al., 2018; Ramond et al., 2018; Scavuzzo et al., 2018; Sharon et al., 2019; Stanescu et al., 2016; Yu et al., 2018, 2019) have shed light on the transcriptional profiles of EPs, the importance of these cells for endocrine cell formation still necessitates a comprehensive temporally resolved and fine-grained mapping of gene expression on the single cell level. The low percentage of EPs within the embryonic pancreas together with the limitation of transgenic reporter lines to reflect accurately the transient expression of the Ngn3 protein have been major obstacles to establishing a detailed roadmap of endocrinogenesis.

Here, we generated a novel reporter mouse line, in which the endogenous Ngn3 is fused to a very bright Venus fluorescence reporter protein (Ngn3-Venus fusion; NVF). The NVF reporter mirrors transient endogenous Ngn3 protein levels and allows specific isolation of EPs in the narrow time window when they actually express Ngn3. Single cell RNA sequencing (scRNA-seq) analysis of EP-enriched pancreatic epithelial cells highlighted the step-wise dynamic changes in gene expression programs from bipotent cells to EP progenitors and precursors, and finally to differentiated endocrine cells as well as ductal and acinar cells. We found a set of 58 previously undescribed EP-signature genes that are specifically expressed in EPs in a similar transient manner as Ngn3. Furthermore, we found that several of these EP-signature genes are differentially expressed in endocrine precursors derived from different developmental stages that defines their lineage restriction towards specific endocrine subtypes. Altogether, our data provide a comprehensive single cell survey of stage-dependent and lineage-specific gene regulation during EP induction, specification and segregation.

We generated an NVF reporter mouse line, in which the Ngn3 translational stop codon was removed and Venus was fused in-frame with Ngn3 using a previously described strategy (Petrezselyova et al., 2015). Heterozygous intercrosses produced offspring with normal Mendelian distribution and homozygous Ngn3^Venus/Venus animals were fertile and indistinguishable from wild-type (WT) littermates (data not shown). To assess whether the NVF reporter reflects the endogenous Ngn3 expression pattern during development, we performed immunostaining of embryonic pancreata from heterozygous Ngn3^Venus/+ mice at E12.5-18.5. Using antibodies against Ngn3 and Venus, we detected a similar nuclear expression pattern of NVF and endogenous Ngn3 in pancreatic sections from different developmental stages (Fig. 1A; Fig. S1A). Plot profile analysis of staining by Venus and two different Ngn3 antibodies also revealed a similar pattern of nuclear signals between NVF and Ngn3 in pancreatic sections (Fig. S1B). Furthermore, using markers of the pancreatic tubular epithelium, such as Sox9, Nkx6-1, Mucin-1 (Muc1), EpCAM and E-cadherin (cadherin 1), we found that NVF⁺ cells exist as both intra- and extra-epithelial cells (Fig. 1B). This data is in line with the idea that at the Ngn3-expressing state, EP progenitors get specified in the pancreatic epithelium and endocrine precursor delaminate into the surrounding mesenchyme to form proto-islet clusters (Gouzi et al., 2011; Sharon et al., 2019). During development, Ngn3 is expressed at low levels in progenitors, reach peak levels in precursors and the levels are reduced upon differentiation of the precursors into the endocrine hormone⁺ lineage (Bechard et al., 2016). To find out whether the NVF protein reflects the transient expression pattern of endogenous Ngn3 protein, we stained E18.5 pancreatic sections and quantified the number of differentiated hormone⁺ endocrine cells that still expressed Ngn3 protein. We observed a subtle and weak NVF signal in the nucleus or cytoplasm of a minor fraction of the hormone⁺ endocrine cells (less than 0.6%) (Fig. 1C,D), confirming that the NVF reporter protein shows a similar transient expression pattern to that of the endogenous Ngn3 protein.

It has been shown that haploinsufficiency of Ngn3 alters endocrine cell development and post-translational modifications of this protein affect its activity (Azzarelli et al., 2017; Krentz et al., 2017; Wang et al., 2010). To confirm that the fusion protein truly represents the activity of endogenous Ngn3 protein, we first checked the expression of several downstream target genes, including Neurod1, Pax4 and Arx. Quantitative PCR (qPCR) analysis revealed no significant differences in the expression levels of these genes in isolated epithelial cells from WT and homozygous NVF pancreata at E15.5 (Fig. S1C). Next, we analyzed endocrine cell composition in WT and NVF reporter embryos. E15.5 and E18.5 pancreata from WT and homozygous NVF embryos showed comparable α-, β-, δ- and PP-cell composition (Fig. 1E,F; Fig. S1D). Furthermore, qPCR data demonstrated comparable expression levels of key endocrine cell markers in isolated islets from adult WT and homozygous NVF mice (Fig. S1E). These data were supported by the detection of normal fasting blood glucose levels in adult homozygous Ngn3^Venus/Venus animals (Fig. S1F) and by comparable levels of glucose-stimulated insulin secretion of isolated WT and Ngn3^Venus/Venus islets (Fig. S1G). Altogether, our analyses indicate that the NVF reporter is transiently expressed in EPs and reflects a spatiotemporal expression pattern mirroring the endogenous Ngn3 protein. Furthermore, the fusion of Venus to Ngn3 does not have a significant impact on TF function during endocrine cell development.

To time-resolve at the single cell level the molecular changes during endocrinogenesis, we performed high-throughput scRNA-seq analysis of pancreatic epithelial cells during the secondary transition (E12.5, E13.5, E14.5 and E15.5). To obtain the rare cell population enriched for EPs (expressing NVF), we first flow-sorted embryonic pancreatic cells using the fluorescence property of the Venus protein in homozygous Ngn3^Venus/Venus pancreata. As we were specifically interested into how epithelial multipotent progenitors become lineage restricted into acinar and bipotent (ductal/endocrine) as well as unipotent progenitors, we isolated from the remaining cells that were negative for NVF (Venus⁻) the pancreatic epithelial fraction by using the epithelial marker EpCAM. This allowed us to deplete non-epithelial cells and highly enrich for epithelial progenitor, ductal and acinar cells. A fraction of these EpCAM⁺/Venus⁻ (NVF⁻) cells were added to the corresponding EP cells (NVF⁺) at each stage (2:1 ratio for NVF⁺:NVF⁻) (Fig. 2A). We used droplet-based scRNA-seq and transcriptionally profiled 10,790 single cells at E12.5, 5042 cells at E13.5, 9633 cells at E14.5 and 10,886 cells at E15.5 (in total 36,351 cells) (Fig. 2B). Of the total cells derived from four different stages, 7260 cells (∼20%) expressed Ngn3. Unsupervised graph-based clustering revealed eight major cell clusters including multipotent pancreatic progenitors (MPCs), tip, trunk, acinar, ductal, EPs, Fev^high and endocrine cells (Fig. 2C; Table S1). These clusters represent distinct pancreatic lineages at different embryonic stages, when the pancreatic epithelium is composed of MPCs. The MPCs undergo step-wise lineage restrictions and segregation into two distinct tip and trunk domains, which further differentiate into acinar and duct/EPs, respectively (Nyeng et al., 2019; Zhou et al., 2007). EPs then differentiate into Fev^high cells that subsequently generate different types of endocrine hormone⁺ cells (Fig. 1D). Fev^high cells have been recently reported as a novel cluster of pancreatic cells (Byrnes et al., 2018). Clusters were annotated based on the expression of well-known marker genes with a described function during pancreas development and/or function, including Dlk1 (MPCs) (Ramond et al., 2018); Cpa1 and Myc (tip) (Zhou et al., 2007); Notch2 (trunk) (Lee et al., 2005); Ptf1a, Cpa1, Cel, Rbpjl (acinar), Sox9, Anxa2 and Bicc1 (ductal) (Lemaire et al., 2016); Ngn3 and Hes6 (EPs) (Ahnfelt-Rønne et al., 2007); Fev, Cck and Neurod1 (Fev^high) (Byrnes et al., 2018); and Rbp4, Pyy and Chgb (endocrine) (Fig. 2E,F). As the expression of Ngn3 is a well-described characteristic of EPs, we annotated Ngn3-expressing cells (except those expressing Fev) in the epithelial populations as EPs (see Materials and Methods). Our data revealed the expression of novel marker genes in each cluster that are uncharacterized or have not been associated with pancreatic cell development and/or function yet, e.g. Mdk and Btf3 (MPCs); Vtn and Jam3 (tip); Cbx3, Hmgn1 and Ybx1 (trunk); Reep5 (acinar); Spp1 (duct); Btbd17 and Gadd45a (EPs); Vwa5b2 and Tox3 (Fev^high); Tmem97 and Fam183b (endocrine) (Fig. 2E,F). Thus, the identification of known markers for different clusters validated our scRNA-seq approach. Furthermore, we identified many novel and functionally uncharacterized genes at different steps of pancreatic and endocrine lineage allocation.

To define the step-wise lineage restriction and allocation and to explore the dynamic heterogeneity of pancreatic progenitors, we performed stage-wise comparison of the molecular signatures of multipotent and bipotent progenitors as well as ductal, endocrine and acinar lineages. We found two subclusters of proliferating and non-proliferating cells for tip, trunk, acinar and ductal populations (Fig. 3A,B; Fig. S2A). MPCs were mainly evident at E12.5 and expressed the multipotent progenitor and acinar marker Cpa1 at high levels, whereas the tip and trunk clusters were prominently found at E12.5 and E13.5. Ductal cells were present throughout all stages and acinar cells appeared mainly at E14.5 and E15.5 (Fig. 3C). EP cells showed different expression levels of Ngn3 as has been previously described (Bechard et al., 2016) and were separated into Ngn3^low progenitors and Ngn3^high precursors (see Materials and Methods).

When looking for the origin of Ngn3^low progenitors, we found that they were scattered within the MPCs, tip, trunk and ductal progenitor clusters, suggesting multiple niches and birthplaces for these progenitors (Fig. 3D,E,G). The contribution of distinct epithelial clusters to Ngn3^low progenitor formation varied at different stages. Ngn3^low progenitors were found within the MPCs only at E12.5. The tip cluster also contained these progenitors from E12.5 to E14.5, but with rapidly decreasing numbers. The Ngn3^low progenitors were detected in the trunk domain mainly until E14.5. The presence of Ngn3^low progenitors from ductal cells was found at all stages and continuously increased from E12.5 to E15.5 (Fig. 3E). Thus, the MPCs and tip cells contain endocrine progenitors until E13.5, the existence of these cells within ductal cells persists until the end of the secondary transition. Moreover, the majority of Ngn3^low progenitors were derived from the ductal population (Fig. 3D,E), supporting the previous notion that this domain is the major epithelial source of endocrine progenitor-generating cells (Solar et al., 2009).

To infer the relationships between early pancreatic lineages, we performed partition-based graph abstraction (PAGA) analysis (Wolf et al., 2019) (see Materials and Methods). PAGA measures subgroup connectivity and draws a graph of possible differentiation paths. Consistent with the common belief (Sznurkowska et al., 2018; Zhou et al., 2007), we found that MPCs are mainly specified into tip and trunk domains. Tip cells then further connected to acinar cells, whereas trunk cells showed strong links to ductal cells. We found strong connections of Ngn3^low progenitors to the tip, trunk and ductal cells, and a weaker connection of these progenitors with MPCs (Fig. 3F). This further supports the finding that all epithelial clusters contribute to Ngn3^low progenitor formation and that a lower number of Ngn3^low progenitors are derived from MPCs compared with other epithelial clusters. Taken together, these findings indicate that lineage restriction from multipotent to bipotent to unipotent fates occurs in a rapid step-wise manner during the secondary transition.

To obtain a dynamic and comprehensive picture of the molecular changes of the EP transcriptome during endocrinogenesis, we compared all Ngn3⁺ EPs (Ngn3^low and Ngn3^high) with all other (Ngn3⁻) cells except acinar cells (Fig. 4A; see Materials and Methods). Differential expression analysis revealed several genes significantly regulated in the Ngn3⁺ EPs (Fig. 4B,C; Table S2) and provided a list of EP-enriched genes. We found increased expression of genes with known function in EPs, such as the endocrine-specification factor Amotl2 (Scavuzzo et al., 2018), the cell-cycle inhibitor Cdkn1a (Miyatsuka et al., 2011), the EP differentiation factor Sox4 (Xu et al., 2015) and the Notch inhibitor Mfng (Svensson et al., 2009) among many others (Fig. 4B,C). In addition, we identified several other upregulated genes that have not been associated with pancreas development before, including the cell-cycle inhibitors Gadd45a and Btg2 (Krentz et al., 2018), the protein phosphatase Ppp1r14a (also known as Cpi17), mitochondrial carrier homolog 1 (Mtch1), Numb-like protein (Numbl), and the Hes1 repressor Hes6 (Fig. 4B,C). Among those, Numbl and Hes6 inhibit Notch signaling highlighting the requirement of efficient suppression of this pathway by multiple mechanisms during endocrine induction.

To identify EP-specific signaling pathways and cellular processes, we analyzed the genes upregulated in Ngn3⁺ cells and identified Notch and Wnt inhibitors, further suggesting that these pathways are negatively modulated in EPs (Fig. 4D). Interestingly, we found increased expression of the Hippo pathway components Amotl2 and Nkd1, which are also antagonists of the Wnt pathway. This highlights that crosstalk between these pathways might be involved in regulation of endocrine cell formation. Furthermore, we found an upregulation of cell-cycle inhibitors in EPs (Fig. 4D), supporting the idea that endocrine differentiation leads to lengthening of the cell cycle. We identified several transcriptional repressors or co-repressors in the EP population, such as the co-repressor Cbfa2t3, the histone-modifying factor Rcor2, the chromatin remolding factor Smarcd2, the endocrine-specification factor Insm1 (Liu et al., 2006; Osipovich et al., 2014) and the Notch repressor Cbfa2t2 (Fig. 4D), that might be important for the suppression of non-EP gene programs.

Next, we investigated whether Ngn3^low progenitors derived from different epithelial clusters or from different developmental stages (E12.5-15.5) are transcriptionally heterogeneous. To this end, using the list of genes enriched in the total EP population (Ngn3^low and Ngn3^high) (Table S2), we first identified those genes that were highly expressed in Ngn3^low progenitors. Then, we compared the expression of these genes in Ngn3^low cells derived from MPCs, tip, trunk and ductal cells (Fig. 4E), as well as from different developmental stages (Fig. 4F). Some genes, such as Spp1, Tmsb10, Mdk, Marcksl1, Cdk4 and Sox4, were differentially expressed in Ngn3^low progenitors; however, expression of these genes only reflect the origin of the Ngn3^low cells. For instance, Ngn3^low progenitors derived from trunk and ductal epithelium expressed high levels of Spp1, which is highly expressed in all trunk and ductal cells. In contrast, the expression of most of the highly enriched genes in Ngn3^low progenitors that could not be directly associated to their origin was not different between Ngn3^low progenitor populations derived from distinct epithelial clusters or from different stages (Fig. 4E,F).

During development, bipotent epithelial cells generate Ngn3^low mitotic progenitors that either reverse to a ductal fate or further generate Ngn3^high precursors (Bechard et al., 2016). Ngn3^high cells then differentiate into Fev^high cells that consequently produce hormone-expressing endocrine cells (Fig. 5A). Consistently, pseudotemporal ordering (Haghverdi et al., 2016) of the cells along this route, which was also inferred by PAGA (Fig. 3F), revealed a continuous differentiation trajectory (see Materials and Methods). Along this pseudotime Ngn3 is transiently expressed as the levels of Ngn3 increased from Ngn3^low to Ngn3^high cells and then decreased again in Fev^high cells (Fig. 5B,C). Next, we investigated whether there are other EP-specific signature genes that follow transient Ngn3 expression in EPs. From the list of genes enriched in EPs we identified 58 genes (including Ngn3) that were expressed in EPs, but not or only at very low level in other pancreatic lineages (see Materials and Methods). Moreover, their expression was transient in EPs, i.e. their expression levels were gradually reduced upon differentiation of EPs into Fev^high cells (Fig. 5B; Fig. S2B). Among these genes were Ppp1r14a, the neuronal determination TF neurogenic differentiation factor 2 (Neurod2) (Gasa et al., 2008), sulfotransferase family cytosolic 2B member 1 (Sult2b1), uroplakin 3b-like protein 1 (Upk3bl), Ly6/neurotoxin-like protein 1 (Lynx1), the cell cycle regulator polo-like kinase 3 (Plk3), G protein γ subunit Gγ13 (Gng13), the Six1 transcriptional co-activator Eya2, regulator of G-protein signaling 16 (Rgs16), semaphorin 3G (Sema3g), the cancer-associated cell-surface antigen six-transmembrane epithelial antigen of the prostate 1 (Steap1), and the uncharacterized gene Gm8773 (Fig. 5C; Fig. S2B). Notably, for most of these genes no function associated with pancreas development has been reported so far. However, genes such as Neurod2, Sult2b1 and Lynx1 are involved in neuronal development and function, which highlights the similarities of the developmental programs required for the formation of EPs and neuronal cells. We also found comparable expression levels of several of these signature genes in E15.5 pancreatic epithelial cells from WT and homozygous NVF pancreata, which further supports proper functioning of the fusion protein (Fig. S2C). Furthermore, scRNA-seq analysis of human fetal pancreatic cells has indicated the expression of some of these genes, including EYA2 and RGS16 in human endocrine progenitors (Ramond et al., 2018).

Next, from the list of genes enriched in the Fev^high population we extracted several genes that followed a similar expression pattern to that of Fev. The expression of these genes started in Ngn3^high precursors, reached maximum levels in Fev^high cells and then decreased again in hormone-expressing endocrine cells. Such Fev^high-enriched genes included Neurod1, von Willebrand factor A domain-containing protein 5B2 (Vwa5b2), the paired box protein Pax4, cholecystokinin (Cck), the cancer suppressor gene Tox3, solute carrier family 35 member D3 (Slc35d3), the Ngn3-regulated transcriptional co-repressor Runx1t1 (Cbfa2t1) (Benitez et al., 2014) and the diabetes-associated protein voltage-dependent calcium channel subunit α₂δ-1 (Cacna2d1) (Mastrolia et al., 2017) and the uncharacterized gene transmembrane protein 185A (BC023829) (Fig. 5D,E; Fig. S2D). Among these, a high expression of RUNX1T1 in FEV⁺ cells has also been reported in human fetal pancreas (Ramond et al., 2018). Moreover, the expression of genes, such as chromogranin-A (Chga) and -B (Chgb), carboxypeptidase E (Cpe), Pax6, the ciliary-localized protein Fam183b (Stauber et al., 2017) and the neuroendocrine protein ProSAAS endopeptidase inhibitor (Pcsk1n) was turned on at the Fev^high stage and persisted in fully differentiated endocrine cells (Fig. S2D,E). Altogether, these data suggest that the development from EPs to fully differentiated endocrine cells occurs in a step-wise manner and involves the expression of unique molecular signatures at each stage (Fig. 5E).

To resolve further the dynamic changes in gene expression during the step-wise endocrine induction, specification and differentiation, we compared the transcriptional profiles of ductal bipotent (as the main source of Ngn3^low progenitors), Ngn3^low, Ngn3^high and Fev^high populations (Fig. S3A; Table S3). We identified a large number of differentially expressed genes between ductal bipotent and Ngn3^low, Ngn3^low and Ngn3^high as well as Ngn3^high and Fev^high cells (Fig. 5F). Ordering the cells along pseudotime supported a continuous developmental trajectory and revealed stage-specific induction or repression of genes involved in signaling pathways, cell-cycle control and organ morphogenesis during EP formation and differentiation. The low expression of Ngn3 in EPs was accompanied by increased expression levels of genes involved in activation of cAMP, Rap1 and Ras/MAPK signaling as well as Notch inhibitors (Fig. 5F; Fig. S3B). Furthermore, compared with ductal bipotent cells, Ngn3^low progenitors expressed decreased levels of genes involved in Wnt signaling activation, suggesting the requirement of Wnt signaling inhibition during EP development. This was supported by upregulation of the secreted TGFβ and Wnt inhibitor Cer1 (Piccolo et al., 1999) and Wnt-deacylation enzyme Notum (non-canonical inhibitor) (Kakugawa et al., 2015) in Ngn3^low progenitors compared with ductal bipotent cells (Fig. 5F). Additionally, alterations in cell adhesion molecules and activation of the Ras/Rap1 pathway in Ngn3^low progenitors likely drive cell morphological changes, suggesting that programs of morphogenesis are initiated at the progenitor level and continue while EPs delaminate from the pancreatic epithelium and differentiate into distinct endocrine cell types. Moreover, the expression of genes involved in cell-cycle exit was initiated in Ngn3^low progenitors and continued in Ngn3^high precursors. The high expression levels of Ngn3 in EPs initiated the expression of the secretory machinery during differentiation towards slowly cycling hormone-expressing cells (Fig. 5F; Fig. S3B). The expression of secretory components was further increased in the Fev^high cells and at this stage cells expressed the hormone secretion machinery (Fig. 5F; Fig. S3B). Interestingly, Fev^high cells highly expressed neurotransmitters and synaptic vesicle genes (Fig. 5F; Fig. S3B), again highlighting the molecular and cellular similarities between endocrine and neuronal cells.

Although all endocrine cells are derived from EPs, at which stage these cells become unipotent to produce distinct endocrine cell types remains largely unknown. It has been shown that mouse endocrine cells are generated in a step-wise manner, in which EPs are consecutively specified towards α-, β-, PP-, δ- and ε-cells during the secondary transition (Johansson et al., 2007). Indeed, the proportions of endocrine subtypes obtained at each time point reflected such a step-wise induction of endocrine cell fates. At E12.5, we only found α-cells, whereas β-, δ- and ε-cells appear between E13.5 and E14.5 (Fig. 6A). Interestingly, a portion of α-cells was also generated at later stages, suggesting that formation of these cells is not restricted to the primary transition (Fig. 6A; Fig. S4B,C). To investigate further the timing of endocrine specification, we estimated the RNA velocity of each cell (see Materials and Methods) (La Manno et al., 2018). RNA velocity is a proxy for the rate of pre-mRNA to mRNA processing and, therefore, is a vector in time that can be used to infer to which state a cell is moving in the high-dimensional transcriptome space. Thus, we can predict towards which endocrine subtype the EPs differentiate at a specific embryonic stage. As expected, we detected a movement of a portion of bipotent cells towards the endocrine lineage and a strong predicted directional flow from Ngn3^low to Ngn3^high to Fev^high cells. Interestingly, EPs and Fev^high cells at E12.5 and E13.5 pointed towards α-cells and to the few detected β-cells, whereas at E14.5 and E15.5 the majority of EPs showed a velocity towards β-cells (Fig. 6B). These data indicate that the α- and β-cell fate is specified in EPs mainly at early and late stages, respectively.

Furthermore, we found stage-dependent differential expression of several genes in differentiated α-cells (Fig. S4A). Genes such as Arx and Slc38a5 showed similar expression levels in early (E12.5) and late (E15.5) differentiated α-cells (Fig. 6C), whereas Sct, Mafb and Meis2 were differentially expressed. We found higher expression levels of Sct in early α-cells and higher expression levels of Mafb and Meis2 in late α-cells (Fig. 6D,E). This differential gene expression of early and late α-cells might be related to the formation of these cells during the first and secondary transition. However, we could not find α-cells with a primary transition signature at E14.5 and E15.5 (Fig. S4B), suggesting that the α-cells derived from early stages either acquire a similar signature to the late-stage α-cells or they have been eliminated during pancreas development.

One of the proposed mechanisms for stage-dependent endocrine specification is a change in the epithelial status and the surrounding niche. In mouse, the formation of endocrine cells coincides with the ramification of epithelial plexus to a single-layer epithelial network that possibly exposes the EPs to different ECM components and defines their fate towards specific endocrine cell types (Bakhti et al., 2019; Bankaitis et al., 2015). However, it is unclear whether EP transcriptional heterogeneity also primes endocrine subtype specification. To this end, we compared the expression levels of the newly described EP-signature genes between EPs derived from different developmental stages. Ngn3^low progenitors derived from different stages showed low expression levels of most of the EP-signature genes and, except for a few genes, they expressed comparable levels of most of these genes (Fig. S4C). In contrast, we found differential expression of several signature genes in Ngn3^high precursors derived from different developmental stages (E12.5-15.5) (Fig. 7A). We discovered that several genes, including Gng13, Steap1 and Gm8773, were expressed more highly in E12.5 Ngn3^high precursors compared with later stages, whereas other genes, including Upk3bl, Sultb21 and Neurod2, showed clearly increased expression in E15.5 Ngn3^high precursors compared with earlier stages (Fig. 7B-D). These data were further supported experimentally by performing qPCR analysis of isolated cells from E13.5 and E15.5 pancreata (Fig. 7E). We speculated that the differential expression of these genes at early and late stages is related to a fate decision of EPs towards distinct endocrine subtypes. In such a scenario, early precursors that express genes such as Gng13 and Gm8773 might be specified towards α-cells, whereas late precursors that express genes such as Upk3bl and Neurod2 might allocate to the β-cell fate. Indeed, PAGA analysis revealed that Ngn3^high precursors expressing high levels of Gng13, Steap1 and Gm8773 are highly connected to α-cells and Ngn3^high precursors expressing high levels of Upk3bl, Sultb2l and Neurod2 are highly connected to β-cells (Fig. 7F). This finding is in line with a recent study showing that Ngn3⁺ cells that express Myt1 are biased toward β-cell fate (Liu et al., 2019). The balance between the expression of Arx and Pax4 is the only currently known mechanism of endocrine cell subtype specification (Collombat et al., 2003, 2009). Although the expression of Pax4 was started at late stages of Ngn3^high precursors, the expression of Arx was initiated only in Fev^high cells (Fig. S4D), suggesting that the Pax4-Arx axis is likely not established in Ngn3^high precursors but rather in the Fev^high cells. The fact that transcriptional heterogeneity of EPs exists in Ngn3^high precursors proposes that these cells might be already determined towards specific endocrine cell fates. Therefore, our findings, together with other previous studies (Bechard et al., 2016; Desgraz and Herrera, 2009), support the idea that EPs are not unipotent progenitors at birth (Ngn3^low), but rather undergo cell-fate determination at the Ngn3^high state.

Improved therapy for diabetic patients could be achieved by islet cell replacement or in vivo regeneration of lost or dysfunctional β-cells. Both approaches rely on a detailed understanding of endocrinogenesis and the molecular programs driving endocrine lineage induction, specification and differentiation. Although gene regulatory networks and signaling cues that orchestrate endocrine cell differentiation have been well studied, the endogenous niche factors, morphogenetic cues and neighboring cell-type interactions priming endocrine cell specification are still only beginning to be understood. Uncovering mechanisms of endocrinogenesis is important for the generation of stem cell-derived islet-like clusters (ILCs) with defined cell composition and improved function (Bakhti et al., 2019).

Here, we generated a novel NVF reporter mouse line, which accurately reflects the spatiotemporal protein expression and regulation of the endogenous Ngn3 protein. Previous Ngn3 reporter mouse lines have been generated either by the replacement of Ngn3 by a fluorescent reporter gene (Lee et al., 2002) to generate heterozygous knock-in/knockout reporter mice or by transgenic expression of fluorescent proteins (FPs) regulated by the Ngn3 promoter (Bechard et al., 2016; Kim et al., 2015; Mellitzer et al., 2004). Although the onset of expression of the transgenic reporter FPs reflects the expression of endogenous Ngn3 protein in these animals, the transient expression of the endogenous Ngn3 protein is not reflected due to the relatively high stability of the FPs. Therefore, analysis of EPs using such reporter mice is affected by the existence of a portion of differentiated endocrine cells that still express FPs but not Ngn3, undermining the efficiency and accuracy of EP labeling/tracking and their isolation. In addition, one of these mouse lines (Lee et al., 2002) can be used only in the heterozygous condition, which impacts endocrine formation and differentiation due to heterozygosity (Wang et al., 2010). In contrast, our NVF reporter mouse line shows no detectable endocrine phenotype even in the homozygous state and the fusion of Venus to Ngn3 results in a transient expression kinetic of NVF that is comparable to that of the endogenous Ngn3 protein. These properties enabled us to isolate a large number of early and late EPs to perform a high-throughput single cell transcriptomic analysis to decipher lineage relationships and molecular changes during endocrinogenesis.

Using scRNA-seq during secondary transition, we identified eight main epithelial clusters and found the existence of Ngn3^low progenitors within MPCs, tip, trunk and ductal clusters. Among these, the presence of Ngn3^low progenitors within the MPCs and the tip cluster, although at lower numbers than in the ductal epithelial cells, was remarkable. Previously, a lineage-tracing study has shown that the pancreatic distal tip domain contains Cpa1⁺ MPCs that at E12.5 are tripotent and are able to give rise to exocrine, duct and endocrine cells, but only generate exocrine cells after E13.5 (Zhou et al., 2007). Moreover, another study has shown that the tip domain cells retain their tripotency until E12.5 and afterwards become restricted towards the acinar fate (Sznurkowska et al., 2018). Lineage-tracing analysis using Ptf1a^CreERTM mice has also indicated contribution of MPCs to all three pancreatic lineages until E14.0, when the number of MPCs significantly declines (Pan et al., 2013). Therefore, our finding supports these previous reports on the contribution of tip cells to EP formation at early stages before tip cells become unipotent, i.e. fate restricted towards the acinar fate. Although our data provide evidence of the existence of endocrine progenitors within the MPCs and tip cluster, they do not show whether EP lineage-primed progenitor cells are physically located within the tip MPC domain. Along this line, a recent study not only reported the contribution of tip cells to the EP pool at E14.5, but also showed the existence of Ptf1a⁺ EPs within the distal domain of the epithelium (Scavuzzo et al., 2018). Therefore, this study may provide the first evidence of the existence of EPs that are physically positioned within the tip domain at early stages before the tip cells become unipotent acinar cells.

We found that the major portion of endocrine progenitors existed within the ductal cluster, supporting a lineage-tracing study using Hnf1β⁺ cells that demonstrated the appearance of EPs from the ductal epithelium until E16.5 (Solar et al., 2009). Furthermore, the appearance of high numbers of Ngn3^low progenitors at E15.5 indicates that the endocrine progenitor pool is not set at very early stages. Instead, it continues to be generated during secondary transition, as has been suggested previously (Bechard et al., 2016). It is important to note that we only detected Ngn3^low progenitors, but not Ngn3^high precursors, within the MPCs, tip, trunk and ductal clusters. Therefore, it still remains unclear whether all cells expressing low levels of Ngn3 eventually differentiate into Ngn3^high precursors or if some of them differentiate back to the respective epithelial cells from which they are derived.

The high number of isolated Ngn3⁺ cells (7260) allowed us to obtain not only a high-confidence EP-enriched gene list (such as Amotl2, Sox4 and Hes6), but also provided a unique set of endocrine signature genes (such as Neurod2, Ppp1r14a, Gng13 and Sult2b1) with a similar expression kinetic as Ngn3 in EPs. To our knowledge, this is the first study reporting the transient expression of several genes other than Ngn3 in EPs. This unique EP-signature offers new markers for the labeling and/or tracking of these cells. Now, the important questions are whether the expression of these genes is controlled by Ngn3 and if they are upstream regulators or downstream effectors of Ngn3 activity. Interestingly, a recent time-resolved single cell study in the intestine identified several of these endocrine-enriched and -signature genes, such as Rcor2, Smarcd2, Sox4, Dll1, Eya2, Mycl and Neurod2 (Gehart et al., 2019). This data set together with our analysis further highlight significant similarities between the developmental programs for pancreatic and intestinal enteroendocrine cells. Remarkably, the Ngn3^low progenitors derived from different epithelial sources or at different developmental stages, did not show transcriptional heterogeneity in regard to the expression levels of EP-enriched and EP-signature genes. These data propose that, although in the Ngn3^low-state progenitors are biased towards an endocrine fate (Bechard et al., 2016), they might not yet be committed towards the endocrine fate and a specific endocrine subtype, leading us to question the unipotency of EPs at birth (Ngn3^low). In contrast, Ngn3^high precursors generated at different developmental stages expressed distinct levels of several EP-signature genes, such as Neurod2, Sult2b1, Gng13 and Gm8773. Moreover, PAGA analysis revealed that the differential expression of these EP-signature genes associates with the fate of Ngn3^high precursors towards α- or β-cells. These findings indicate that Ngn3^high precursors are not a homogenous population on the transcriptional level and suggest that endocrine subtype specification likely occurs in the Ngn3^high precursor state. This conclusion supports the previous in vivo clonal study that has shown the unipotency of EPs expressing high levels of Ngn3 (Desgraz and Herrera, 2009). Future work should address the function of those differentially expressed signature genes in early and late Ngn3^high precursors to provide deeper insights into the mechanisms underlying endocrine cell allocation. Such functional analysis might uncover the upstream cues defining the fate of endocrine cells through the establishment of the Arx-Pax4 axis. In addition, profound analysis of potential receptor-ligand interactions between EPs and endocrine subtypes could reveal unknown (feedback) cues from differentiated endocrine cells that further regulate induction of specific endocrine subtypes. Moreover, how many of those EP-signature genes are also expressed in human endocrine progenitors-precursors and whether they have a conserved function in endocrinogenesis deserves more attention. Finally, integration of real experimental time into pseudotime and trajectory inference can provide further insights into the dynamics of endocrine differentiation and the evolution of Ngn3^high precursors and their descendants over time (Fischer et al., 2019; Schiebinger et al., 2019). We anticipate that further development of computational methods to include large asynchronicity as found here will enable investigators to track back the differentiated endocrine subtypes to their origin and further elucidate their differentiation trajectories.

The Fev^high population has been recently reported as a novel endocrine cell state during pancreas development in mouse and human (Byrnes et al., 2018; Krentz et al., 2018; Ramond et al., 2018). These cells express low levels of Ngn3 (declining after the Ngn3^high precursor state), but also several early endocrine TFs and marker genes, such as Pax4, Neurod1 and Chga/b. Lineage-tracing analyses have indicated that these cells are derived from Ngn3⁺ EPs and differentiate into endocrine cells (Byrnes et al., 2018). This finding is supported by our PAGA and RNA velocity analysis indicating that differentiation of EPs into endocrine cells passes through a Fev⁺ stage. Indeed, most of EP-signature genes, such as Ppp1r14a, Sult2b1 and Sema3g, are still expressed (though with lower levels compared with that of EPs) in Fev^high cells, whereas the expression of many endocrine-specific genes, such as Chga/b, Pax6 and Cpe, is initiated and increasing in these cells. This observation suggests that although each stage of endocrinogenesis involves the expression of a set of highly specific genes, the differentiation of EPs into endocrine cells does not occur through a rapid transcriptional switch but rather passes through gradual alterations in gene expression profiles (Fig. 5E). In addition, we found that the transcriptional heterogeneity of Ngn3^high EPs from different developmental stages is also partially reflected in Fev^high cells (data not shown), suggesting that similar to Ngn3^high precursors, Fev^high cells are also unipotent cells with a defined fate towards a specific endocrine subtype.

In summary, we provide here a comprehensive high-resolution single cell gene expression profile of distinct pancreatic and endocrine progenitors and lineages during secondary transition. This information can be used as a blueprint for the generation of ILCs from stem cells. We identified the formation of Ngn3^low progenitors from the heterogeneous populations of tip, trunk and ductal cells, suggesting that only a subset of specified pancreatic progenitors can be primed towards EPs. Defining this subpopulation and the involved pathways will help to expand adequate numbers of endocrine progenitors. Furthermore, our analysis indicated step-wise activation/suppression of different signaling pathways such as cAMP and Wnt signaling, as well as changes in morphogenesis programs during the induction and differentiation of EPs that will aid the efficient production of differentiated endocrine subtypes in vitro. Finally, we identified 58 transiently expressed endocrine signature genes that subdivide Ngn3^high precursors into distinct subpopulations. Such markers and factors might help to discriminate EP subpopulations and allow differentiation towards a specific endocrine cell fate thus improving current protocols to generate ILCs with defined cell-type composition. In addition, understanding how EPs differentiate towards non-β- and β-cells might uncover the molecular programs required for redifferentiation of dedifferentiated β-cells or transdifferentiation of non-β- into β-cells for regenerative diabetes therapy.

Mice were kept at the central facilities at Helmholtz Zentrum München German Research Center of Environmental Health in accordance with the German animal welfare legislation and acknowledged guidelines of the Society of Laboratory Animals (GV-SOLAS) and of the Federation of Laboratory Animal Science Associations (FELASA). Post-mortem examination of organs was not subject to regulatory authorization.

The Ngn3-Venus fusion (NVF) construct was generated as described previously (Petrezselyova et al., 2015). Mouse IDG3.2 embryonic stem cells were electroporated with pCAG195 Cas9v2D10A-bpA (encoding for the nickase Cas9), pbs-U6-chimericRNA and linearized targeting vector. Using geneticin treatment, positive Neo-resistant clones were selected and the targeting of Ngn3-Venus was assessed by PCR genotyping. Next, the NVF-positive clones were aggregated with CD1 morulae and the resulting chimeric mice passed the NVF allele through their germline cells. Finally, the Neo cassette flanked by LoxP sites was deleted by intercrossing with the ROSACre mouse line. The excision of the Neo cassette was confirmed by genotyping PCR. The mouse colony was maintained with a mixed background (C57BL/6J×129/SvJ). We named this mouse line the Ngn3-Venus fusion (NVF) reporter line. The mouse line will be available to the research community upon direct request to H.L.

Embryonic pancreata from NVF homozygous mice were dissected and pooled together for each stage as follows: 29 pancreata from E12.5, 35 pancreata from E13.5, ten pancreata from E14.5 and seven pancreata from E15.5. Next, pancreata were kept in 0.25% Trypsin for 5 min on ice and then incubated at 37°C for 10 min. The single cell samples were then centrifuged at 1700 rpm (290 g) for 5 min at 4°C. The supernatant was removed and cells were counted. 5 µl anti-mouse CD326 (EpCAM) PE (eBioscience, 12-5791-81) and rat IgG2a K isotype control (eBioscience, 12-4321-42) were used for 1×10⁶ cells in 100 µl total volume. After staining for 30 min at 4°C, the cells were stained with DAPI to detect dead cells. The samples were then washed twice and resuspended in FACS buffer (PBS, 1% BSA, 0.5 mM EDTA) and loaded for FACS sorting. The gating strategy was as follows: main population>single cells>living cells (DAPI negative)>Ngn3⁺ (FITC⁺) and EpCAM⁺ (PE⁺) cells. The collected cells were enriched for Ngn3⁺ (∼2/3) and EpCam⁺ cells (∼1/3) in FACS buffer containing 1% BSA and 0.1 mM EDTA. Finally, cells were counted and the number of dead cells was identified by Trypan Blue staining. If less than 20% of the cells were dead, cells were processed for single cell RNA seq.

Single cell libraries were generated using the Chromium Single Cell 3′ library and gel bead kit v2 (PN #120237) from 10x Genomics. Briefly, to reach a target cell number of 10,000 cells per sample 16,000 cells per sample were loaded onto a channel of the 10x chip to produce Gel Bead-in-Emulsions (GEMs). This underwent reverse transcription to barcode RNA before cleanup and cDNA amplification followed by enzymatic fragmentation and 5′adaptor and sample index attachment. Libraries were sequenced on the HiSeq4000 (Illumina) with 150 bp paired-end sequencing of read2.

qPCR was performed using TaqMan probes and 25 ng cDNA per reaction. Each reaction consisted of 4.5 µl cDNA in nuclease-free water, 5 µl TaqMan Advanced master mix (Life Technologies) and 0.5 µl TaqMan probe. The qPCR was performed using Viia7 (Thermo Fisher Scientific). C_t-values were normalized among samples, transformed to linear expression values, normalized to reference genes and to control samples, i.e. relative expression (gene)=(2^{Ct(mean genes)−Ct(gene)})/(2^{Ct(mean references)−Ct(reference)}), and normalized expression (gene)=relative expression (gene)/relative expression_control (gene).

See supplementary Materials and Methods for details of the probes used.

Embryonic pancreata were dissected and fixed in 4% paraformaldehyde (PFA) in PBS overnight at 4°C. The tissues were equilibrated in 10% and 30% sucrose-PBS solutions at RT (2 h each solution) followed by 1:1 tissue-freezing medium:30% sucrose overnight at 4°C. Afterwards, they were embedded in cryoblocks using tissue-freezing medium (Leica, 14020108926) and sections of 20 μm were cut. Next, samples were permeabilized (0.1% Triton X-100, 0.1 M glycin) for 15 min and incubated in blocking solution (10% fetal calf serum, 3% donkey serum, 0.1% BSA and 0.1% Tween-20 in PBS) for 1 h at room temperature (RT). Primary antibodies were diluted in blocking solution and incubated with the samples overnight at 4°C. After washing with PBS, samples were stained with secondary antibodies diluted in the blocking solution for 3-5 h at RT. The samples were then incubated with 4′, 6-diamidin-2-phenylindol, followed by three washes with PBS and embedded in commercial medium ProLong Gold (Life Technologies). All images were obtained with a DMI 6000Leica microscope using LAS AF software. Images were analyzed using LAS AF and ImageJ software programs. See supplementary Materials and Methods for details of the antibodies used.

The isolation of islets was performed using collagenase P (Roche) digestion of adult pancreas as described previously (Bastidas-Ponce et al., 2017b). See supplementary Materials and Methods for further details.

For GSIS analysis, the isolated islets were cultured overnight before transferring to a 96-well plate containing modified Krebs Ringer phosphate Hepes (KRPH) buffer with 1 mM glucose for 1 h. Different glucose concentrations (2.8 and 16.8 mM) were added to the islets (1 h for each). The supernatants were collected and used for insulin measurement. At the end, the islets were lysed for DNA measurements. Insulin concentrations were measured and quantified using an ultrasensitive insulin ELISA kit (CristalChem). The analysis was performed using a standard curve, and the data were normalized to the DNA content.

Mice were maintained in standard conditions, and were starved 4 h before the measurements. Blood glucose values were determined from venous blood using an automatic glucose monitor (Glucometer Elite, Bayer).

The following data processing and computational analyses of single cell data were used: preprocessing of droplet-based single cell RNA-seq data, low dimensional embedding, visualization and clustering, marker gene identification and subtype characterization, cell cycle classification, reconstruction of lineage relationships and differentiation trajectories, directionality of endocrine differentiation using RNA velocity estimation and software specifications. See supplementary Materials and Methods for further details of each analysis.

Data are presented as mean±s.d. Two-tailed unpaired t-tests were used to analyze the data.

We thank Jessica Jaki, Bianca Vogel, Kerstin Diemer and Julia Beckenbauer for technical support. We thank Maike Sander and Helena Edlund for providing Ngn3 antibodies.