Single-Cell Heterogeneity Analysis and CRISPR Screen Identify Key β-Cell-Specific Disease Genes

Drop-Seq Analysis of Human Islet Samples

We prepared Drop-Seq libraries with fresh human islet samples from 6 healthy (3 overweighed with BMI >30) and 3 T2D donors (2 overweighed). In total, we obtained transcriptome data from 39,905 single cells (1,206–9,409 cells from each donor, Figure 1A) and used a very stringent clustering-based analysis pipeline to determine the types of 28,026 “clean” cells without ambiguity (Figure S1; Data S1). When projecting the cells to a two-dimensional t-distributed stochastic neighbor embedding (tSNE) plot, we observed a clear distinction between endocrine cells and a few non-endocrine cell types, mainly pancreatic ductal cells (PDCs) marked by several keratin genes (KRTs), and pancreatic stellate cells (PSCs) marked by collagen genes (Figures 1B and 1C). We observed very few acinar cells marked by REG1A and PRSS1 genes, which were identified as PCA outliers but failed to form a distinct cluster in t-SNE due to the scarcity (n = 108, Figures 2A-2D). We further performed a second-round unsupervised clustering with the endocrine cells and distinguished four major endocrine clusters, which are recognized as α, β, δ, and PP cells based on the enrichment of corresponding marker genes (Figures 1D and 1E). We could not observe a distinct cluster of ε cells in tSNE due to the extreme scarcity of this cell type in our samples: only 13 of the 28,026 “clean” cells express the ε cell hormone gene GHRL (Figures 2A-2D). Taken together, all of the samples contain 10%–20% non-endocrine cells (Figure 1F), consistent with an estimated 80%−90% islet purity, and ~90% of endocrine cells in every donor are α or β cells (Figure 1F).

Gene Signatures of Non-endocrine Cell Types

We first used a negative binomial model to define the non-endocrine cell marker genes (STAR Methods), including a number of transcription factors (TFs) that may function as master cell type regulators (Figures 2E and 2F; a complete gene list is included in Data S2). As expected, PSCs express collagen genes; ductal cells express keratins, CFTR, and TSPAN8; and exocrine acinar cells express REG1A, REG1B, REG3A, and a number of digestive enzymes (Figures 2E and 2G). Interestingly, we also observed duct-exclusive expression of inflammation genes such as CCL2, CXCL2, MMP7, and DEFB1 (Figures 2E and 2G). It is not clear whether these duct-expressed inflammatory genes may contribute to disease initiation or development, although the elevated levels of cytokines in pancreatitis and pancreatic cancer is well documented.

We observed two clusters of PSCs representing activated and quiescent populations; ~70% of PSCs are activated (Figure 1B). Both populations contain cells from multiple donors (comparing Figures 1B with S2A). It is known that under quiescent state, the PSCs contain abundant lipid droplets in the cytoplasm, and pancreatic injury or other pathogenic insults can activate PSCs, leading to increased migration and proliferation, more production of extracellular matrix (ECM) proteins, metalloproteinase (MMPs), and tissue inhibitors of metalloproteinase (TIMPs) (Omary et al., 2007). Indeed, we found that the activated PSCs expressed more collagen, TIMPs, and stem cell or cancer signature genes but lost the expression of lipid-processing or adipogenesis genes such as FABP4 and ADIRF (Figures 2H and 2I).

Gene Signatures of Endocrine Cell Types

We next defined hundreds of genes specifically expressed in various endocrine cell types (Figure 3A; Data S2). As expected, β-cell-specific genes are associated with insulin secretion and diabetes (Figure 3B). We identified a number of endocrine-specific transcription factors with dramatic specificity between different endocrine cell types (Figure 3C). For example, PAX6, NEUROD1, and INSM1 are expressed in all four, and MAFB is expressed in three endocrine cell types. Our top list includes many classically known cell identity master transcription factors, such as PDX1, MAFA, and NKX6-1 in β cells, ARX and IRX2 in α cells, and HHEX in δ cells (Figure 3C). Consistent with a few other single-cell results (Baron et al., 2016; Lawlor et al., 2017; Segerstolpe et al., 2016), we also observed dual specificity of PDX1 (β and δ) and ARX (α and PP). Some signaling-dependent TFs are also cell-type-specific, such as ESR1 (β cell) and EGR2, EGR3, and STAT4 (PP cell) (Figure 3C). These signal-dependent TFs may have important physiological functions. For example, studies in both human and mouse have shown that the activation of estrogen receptor α (ESR1) protects β cells from apoptosis and preserves functional β cell mass in diabetes (Tiano and Mauvais-Jarvis, 2012).

Diabetes and obesity are two complex multigene disorders contributed by many tissues. We collected over 1,000 GWAS risk genes (Data S3) and reasoned that knowing their expression pattern in endocrine cells, especially β cells, would be helpful for the understanding the disease etiology. We found 163 GWAS risk genes specifically expressed in one or more endocrine cell types (Figure 3D; Data S3) including some well-known β cell genes RNF6, PDX1, SLC2A2, MEG3, and DLK1. Interestingly, a number of risk genes are δ-cell-specific, including HHEX, LEPR, ERBB4, etc., suggesting a particularly important role of δ cells to diabetes or obesity (Figure 3D).

Last, we examined the endocrine-specific cell surface genes due to their potential as cellular markers, or targets for pharmacological intervention (Figure 3E; Data S3). Interestingly, we noticed that δ cells specifically express several important hormone receptor genes, including LEPR (receptor of leptin), GHSR (receptor of Ghrelin), GHR (receptor of growth hormone), ERBB4 (receptor of EGF), and two GABA receptors (GABRG2 and GABRA1) (Figure 3E). Receptor activity is actually one of the top functional categories for the δ-cell-specific genes (Figure 3B). These results suggested a key regulatory role of δ cells by integrating multiple cell signaling (Segerstolpe et al., 2016).

Identifying α or β Cell Subpopulations

Previous studies, including several recent single-cell RNA-seq analyses, have investigated β cell heterogeneity but reached discrepant conclusions regarding the existence of β cell subpopulations, partly due to the low cell numbers analyzed (Gutierrez et al., 2017). We posited that Drop-Seq would be more sensitive in resolving distinct cell populations with greatly improved throughput. In this study, we pooled α or β cells from all donors and used principle component analysis (PCA) to identify cell subpopulations as outliers (Figures S2B and S2C) then projected the cells onto t-SNE plots to confirm their presence in multiple donors (Figures S2D and S2E). Marker genes for these subpopulations are listed in Data S2.

The most obvious α cell subpopulation consists of a small number of proliferating α cell expressing TOPA1, CENPF, and AURKB (Figures 2F and S2B). Importantly, proliferating α cell is a reproducible population that can be found from multiple donors including H1, H3, H4, H5, H6, and T2D3 (Figure S2E, cells in dark green). Notably, the proliferating α cells were also reported in another recent single-cell study but with low sensitivity due to limited cell throughput (Segerstolpe et al., 2016). Our study has therefore confirmed the existence of this rare α cell population, which may also play an important role in β cell replenishment (Thorel et al., 2010).

The biggest subpopulations are α-HS and β-HS cells expressing the same set of heat shock genes including DNAJA1, DNAJB1, HSPAs, and HSPBs, etc. (Figures 1E, 2F, 2G, and S2B-S2E). However, both α-HS and β-HS cells are mainly found in donor H3 and T2D3 (Figure S2E). Many types of stresses can activate heat shock genes (HSPs). Reduced expression of HSPs is associated with diabetes, and upregulation of HSPs may provide a cytoprotective effect to β cells (Hooper and Hooper, 2009). Interestingly, a recent single-cell study also observed a correlation between stress gene expression and aging (Enge et al., 2017). Last, we also identified two other minor cell clusters (α-KCNQ1OT1 and β-KCNQ1OT1), both of which were exclusive to donor H5 (Figure S2). More donors are necessary to evaluate the physiological relevance of these individual-specific subpopulations.

Importantly, we did not find any of these subpopulations correlate with the BMI or T2D status of the donors, leading to a conclusion that disease-associated effect on α or β cells does not create distinctive cell subpopulations. Therefore, we excluded the minor cell populations in the following disease association analysis and focused on the continuous changes of α or β transcriptome associated with disease status.

Obesity and T2D Cause Different Single-Cell Transcriptome Heterogeneity

We next investigated the molecular differences between cells from healthy, obese, and T2D donors. A few recent single-cell studies have explored this question (Lawlor et al., 2017; Segerstolpe et al., 2016; Xin et al., 2016). These studies used conventional statistical models comparing cells from normal and patients and looked for differentially expressed genes correlated with disease states. However, these models do not account for cellular heterogeneity and treat all cells from the same donor equally, which may result in low sensitivity for the following reasons: (1) disease-causing cells can be too rare to be detected when a majority of cells are normal (Figure 4A); and (2) the difference between normal and patient cells can be small and masked by large individual variance. In both scenarios, many human samples become necessary to increase the statistical power.

Here, we propose a strategy to improve the sensitivity to identify disease signature genes by dissecting disease-associated single-cell heterogeneity. Our strategy is based on a fundamentally different concept assuming the presence of disease relevant cellular heterogeneity from the same donor. Specifically, obesity or T2D may cause different transcriptome heterogeneity at single-cell level. Theoretically, by ranking and directly comparing cells at normal, transitional, and disease states, it is possible to improve the sensitivity even if very few donors are available (Figure 4A). Figure 4B shows scatterplots of all the β cells in the space of top three principle components (PCs). Remarkably, although cells from different donors do not segregate in the PCA plots, there is a clear continuous shift correlated with BMI or T2D status. Importantly, the two disease conditions show different trajectories, indicating different effects on the transcriptome.

RePACT: A Sensitive Approach to Identify Disease-Relevant Gene Signatures from Single-Cell Data

We therefore developed RePACT (regressing principle components for the assembly of continuous trajectory), to identify genes associated with disease relevant cellular heterogeneity. In this algorithm, we first performed PCA to reduce the dimension of transcriptome data. Next, we used regression analysis to draw two optimal trajectory lines reflecting the obesity- or T2D-relevant variation. Figure 4B is an example when we only used the top three PCs to plot the trajectory lines. In this study, we actually used top 10 PCs as predictors (STAR Methods). The numeric projection of each cell on the obesity or T2D trajectory (BMI index and T2D index) served as a measurement of the degree to which the cell has transformed during disease development. We then binned the cells into a number of pseudostates according to the index values (Figure 4C). By comparing cells from different pseudo-states, RePACT greatly improved the statistical power to identify gene signatures for obesity or T2D status in α or β cell (Figures 4D and S3A-S3C; STAR Methods). To evaluate the robustness of RePACT method, we called top 200 T2D signature genes in α or β cells using data from all three T2D donors and then compared them to the results when only two T2D donors were included (dropout analyses). We found that 91%–95% of the T2D signature genes can be recovered from at least two of the three dropout analyses (Figures S3D-S3E), proving that the signature genes identified from RePACT were highly reproducible.

RePACT Identifies Common and Specific β Cell Gene Signatures Associated with Obesity and T2D

We identified 1,368 T2D trajectory genes and 1,188 obesity trajectory genes in β cells (Figure 4D; Data S4). It is known that obesity increases diabetes risk, and as expected, the two signatures shared many common genes (Figure 4E). For example, GAPDH, IAPP, SPP1, and CPE are downregulated in both obesity and T2D trajectories (Figure 4D). GAPDH is a key enzyme for glycolysis and glucose metabolism. IAPP, SPP1, and CPE are all well-known diabetes or obesity risk genes. Particularly, IAPP forms islet amyloid, which is linked to cellular toxicity in T2D (Clark and Nilsson, 2004); SPP1 mediates the obesity-induced macrophage infiltration into adipose tissue and insulin resistance (Nomiyama et al., 2007); and CPE is important for proinsulin processing and its mutation leads to obesity and hyperproinsuline-mia in mouse (Naggert et al., 1995). With the single-cell trajectory data, we also determined the transcription dynamics of 149 GWAS risk genes (Figure S4A), which may help the dissection of etiology of obesity or T2D.

Most interestingly, we found many genes specific to one trajectory but not the other, including genes with opposite trends in the two trajectories (Figure 4E). The best example is probably the insulin gene (INS) itself, which is upregulated in obesity (consistent with the hyperinsulinemia in obesity) but downregulated in T2D (consistent with the β cell dysfunction in T2D) (Marchetti et al., 2008; Templeman et al., 2017). On the contrary, two ferritin genes (FTL and FTH1) are downregulated in obesity but upregulated in T2D (Figure 4D). Ferritin is the major intracellular iron storage protein. Clinical results have shown that low-serum iron concentration is associated with obesity (Lecube et al., 2006; Nead et al., 2004), while iron overload is a risk factor for T2D (Simcox and McClain, 2013).

Another interesting observation is that although only 15 obesity-upregulated genes are downregulated in T2D (including INS), many more (99) obesity-downregulated genes shift their trends in to upregulation in T2D (Figure 4E). We therefore performed a gene set enrichment analysis (GSEA) on genes with common or opposite trends in obesity or T2D (Figure 4F). For example, consistent with the connection to hypoxia in both T2D and BMI (Ye, 2009), we found upregulation of hypoxia, glycolysis, and HIF1A target genes, but downregulation of aerobic respiration pathways in both trajectories (Figure 4F). Surprisingly, the proteasome genes or pathway are elevated in T2D but downregulated in obesity trajectory (Figure 4F). Polymorphisms at proteasome genes have been associated to obesity (Kupca et al., 2013), and proteasome activity is weakened in obese liver that can induce endoplasmic reticulum (ER) stress and insulin resistance (Otoda et al., 2013). However, a recent study showed that in β cells, the inhibition of proteasome activity improves insulin production (Weisberg et al., 2016). Our results therefore suggest an important mechanism that is differentially presented in obesity and T2D governing β cell function.

A Genome-Wide Gene Deletion Analysis to Identify Insulin Regulators

The primary functions of β cell are to produce, store, and secrete insulin. Although we have identified a number of signature genes in β cells with Drop-Seq, the remaining question is whether these genes affect β cell functions. To address this problem, we decided to perform an unbiased genome-wide CRISPR screen to identify insulin regulators. We choose to use the MIN6 insulinoma cell because it is one of the most robust β cell lines for the study of glucose sensing and insulin secretion (Skelin et al., 2010). Although MIN6 is a mouse cell line, and the functions of screen hits eventually need to be test in native human islets, we reasoned that most of the mouse insulin-regulatory genes should have conserved functions with their human orthologs. Therefore, the genome-wide mouse CRISPR data should still serve the purpose of providing an additional layer of information about gene functions.

Briefly, we packaged a lentiviral library using the mouse GeCKOv2 single guide RNA (sgRNA) library that contains ~130 K sgRNAs (6 sgRNAs per protein-coding gene in the mouse genome), and transduced >100 million MIN6 cells at low MOI so that most cells were infected with only one virus. After 7 days of puromycin selection expansion of transduced cells, we fixed and permeabilized the cells and sorted ~200 million cells based on the intracellular insulin protein level (Figure 5A; STAR Methods). Flow cytometry analysis demonstrated an increased variation of the insulin signal among the infected cells, suggesting that the screen has successfully targeted both positive and negative insulin regulators (Figure 5B). We collected top (Ins^Hi) and bottom (Ins^Lo) 10% cells and amplified the sgRNA sequences from integrated viral DNA for next generation sequencing. By comparing the abundance of each sgRNA in the two cell populations, we can determine the effect of target genes on insulin level in the cells.

We performed two independent screens and called hit sgRNAs only if they are reproducible between two replicates. Based on their significance and reproducibility, (Figure 5C; STAR Methods), we classified hit sgRNAs into five tiers (tier one sgRNAs have p < 0.001 in both replicates, false discover rate [FDR] = 0.1) and called 373 hit genes with at least one tier-one sgRNA. The design of sgRNA redundancy (6 sgRNAs per gene) also allowed us to evaluate the off-target risk for all hit genes from the screen. We found that most tier-one hit genes (223 of 373, or 59.8%) are supported by two or more sgRNAs (Data S5). However, the more likely reason for a sgRNA to fail the significance test is its under-representation in the GeCKO sgRNA library: the concentration of most abundant sgRNAs can be over 100-fold higher than the low-abundance ones (Figure S5E). Therefore, in this study, we choose to report all tier 1 genes, but also divide tier 1 hit genes into sub tiers to indicate if they are supported by multiple sgRNAs: 1A (119 genes, ≥3 supporting sgRNAs, FDR = 0.001), 1B (104 genes, 2 supporting sgRNAs, FDR = 0.038), and 1C (150 genes, 1 supporting sgRNA, FDR = 0.28). Note that there is a significant off-target risk for Tier 1C genes because they only have one supporting sgRNA (STAR Methods). Proteins from the same complexes, such as Nelfa-Nelfb-Nelfe and Mau2-Nipbl, were often identified together, suggesting that our CRISPR screen is sensitive and robust (Figure 5C). We called a gene “positive regulator” if its’ sgRNA was enriched in Ins^Lo population, meaning that the gene increased the intracellular insulin amount. As expected, the strongest positive insulin regulator was Ins2, the mouse ortholog of human insulin gene INS. Conversely, we also called “negative regulators” from sgRNAs enriched in the Ins^Hi population (Figure 5C).

Multiple mechanisms may regulate the intracellular insulin levels. Among the top hits (Figure 5C), we found transcriptional regulators (such as Glis3 and Nelfs), translational regulators (such as Eef2), and post-translational regulators (such as unfolded protein response regulator gene Xbp1) (Lee et al., 2011). Importantly, our screen also identified regulators of insulin secretion, especially genes involved in glucose induced insulin secretion (GSIS). For example, Glp1r (Trujillo and Nuffer, 2014), a known GSIS regulator, was identified as one of the top negative regulators (Figure 5C). Notably, although 74 of 373 insulin regulators are also “fitness genes” (leading to slow growth upon deletion) (Hart et al., 2015), the distribution of fitness genes among positive (17 of 117, 15%) or negative insulin regulators (57 of 256, 22%) is not biased (p = 0.094, Fisher’s exact test) suggesting that slow-growing cells do not cause a bias in our screen.

To validate the hits from the screen, we also infected MIN6 cells with individually cloned sgRNAs and used flow cytometry to measure the intracellular insulin levels (STAR Methods). We verified 13 out of 20 (65%) positive regulators and 30 out of 33 (91%) negative regulators (Figure 5D). We further test 20 of the verified genes with newly designed sgRNAs and re-validated 19 of them, indicating a low off-target rate (Figure S5F). It is worth noting that our genome-wide CRISPR analysis tends to report higher fold changes than individual experiments. This is mostly likely because in the individual experiments, a significant fraction of cells did not obtain loss-of-function mutations, whereas the selection of Ins^Hi and Ins^Lo populations in the screen enriches the cells with successful gene knockout events. As a result, the pooled screen had a higher sensitivity, which may also explain why some CRISPR hits were not validated in the individual experiments.

Integrative Analysis of Single-Cell and CRISPR Data Revealed Disease-Relevant Insulin-Regulating Modules

We performed STRING analysis (Szklarczyk et al., 2015) on the 373 tier-one insulin regulators and drew all the functional or physical associations between these genes (Figure 6; STAR Methods). Compared to Drop-Seq data, we identified 100 insulin regulators up- or downregulated in the diabetes or obesity trajectory, including Cdkn2a, Cox7c, Glis3, and Xbp1 (Figure 6 and S4B). Therefore, these transcriptional changes may play a causal role in the disease development. We also identified 40 insulin regulators specifically expressed in β cell including Glp1r. These genes may have important selective functions in β cells (Figure 6; Data S5). Seventeen β-cell-specific insulin regulators also change their expression in diabetes or obesity trajectories (Atp6v1h, Bcat2, Cct8, Cdkn2a, Cops6, Dhps, Eif3e, Gnas, Ins2, Lmo2, Mrpl22, Nelfcd, Rph3al, Rtcb, Srp9, Taf7, and Xbp1) (Figure 6).

The network analysis also expanded the repertoire of disease relevant insulin regulators. Our analysis highlighted nine notable modules of insulin regulators; each module contains at least one diabetes or obesity signature genes in β cell, and all of them represent multi-protein complexes with cooperative functions (Figure 6). The identification of a whole protein complex from the CRISPR screen is strong evidence that the complex is key for β cell function. Alteration of signature genes in these complexes may influence their normal functions and contribute to the disease development. For example, the largest module is the protein translation and processing network. Most of the genes in this module are positive insulin regulators, presumably because they positively regulate insulin protein production and maturation. Negative regulators in this module, such as PP2A protein phosphatase members Ppp2r1a and Ppp2ca, may function as inhibitors of protein translation or processing.

Six protein modules appear to regulate insulin release (Figure 6): (1) GPCR signaling complex (GPCR signal such as Glp1r amplifies insulin secretion); (2) mitochondria ATP production module (important for ATP production, therefore ATP-dependent insulin release); (3) synaptic vesicle exocytosis module (insulin release); (4) COP9 signalosome (protein ubiquitination and ubiquitin-dependent endocytosis); (5) LAMTOR Ragulator complex (late endosome or lysosome scaffold proteins); and (6) chaperonins (promoting insulin secretion by assisting protein folding). Interestingly, nearly all genes in these modules are negative insulin regulators, i.e., once these genes are deleted, insulin accumulates in the cells due to the impaired secretion. The only positive regulator in these modules is Sytl4, a known inhibitor of exocytosis. Consistently, another exocytosis inhibitor Stxbp5 is also a positive insulin regulator (Figure 6). It is interesting that the endocytic modules are also negative insulin regulators, because this suggests that insulin granules require not only the exocytosis process for secretion, but also the endocytosis process to ensure vesicle recycle or membrane recapture (MacDonald and Rorsman, 2007). Misregulation of either of these two trafficking pathways will lead to the failure of sustainable insulin release.

Our CRISPR screen also identified two previously unrecognized insulin regulators modules: the cohesin function module and the NuA4/Tip60 histone acetyltransferase (HAT) complex module (Figure 6). The β cell functions of these two complexes are unknown, so we performed additional analyses to gain further insights into the underlying mechanisms (Figure 7).

Mau2-Nipbl Cohesin Loading Complex Regulates Insulin Gene Transcription

Cohesin is a multifunction protein complex. During S phase and M phase, cohesin tethers two sister chromatids together and is essential for proper chromosome segregation in mitosis (Peters et al., 2008). In the G1 phase, cohesin is important for long-range chromosome interactions at promoters, enhancers, and CCCTC binding factor (CTCF)-occupied insulator elements. Our CRISPR screen identified two cohesin proteins, Smc3 and Pds5b, as negative regulators, but intriguingly, cohesin loading factors Mau2 and Nipbl were identified as positive insulin regulators (Figure 6). In the mammalian genome, Nipbl mainly binds to promoters and enhancers, while the strongest peaks of cohesin are at intergenic CTCF sites (Busslinger et al., 2017; Kagey et al., 2010; Zuin et al., 2014). Therefore, the cohesin loading complex may play a direct role regulating gene transcription independent of the CTCF-bound cohesin or mitotic cohesin (Zuin et al., 2014). Consistent with these findings, we found that Mau2 and Nipbl bind to Ins2 promoter, and the transcription level of Ins2, rather than Ins1, was lower upon Mau2 or Nipbl deletion (Figures 7A and 7B). Our single-cell data demonstrated that in human islets, MAU2 is expressed at higher levels in β cells than in other cell types (Figure 7C). Interestingly, from the 171 human pancreas samples collected by the Genotype-Tissue Expression (GTEx) project (GTEx Consortium, 2015), we found a significant positive correlation between the transcription of MAU2 and INS gene (Figure 7D). Taken together, our results suggested a β-cell-specific function of cohesin loading factors in regulating insulin gene transcription. These results may also shed light on the mechanism of Cornelia de Lange syndrome (CdLS), which has been linked to protein mutations in NIPBL and cohesin complex (Liu and Baynam, 2010; Liu et al., 2009).

NuA4/Tip60 HAT Complex Regulates Insulin Secretion

We were interested in the β cell function of NuA4/Tip60 HAT complex module because Kat5 and Dmap1 were the strongest hits from our CRISPR screen (Figure 5C), suggesting a direct insulin regulatory function. The best-characterized function of the NuA4/Tip60 complex is transcriptional activation by acetylation of histones H4 and H2A. However, this trans-activity cannot explain the increase of intracellular insulin level upon deletion of NuA4/Tip60 complex members Kat5, Dmap1, and Brd8 (Figures 5D and 7F). We therefore tested the possibility that the NuA4/Tip60 complex may control insulin secretion. First, we observed compromised GSIS when treating either MIN6 cells or primary human islets with NU9056, a KAT5-specific acetyl-transferase inhibitor (Coffey et al., 2012). Acetate can enhance GSIS, presumably by increasing the cellular level of acetyl-CoA, the substrate of all acetyl-transferases including KAT5 (Figure 7E) (Shimazu et al., 2010). Furthermore, deletion of NuA4/Tip60 complex proteins (Kat5, Dmap1, and Brd8) as well as known insulin secretion regulators (Glp1r, Gnas, and Stxbp1) all led to accumulation of intracellular insulin, lower baseline insulin secretion, and GSIS defect (Figures 7F-7H). As a control, deletion of positive regulator Mau2 and Nipbl did not affect glucose responses. Interestingly, the baseline insulin secretion from Mau2 and Nipbl knockout cells increased (Figure 7G); the reason remains elusive because knocking out these two proteins may cause widespread effects on gene transcription, genome architecture, and cell cycle. Taken together, these data strongly suggest that the acetyltransferase activity of NuA4/Tip60 complex is key for insulin secretion from β cell. It will be interesting to elucidate the downstream targets of NuA4/Tip60 complex responsible for insulin release.

Higher Throughput Single-Cell Analysis

We have performed an in-depth analysis of 39,905 single islet cell transcriptome from nine human donors. To our knowledge, this is the largest single-cell transcriptome dataset in primary human islets. The data of identified cell types and their signature genes will serve as a valuable resource to understand islet cell differentiation and disease development. Due to improved throughput, we achieved high sensitivity and robustness to detect rare cell subpopulations, such as the proliferating α cells. We also observed a few rare cell subpopulations that are relatively individual-specific; more donors are necessary to obtain a comprehensive picture of endocrine cell heterogeneity.

Highly Sensitive RePACT Algorithm Identifies Disease-Specific Cellular Changes

Obesity and T2D are two highly related diseases that share many common characteristics. Understanding the commonality and differences between the two diseases can be challenging because many patients have both diseases. To address this challenge, we developed RePACT, a general single-cell analysis algorithm to identify disease relevant genes with high sensitivity. Similar to several published single-cell analysis methods in differentiation systems (Trapnell et al., 2014), the key step of RePACT is to determine a series of pseudo states best reflecting disease-relevant variance. By doing this, RePACT can discern a continuous disease trajectory that is too subtle to be detected from bulk analysis (Figure 4A). In this study, with merely nine donors, RePACT defined different trajectories for obesity and T2D and further identified numerous common and specific signature genes in b cells. We believe that the use of RePACT in future single-cell studies with more donors will significantly improve the sensitivity and robustness. The RePACT strategy is applicable to the studies of many other human diseases, especially when the sample availability is an issue.

Integration of Genome-wide CRISPR Analyses to Reveal Disease-Contributing Genes and Pathways

Recent development of CRISPR-based genome editing has enabled pooled genome-wide screens in mammalian cells (Gilbert et al., 2014; Konermann et al., 2015; Parnas et al., 2015; Shalem et al., 2014; Wang et al., 2014). In this study, in order to reveal the functions of β cell signature genes from Drop-Seq data, we also performed an unbiased CRISPR screen for insulin regulators. Simply using resting-state intracellular insulin level as a readout, our screen identified not only transcriptional and translational regulators, but also numerous genes that control energy production, protein folding, vesicle trafficking, etc., suggesting that insulin production and release is controlled by a highly complex network. Importantly, the integrative analysis of single-cell transcriptome and CRISPR screen data highlighted potential causal genes in diabetes or obesity disease development, among which we have identified the Mau2-Nipbl cohesin loading complex and the NuA4/Tip60 HAT complex as two previously unrecognized insulin regulating protein modules. We further revealed that the Mau2-Nipbl complex regulates insulin gene transcription, and a surprising role of the NuA4/Tip60 HAT complex in regulating insulin protein release. Our study provides a general strategy for systematically characterizing diabetes genes in pancreatic islets, as well as disease genes in other complex tissues.

CONTACT FOR REAGENT AND RESOURCE SHARING

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Yan Li (yxl1379@case.edu).

EXPERIMENTAL MODEL AND SUBJECT DETAILS

METHOD DETAILS

QUANTIFICATION AND STATISTICAL ANALYSIS

DATA AND SOFTWARE AVAILABILITY

Raw data files of Drop-seq and Chip-seq are accessible at GEO: GSE101207; RePACT package is available at https://github.com/chenweng1991/RePACT.