Lineage dynamics of murine pancreatic development at single-cell resolution

Abstract

Organogenesis requires the complex interactions of multiple cell lineages that coordinate their expansion, differentiation, and maturation over time. Here, we profile the cell types within the epithelial and mesenchymal compartments of the murine pancreas across developmental time using a combination of single-cell RNA sequencing, immunofluorescence, in situ hybridization, and genetic lineage tracing. We identify previously underappreciated cellular heterogeneity of the developing mesenchyme and reconstruct potential lineage relationships among the pancreatic mesothelium and mesenchymal cell types. Within the epithelium, we find a previously undescribed endocrine progenitor population, as well as an analogous population in both human fetal tissue and human embryonic stem cells differentiating toward a pancreatic beta cell fate. Further, we identify candidate transcriptional regulators along the differentiation trajectory of this population toward the alpha or beta cell lineages. This work establishes a roadmap of pancreatic development and demonstrates the broad utility of this approach for understanding lineage dynamics in developing organs.

Fig. 1

Single-cell sequencing identifies broad patterns of cellular heterogeneity in E14.5 murine pancreas. a Overview of murine pancreatic development. b Schematic of experimental approach. c t-Distributed stochastic neighbor embedding (t-SNE) visualization of populations from pooled E14.5 mouse pancreata (n = 14). Each dot represents the transcriptome of a single cell, color-coded according to its cellular identity (epithelial, mesenchymal, or immune/vascular). Each cell compartment contains multiple subpopulations, represented by varying degrees of color shading. d Established marker genes identify epithelial cells (Cdh1+), endocrine cells (Chga+), mesenchymal cells (Vim+ and Col3a1+), endothelial cells (Pecam1+), and immune cells (Rac2+). e Heatmap depicting greater than two-fold differentially expressed genes in each cluster compared to all other clusters. Cells are represented in columns, and genes in rows. Specific genes used to annotate clusters are indicated to the right of the heatmap

Fig. 2

Identification of multiple uncharacterized mesenchymal populations. a t-SNE visualization of subclustered E14.5 mesenchymal clusters (from n = 14 pancreata). b Density plot depicting Pearson’s correlation values (depicted in heatmap in Supplementary Fig. 3b) within the epithelial and mesenchymal populations based on average gene expression in each cluster. c Dot plot of top differentially expressed markers of each mesenchymal population. Bars are color-coded by cluster identity in a. The gray bar represents pan-mesenchymal markers. The size of each dot represents the proportion of cells within a given population that expresses the gene; the intensity of color indicates the average level of expression. d Pathway analysis of genes greater than 2-fold differentially expressed by cells in clusters 1, 2, 4, and 5. e Expression of genes marking clusters 1 (Cav1), 2 (Stmn2), 4 (Cxcl12), and 5 (Barx1) in all E14.5 mesenchymal cells. Color intensity indicates level of expression. f–h Multiplexed fluorescent ISH combined with Epcam IF validates clusters 2 and 5 (f) and cluster 1 (g-h). Epcam marks pancreatic epithelium. In (f), Barx1 + cells (red arrows, cluster 5) are distinct from Stmn2+ cells (green arrows, cluster 2), validating the single-cell data. In (g), Cav1+ cells (red arrows, cluster 1) are distinct from Stmn2+ cells (green arrows, cluster 2). In (h), Barx1+ cells that do not express Cav1 (red arrows) represent cluster 5, whereas Barx1+/Cav1+ cells (yellow arrows) represent cluster 1. Cav1+ cells that do not express Barx1 are also identified (green arrows), likely representing endothelial cells. Scale bar represents 50 µm in f–h

Fig. 3

Mesothelial cells are dynamic over developmental time and are predicted to give rise to vascular smooth muscle populations. a t-SNE visualization of merged mesenchymal clusters from E12.5 (n = 18 pancreata), E14.5 (n = 14 pancreata for batch 1; n = 11 for batch 2), and E17.5 (n = 8 pancreata) tissue. Mesenchymal clusters were identified at each timepoint, subclustered, merged together, and reanalyzed. Cells are colored by cluster or timepoint. Dotted circle highlights timepoint-segregated mesothelial clusters. b Dot plot of top differentially expressed genes in timepoint-specific mesothelial clusters (clusters 1, 11, and 17). Size of the dot represents proportion of the population that expresses each specified marker. Color indicates level of expression. c ISH for Pitx2 and Msln in E12.5 and E17.5 pancreata. Pitx2 expression was detected in E12.5, but not E17.5 mesothelium, whereas Msln was detected in E17.5, but not E12.5 mesothelium. Vimentin (Vim) IF staining depicts pancreatic mesenchyme. Dotted line indicates tissue boundary. Yellow arrows identify Pitx2+ mesothelial cells. Red arrows identify Msln+ mesothelial cells. Scale bar represents 50 µm. d Expression levels of VSM-related genes in merged mesenchymal clusters. Color intensity indicates level of expression. e Pseudotime ordering of mesothelial and VSM-related merged mesenchymal clusters. Colors correspond to t-SNE in a. All clusters are individually plotted in Supplementary Fig. 3j. f Cluster proportions over pseudotime. Pseudotime was binned into ten groups and the proportion of each cluster within that bin of pseudotime was calculated. g Model of lineage relationships among mesothelial and VSM-related mesenchymal populations based on pseudotime ordering in e

Fig. 4

Identification of epithelial cell populations in E14.5 mouse pancreas. a t-SNE visualization of epithelial groups only, as defined in Fig. 1. b Dot plot depicting known and uncharacterized markers of epithelial populations, as well as markers specific to the Fev^Hi population. Size of the dot represents proportion of the population that expresses each specified marker. Color indicates level of expression. c Expression of Fev and Ngn3 within epithelial cells. Color indicates level of expression. d Gene expression comparison between the Ngn3+ and Fev^Hi population. Genes greater than 2-fold differentially expressed are highlighted in dark blue (higher in Fev^Hi cells) or light blue (higher in Ngn3+ cells). e Pathway analysis of genes greater than 2-fold differentially expressed in Ngn3+ and Fev^Hi populations. f t-SNE visualization of the 661 cells of the endocrine lineage (Ngn3+, Fev^Hi, alpha, beta, and epsilon populations). g Pseudotime ordering of Ngn3+, Fev+/Pax4+, Fev^Hi, alpha, and beta cell populations place Fev+ cells between Ngn3+ and hormone+ populations

Fig. 5

Fev^Hi cells are endocrine progenitors. a In situ hybridization (ISH) for Ngn3, Fev, and Isl1 in lineage-traced Ngn3-Cre; Rosa26^mTmG E14.5 pancreata where Ngn3-lineage-traced cells are mGFP+. Gray arrowheads identify Ngn3+ cells, presumably not yet Ngn3-lineage labeled due to the transient nature of Ngn3 expression and the delay of Cre-mediated recombination that permits expression of mGFP. Blue arrowheads identify Ngn3+/Fev+ cells that are Ngn3-lineage-traced. Yellow arrowheads identify Ngn3-lineage-traced cells that are Fev+, but do not express Ngn3 or Isl1. Purple arrowheads identify Fev+/Isl1+ cells that are Ngn3-lineage-traced. Magenta arrowheads identify Isl1+ cells that are Ngn3-lineage-traced. b–c Dual ISH/immunofluorescence (IF) for NGN3 and FEV mRNA and CHGA protein in human fetal pancreas at 23 weeks of gestation (n = 1 pancreas). Gray arrowheads identify NGN3+ cells. Yellow arrowheads identify FEV+ cells. Purple arrowheads identify FEV+/CHGA+ cells. Magenta arrowheads identify CHGA+ cells. d Multiplexed fluorescent ISH for NGN3, FEV, and ISL1 mRNA in hESC-derived endocrine progenitor cells. Blue arrowheads identify NGN3+/FEV+ cells. Yellow arrowheads identify FEV+ cells. Purple arrowheads identify FEV+/ISL1+ cells. e Quantification of each population detected in Ngn3-lineage-traced pancreata as a percentage of Ngn3-lineage-traced cells (n = 464 cells, 6 pancreata). Data are represented as mean + standard deviation (SD). f Quantification of each population detected in hESC-derived progenitor cells as a percentage of total stained cells (n = 418 cells, 3 clusters representing technical replicates from one hESC differentiation). Data are represented as mean + SD. g Proposed model for the derivation of Fev^Hi endocrine cells from Ngn3+ cells, and their differentiation into hormone+/Fev^Lo endocrine cells. Colors of arrowheads and bars in a–f correspond to cell identity in g. a, d Scale bar: 10 µm. b, c Scale bar: 20 µm. h t-SNE visualization of v2 merged endocrine timecourse (E12.5, E14.5, aggregated E17.5). Clusters are annotated based on correlation with v1 dataset or top differentially expressed genes. i Timepoint labels for v2 merged endocrine timecourse data. t-SNE is the same as h. j Cell type proportions at each timepoint, calculated from the clusters depicted in h

Fig. 6

Differentiated hormone+ endocrine cells transit through a Fev-expressing stage during pancreatic development. a–e Dual IF (for membrane GFP) and fluorescent ISH for hormones in Fev-Cre; ROSA26^mTmG lineage-traced animals at E14.5. n = 46 cells of 4 pancreata for Ins1 (100% labeled-lineage); n = 103 of 4 pancreata cells for Gcg (100% lineage-labeled); n = 6 cells of 2 pancreata for Sst (100% lineage-labeled); n = 26 cells of 2 pancreata for Ghrl/Gcg (23.2% lineage-labeled); n = 71 cells of 8 pancreata for Ppy (90.1% lineage labeled). Scale bar represents 10 µm. f Schematic of E14.5 Fev-Cre; ROSA26^mTmG FACS sorting and single-cell RNA sequencing. g Representative FACS plots of sorted single, live GFP+ and TdTomato+/GFP− cells from dissociated pancreata used for single-cell sequencing. h t-SNE visualization of endocrine cells in Fev-lineage-traced E14.5 mouse pancreata (n = 3). i Expression of major markers of endocrine cell types. Color indicates level of expression, except for the eGFP plot, which indicates presence or absence of eGFP counts

Fig. 7

Identification of candidate regulators of beta and alpha cell fate decisions. a Pseudotime ordering of the endocrine cells at E14.5 depicted in Fig. 6h yields a bifurcated tree in which the two main branches terminate in cells that highly express Ins1 (beta cell branch) or Gcg (alpha cell branch). b Heatmap depicting the expression of genes along each branch, in pseudotime. An independent expression pattern is calculated across the entire pseudotime trajectory for each branch. Therefore, the portion of the trajectory before the branch point is displayed for each branch separately. Genes are clustered based on expression pattern across pseudotime; selected genes with differential expression along the branches are highlighted to the right. c Gene expression plots depicting the kinetic trends along each branch. d–e Multiplexed fluorescent ISH for Fev, Gng12, and Islet1 (d) or Fev, Peg10, and Islet1 (e) in lineage-traced E14.5 Ngn3-Cre; ROSA26^mTmG pancreas. Arrowheads identify lineage-traced Fev+/Islet1− cells with Gng12 (d, teal gradient arrowheads) or Peg10 (e, indigo gradient arrowheads) expression. f Multiplexed fluorescent ISH for Fev, Gng12, and Ins1. Teal arrowheads identify lineage-traced Ins1+ beta cells that express Gng12. g Multiplexed fluorescent ISH for Fev, Peg10, and Gcg. Indigo arrowheads identify lineage-traced Gcg+ alpha cells that express Peg10. h Model for Fev^Hi (yellow) cell differentiation into distinct alpha or beta cells. Peg10 and Gng12 expression in Fev^Hi cells may represent progenitors pre-fated towards the alpha and beta lineages, respectively, during endocrine lineage allocation. d–g Scale bars represent 10 µm. Blue staining represents DAPI-labeled nuclei. Colors of arrowheads match colors of cells represented in h