Single-cell transcriptomes reveal characteristic features of human pancreatic islet cell types

Abstract

Pancreatic islets of Langerhans contain several specialized endocrine cell types, which are commonly identified by the expression of single marker genes. However, the established marker genes cannot capture the complete spectrum of cellular heterogeneity in human pancreatic islets, and existing bulk transcriptome datasets provide averages across several cell populations. To dissect the cellular composition of the human pancreatic islet and to establish transcriptomes for all major cell types, we performed single-cell RNA sequencing on 70 cells sorted from human primary tissue. We used this dataset to validate previously described marker genes at the single-cell level and to identify specifically expressed transcription factors for all islet cell subtypes. All data are available for browsing and download, thus establishing a useful resource of single-cell expression profiles for endocrine cells in human pancreatic islets.

Keywords: alpha cells; beta cells; diabetes; marker genes; single‐cell RNA‐seq.

Figure 1. Single‐cell transcriptomes recapitulate the major pancreatic cell types

1. Workflow for obtaining and analyzing single‐cell RNA‐seq data from human pancreatic islets.

2. Iterative PCA/biplot‐based approach for the identification of cell types and cell type‐defining transcripts from single‐cell RNA‐seq data.

3. Expression (scaled RPKM values) of cell type‐defining genes as identified in (B) across all single cells. Transcripts and single cells are grouped by cell type as identified in (B).

4. Display of transcriptional similarity between all single cells by MDS. The coloring scheme is based on the cell types as identified in (B).

5. Relative expression (scaled RPKM value) of canonical marker genes for the 6 identified pancreatic islet cell populations represented by bubble size and projected onto the MDS profile.

Figure EV1. Statistical analysis of single‐cell RNA‐seq data

1. Number of detected transcripts and total aligned reads for each single cell. The red line denotes 500 transcripts, below which samples were excluded from the analysis.

2. Scatter plots displaying the correlation between the number of input ERCC RNA molecules and measured RPKM values in four representative single cells. The Pearson correlation (r) is noted in the upper left corner.

3. Scatter plots correlating raw read counts (left) and RPKM normalized expression values (right) with transcript length. All adequately covered transcripts (> 25 reads) of all 64 samples were included in this analysis. The observed negative correlation after RPKM normalization lies in the range of what had been reported previously 19.

4. Raw counts and RPKM normalized expression values for 10 groups of ERCC spike‐in controls. The amount of molecules spiked into the sequencing reaction is constant within one group, whereas the length of the transcripts varies considerably. Within each group, equal expression values across different transcript lengths thereby confirm that RPKM normalization does not systematically penalize longer transcripts in the assessed ERCC transcript length range of ˜200 to 2,000 bp. ERCC spike‐in controls covered by more than 25 reads are indicated in blue and those with ≤ 25 reads in red.

5. MDS displaying transcriptional similarity between corresponding cell types of a published dataset 20 (prefix: ext) and the current dataset.

Figure EV2. Assessing cross‐contamination between alpha and beta cells

1. Scatter plot displaying single alpha and beta cells, 500‐cell islet samples, as well as bulk islet and beta cell samples from published datasets according to their weighted mean of scaled expression values in alpha and beta cell‐specific profile genes. The three selected profile cells for each cell type are indicated by their sample ID.

2. Pure and mixed expression profiles consisting of 233 alpha cell‐specific genes and 252 beta cell‐specific genes. Alpha and beta cell‐specific profiles are calculated based on the expression values of the three selected profile cells only, while profile genes were selected based on all single cells classified as alpha or beta cells, which is why the expression gradients in the mix profiles do not always follow the same direction.

3. Profile correlation curves for each individual sample. The maximum of each curve defines the maximum variance that can be explained (y‐axis) by the corresponding mix profile (x‐axis) providing a measure for the composition of the respective sample.

4. Diagram explaining the transition from profile correlation curves to sample composition estimates. The profile composition that explains most variance is linearly scaled to the maximum variance explained.

5. Sample composition estimates for each assessed sample. The differences between the 500‐cell islet samples and bulk samples might be explained by technical effects that enrich for alpha cells during islet cultivation, disassociation, and FACS purification.

Figure 2. Expression of cell type‐specific transcription factors at single‐cell resolution

1. A

   Merged UCSC Genome Browser tracks for the PDX1 and ARX loci. The respective tracks for all single cells are presented in Appendix Fig S2.

2. B

   Relative expression (scaled RPKM value) of important transcription factors represented by bubble size and projected onto the MDS profile.

3. C–F

   Cell type‐specific expression of pan‐endocrine (C), alpha cell (D), beta cell (E), and PP cell (F) transcription factors (red bar: mean expression). The statistical significance of the differential gene expression is presented in Appendix Fig S6.

Figure EV3. Assessing pairwise correlation of endocrine marker genes

Correlation matrix displaying all genes (y‐axis) that are highly correlated (r > 0.9) with at least one of the endocrine marker genes (x‐axis). Different transcripts of the same gene are indicated by “Tx”.

Figure 3. Single‐cell transcriptomes reveal unique features of human islets

1. A

   Heatmap displaying the P‐values obtained by pairwise Gene Set Enrichment Analysis (GSEA) for the REST‐binding motif.

2. B

   Relative expression (scaled RPKM value) of genes contained in the REST‐binding motif gene set.

3. C

   Merged UCSC Genome Browser tracks for REST. The respective tracks for all single cells are presented in Appendix Fig S3.

4. D

   Expression of DLK1 and GC in human islet cell types (red bar: mean expression). The statistical significance of the differential gene expression is presented in Appendix Fig S6.

5. E–H

   Co‐staining of DLK1 (E, F) or GC (G, H) with insulin and with glucagon in representative human (E, G) and mouse (F, H) islets.

Figure EV4. Specific expression of selected marker genes

Relative expression (scaled RPKM value) of interesting genes across all single cells represented by bubble size and projected onto the MDS profile as displayed in Fig 1D.

Figure EV5. Assessing cell type specificity of genes identified in diabetes‐related GWAS

1. Cell type specificity for genes reported in diabetes‐related GWAS. Each gene reported in a diabetes‐related GWAS (search for “Diabetes” on GWAS Catalog) was assigned to the pancreatic cell type in which it was found to be most specifically expressed. Ranking was performed as described in Appendix Fig S7B and Dataset EV6.

2. Heatmap showing mean expression values for the most cell type‐specific diabetes‐associated GWAS genes in the different here identified human islet cell types. Specifically, only genes with a specificity rank lower than 500 (dashed line in panel A) are listed, and genes with equal expression in multiple cell types are not shown. The numbers in the colored boxes indicate the number of studies in which the respective gene has been reported. Heatmaps are colored by mean ln(RPKM), the mean of the natural logarithm of the RPKM values across all cells of the respective cell type.