Single-cell transcriptomes identify human islet cell signatures and reveal cell-type–specific expression changes in type 2 diabetes

1Pancreatic islets of Langerhans are clusters of at least four different hormone-secreting endocrine cell types that elicit coordinated—but distinct—responses to maintain glucose homeostasis. As such, they are central to diabetes pathophysiology. On average, human islets consist mostly of beta (54%), alpha (35%), and delta (11%) cells; up to a few percent gamma/pancreatic polypeptide (PP) cells; and very few epsilon cells. Human islet composition is neither uniform nor static but varies between individuals and across regions of the pancreas. Cellular heterogeneity complicates molecular studies of whole human islets and may mask important role(s) for less common cells in the population. Moreover, it complicates attempts to identify epigenetic and transcriptional signatures distinguishing diabetic from nondiabetic (ND) islets, leading to inconsistent reports of genes and pathways affected. Conventional sorting and enrichment techniques are unable to specifically purify each human islet cell type, thus a precise understanding of the transcriptional repertoire governing each cell type's identity and function is lacking. Identifying the cell-type–specific expression programs that contribute to islet dysfunction and type 2 diabetes (T2D) should reveal novel targets and approaches to prevent, monitor, and treat T2D.

2In this study, we sought to decipher the transcriptional repertoire of each islet cell type in an agnostic and precise manner by capturing and profiling pancreatic single cells from ND and T2D individuals. From these profiles, we identified transcripts uniquely important for each islet cell type's identity and function. Finally, we compared T2D and ND individuals to identify islet cell-type–specific expression changes that were otherwise masked by islet cellular heterogeneity. The insights and data from this study provide an important foundation to guide future genomics-based interrogation of islet dysfunction and diabetes.

3Pancreatic islets (>85% purity and >90% viability) were obtained from eight human cadaveric organ donors (five ND, three T2D) (Fig. 1A; Supplemental Table S1). Each islet sample was processed to generate single-cell RNA-seq libraries (Fig. 1A; single cell) and paired bulk RNA-seq libraries at three different stages of islet processing (Fig. 1A; baseline, intact, and dissociated). All RNA-seq methods employed SMARTer chemistry (Methods), and bulk islet cDNA libraries were sequenced to an average approximate depth of 34 million reads (Supplemental Table S2). Baseline, intact, and dissociated transcriptomes from each person were highly correlated (Supplemental Fig. S1). Transcriptomes clustered by donor and not by processing condition or incubation time (Fig. 1B), strongly suggesting that islet processing did not significantly alter islet transcriptomes.

4A total of 1050 islet cells (622 ND and 428 T2D) were captured on 11 Fluidigm C1 chips. cDNA libraries were constructed from captured cells and barcoded, fragmented, pooled, and sequenced to an average depth of 3 million reads (Supplemental Table S2). Two separate library preparations from the same amplified cDNA for each of 83 single cells demonstrated remarkable correlation, suggesting minimal batch effects resulting from the cDNA processing and sequencing steps. Resequenced samples are highlighted in Supplemental Table S2 but were not included in subsequent analyses. Transcript coverage is indicated in Supplemental Figure S2. Approximately 81% (21,484/26,616) of protein-coding genes and long intergenic noncoding RNAs (lincRNAs) were detected in at least one cell from the collection. On average, each single cell expressed 5944 genes (Fig. 1C). Cells expressing less than 3500 genes (n = 72) also exhibited high mitochondrial alignment rates and other reported transcriptional metrics of cell death and/or poor quality and were removed from subsequent analyses (Fig. 1C).

5We next assessed the extent to which the remaining 978 single-cell transcriptomes represent the expression patterns observed in intact islets. Single-cell transcriptome ensembles from each person were highly correlated (Pearson's R2 ranged from 0.91–0.98) (Supplemental Fig. S3), regardless of disease state. Pearson's R2 values between individuals’ single-cell ensembles and corresponding “bulk” transcriptomes ranged from 0.75–0.86 (Supplemental Fig. S4) and did not differ substantially between ND (R2 = 0.87) and T2D (R2 = 0.85) samples (Supplemental Fig. S5). Overall, ensemble/aggregate single-cell transcriptome profiles correlated well with those of pooled bulk islet transcriptomes from all individuals (Fig. 1D, R2 = 0.87). These results suggest that the islet single-cell transcriptomes are high quality and effectively reflect bulk islet transcriptomes.

6Five islet endocrine cell types have been assigned based on exclusive and robust expression of the peptide hormone genes INS (beta), GCG (alpha), SST (delta), PPY (PP/gamma), and GHRL (epsilon). The pancreas also contains three exocrine cell types—acinar, stellate, and ductal—that critically support digestion through synthesis and transport of digestive enzymes. Each also has been identified by specific marker gene expression, including the serine peptidase gene PRSS1 (acinar), the extracellular matrix protein gene COL1A1 (stellate), and the structural keratin gene KRT19 (ductal). We used these marker genes to determine the representation of each islet cell type among our 978 profiled single cells.

7Density plots (Fig. 2A) revealed bimodal expression of each marker gene across the population of single cells. Therefore, we employed Gaussian mixture modeling (GMM) to classify the cells unambiguously (Fig. 2B). Approximate log2 counts per million (CPM) thresholds for each marker gene used to classify cell types are provided in Supplemental Table S3. This approach identified 617 single cells (∼63%) from T2D and ND islets expressing a single marker gene representative of each major endocrine and exocrine cell type, examples of which are shown in Figure 2C. This included 239 alpha, 264 beta, 25 delta, and 18 PP/gamma cells (Table 1); the proportions of each cell type are in the ranges previously reported. Only one cell expressing high levels (log2CPM > 15) of GHRL was identified, which we presume to be an exceedingly rare epsilon cell. Additionally, we obtained 19 stellate, 24 acinar, and 27 ductal cells (Table 1), presumably from exocrine contamination of the islet cell preparations. Only 21 cells (∼2%) expressed none of the specified marker genes (Table 1). Approximately one-third (340/978) of cells expressed more than one marker gene; these were removed from subsequent analysis due to concerns that these represent two vertically stacked cells in a given capture site (for details, see Methods). Similar ratios of potential stacked cells have been reported in other studies using the Fluidigm C1 platform to capture mouse and human islet cells. Collectively, these single-cell data capture transcriptome profiles representing each of the major and minor pancreatic endocrine and exocrine cell types. Genome Browser tracks representing aggregate single-cell expression for each islet cell type have been generated using HOMER and are made available (see Data Access) to facilitate their use and investigation by members of the islet biology and diabetes research communities.

8To determine if and how islet cell transcriptomes cluster, we completed unsupervised dimensionality reduction via t-distributed stochastic neighbor embedding (t-SNE) on 380 ND single-cell samples (excluding “multiple” labeled samples). t-SNE assembled single-cell transcriptomes into discrete clusters based upon 1824 highly expressed genes (see Methods; Supplemental Table S4); GMM-based marker gene analysis revealed that each cluster corresponded to a distinct endocrine and exocrine cell type (Fig. 3A; Supplemental Fig. S6). Unsupervised hierarchical clustering also grouped single-cell transcriptomes into discrete cell types (Fig. 3B). Despite being obtained from different individuals, 161/168 beta, 128/138 alpha, 15/16 delta, and 12/12 PP/gamma cell transcriptomes clustered into the same dendrogram branches, strongly suggesting that cell type encodes the greatest variation in the data. Exocrine cells and those expressing none of the specified marker genes (“none”) clustered separately from the endocrine cell types. Importantly, this clustering was driven by neither sequencing depth (Supplemental Fig. S7B) nor expression of classic marker genes (INS, GCG, SST, PPY, GHRL, COL1A1, PRSS1, and KRT19), as cells continued to cluster into discrete cell types even when all marker genes were removed from the expression data sets (Supplemental Figs. S7C, S8). Recent studies have reported heterogeneity among beta cells. Specifically, Dorrell et al. characterized four subpopulations of human beta cells based on differing ST8SIA1 and CD9 expression. Similarly, distinguished two populations of proliferating (Fltp+) and mature (Fltp−) mouse beta cells. We did not find evidence of beta cell subpopulations (Supplemental Fig. S9), nor did we identify numerous proliferating cells (Supplemental Table S5). T2D single-cell transcriptomes (n = 258) also demonstrated clear clustering by cell type in unsupervised analyses (Supplemental Figs. S10–S14) based on 1908 highly expressed genes (Supplemental Table S4). Thus, each endocrine and exocrine pancreatic cell type exhibits a complex characteristic expression signature that uniquely identifies it.

9To identify gene signatures distinguishing each islet cell type, we completed a series of pairwise differential expression analyses (Supplemental Table S6) between each cell type (see Methods). After intersecting the results from each pairwise comparison, we identified a conservative collection of 154 islet endocrine cell-type “signature” genes (61 beta, 51 alpha, 17 delta, 25 gamma), as well as 202 exocrine genes (109 stellate, 31 acinar, 62 ductal) at 5% false-discovery rate (FDR) (Fig. 3C; Supplemental Table S7). Two genes exhibited overlap between the endocrine and exocrine signature lists: FAP (alpha and stellate cell overlap) and TNS1 (beta and stellate cell overlap). Gene set enrichment analysis (GSEA) identified enrichment (FDR-adjusted P-value <0.05) of insulin signaling, oxidative phosphorylation, maturity-onset diabetes of the young (MODY), and glycolysis/gluconeogenesis KEGG pathways in beta cells relative to the other endocrine cells (Supplemental Table S8).

10Signature genes included previously reported beta-specific genes like NKX6-1, DLK1, and ADCYAP1 (Fig. 3C, right) and alpha cell–specific genes like IRX2, LOXL4, and DPP4, a cell surface receptor and diabetes drug target. Among delta cell signature genes, we detected exclusive expression of HHEX, a transcription factor reported to govern delta cell identity and function and linked to T2D GWAS. Delta cells also specifically expressed BCHE, which encodes butyrylcholinesterase. BCHE catalyzes the breakdown of acetylcholine and ghrelin, thus providing a mechanism for delta cells to exert local inhibition of islet-influencing endocrine signals. PP/gamma cell–specific transcriptomes included CTD-2008P7.8, a lincRNA of unknown function; CNTNAP5, a member of the neurexin family of cell adhesion molecules; and ID4, which encodes an inhibitor of DNA-binding protein. In addition to DPP4, we detected 30 islet signature genes whose proteins SWISSPROT predicts to localize to the cell surface (Supplemental Table S9). DPP4 antibodies have recently been used to isolate purer alpha cell populations from islets. Thus, antibodies against the other candidate cell-type–specific surface markers we have identified may be useful to purify other islet cell types.

11Cell sorting and enrichment methods such as fluorescence-activated cell sorting (FACS) have been used to identify characteristic alpha and beta cell genes. However, expression of SST or PPY in the reported alpha and beta cell gene sets suggests the presence of the less abundant delta and PP/gamma islet cell types in the enriched cell preparations. To distinguish genes exhibiting alpha- and beta-specific gene expression from those expressed also in delta and PP/gamma cells, we investigated the expression of previously reported alpha- and beta-specific genes (Supplemental Table S10; Supplemental Fig. S15) in our ND endocrine single-cell transcriptomes. Only 115/1683 previously reported beta-specific genes were expressed greater than fourfold higher in beta cells relative to the other endocrine cells (FDR < 0.05; one-way ANOVA followed by Tukey's honest significant difference [THSD]) (Fig. 3D). Similarly, 75/1853 reported alpha-specific genes were alpha cell enriched (Fig. 3E). Several genes previously reported to be enriched in the major islet cell types, such as MAFA, SLC2A2, SIX3, and DLK1 in beta cells and IRX2, DPP4, and ADORA2A in alpha cells, were confirmed to be signature genes. Surprisingly, we found that 37 and 33 reported beta- and alpha-specific genes were also expressed in delta and PP/gamma cells, respectively (Fig. 3F; Supplemental Table S10). Notable examples included beta and delta cell expression of the congenital hyperinsulinemia (CHI) gene HADH and alpha and PP/gamma cell expression of the ARX transcription factor. HADH is typically associated with beta cell expression and, when mutated, leads to insulin hypersecretion and CHI; these data implicate the delta cell in the molecular genetics of CHI. Misexpression of ARX has been shown to convey both alpha and PP/gamma cell features to cells, suggesting that its expression in each cell type is important for identity and function.

12Perturbations in genes involved in glucose sensing and proper maintenance of blood glucose levels contribute to T2D pathophysiology. Beta cells regulate blood glucose through the secretion of insulin and are thus exquisitely sensitive to blood glucose levels. Glucose-stimulated insulin secretion (GSIS) is linked to universal basic pathways of cellular metabolism in beta cells. To gain insight into beta cell-type–specific transcriptomic features associated with GSIS, namely, glucose uptake and glycolysis, we examined the expression of relevant genes in our islet single-cell transcriptomes (Fig. 4A).

13GSIS pathway genes associated with glucose sensing and uptake displayed highly beta cell–specific expression, including SLC2A2, which encodes the glucose transporter GLUT2; G6PC2, which encodes a subunit of glucose-6-phosphatase; and PFKFB2, which encodes an enzyme involved in regulation of glycolysis (Fig. 4A;; ). While expressed in all cell types, the enzyme, ALDOA1, immediately downstream from PFK1 and associated with the glycerol phosphate (GP) shuttle, is enriched in beta cells, perhaps reflecting an additional point of GSIS control. Protein-coding genes for five subsequent glycolytic enzymatic steps from glyceraldehyde-3-phosphate to pyruvate were not significantly differentially expressed between cell types. Beta cells are known to be limited in their ability to produce lactate from pyruvate; this is reflected by high LDHB/LDHA ratios that favor the lactate to pyruvate flux in beta cells.

14The glycerol-3-phosphate shuttle allows NAD+ regeneration in the cytosol to sustain glycolytic flux essential for GSIS. Cytoplasmic NAD+ generation has been shown to be essential for GSIS. Both components of the glycerol-3-phosphate shuttle, cytoplasmic GPD1 and mitochondrial GPD2, were expressed in beta cells, with the former representing a beta cell signature gene (Fig. 4A). Additionally, we identified the mitochondrial solute transporter SLC25A34 as beta cell specific. While its transport specificities have yet to be determined, the closest yeast ortholog of SLC25A34, Oac1p/YKL120w, is thought to import oxaloacetate into the mitochondria. This is particularly intriguing considering our data and others show the complete absence of pyruvate carboxylase (PC) expression in human beta cells, despite the essential role PC is known to play in rodent GSIS through mitochondrial production of oxaloacetate. We hypothesize that SLC25A34 may provide an alternate, cytoplasmic source for mitochondrial oxaloacetate in the human beta cell.

15Single-cell profiling also allowed us to interrogate the transcriptional repertoire of less abundant delta and PP/gamma cell types, which have been elusive in both whole islet and sorted islet studies. While it is difficult to determine epsilon cell expression signatures with one ghrelin-positive cell, our ND data set includes 16 delta cells and 12 PP/gamma cells. Among the top 100 differentially expressed (FDR < 5%) genes in delta versus other islet endocrine cells are receptors for the appetite-regulating hormones leptin (LEPR) and ghrelin (GHSR), the growth factor neuregulin 4 (ERBB4), and the neurotransmitter dopamine (DRD2) (Fig. 4B). GHSR has recently been shown to be specifically expressed and functional in both human and mouse delta cells, reducing GSIS in human and mouse beta cells when induced. LEPR, DRD2, and ERBB4 expression is specific to human delta cells. In situ analyses (ViewRNA, Affymetrix) detected coexpression of LEPR in 79/102 (77%) of SST-expressing cells (Fig. 4D, arrowheads) in ND islets, confirming the delta cell–specific expression detected in Fluidigm C1 profiling. Thus, our data suggest intriguing roles for islet delta cells in the integration of metabolic signals via leptin, ghrelin, and dopamine signaling pathways.

16PP/gamma, along with epsilon cells, are among the least studied islet cell types due to their scarcity in islets. Recent studies show that PP/gamma cells are crucial regulators of energy homeostasis. In response to food intake, these cells secrete the anorexigenic hormone PPY to facilitate vagal stimulation of neuropeptide Y receptors in the hypothalamus and induce satiety. Our data suggest interesting parallels in expression between PP/gamma cells and serotonergic neurons, a group of neurons that influence various cognitive and physiological processes including anxiety, mood, sleep, and satiety. We report expression of FEV, a serotonergic transcription factor and necessary driver of neuronal maturation previously reported in mouse beta cells, in PP/gamma cells (average log2CPM of 2.172). Interestingly, FEV has also been implicated in beta cell differentiation, and Fev−/− mice exhibit insulin production, insulin secretion, and glucose clearance defects. Other related signature genes in PP/gamma cells include TPH1, encoding a tryptophan hydroxylase essential for the initial catalysis of serotonin, and SLC6A4, a serotonin reuptake transporter. Serotonin colocalizes with insulin in beta cells and promotes GSIS. Mice lacking TPH1 are diabetic and exhibit impaired insulin secretion due to a lack of pancreatic serotonin. Elevated FEV, TPH1, and SLC6A4 expression suggests PP/gamma cells share a suite of characteristic genes with serotonergic neurons that, in the pancreas, integrate central and peripheral hunger and satiety cues. We also observed high PP/gamma expression of muscarinic acetylcholine receptor M3, CHRM3, which stimulates exocrine pancreatic amylase, insulin secretion, and smooth muscle contraction and gastric emptying. These data implicate the less abundant delta and PP/gamma cell types as critical for islet function via the integration of systemic cues and warrant further studies to elucidate the function and health of these cells in normal and diabetogenic conditions.

17We next sought to understand the cell type(s) involved in rare forms of diabetes, including transient/permanent neonatal diabetes (T/PNDM), CHI and MODY, as well as more common forms of islet dysfunction and diabetes (T1D/T2D). Monogenic diabetic disorders, including CHI, MODY, and neonatal diabetes, are characterized by mutations in a single gene, often resulting in beta cell dysfunction and death. Five monogenic diabetes risk genes (Supplemental Table S11;;; ), were enriched in beta cells (i.e., greater than fourfold change in expression in specific islet cell type relative to other endocrine cells), including glucose transporter SLC2A2 (data not shown), beta cell maturation transcription factor PDX1, and the sulfonylurea drug target ABCC8 (Fig. 4C). PDX1 expression has been reported in human and mouse beta and delta cells. Despite the modest number of delta cells sampled, our data also suggest moderate PDX1 expression in delta cells (four of 16 delta cells with expression ≥16 CPM). Robust expression of HADH in both beta and delta cells (Fig. 4C) was confirmed by in situ (View RNA) analyses (Fig. 4E,F). Approximately 386/457 cells (84%) in HADH and INS labeled sections coexpressed both markers (shown in Fig. 4E). Adjacent SST/HADH colabeling yielded an approximately equal proportion (n = 255/306; 83%) of SST-negative/HADH-positive cells. Finally, 43/457 (9%) cells were INS negative/HADH positive, and 41/306 (13%) cells coexpressed SST and HADH (shown in Fig. 4F) in the respective serial sections. Another CHI-associated gene, UCP2, which was reported to be highly expressed in human beta cells and to suppress insulin secretion, was enriched in delta cells (Fig. 4C). Delta cell expression of monogenic diabetes genes thus implicate this cell type in the molecular genetics of rare islet dysfunction and diabetes disorders, particularly CHI.

18We also investigated cell type expression patterns of 536 islet expression quantitative trait loci (eQTL) target genes. The majority of these genes (n = 309; Supplemental Table S11) were lowly expressed in both the endocrine islet single-cell transcriptomes and in the paired bulk islet transcriptomes (Supplemental Fig. S16A). One hundred fifty-nine additional genes did not exhibit a greater than or equal to fourfold expression change in any endocrine islet cell type. Of the remaining 68 eQTL genes, 54, 46, 51, and 43 were expressed in beta, alpha, delta, and PP/gamma cells, respectively. Surprisingly, beta and delta cells possessed the highest numbers of cell-type–specific eQTL genes (Supplemental Table S11).

19Genome-wide association studies (GWAS) have identified more than 100 loci associated with T2D and related quantitative traits. Because GWAS identify genetic variants associated with a disease, but not the specific gene(s) affected, we took two approaches to assess cell-type expression of patterns of putative GWAS genes. First, we compiled and examined a list of 197 reported putative T1D and T2D GWAS genes. Of these genes, 37 were expressed in beta, 24 in alpha, 28 in delta, and 22 in PP/gamma cells (Supplemental Table S11). Similarly, genes that were cell-type specific were expressed at higher levels in ND bulk intact islets compared with those genes without cell-type specificity (Supplemental Fig. S16B). Ten genes were uniquely expressed in beta cells, including MEG3, a type 1 diabetes (T1D)–associated lincRNA with reported expression in mouse beta cells and potential tumor suppressor activity, and IAPP, whose protein product, when aggregated, possesses cytotoxic properties that may contribute to beta cell death and dysfunction in T2D. We also identified five putative T2D GWAS genes (including HHEX) to be uniquely expressed in delta cells. To conduct a more liberal analysis of putative GWAS genes, we identified all single-nucleotide polymorphisms (SNPs) associated with polygenic diabetes and related traits from the GWAS catalog (https://www.ebi.ac.uk/gwas/). For each reported SNP associated with T2D, T1D, fasting insulin, fasting glucose, and proinsulin, we examined the expression of all genes overlapping within one megabase of the chromosomal locus and identified 263 genes with cell-type–specific expression (Supplemental Table S12). Together, our observations of cell-type–specific expression of eQTL and monogenic and common (T2D GWAS) diabetes genes both confirm beta cell–specific expression of multiple diabetes-associated genes (MEG3, DLK1, SLC2A2, etc.) and implicate other cell types in the molecular genetic pathogenesis of diabetes. In light of recent studies and our data, which suggest that delta cells may be critical regulators of glucose homeostasis and islet function, this provides a new avenue for investigation of T2D pathogenesis, as well as potentially new therapeutic targets and treatment options.

20Finally, we compared single-cell transcriptome profiles from T2D and ND donors to identify differentially regulated genes and obtain greater insight into the molecular genetic pathogenesis of diabetes. After unsupervised hierarchical clustering (Fig. 5A) and t-SNE analysis (Supplemental Figs. S17, S18) using 2754 of the most highly expressed genes (Supplemental Table S4), we observed that transcriptomes clustered by cell type regardless of disease state. As previously observed, clustering was not driven by marker gene expression (Supplemental Figs. S19, S20). For regions of the dendrogram (Fig. 5A) where samples appeared to cluster by disease state, we found that islet donor identity was an underlying factor that reflected sample subclustering (Supplemental Fig. S21). We obtained fewer beta cells among the T2D islet cells sampled compared with ND samples (Fig. 5B). However, observed differences in T2D and ND single-cell proportions did not differ significantly from expected cell-type proportions (Fig. 5B, χ2 P-value = 0.2733), and none of the islets from these newly diagnosed T2D individuals exhibited as significant a decrease as previously reported.

21Recent studies have reported features of beta cell de-differentiation under diabetogenic and stress conditions. However, we did not identify significant shifts in islet cell populations, increases in number of hormone-negative “none” cells, or appearances of new or more abundant populations of cells in T2D islets that clustered distinctly from the known islet cell types in this study. Moreover, expression of reported de-differentiation genes including FOXO1, NANOG, and POU5F1 did not differ significantly between T2D and ND islet cell types nor the paired bulk intact islet preparations (Supplemental Fig. S22). Finally, other de-differentiation markers such as NEUROG3 and MYCL were not detected in our single-cell or bulk intact islet data. Thus, our analysis did not identify transcriptional evidence of de-differentiated cells in T2D islets.

22Comparison of islet cell-type transcriptomes (e.g., T2D beta vs. ND beta) did, however, identify 410 genes that were differentially expressed (FDR < 5%) between T2D and ND donors (Supplemental Table S6) beta, (Fig. 5C, n = 248), alpha (Fig. 5D, n = 138), and delta cells (Fig. 5E, n = 24). We also identified differentially expressed genes in acinar (n = 74), ductal (n = 35), and stellate (n = 28) exocrine cell types (Supplemental Fig. S23; Supplemental Table S6). T2D beta cells exhibited a 1.4-fold decreased INS expression compared with ND beta cells (Fig. 5C). STX1A was significantly reduced (log2FC −1.5178) in T2D beta cells, consistent with reported decreases in STX1A protein levels in T2D beta cells. STX1A combines with SNAP-25 and VAMP2 to form a tertiary SNARE protein complex important for insulin secretion in beta cells, and STX1A inhibition drastically reduces GSIS and exocytosis. Additionally, we detected elevated DLK1 expression in T2D beta cells (log2FC 2.010), which has been implicated in T1D/T2D GWAS and is part of a dysregulated locus in T2D islets. Dlk1−/− mice exhibit increased glucose sensitivity and insulin secretion, and high levels of serum DLK1 have been associated with insulin resistance in both rodents and humans. Immunofluorescence indicates that DLK1 is beta cell specific in human but not mouse islets, and FACS-enriched mouse beta cells show low expression of Dlk1 in comparison to other sorted islet alpha and delta cells, potentially implicating a unique role of this gene in human T2D progression. These findings suggest that perturbations in STX1A and DLK1 expression may contribute to the beta cell dysfunction and impaired insulin secretion that is commonly observed in T2D pathogenesis.

23Decreased beta cell function and mass are hallmarks of T2D pathophysiology. Our analyses suggest that transcriptional changes in nonbeta cells may also contribute to T2D pathogenesis. Specifically, we highlight increased expression of fatty acid translocase gene CD36 (log2FC 2.296), as well as decreased expression of the guanine deaminase gene, GDA (log2FC −1.062), in T2D alpha cells. Soluble CD36 is a biomarker of T2D and diabetic nephropathy and coordinates activation of the NLRP3 inflammasome, leading to proinflammatory cytokine release and reduced insulin secretion. Within T2D delta cell transcriptomes, we note increased LAPTM4B expression (log2FC 2.871) and drastically reduced RCOR1 expression (log2FC −10.128). The underlying biological significance of these differentially regulated genes remains unclear and thus requires further investigation of their roles in nonbeta cell types and T2D pathology. We also compared the transcriptional differences between T2D and ND endocrine cells without first segregating them into islet cell types (334 ND and 212 T2D single-cell profiles). Approximately 66% of beta cell–specific (n = 165/248), 50% of alpha cell–specific (n = 67/138), and >90% of delta–specific (n = 23/24) changes in gene expression were missed when cell types were not defined and specifically compared (Fig. 5F). The decreased heterogeneity in the transcriptional profiles of cell-type–specific comparisons provides increased power to detect the transcriptomic differences and argues the importance of single-cell analysis in understanding the molecular basis of T2D.

24Recent islet single-cell studies emerged while this study was under review. We therefore sought to validate our observed cell-type–specific differences in T2D islets using these independent data sets. We found that 54/77 genes and 32/171 were also significantly up- and down-regulated, respectively, in T2D beta cells in these studies (P < 0.05, two-sided Wilcoxon rank-sum test) (Supplemental Fig. S24A,B; Supplemental Table S13). Notably, DLK1 consistently exhibited approximately fourfold induction in T2D beta cells in each study (Supplemental Fig. S24C,D) Similarly, 39/60 and 14/78 genes were significantly induced or repressed, respectively, in T2D alpha cells (Supplemental Fig. S24E,F). This included approximately twofold CD36 induction in each study (Supplemental Fig. S24G,H). Validation rates for delta cells was notably lower, likely due to the relatively few cells profiled for comparison. However, we did note a significant increase (log2FC 1.203) in LAPTM4B in T2D delta cells from, consistent with our data.

25In this study, we completed transcriptome profiling and analysis of 638 single islet cells from ND and T2D individuals. Single-cell RNA-seq protocols are often limited by their capture efficiency due to the fact that a limited proportion of each cell's total transcripts is represented in the final sequencing library. Additionally, these approaches have difficulty detecting expression and changes in expression of low abundance transcripts. Despite these limitations, we observed a strong correlation between the transcriptomes of paired bulk islets and single cells, indicating these are high-quality and representative data sets. Based on single-cell transcriptome profiles, we have identified cells of each endocrine (alpha, beta, delta, PP/gamma, epsilon) and exocrine (stellate, ductal, acinar) type in the pancreas in an agnostic and data-driven manner.

26This approach has defined expression signatures of each cell type with single-cell precision. Cell-type–specific expression patterns in our data such as MAFA in beta cells and IRX2 in alpha cells are concordant with and extend those generated on a smaller set of cells and an independent platform. Notably, our approach also uncovered important instances of shared expression between these cell types and the less common delta and PP/gamma islet populations, including genes mutated in CHI (HADH) and transcription factors regulating cell fate/identity (ARX). HADH encodes hydroxyacyl-CoA dehydrogenase, an important enzyme and negative regulator of glutamate dehydrogenase (GDH) and insulin secretion. Expression of HADH in islets has been shown to be beta cell specific, and indeed, knockdown of HADH in rat 832/13 beta cells increases insulin secretion. Surprisingly, our combined transcriptomic analyses and in situ (ViewRNA) validation of HADH revealed shared expression in beta and delta cells. These findings suggest delta cell dysfunction, in addition to beta cell dysfunction, as potential contributing factors to the development of monogenic diabetic disorders.

27Most importantly, analysis of the delta and PP/gamma islet cell transcriptomes revealed cell-type–specific expression of multiple genes that suggest important roles for these cells in islet physiology and the molecular genetics of islet dysfunction in rare (e.g., PNDM, TNDM, and MODY) and common (e.g., T2D) forms of diabetes. The novel transcriptome signatures uncovered for human delta and PP/gamma cells includes genes that strongly suggest important roles for each cell type in sensing and integrating specific systemic cues to govern islet (dys)function. This clearly warrants additional work to better understand the regulation and function of these cells in normal and diabetic states. New cell surface markers identified for each of these cell types could be used to specifically enrich and purify these populations for detailed functional analysis.

28Finally, by comparing single-cell transcriptomes from T2D and ND islets, we were able to study quantitative changes in cell populations and cell-type–specific expression in T2D pathogenesis. Although not reaching statistical significance, we did observe a trend of decreased beta cells in T2D islets versus ND islets. This difference was not as pronounced as in previous reports, possibly due to the relatively modest number of cells sampled per individual. Alternatively, as most of the T2D islet single-cell transcriptomes came from newly diagnosed individuals, this difference may also reflect the shorter duration or decreased severity of T2D in these samples compared with other studies. Previous studies suggested that beta cell de-differentiation may underlie beta cell loss in T2D. However, a subsequent study comparing human islets from 14 T2D and 13 ND individuals did not identify clear evidence of this phenomenon. Similarly, our data do not provide transcriptome-based evidence of trans-differentiation or de-differentiation phenomena in T2D islets. We observed neither the appearance of new or distinct subpopulations among the T2D islet single cells nor significant changes of reported de-differentiation genes between T2D and ND cell types (e.g., T2D beta cells vs. ND beta cells), although it is possible that de-differentiated cells were simply not captured in our study. Overall, we identify 248, 138, and 24 genes exhibiting differential expression in T2D versus ND beta, alpha, and delta cells, respectively. Consistent with Simpson's paradox, approximately half of these genes in each major islet cell type (64% beta, 45% alpha) and ∼90% of these in the less abundant delta cells were not detected in whole islet or single-cell islet transcriptomes when they were not stratified by cell type. Each of these differentially regulated genes may represent important new candidate genes in T2D pathogenesis and therapeutic targeting.

29Islets were procured from ProdoLabs or the Integrated Islet Distribution Program (IIDP) and shipped in PIM(T) media (ProdoLabs) overnight on cold packs. Upon arrival, islets were washed and transferred into PIM(S) media with PIM(G) and PIM(ABS) supplements according to the manufacturer's instructions (ProdoLabs) and incubated at 37°C in a 5% CO2 tissue culture incubator. Twenty-four hours after transfer, approximately 500 islet equivalents (IEq) were aliquoted and centrifuged at 180g for 3 min at room temperature (RT). One aliquot (100 IEq) was immediately flash frozen (Fig. 1A, baseline), one (200 IEq) was resuspended in 1 mL Prodo-media (Fig. 1A, intact), and one (200 IEq) was resuspended in 1 mL Accutase (Innovative Cell Technologies) (Fig. 1A, dissociated and single cell) and incubated for 10 min in a 37°C water bath, with pipetting every 2 min. Accutase-dissociated cells were filtered through a prewet cell strainer (BD) to collect single cells, rinsed with 9 mL prewarmed CMRL + 10%FBS media to stop the reaction, and centrifuged at 180g for 3 min at RT. Dissociated cells were resuspended in 300 µL CMRL + 10%FBS media. Cell size, number, and viability were assessed using Countess II FL (Thermo Fisher Scientific), and the cell suspension was diluted to a final concentration of 300 cells/µL. Total processing and handling time for each islet was ≤60 min.

30After counting, cells were diluted to a final concentration range of 250–400 cells/μL and 5 µL loaded onto each C1 integrated fluidic circuit (IFC; 10- to 17-μm chip) for cell capture on the C1 single-cell Autoprep system. For each islet preparation, up to two microfluidic chips were used. After capture, cells were imaged within each capture nest with an EVOS FL auto microscope (Life Technologies). IFCs were subsequently loaded with additional reagents for subsequent cell lysis; SMARTer v1- based (Clontech), olio-(dT)-primed reverse transcription; template switching for second-strand priming; and amplification of cDNA on the C1 System. Qualitative and quantitative analysis of all single-cell cDNA products was performed on a 96 capillary fragment analyzer (Advanced Analytical). Only cell singlets, as determined by imaging, with adequate cDNA yield and quality were processed for subsequent sequencing. Fragmentation and tagmentation of cDNA was done with Nextera XT reagent (Illumina) using dual indices to prepare single-cell multiplexed libraries.

31Bulk cells were pelleted and RNA purified using the PicoPure RNA isolation kit (Life Technologies). All RNA-seq libraries from bulk-sample RNA were generated with the same SMARTer v1 chemistry (Clontech) as for the C1 single-cell data largely following the manufacturer's instructions. Unlike the C1 workflow, after first-strand DNA synthesis, cDNA was purified using Agencourt AMPure beads (Beckman Coulter). cDNA was subsequently amplified through 12 PCR cycles. The cDNA yield and fragment size were measured on a 2100 Bioanalyzer (Agilent). For sequencing library preparation, amplified cDNA was sheared using a Covaris LE220 system to obtain fragments of ∼200 bp. The fragmented cDNA was prepared for sequencing using the NEBNext DNA library prep kit for Illumina sequencing (New England BioLabs).

32All sequencing was performed on a NextSeq500 (Illumina) using the 75-cycle high-output chip. RNA-seq reads were subjected to quality control using custom scripts developed at the computational sciences group at The Jackson Laboratory. Briefly, reads with >30% of bases with quality scores less than 30 were removed from the analysis, and samples with >50% of the low-quality reads were removed from the cohort. Trimmed reads were mapped to human transcriptome (GRCh37, Ensembl v70) using Bowtie 2, and expression levels of all genes were estimated using RSEM. Transcript per million (TPM) values as defined by RSEM were added a value of one prior to log2 transformation to avoid zeros. GRCh37 was selected for mapping to facilitate integration and comparative analyses with existing islet data sets (e.g.,;; ) and ENCODE and NIH Roadmap data by members of the islet biology, diabetes, and functional genomics communities. The observation of expected “positive control” genes for each cell type strongly suggested that mapping to GRCh37 instead of GRCh38 did not mask or alter any important conclusions that could be drawn from the data.

33We used 26,616 protein-coding genes and lincRNAs from the GRCh37, Ensembl v70 build in our study. Genes with expression five or more TPMs in a sample were defined as expressed. Seventy-two single-cell samples that expressed fewer than 3500 genes according to these criteria were removed from downstream analysis.

34GMM of islet marker genes was performed on a per gene basis using the R-package mclust_5.2. Each single-cell sample was classified as a specific pancreatic cell type if and only if a single gene from the selected marker gene list—INS (beta), GCG (alpha), SST (delta), PPY (PP/gamma), KRT19 (ductal), PRSS1 (acinar), and COL1A1 (stellate)—was expressed in the sample and none of the other marker genes were expressed. Cells expressing no marker genes were labeled as “none,” and those expressing >1 marker gene were labeled as “multiple.” Fluidigm released a white paper report detailing the potential for single cells to “z-stack” in up to 30% of capture nests on the medium (10–17 µm) Fluidigm C1 chip (http://info.fluidigm.com/rs/673-MRG-416/images/C1-Med-96-IFC-Redesign_wp_101-3328B1_FINAL.pdf). DAPI staining identified z-stacked islet cell doublets in 10% and 30% of capture nests from two additional C1 single-cell captures. Because the proportion of “multiple” labeled cells approximately equaled that of z-stacked doublets, we discarded these cells (n = 340) from subsequent analyses.

35Barnes-Hut variant of t-SNE was implemented using the R-package Rtsne_0.10 (https://github.com/jkrijthe/Rtsne). ND single-cell transcriptomes were reduced to two dimensions in an unsupervised manner using genes with log2 CPM values greater than 10.5 in at least one sample. Similar analyses were repeated using only the T2D single-cell data and the combined single-cell data. Hierarchical clustering of cell transcriptomes was performed using Euclidean distance, Ward.D2 linkage, and the same gene selection criteria. Resultant “fan” dendrogram images were produced using the R-packages dendextend_1.1.8 and ape_3.5. Bulk islet transcriptomes were clustered using the same criteria.

36Differential expression analyses were performed using the Bioconductor package edgeR_3.14.1. Gender was used as a blocking factor to account for variability between male and female patient islets. In each comparison, protein-coding genes and lincRNAs with 20 or fewer counts in at least 20% of either cell type population being compared or at least a minimum of three cells were used. Differentially expressed genes with FDR < 5% were regarded as significant results. Endocrine cell signature genes were identified by first performing the above differential expression analysis procedure between each endocrine cell type (e.g., beta vs. alpha, beta vs. delta, and beta vs. PP/gamma). Afterward, the intersection of these results was performed to identify genes exclusively enriched in the cell type. Exocrine cell signature genes were identified after pairwise comparisons between each respective exocrine cell type (e.g., acinar vs. stellate, acinar vs. ductal). Comparisons between T2D and ND single-cell transcriptomes were performed for the same cell types (e.g., T2D beta cells vs. ND aeta cells) to identify cell-type–specific differences in gene expression between T2D and ND states.

37For each collection of diabetes-associated and eQTL genes examined, one-way analysis of variance (ANOVA) was used to identify statistically significant differences (FDR > 5%) in islet cell-type gene expression. Following this, we performed a post hoc analysis using a THSD test to determine genes with cell-type–specific expression patterns (fold change > 4). Genes were classified as “pan-islet” if they had an average log2(CPM) expression greater than four in all islet cell types. Genes that were not enriched in a cell type or pan-islet were classified as “lowly expressed” (average log2(CPM) < 2 in all cell types), and the remaining genes were classified as “less than fourfold change.” This same methodology was used to characterize expression of the previously reported alpha- and beta-specific genes from,,, and. Similar methods were used to characterize expression patterns of genes nearby diabetes-related GWAS SNPs (downloaded from the GWAS Catalog, https://www.ebi.ac.uk/gwas/, and available in Supplemental Table S12). Protein-coding RNAs and lincRNAs that overlapped within one megabase upstream of and downstream from the diabetes-associated SNPs were analyzed.

38The Bioconductor package gage_2.22.0, was used with default parameters to identify enrichment (FDR < 5%) of human Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways in each of the ND islet cell transcriptomes. Enriched pathways were determined by comparing each cell-type transcriptome against the other aggregate islet cell-type transcriptomes (e.g., beta vs. alpha, delta, and PP/gamma).

39RNA transcripts were visualized in OCT-embedded pancreatic islet sections from two ND donors (P3 and P4) using QuantiGene ViewRNA ISH cell assay kit (catalog no. QVC0001, Affymetrix). Human INS ViewRNA type 6 probe (Catalog no. VA6-13248-06), SST ViewRNA type 6 probe (catalog no. VA6-17244-06), LEPR ViewRNA type 1 probe (catalog no. VA1-15221-06), and HADH ViewRNA type 1 probe (catalog no. VA1-12106-06) were purchased from Affymetrix (Santa Clara). The assay was performed according to the cell-based ViewRNA assay protocol with a 15-min formaldehyde fixation and a 10-min protease treatment (dilution factor 1:4000). ViewRNA probes were detected at 550 nm (Cy3) and 650 nm (Cy5) using a Leica TSC SP8 confocal microscope at 63× magnification.

40previously defined four beta cell subpopulations with differing expression of 59 genes. With this gene set, we attempted to validate the presence of these four subpopulations via unsupervised t-SNE and hierarchical clustering of all ND beta cell transcriptomes (n = 168). Similarly, characterized proliferative (Fltp+/FVR+) and mature (Fltp−/FVR−) mouse beta cells that showed differential expression of 996 transcripts. By using the Mouse Genome Informatics (MGI; http://www.informatics.jax.org) database, these 996 transcripts corresponded to 691 human orthologs that were detected in our data set. Beta cell transcriptomes were clustered using these orthologs to attempt to identify mature and proliferating subpopulations. Finally, we used the functions available in scran_1.04 (http://bioconductor.org/packages/release/bioc/html/scran.html) to computationally assign single-cell samples into cell cycle phases (G1, G2/M, or S phase) and investigate whether our data set contained proliferating islet cells.

41Raw sequence data from this study have been submitted to the databases of NCBI Sequence Read Archive (SRA; http://www.ncbi.nlm.nih.gov/sra) under accession number SRP075970 and BioProject (http://www.ncbi.nlm.nih.gov/bioproject/) under accession number PRJNA323853. Processed data sets from this study have been submitted to Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) under accession number GSE86473. UCSC Genome Browser tracks of aggregate ND islet single-cell-type transcriptomes are available at http://genome.ucsc.edu/ and may be accessed with the user name “lawlorn” and session name “Islet_Single_Cell_Type_Transcriptomes.” The source code used to produce the figures and tables in this paper is available in the Supplemental_Methods_Source_Code folder as suggested by.