Single-cell resolution analysis of the human pancreatic ductal progenitor cell niche.

We have described multipotent progenitor-like cells within the major pancreatic ducts (MPDs) of the human pancreas. They express PDX1, its surrogate surface marker P2RY1, and the bone morphogenetic protein (BMP) receptor 1A (BMPR1A)/activin-like kinase 3 (ALK3), but not carbonic anhydrase II (CAII). Here we report the single-cell RNA sequencing (scRNA-seq) of ALK3bright+-sorted ductal cells, a fraction that harbors BMP-responsive progenitor-like cells. Our analysis unveiled the existence of multiple subpopulations along two major axes, one that encompasses a gradient of ductal cell differentiation stages, and another featuring cells with transitional phenotypes toward acinar tissue. A third potential ducto-endocrine axis is revealed upon integration of the ALK3bright+ dataset with a single-cell whole-pancreas transcriptome. When transplanted into immunodeficient mice, P2RY1+/ALK3bright+ populations (enriched in PDX1+/ALK3+/CAII- cells) differentiate into all pancreatic lineages, including functional β-cells. This process is accelerated when hosts are treated systemically with an ALK3 agonist. We found PDX1+/ALK3+/CAII- progenitor-like cells in the MPDs of types 1 and 2 diabetes donors, regardless of the duration of the disease. Our findings open the door to the pharmacological activation of progenitor cells in situ.

The existence of progenitor-like cells in the human pancreatic duct tree, despite conflicting findings in mice (1, 2), is supported by experimental evidence (reviewed in ref. 3). We have recently characterized one such cell population in human major pancreatic ducts (MPDs) and pancreatic duct glands (PDGs) (4, 5). These cells express the ductal/β-cell marker pancreatic and duodenal homeobox-1 gene (PDX1) and the bone morphogenetic protein (BMP) receptor 1A (BMPR1A), also known as activin-like kinase-3 (ALK3). Antibodies against ALK3 and the purinergic receptor P2Y1 (P2RY1), which we have established as a surrogate surface marker for PDX1 (4), can be used to sort live progenitor-like cells that expand in vitro in the presence of BMP-7 and differentiate along all three adult pancreatic lineages (endocrine, acinar, and ductal) upon removal of this factor (4).

This marker combination (P2RY1+/ALK3bright+) specifically excludes islets, which are ALK3− (4). Another salient feature of these cells is that, notwithstanding their location in the MPD epithelium, they lack protein expression of carbonic anhydrase II (CAII), a marker previously thought to be pan-ductal (6). Thus, PDX1+/ALK3+/CAII− cells, which we hypothesized to be progenitors, are readily detected in human MPDs and PDGs, often intercalated with regular ductal PDX1+/ALK3+/CAII+ cells (4).

In this report, we conducted single-cell RNA-sequencing (scRNA-seq) analysis of sorted human pancreatic ALK3bright+ cells, which offers a “high-magnification” view of all of the cell types within the ductal tree. This analysis shows the existence of multiple ductal subpopulations that run the gamut from mature functional ductal cells to undifferentiated/transitional cells, including progenitor-like cells with a transcriptional profile consistent with stress-induced de-differentiation. P2RY1+/ALK3bright+ sorting enriches in CAII− cells from such progenitor-like clusters. We demonstrate that sorted P2RY1+/ALK3bright+ cells engraft in immunocompromised mice, where they spontaneously mature into multiple pancreatic endocrine and exocrine cell types. This process was accelerated by the systemic administration of THR-123, a cyclic peptide with ALK3 agonism activity. Next, we sought to explore whether PDX1+/ALK3+/CAII− existed in samples from donors with type 1 and 2 diabetes (T1D/T2D). We confirm that progenitor-like cells exist in the MPD/PDG epithelium of all patients examined, regardless of the duration of the disease. Our findings that human P2RY1+/ALK3bright+ cells can be pharmacologically stimulated in vivo, combined with the observation that similar cells may be intact in T1D/T2D patients, have potential therapeutic implications.

Results

scRNA-Seq Analysis of ALK3+ Cells Reveals Ductal Heterogeneity and Suggests Transitional Phenotypes.

We have previously determined that the P2RY1+/ALK3bright+ fraction of the human exocrine pancreas is enriched in PDX1+/ALK3+/CAII− progenitor-like cells that can be cultured in defined conditions, expand in the presence of BMP-7, and exhibit multipotency when BMP-7 is removed. To further analyze in depth the composition of the human ALK3+ ductal compartment, we conducted scRNA-seq of sorted ALK3bright+ cells from human nonendocrine pancreatic tissue obtained after islet isolations at the Diabetes Research Institute (two donors) and the University of Pittsburgh (one donor) (). Our rationale was that, by analyzing the ALK3bright+ population at large, we would be able to study the PDX1(P2RY1)+/ALK3+ cells (a subpopulation thereof) in their broader context within the ductal tree. Fig. 1 presents a scheme of the experimental plan. The expected outcome of this sorting strategy is the selection of epithelial cells from the MPDs and PDGs with the exclusion of a significant percentage of ALK3dim+ cells (likely from small ducts), as determined in Qadir et al. (4). We also depleted contaminating mesenchymal cells by CD90− selection, further enriching in epithelial cells.

Fig. 1.

scRNA-seq identifies cellular heterogeneity across human ALK3bright+ pancreatic ductal cells. (A) Scheme of the experimental approach employed to isolate ALK3bright+ cells and derive single-cell libraries using the 10X Genomics microfluidics system. Data generated from sequenced single-cell libraries are analyzed to dissect cellular substructure. (B) UMAP plots of ductal populations from individual donor pancreas (n = 3, 4,878 cells). Each dot represents the transcriptome of a single ALK3bright+ cell, with color coding defining clusters of cells having similar transcriptional identities. (C) Gene-expression heatmap of the top 20 most DE genes in each cluster compared to all other clusters. Genes are represented in rows and cell clusters in columns. Select genes are indicated to the left. (D) Single-cell gene expression of known markers of ductal (KRT19 and SOX9), epithelial (CDH1), acinar (AMY2B), and endocrine (ISL1, INS, GCG, SST, PPY) cells. Scale bars represent z-test–normalized gene expression in C and gene counts in D.

Approximately 2,000 cells per preparation were sequenced at 250,000 reads per cell (Genewiz). Only single cells with mitochondrial RNA less than 20% of total RNA were included (). Using principal component analysis and clustering with the R v.3.5.3 package Seurat 3.1.1 (), an uniform manifold approximation and projection (UMAP) plot of the combined three preparations (4,878 cells) was generated (Fig. 1 and ). UMAP is a novel nonlinear dimensionality reduction algorithm for the analysis of high-dimensional data (7), which has proven superior to t-distributed stochastic neighborhood embedding for resolving subtle differences in cell populations (8–11). A scatterplot of differentially expressed (DE) genes and a heatmap with the top 20 DE genes/cluster (vs. all of the other clusters) are presented in and Fig. 1C, respectively. Clustree v0.4.1 (12) was used to choose the highest clustering resolution without cluster destabilization or mixing (). present the individual contribution of each donor to the combined dataset. The overall epithelial/ductal identity of the sorted cells was confirmed by the widespread expression of the ductal markers KRT19 and SOX9, the epithelial marker CDH1, and the absence of any meaningful expression of the archetypal acinar marker AMY2B and the endocrine markers ISL1, insulin (INS), glucagon (GCG), somatostatin (SST), and PPY (Fig. 1).

The only exception was the smaller cluster 7, which was largely KRT19 and SOX9. The analysis of DE genes and gene ontology (GO)-term pathways (Fig. 2 and Datasets S1 and S2) suggests an immune cell identity. BMP signaling regulates proliferation/differentiation in the immune system (13), and ALK2/ALK3 are functional BMP receptors in macrophages. Therefore, it is likely that ALK3+ sorting resulted in the carryover of this small immune-related subpopulation.

Fig. 2.

Identification of multiple subpopulations of ALK3bright+ ductal cells in the human pancreas. (A) UMAP plot mapping multiple subclusters across the single-cell dataset (n = 3). Each cluster is defined by a specific color and a representative DE gene. The dotted line divides ducto-acinar or “transitional to acinar” clusters (Upper) from “productal” clusters (Lower). (B) Gene-expression dot plot of select DE genes in each cluster. Rows depict clusters, while columns depict genes. The intensity of any given point indicates average expression, while its size represents the proportion of cells expressing a particular gene. Dotted regions identify select genes for a particular cluster. The scale bar represents average expression in gene counts. (C) UMAP expression plots of representative DE genes for clusters 3 (small ducts; KRTAP2-3, AKAP12); 2 (activated/migrating progenitor cells; TFF1, IGFBP3); 1 (stress/harboring progenitor-like cells; CRP, SPP1); 4 (centroacinar; WSB1, XIST); 5 (transitional to acinar 1; CEL, OLFM4); 6 (transitional to acinar 2; SYCN, CPA2); and 7 (immune cells; SRGN, IDO1). (D) GO-driven pathway analysis of DE genes with a ≥twofold change across clusters 1 to 7. (E) Violin plots for the expression of acinar (Left) and ductal genes (Right). Expression in each cell is shown along with the probability density of gene expression, denoted by the shape of the plot.

Fig. 2 shows the predicted identity of each cluster according to the above analysis. Representative top DE genes and GO pathways are indicated in Fig. 2 and , as well as Datasets S1 and S2. For identification purposes, we selected those whose expression is most specific to any given cluster vs. the rest. We detected two major groups of cells: Productal (clusters 1 to 4 in Fig. 2 , Lower) and ducto-acinar (clusters 5 to 6 in Fig. 2 , Upper). The distinction between ducto-acinar and productal is based on the DE of numerous acinar genes in clusters 5 to 6 (Fig. 2 and )

Cluster 1 is represented by osteopontin/secreted phosphoprotein 1 (SPP1). This gene was previously identified as a marker of undifferentiated pancreatic progenitors in mice (14), and also by scRNA-seq in the developing mouse pancreas as a marker of “proliferative ducts” (15). Cluster 1 also has the highest PDX1 and inhibition of differentiation (ID)-1, -2, and -4 expression in the entire dataset (). High expression of ID1-4 is typical of progenitors, where they sequester basic helix–loop–helix transcription factors, leading to the inhibition of lineage-specific and cell cycle-inhibitory genes (16). ID2, in particular, has been shown to regulate BMP-dependent pancreatic progenitor expansion in mice (17). Elevated ID expression was also a hallmark of the P2RY1+/ALK3bright+ cells previously described by our team as progenitor-like, especially after the addition of BMP-7 (4). In line with the hypothesis that cluster 1 may harbor the multipotent progenitors therein described, among the most DE genes we also found HES1, FOS, and JUN (Datasets S1 and S2). HES1 is an inhibitor of differentiation up-regulated by Notch signaling, which is associated with progenitor proliferation during pancreatic development (18). FOS and JUN are the two components of the activator protein 1 early response transcription factor, which is essential for the proliferation of ductal epithelial cells in the context of pancreatic cancer. While there is no evidence that the donors used in this study had any malignancy, we conducted a search for DE pancreatic cancer pathways (19) to rule out any such conditions, which may have confounded the interpretation of our data. This search was negative, further reinforcing the notion that this cluster may harbor progenitor cells instead. In fact, also among the top DE genes in this cluster is TFPI2, a pancreatic tumor suppressor (20).

Inflammation/stress responses were also recurrent themes from the DE and GO pathway analyses conducted for cluster 1. In fact, C-reactive protein (CRP) was the top DE gene. The most likely source of stress in the ducts is protease secretion. Several serpins [a superfamily of proteins that protect exocrine cells against their own secreted proteases (21)] are indeed among the top DE genes in this cluster (Datasets S1 and S2). Therefore, the above progenitor-like signature is also consistent with potential de-differentiation processes associated with stress. Although our data do not offer direct confirmation of this hypothesis, stress-mediated de-differentiation has been extensively reported in the ductal tree (reviewed in ref. 3).

To validate gene-expression findings, we examined the immunostaining of select DE expressed markers for each cluster in the Human Protein Atlas, a comprehensive resource of more than 24,000 antibodies and 13 million immunohistochemistry images (22, 23). Since differential expression is not exclusive expression, we often found that markers with a positive signal were not restricted to specific pancreatic compartments. Due to this, we focused on the regions that exhibited the most defined staining for any given DE gene. For some of the top DE genes in cluster 1 (FOS, TSPAN8, SPP1) (Datasets S1 and S2), a strong signal was commonly observed within larger ductal structures containing a fibromuscular wall and in connection with PDGs (). The presence of Mucin 20 among the DE genes further confirms the potential link of these cells to PDGs, which have long been proposed to be niches of pancreatic progenitors (24–27).

DE GO pathways in cluster 2 (trefoil factor 1, TFF1) reflect active extracellular matrix remodeling, as evidenced by the presence of multiple metalloproteinases (MMP1, MMP7, MMP10) among the top DE genes. This list is headed by TFF1, with TFF2 close behind. The concerted action of these two genes mediates the migration of cells out of PDGs following inflammatory injury (26). Epithelial cell migration, alongside extracellular matrix remodeling are, in fact, the top GO pathways differentially up-regulated for this cluster. Another interesting DE gene is c19orf33, which encodes the hepatocyte growth factor activator inhibitor type 2-related small peptide (H2RSP). This protein is expressed in the epithelium of the gastrointestinal tract, where it is up-regulated in response to injury. The nuclear translocation of H2RSP has been suggested to have a role in the transition from proliferation to differentiation (28). Progenitor markers PDX1, and ID1-4 are also relatively elevated (). HMGA1 [a transcription factor that mediates motility and has been associated with pancreatic adenocarcinoma invasiveness (29, 30)] and S100A6 [calcium-binding protein also associated with pancreatic cancer cell motility (31)] are two other DE genes in this cluster that further suggest a migratory phenotype. As before, a search for DE pancreatic cancer pathways (19) was negative. Immunostaining for these markers shows that the cells with the strongest signal are either in major ductal/PDG structures or seemingly delaminating from them (). Single cells positive for these markers can be detected throughout the parenchyma, including acini and islets.

Cluster 3 (AKAP12) exhibits the DE of multiple cytokeratins (KRT17, KRT81, KRT23, KRT7, KRTAP2-3), which are markers of epithelial differentiation (32). GO analysis shows the differential elevation of homotypic cell–cell adhesion, positive regulation of cell differentiation, and negative regulation of growth pathways. Taken together, these features suggest terminally differentiated ducts. Immunostaining of some of the top DE genes (KRT17, KRT23, ANXA1, ANXA3, LAMA3) () confirms widespread staining in small ducts, with larger ones sometimes displaying intermittent signal patterns.

GO pathway analysis of cluster 4 (WSB1) reveals cells with an active transcriptional profile and the up-regulation of pathways, such as negative regulation of differentiation, Notch, and Wnt signaling. Interestingly, the immunostaining of some of the top DE genes for this cluster suggests a centroacinar location (). Centroacinar cells are difficult to study because of the lack of markers to distinguish them from the adjacent intercalated small ductal cells (33). Studies in mice and zebrafish suggest that these cells express high levels of aldehyde dehydrogenases (ALDH), which have been proposed as markers for populations enriched in progenitors (34). However, this has not been confirmed in humans, where centroacinar ALDH activity has only been detected recently in children and during pathological processes (35). We did not observe DE of ALDHs in this cluster, and the Human Protein Atlas does not show differential centroacinar expression of these markers. However, ALK3bright+ sorting may have excluded many centroacinar cells, as ALK3 expression is brightest in major ducts but tends to be dimmer in smaller ducts (4). Therefore, it is possible that our dataset comprises only a fraction of the overall centroacinar population. Still, we base the hypothetical assignment of centroacinar identity to this cluster on several observations. First, its position at the crossroads between productal and ducto-acinar axes is consistent with cell types known to be at the interface between acinar and ductal tissues. Second, as indicated above, the immunostaining profile of some DE markers (), including epithelium-specific transcription factor-3 (ELF3), vacuole membrane protein-1 (VPMP1), or WD repeat and SOCS box containing-1 (WSB1) suggests a centroacinar location. Many studies suggest a direct involvement of centroacinar cells in neoplastic transformation trough centroacinar-to-ductal transdifferentiation (36). In this context, mutated forms of the above DE markers ELF3 (37), VMP1 (38), and WSB1 (39) have all been associated to neoplastic transformation of the pancreas. Finally, Notch signaling, which has been characteristically associated to centroacinar cells (40), is one of the top DE pathways in GO analysis. Taken together, the above evidence is consistent with the hypothesis that cells from cluster 4 are centroacinar.

The ducto-acinar clusters (5, 6), while predominantly ductal in nature (they are mostly negative for the archetypal acinar marker AMY2B and positive for the ductal marker KRT19), also feature the differential expression of acinar-specific genes, such as CEL, REG3A, REG1A, or CPA1 for cluster 5 (OLFM4+) and CLPS, CELA2A/3A/3B, or PGA5 for cluster 6 (CPA2) (Fig. 2 and ). The latter already shows up-regulation of typical acinar cellular processes, such as zymogen activation (41) and the metabolization of glutathione (42). Immunostaining of most of these DE genes for clusters 5 and 6 is widespread throughout the pancreas and difficult to pinpoint to specific locations, since proteases are secreted by the acinar tissue into the ducts. However, some of these makers exhibit a distinct centroacinar position ().

Taken together, our data suggest that the ALK3bright+-sorted population harbors an assortment of ductal cell types, ranging from cells with progenitor-like features to mature ducts, as well as subpopulations that exhibit some acinar-like expression patterns.

Pseudotemporal Analysis Suggests Two Transitional Axes in the ALK3bright+ Ductal Compartment.

Given the indications of potential differentiation gradients between ductal and acinar cells in clusters 4-5-6, as well as between different ductal cell types in clusters 1-2-3, we used Monocle v2.10.1 to arrange the cells in an unbiased pseudotemporal manner; that is, following a “developmental/transitional” path according to their transcriptomic similarity. This strategy has been used to determine lineage dynamics during pancreatic development (43). The results of this analysis after the bioinformatics removal of the immune cluster 7 are shown in Fig. 3 . The top “ducto-acinar” arm of the graphic in Fig. 3 indicates that there is a pseudotemporal continuum, with a progressive acquisition of acinarity following the sequence from clusters 4 (centroacinar) to 6 (transitional to acinar 2). Such inference, while solely based on bioinformatics tools, strongly suggests the existence of cells transitioning from ductal to acinar phenotypes (4-5-6) or vice versa (6-5-4). Similarly, the bottom “productal” arms are also arranged in a manner that could be interpreted as progenitor cells (cluster 1) progressively becoming, upon activation and migration, cells of cluster 2. These cells lie at the crossroads between two differentiation paths: One toward centroacinar-acinar phenotypes (2-4-5-6) and another toward differentiated ducts (2-3). The widespread distribution of the cells of cluster 1 across the entirety of the productal arms suggests that these cells may give rise to (or result from the de-differentiation of) all mature ductal cell types.

Fig. 3.

Ductal cells show transitional dynamism in pseudotime. (A) Monocle-generated plots presenting pseudotime ordering and differentiation trajectory of cell clusters. Color coding is as shown in Fig. 2. Clusters 4-5-6 are arranged in a manner that suggests ductal-acinar transition (bottom left arm), while activated/migrating progenitor-like cells (cluster 2) are saddled between the productal and ducto-acinar clusters. (B) Monocle-generated plots showing pseudotime-ordered expression of selected marker genes: TFF1 (cluster 2), SPP1 (cluster 1), CPA2 (cluster 6), and AKAP12 (cluster 3). The color of a cell shows its cluster identity (as in Fig. 2). Lines denote relative expression of each marker in pseudotime. (C) UMAP plots and corresponding violin plots showing relative gene expression for markers previously associated with ductal progenitors. Each dot in the violin plots represents the expression of a single cell, while the shape of the plot denotes probability density. Scale bars denote average expression in gene counts.

The notion that adult ducts harbor progenitor cells has been extensively studied, and the literature describes numerous markers for such cells, both in mice and humans. Our analysis of the ductal progenitor cell niche allowed us to map some of these against the UMAP plot. As shown in Fig. 3, several have a widespread distribution (i.e., are found in both productal and ducto-acinar axes). That is the case of CTNND1 (44) and CD24 (45). While the expression of SOX9 (27, 46) is weaker in clusters 4 and 5, it is robustly expressed in clusters that are in both the productal axis (1–3) as well as in the proacinar cluster 6. Other markers, such as CEACAM6 (47), Prominin-1/CD133 (48), and F3/CD142 (49, 50), are detected preferentially in the productal axis. Interestingly, GP2 (51, 52) is expressed almost exclusively in the ducto-acinar axis (clusters 5 and 6), which is consistent with the recent discovery that this gene marks progenitor-like de-differentiated acinar cells (49). Taken together, our data confirm the existence of two major lineage differentiation axes in the adult ALK3bright+ ductal compartment.

Alignment of the ALK3bright+ Dataset Against a Single-Cell Resolution Transcriptome of the Whole Pancreas Confirms the Ducto-Acinar Transitional Axis and Suggests an Association between Productal and Endocrine Compartments.

To ascertain the position of sorted ALK3bright+ populations in the larger context of the pancreas, we plotted our dataset against two other datasets previously reported for human islets (47, 53), which are known to contain residual amounts of ductal and acinar cells and therefore offer a rudimentary atlas of most cell types in the organ (53–55). Fig. 4 and present the combined UMAP plot, cluster resolution, heatmap, and dot-plot of select DE genes of the ALK3bright+ clusters vs. the whole-pancreas integrated dataset, respectively. The cluster color coding of our original dataset has been preserved in the integrated graphic. The fact that the original clusters are largely maintained following the integration with a whole-pancreas dataset further confirms their distinctive identity. The only difference is that, upon integration, former cluster 2 (activated/migrating progenitors) is split into two subclusters, now identified by two different shades of green, and cluster 4 (centroacinar cells) disappears. Representative genes for each cluster are indicated in Fig. 4 (heatmap in ). The contribution of the ALK3bright+ dataset (teal) is shown against the scRNAseq “map” of the whole pancreas (coral) in the circular inset of Fig. 4. This analysis confirms that ALK3bright+ cells fall under any of the two categories described above: Productal (bottom area in Fig. 4) or ducto-acinar (top right area in Fig. 4). The latter appears to seamlessly merge with the population defined as “acinar” in the whole-pancreas dataset. shows that the whole-pancreas–derived acinar cells are more acinar than the ALK3bright+-sorted “transitional to acinar” clusters (higher CPA1/AMY2A and expression of AMY2B, which is undetectable in the ALK3bright+ dataset). This, coupled with the expression of KRT19 in those clusters, further confirms their transitional nature. Importantly, no cells from the ALK3bright+ dataset correspond to INS-expressing cells (i.e., all INS+ cells derive from the whole-pancreas dataset) ().

Fig. 4.

Alignment of the ALK3bright+ clusters against a scRNA-seq dataset of the whole pancreas confirms the ducto-acinar transitional axis and suggests potential association of the productal and endocrine compartments. (A) UMAP plot showing cellular distribution when the ALK3bright+ dataset (n = 3, 4,878 cells) is mapped against two combined islet datasets claimed to represent the whole pancreas (3,277 cells). Clustering of pancreatic subtypes reflects differential gene expression and was validated by visualizing gene expression in UMAP plots (left side: endocrine; top: endothelial, immune, mesenchymal and acinar; right side: ductal). Scale bars for each marker denote average expression in gene counts. (Inset) UMAP plot showing distribution of ALK3bright+ cells (teal color) vs. the whole-pancreas dataset (coral color). (B) Single-cell arrangement based on pseudotime calculations of the integrated whole-pancreas and ALK3bright+ datasets. Cell ordering analysis results in four main arms: I (productal), II (ducto-acinar), III (ducto-endocrine), and IV (immune/mesenchymal). (C) Pseudotime-ordered ductal and acinar components of the integrated whole-pancreas and ALK3bright+ cell datasets. (Inset) UMAP plot showing the distribution of clusters based on cellular identity. (D) Pseudotime-ordered ductal and endocrine components of the integrated whole-pancreas and ALK3bright+ cell datasets. (Inset) UMAP plot showing the distribution of clusters based on cellular identity.

The Monocle arrangement of the integrated dataset is depicted in Fig. 4. Comparing this to the ALK3bright+ dataset alone, cells of cluster 2 (activated/migrating ductal progenitors) remain at the crossroads between the productal axis (Fig. 4, box I), the ducto-acinar arm (Fig. 4, box II) and what appears to be a ducto-endocrine branch (Fig. 4, box III). Box IV in Fig. 4 indicates nonepithelial cells (mesenchymal, endothelial, immune). This pseudotemporal arrangement suggests that the differentiation axes previously identified within ALK3-sorted ductal cells are not only conserved upon integration with the whole-pancreas dataset, but may also span the endocrine compartment as well. This would be consistent with the observation of a spatial continuity between ductal and β-cell clusters (56), as well as with the historical evidence of ductal-to-endocrine β-cell neogenesis (reviewed in ref. 3).

Fig. 4 show additional pseudotemporal analyses focusing specifically on ductal+acinar cells (i.e., excluding endocrine and immune/mesenchymal clusters) and ductal+endocrine cells (i.e., excluding acinar and immune/mesenchymal clusters). The first (Fig. 4) is similar to that previously shown in Fig. 3 for the sorted ALK3bright+ cells, which is not surprising, since most cells derive from that dataset. The analysis of the ductal+endocrine clusters (Fig. 4) further supports the existence of a differentiation gradient from progenitors (cluster 1, blue) to activated/migrating progenitors (cluster 2, green), which then bifurcates into GCG+ (top arm in Fig. 4) and INS+ (bottom arm in Fig. 4) arms.

Regular UMAP plots are bidimensional reductions of multidimensional data. While 2D renderings sometimes distort the relative position of clusters, this was not the case in the integrated dataset, whose 3D analysis showed maintenance of the previously determined spatial cluster distribution (Movie S1).

Sorted Human P2RY1+/ALK3bright+ Cells Exhibit Multilineage Differentiation Potential In Vivo, and This Outcome Is Enhanced upon Systemic Administration of a BMPR Agonist.

We have previously shown that sorted ALK3bright+/P2RY1+ cells exhibit a progenitor-like behavior following exposure to BMP receptor agonists in vitro (4). The replication of such outcomes in vivo may set the stage for therapies aimed at pharmacologically restoring islet mass. In an attempt to correlate our ALK3bright+ scRNA-seq findings and the previous analyses on a subset of such cells (ALK3bright+/P2RY1+), we observed that PDX1 (the marker for which P2RY1 is a surrogate) is almost exclusively expressed in cluster 1 (Fig. 5 , Left). We hypothesized earlier that this cluster harbors progenitor-like cells, characterized by high PDX1 and ID1-2 and -4 expression. Sorted ALK3bright+/P2RY1+ cells are highly depleted of BMP-7–unresponsive CAII+ cells (4) (). Intriguingly, cluster 1 had a significant number of CAII+ cells (Fig. 5 , Center). However, we found that a vast majority of the cells with high PDX1 expression did not coexpress CAII (Fig. 5 , Right). Taken together, these observations suggest that 1) only a subset of cells from cluster 1 is progenitor-like and 2) ALK3bright+/P2RY1+ sorting enriches for the progenitor-like cells (PDX1+/ALK3+/CAII—), which are preferentially found within cluster 1.

Fig. 5.

Sorted P2RY1+/ALK3bright+ cells display multilineage differentiation potential upon transplantation in immunodeficient mice. (A) The ALK3bright+/P2RY1+ fraction is enriched in PDX1+/ALK3+/CAII− cells, which map chiefly in cluster 1, as the expression of PDX1 (the marker for which P2RY1 is a surface surrogate; green) and that of CAII (red) are mutually exclusive. (B) Experimental outline for the transplantation of FACS-sorted P2RY1+/ALK3bright+ cells under the kidney capsule of immunodeficient nu/nu mice. (C) Circulating peripheral hC-pep measured at day 0 and weeks 2 and 4. (D) Circulating peripheral hC-pep measured during an IPGTT at 0 min and 60 min after glucose challenge. Each of the first three charts in C and D represents an individual experiment, with the sex of the recipient mice indicated on top. The fourth chart, in a larger size, is the combination of all of the experiments. *P < 0.05, **P < 0.01, ***P < 0.001. Error bars: SD. (E, Left) Quantification of the percentage of endocrine (INS, GCG, SST), acinar (AMY), and ductal (CK19) DAPI+ cells in the grafts of n = 3 transplanted, THR-123–treated animals. Only transplanted pancreatic cells (defined as expressing any of the above markers plus DAPI) were counted in the graft region, and values represent percentages vs. the total counted cells. (Right) Relative estimated percentages of engrafted cells according to their general lineage (endocrine, acinar, ductal). Error bars: SD. (F) Representative confocal immunofluorescence imaging of z-stack maximal projections showing expression of select markers of the grafted cells in the kidney subcapsular space of recipient mice. (i) Example of disorganized acinar (AMY, green) and ductal (CAII, red) tissue. (Scale bar, 100 µm.) (ii) Example of well-organized engrafted tissue showing fully formed acini (AMY, green), typical pancreatic epithelial architecture (E-cadherin, ECAD, red), and a small cluster of INS-producing cells (INS, gray). Details are shown in (iii). [Scale bars, 100 µm in (ii) and 50 µm in (iii).] (iv) Primitive ductal structures are shown in magenta (E-CAD, red + cytokeratin 19, CK19, white). (Scale bar, 25 µm.) (v) Endocrine clusters surrounded by ductal structures: Chromogranin A (CHGA, green); PDX1 (red); CK19, gray; and DAPI (blue). (Scale bar, 50 µm.) (vi) Example of endocrine cluster showing multiple single-hormone expressing cells: GCG (green); INS (red), SST (gray); and DAPI (blue). (Scale bar, 25 µm.) Merged images are depicted following those of single channels (n = 3 biological replicates with 6 mice per group).

Based on the above, we transplanted ALK3bright+/P2RY1+ cells to test the prediction that this subset of cells that maps in the region of cluster 1 would act as progenitors in vivo. We reasoned that this approach would be preferable to the alternative of transplanting cells sorted using yet-undefined surface markers specific for this cluster (which may or may not yield separation between CAII+ and CAII− cells). These experiments were conceived as proof-of-principle that human progenitor-like cells can differentiate into adult phenotypes in vivo, and not to model any potential therapeutic intervention. Therefore, we only transplanted subtherapeutic cell dosages in nondiabetic animals.

Transplantation of sorted P2RY1+/ALK3bright+ cells was done under the kidney capsule of immunodeficient (nu/nu) mice. These cells are hormone-negative (4) (Figs. 1 and 4 and ). The strategy is schematized in Fig. 5. Five- to 6-wk-old nu/nu mice received 5 × 105 P2RY1+/ALK3bright+ cells sorted from n = 3 donors in three independent experiments (demographics in ). Six mice were transplanted in each of the three experiments. Males were used in Exps. 1 and 3, whereas in Exp. 2 only females were used.

Although the primary endpoint was to test the hypothesis that sorted progenitor-like cells spontaneously differentiate upon transplantation, we also wanted to test whether the process could be enhanced by treating the recipient animals with a BMP receptor agonist. Thus, half the mice in each experiment (three) also received daily intraperitoneal injections (10 mg/kg of body weight) of THR-123, a small cyclic peptide with ALK3 agonism activity (5, 57). The treatment was for 10 d, starting 24 h after the transplantation. The other three animals per experiment received a similar regimen of saline injections.

Circulating human C-peptide (hC-pep) values were measured in both groups at 2 and 4 wk after transplantation. As shown in Fig. 5, there was measurable hC-pep in all three experiments at both time points, indirectly establishing the formation of INS-secreting cells from the engrafted P2RY1+/ALK3bright+ cells. Average basal hC-pep in plasma at 4 wk was ∼0.1 ng/mL. This figure is in line with hC-pep values reported at any time point prior to week 8 postprocedure when human embryonic stem cell-derived stage 4 undifferentiated progenitors are transplanted (58), despite a 10× lower cell dosage per mouse. Of note, the concentration of hC-pep in blood was higher in the THR-123 group than in saline-injected mice (P < 0.05 at week 2 and P < 0.001 at week 4). An intraperitoneal glucose tolerance test (IPGTT) conducted at day 15 revealed that engrafted cells in both groups release INS in response to glucose stimulation (Fig. 5), further supporting that INS-expressing cells are functional neogenic β-like cells. Again, while this was observed in both groups, responses were significantly stronger in those animals treated with THR-123. No statistically significant differences were observed between males and females. All animals were killed at day 16, and immunofluorescence of the grafts was performed. The quantification of the percentage of endocrine (INS+, GCG+, SST+), acinar (AMY+), and ductal (CK19+) cells was done by sampling transplanted kidney sections of n = 3 THR-123–treated (Fig. 5, red columns) and n = 3 saline-treated (Fig. 5, black columns) mice and counting pancreatic cells within the graft by FIJI/ImageJ (Fig. 5). We observed a 2.7-fold reduction in the overall number of engrafted cells in sampled sections of controls vs. THR-123–treated mice (P = 0.025). In particular, there was a 2.5-fold (P = 0.08, not significant, n.s.) increase in the number of INS+ cells detected in the THR-123 group vs. controls. These observations are aligned with our previous determination that BMP-7 increased the mitotic rate of sorted ALK3bright+/P2RY1+ cells in vitro by ∼2.5-fold (4). However, there were no statistical differences between both groups in terms of relative percentages of specific cell types (excluding SST-expressing cells, undetected in controls). Therefore, differences in circulating C-peptide are not due to differentiation biases, but probably only reflect a differential expansion of the transplanted undifferentiated cells when animals were treated with THR-123 or saline. Of note, these values represent only estimations based on the sampling of a small number of animals.

Fig. 5 shows representative images of pancreatic cell types in the grafts of THR-123–treated mice. Acinar tissue was detected in discrete regions either in a disorganized fashion, intermingled with ductal cells (Fig. 5 , i), or in well-developed acinar structures (shown in Fig. 5 , ii with an endocrine cluster) with typical epithelial morphology (Fig. 5 , iii) and small ductules (Fig. 5 , iv). Although CK19 is a marker already found in the sorted fraction (KRT19) (Fig. 1), the fact that we observed CAII+ cells when P2RY1+/ALK3bright+ cells are nearly 100% CAII− (4) (Fig. 5 and ) suggests that they have arisen by differentiation. AMY+ cells may have also differentiated along the ducto-acinar axis. Our data do not offer evidence of draining into ducts, and we cannot discard some autodigestion. However, pancreatic amylase is allosterically activated by chloride in the duodenum (59, 60). Other pancreatic enzymes also require the duodenum/intestinal mucosa microenvironment (absent in grafts) for activation. Therefore, we speculate that zymogen activation levels in the graft may be insufficient for extensive autodigestion.

Endocrine clusters were small (8 to 15 cells) and contained monohormonal cells expressing either INS, GCG, or SST (Fig. 5 , v and vi). Considering that the transplanted populations are devoid of endocrine cells (4) (Figs. 1 and 4 and ) and that the circulating hC-pep levels increase with time (Fig. 5), we hypothesize that these clusters also arise by differentiation. Of note, while we did not transplant the negative fraction as a negative control, our previous experiments strongly suggest that it is strongly depleted in mature pancreatic markers, and completely devoid of progenitor-like cells (4). In summary, sorted progenitor-like cells transplanted into immunodeficient animals differentiate along all pancreatic lineages.

PDX1+/ALK3+/CAII− Progenitor-like Cells Exist in Pancreata of T1D/T2D Donors Regardless of the Duration of the Disease.

If progenitor-like cells were unaffected by T1D autoimmunity, such discovery may open the door to therapies aimed at their in situ activation for regenerative purposes. We have determined that the P2RY1+/ALK3bright+ fraction is highly enriched in cells with a PDX1+/ALK3+/CAII− phenotype (4) (), which is quantifiable by immunofluorescence in human samples. To determine if such cells are present in ducts after disease onset, we examined the expression of PDX1, ALK3, and CAII in short- (0 to 5 y), medium- (6 to 10 y), and long- (>10 y) term T1D donor samples. Nondiabetic controls were used for comparison; n = 3 donors per group, plus one T2D donor, were selected for analysis (demographics in ). For each donor, three histological sections containing MPDs (characterized by a fibromuscular layer) were used for each of the three anatomical regions of the pancreas (head, body, and tail).

Immunofluorescence of PDX1 in paraffin-embedded nPOD (Network for Pancreatic Organ Donors with Diabetes) pancreatic samples is challenging and has not been reported thus far. Therefore, we optimized the protocol using control and T1D sections that encompassed both islets and ducts in the same field. We found that PDX1 signal was very abundant not only in islets (as expected), but also throughout the entirety of the ductal tree (Fig. 6 and ). shows gradual loss of INS, PDX1, and P2RY1 in islets as T1D progresses. However, PDX1 signal remains strong in ductal structures regardless of the stage of the disease. This protocol is described in .

Fig. 6.

PDX1+/ALK3+/CAII− cells are present in the large ducts of T1D/T2D donors. (A) Representative confocal immunofluorescence imaging of a z-stack maximal projection showing expression of ALK3 (red), CAII (green), PDX1 (gray), and DAPI (blue) from human donors. (Left) Nondiabetic (Upper) and short-term T1D (Lower). (Right) Medium-term (Upper) and long-term (Lower) T1D. White arrowheads indicate cells with a PDX1+/ALK3+/CAII− phenotype. (B) Quantification of PDX1+/ALK3+/CAII− cells calculated as a percentage of overall PDX1+/ALK3+ ductal cells in large human pancreatic ducts across the head, body, and tail of the pancreas, for control, short-term, medium-term, and long-term diabetic donors. Error bars: SE. No statistically significant differences were detected between groups. (C) Representative confocal immunofluorescence imaging of a z-stack maximal projection showing the presence of PDX1+/ALK3+/CAII− cells in a T2D donor. (Scale bars, 200 µm for low-magnification and 50 µm for high-magnification images.)

We also detected nuclear PDX1 in mature ductal cells that also expressed CAII (Fig. 6 ). While scRNA-seq indicated that PDX1 expression is restricted to a small proportion of ALK3+/CAII− cells, predominantly within cluster 1, this analysis was done on sorted ALK3bright+ cells, which represent only a fraction of all ALK3+ cells in the ductal tree. In fact, this analysis confirmed our previous observation (4) that PDX1+/CAII+ are generally ALK3dim+, whereas PDX1+/CAII− cells are invariably ALK3bright+ (Fig. 6). This also explains why the ALK3bright+/P2RY1+ sorting conducted for transplantation results in the depletion of most CAII+ cells (4) (Fig. 5 and ).

PDX1+/ALK3+/CAII− ductal cells, often in an intercalated pattern with PDX1+/ALK3+/CAII cells, were detected in head, body, and tail of all donors, regardless of the duration of the disease (Fig. 6). The percentage of PDX1+/ALK3+ cells vs. DAPI, and that of CAII− and CAII+ cells within PDX1+/ALK3+ cells, was calculated by FIJI/ImageJ for the epithelial lining of all of the major ducts found in each section (Fig. 6). The overall percentage of PDX1+/ALK3+/CAII− cells ranged between 5% and 12% of all PDX1+/ALK3+ ductal cells examined (the remaining being CAII+). There were no statistically significant differences between control, short-, medium-, or long-term samples. No meaningful differences were observed either when comparing the different regions of the pancreas within each donor. These cells were also detected in a sample from a T2D donor (Fig. 6), but no quantification was attempted. Taken together, these results indicate that cells with a progenitor-like immunofluorescence signature (PDX1+/ALK3+/CAII−) exist in the ductal tree of T1D patients regardless of the duration of the disease. It remains to be established whether these cells exhibit the same multipotency previously established for nondiabetic donors.

Discussion

We present the high-resolution characterization of a human pancreatic ductal progenitor cell niche. All previous analyses of the human pancreas by scRNA-seq have been performed on isolated islets (53–56, 61), which contain remnants of the acini and ducts. However, the number of sequenced ductal cells was invariably low, and probably depleted of the major ductal cells that are the subject of our study. Furthermore, the ductal identity of any given cell was determined indirectly by bioinformatics a posteriori. Our study design circumvents these limitations by focusing on a specific lineage of the pancreas and prevalidating its identity by cell sorting.

scRNA-seq of a sorted ductal population previously shown to harbor cells with progenitor-like characteristics (4) revealed an unexpected degree of heterogeneity, with cells spanning a broad range of differentiation stages across two major axes, one encompassing mature ductal cell types and the other showing progressive acinarity. The existence of a gradient of expression of acinar genes in clusters 5 and 6 rules out the possibility of ambient contamination of acinar RNA during the procedure (which would be equally detected in all clusters), and suggests an active transition between the cells of these clusters. Acinar-to-ductal metaplasia is a phenomenon by which, in response to inflammation or stress, acinar cells de-differentiate into progenitor-like (metaplastic) ductal cells. Such cells are then endowed with the capacity to regenerate the damaged organ (62–64). This process only becomes irreversible when cells acquire oncogenic Kras mutations or are subjected to persistent/aberrant signaling, which prevent redifferentiation and initiate further progression into pancreatic adenocarcinoma (65). We believe that the ducto-acinar ALK3bright+ axis described herein represents cells in the process of transitioning between these two lineages, which suggests that this reversible transdifferentiation may be more common than previously thought. Consistent with this hypothesis is the observation that GP2, an acinar-to-ductal transitional marker (49), is almost exclusively detected in the ducto-acinar clusters 5 and 6 (Fig. 3).

Cluster 1 cells displayed strong inflammation/stress responses. While dissociation and sorting may induce stress, it is doubtful that this is the reason behind our observation. If that were the case, we would also probably observe it in the other clusters, which have also undergone the same process. Stress in the pancreatic microenvironment is a well-studied phenomenon. The exocrine pancreas is the organ with the highest protein synthesis, producing liters of pancreatic juices daily. Leaving aside the potential risk of intracellular activation of digestive enzymes (the main cause of pancreatitis), protein synthesis itself is a key stressor of the exocrine pancreas at all times (66). Stress has been associated with de-differentiation and acquisition of progenitor-like characteristics in ducts (67–69), but whether the stress signature of cells in cluster 1 is linked to de-differentiation is merely presented as a hypothesis.

This cluster also had the highest levels of PDX1 and ID1-2, and -4, known housekeepers of progenitor-like phenotypes (4, 17). Of note, while P2RY1+ cells were also found throughout the dataset, their scarcity made it difficult to analyze them. This is because P2RY1 expression is very low. As shown in , the relative level of expression (2−ΔCt) of P2RY1 was 5-fold lower than that of PDX1 and 60-fold lower than that of KRT19. Our findings are confirmed by the entry for pancreas in the Human Protein Atlas, which shows an average of transcripts per million across 248 samples to be 0.8 (i.e., there is a probability of finding 0.8 transcripts of P2RY1 per one million reads). As reference, PDX1 has 6.2 and INS has 3535.7. As recently reviewed in ref. 70, scRNA-seq cannot be used to draw meaningful conclusions from low-expressing genes.

Interestingly, albeit relatively widespread throughout all ALK3bright+ clusters, the expression of the progenitor marker CTNND1 (p120ctn) was highest in cluster 1. During development, cells with high p120ctn expression define “trunk” (major) ductal populations that give rise to both ducts and endocrine cells, whereas p120ctnlow cells, typically located in the tip of the ducts, only differentiate into acini (44). Similarly, cluster 1 had the highest expression of CEACAM6, which was independently predicted to be a progenitor marker in the pancreas (47); as well as PROM1, a marker of tripotent colony-forming progenitor cells in the adult murine pancreas (48). Enrichment in progenitor-like cells from cluster 1 was predicted to occur by P2RY1+/ALK3bright+ sorting, since P2RY1 (PDX1) selection mostly excludes the nonprogenitor CAII+ fraction therefrom.

We did not expect to detect ductal-to-endocrine transitional cells in our dataset due to the very design of the experiment, where endocrine cells (ALK3−) are specifically excluded by gating on the ALK3bright+ fraction. Our study is based on a static snapshot of the pancreas, and therefore conclusions about potential transitions are necessarily speculative and based on bioinformatics predictions. Even so, the integration with the whole-pancreas datasets, and especially the pseudotemporal arrangement analysis, strongly suggests the existence of a ducto-endocrine axis. In this particular setting, where none of the donors had diabetes and sources of stress appear to be largely physiological, it is plausible that there may be a basal level of endocrine cell formation from de-differentiated ductal cells, even if the full activation of ductal-mediated regenerative mechanisms may require more extensive damage (3, 71).

Transplantation is stressful for dissociated pancreatic epithelial cells (4). In this context, we further hypothesized that xenotransplanting sorted P2RY1+/ALK3bright+cells may offer us a glimpse of their potency under stress in vivo. We show here that these cells, which are hormone-negative at the time of transplantation, can differentiate into monohormonal GCG-, SST-, and INS-expressing cells, as well as cells of the acinar and ductal lineages. While some regions of the graft did not exhibit mature histological patterns, the organization of other regions in acini, small endocrine clusters, and ductules was remarkable, considering that only single cells sorted for ALK3 and P2RY1 (none of which is expressed in human acinar tissue) were transplanted. Our ultimate goal is not transplantation, but rather the activation of progenitor cells in situ. As indicated earlier, these experiments were only designed to establish proof-of-concept of in vivo multilineage differentiation, and not to revert diabetes in a disease model. Therefore, we did not conduct any comparison with islets or extensive functional studies. Still, we show that engrafted INS+ cells secreted INS in response to glucose stimulation in vivo, which is the hallmark of functional β-like cells. This process was enhanced when the animals were treated systemically with THR-123, an ALK3 receptor agonist (57). This finding, coupled with our discovery that PDX1+/ALK3+/CAII− cells (a marker combination signature equivalent to P2RY1+/ALK3bright+) (4) are present in T1D/T2D donors, suggests potential therapeutic avenues to induce β-cell regeneration in situ. To determine if that is a possibility, future studies will explore whether cells sorted from T1D/T2D donors also exhibit multipotency.

The existence of progenitor cells with the capacity to regenerate the endocrine and exocrine compartments in the adult human pancreas remains a contentious issue. However, the recent emergence of high-resolution analytical tools is bringing much needed nuance to the debate. Lineage barriers are now blurrier than ever, a discovery that calls into question the validity of many earlier lineage-tracing conclusions. The high degree of developmental heterogeneity supports the dynamic fate plasticity and the constant state of flux of the ductal compartment. If harnessed, these findings could be the basis for the development of pharmacological therapies for endocrine regeneration.

Materials and Methods

Human pancreata were processed as detailed in the . Methods for FACS of cells for scRNA-seq analysis and transplantation are also expanded therein. Animal experiments were conducted under the oversight of the University of Miami Institutional Animal Care and Use Committee. Procedures for transplantation and metabolic analyses of in vivo differentiation as well as immunofluorescence analyses are detailed in the . Single-cell capture, library preparation, and sequencing were performed at Genewiz (South Plainfield, NJ). Details on library preparation, computational analyses, and statistics are provided in .

Data Availability.

Sequencing data have been deposited in the Gene Expression Omnibus (GEO), accession no. GSE131886. Previously published scRNA-seq data reanalyzed here are available at under GEO accession nos. GSE81076/GSE85241. Source data for figures are provided in the reagent table. Code/Rscript files for the analyses reported are available in the GitHub repository https://github.com/JDBLab/Pancreas_ductal_scRNAseq (72).