Single-cell RNA sequencing of murine islets shows high cellular complexity at all stages of autoimmune diabetes

Type 1 diabetes (T1D) is a chronic progressive autoimmune disease that results in β-cell demise caused by autoreactive CD4 and CD8 T cells. Irremediable control of the disease leads to permanent insulin dependence. Therefore, it is crucial to understand the immunological events during initiation and progression of T1D. The information from such studies may facilitate the development of targeted therapies for restricting the autoimmune process. Examination of patients with T1D has led to the identification of clinical stages in disease progression. Early stages can be identified by the presence of autoantibodies to β-cell proteins (Pihoker et al., 2005; Pietropaolo et al., 2012; Ziegler et al., 2013), particularly to those reactive with native insulin. Insulin is likely the primary autoantigen (Nakayama et al., 2005; Zhang et al., 2008). Autoantibodies to other β-cell proteins follow, showing an expansion of the autoreactivity that eventually progresses to overt clinical disease (Ziegler et al., 2013). Limited studies have been performed on human islets during the development of T1D. Histopathologic changes showed heterogeneity in the degree of islet inflammation, but a strict correlation between immune cell composition and activation status with stages of disease progression is difficult to establish using human specimens, in which the time of initiation cannot be estimated.

Recent studies applied advanced imaging mass cytometry techniques to compare the human pancreas in T1D patients and in control individuals. Changes were found in islet structure and size, vascularization, and levels of various immune cells, as well as progressive cell loss preceded by infiltration of cytotoxic and helper T cells in the pancreatic islets (Damond et al., 2019; Wang et al., 2019). However, the analysis of stages during T1D progression is complicated for multiple reasons, particularly the late manifestation of the disease and limited availability of samples. Most human studies include recent-onset donors, in whom the majority of β cells are already affected: the key information regarding the autoimmune attack preceding β-cell death is lacking.

The nonobese diabetic (NOD) autoimmune mouse model allows us to obtain islets and examine at greater depth their cellular contents during all stages of the autoimmune reactivity. The disease process in the NOD mouse is chronic, a feature similar to human T1D, but more compact in time, as would be expected, such that the different stages can be precisely identified. Because it is an inbred strain with a predictable outcome (>80% penetrance of diabetes in our colony), the staging is less variable among the mice. Overall, such information may provide insight for interpreting human disease. By conventional flow cytometry analysis and immunoimaging, we and others have identified the cellular landscape of the islets in NOD mice and found that its progression can be parsed into several discrete phases (Jansen et al., 1994; Anderson and Bluestone, 2005; Carrero et al., 2013). We detected the first T cells in the islets at ∼3 wk, as soon as mice were weaned. (Our mouse colony is representative; diabetes develops by 18–20 wk with an incidence of >80%.) T cells are found always in contact with the resident macrophages, much before the development of dysglycemia (Carrero et al., 2013; Mohan et al., 2017). Macrophages are the normal resident myeloid cells of islets, but dendritic cells (DCs) enter at the start of the diabetic process (Anderson and Bluestone, 2005; Jansen et al., 1994; Ferris et al., 2014; Klementowicz et al., 2017). Both present peptides derived from insulin (Ferris et al., 2014). Insulin is the main autoantigen driving this process (Wegmann et al., 1994; Nakayama et al., 2005). A recent study delineated important features of insulin-specific CD4 T cells based on selected gene expression profiles using single-cell quantitative PCR (qPCR), showing plasticity of T cells in islets at the single-cell level for the first time (Gioia et al., 2019). Given that T cells specific to β-cell products make islets receptive for nonspecific T cells to enter (Calderon et al., 2011a), and given the complex interactions between innate and adaptive immune cells that occur in islets, a comprehensive study to analyze islets at both cellular and molecular levels during disease progression from the very early stage is needed.

Our laboratory performed a chronologic examination of transcriptional signatures of NOD islets during diabetes development, showing progressive changes in the expression of genes linked to inflammatory responses (Carrero et al., 2013). Several reports from other laboratories in which examination of islets was performed at selected time points agree in general with our conclusions (Miyazaki et al., 1985; Jansen et al., 1994; Magnuson et al., 2015). Whole-islet transcriptional analysis of highly diverse infiltrating cells, however, is difficult to relate to particular sets of lymphocytes and innate cells. Single-cell RNA sequencing (sc-RNASeq) is a powerful technology enabling unbiased identification of new cell subsets and their activation programs. Here, we used this tool to understand the diversity of immune cell subpopulations and their transcriptional heterogeneity at three critical time periods. T cells show a complex cellular composition from the very early stage, including subsets with memory and cytotoxic activity harboring finely tuned regulatory mechanisms. Considering the importance of the innate myeloid cells in the initiation of autoimmunity, we placed a major emphasis on the analysis of DCs and islet-resident macrophages.

We performed sc-RNASeq on immune cells, endothelial cells, and mesenchymal cells isolated from pancreatic islets (Fig. 1 A). The leukocytes (CD45⁺) and vascular (CD45⁻CD31⁺) and islet mesenchymal cells (CD45⁻CD90⁺; Chase et al., 2007; Carlotti et al., 2010) were subjected to 10X Genomics pipeline barcoding, library preparation, and sequencing (Zheng et al., 2017; Fig. 1 A and Fig. S1 A). We chose three time points, reflecting the major steps in autoimmune development: 4, 8, and 15 wk of age. These three stages corresponded to a progressive increase of intra-islet leukocyte infiltration (Fig. 1 B). 4 wk is about the first time that the first infiltrating T cells can be identified; 8 wk represents a time when leukocyte infiltration is prominent in most islets, still with no evidence of dysglycemia; and 15 wk is just before the time when clinical diabetes becomes evident. Thus we cover the early stage, the apparent “control” stage, and the advanced prediabetic period (Jansen et al., 1994; Carrero et al., 2013).

We identified the major immune and nonimmune cell populations (Fig. 1 C and Fig. S1 B) using a graph-based clustering approach and t-distributed stochastic neighbor embedding (t-SNE) dimensionality reduction implemented in Seurat package (Butler et al., 2018; Stuart et al., 2019). Cell identity was assigned based on the expression of hallmark genes (Fig. 1 D and Fig. S1 C), confirmed in the ImmGen database to be the appropriate markers (Fig. S1 D). In total, we obtained 16,702 leukocytes and 28,414 nonimmune cells at the three combined time points. In this report, we focus on leukocytes.

The major populations of CD45⁺ cells included T cells, B cells, macrophages, conventional DCs (cDCs), and plasmacytoid DCs (pDCs; Fig. 1 C). Among these, macrophages are the only cells resident in islets since birth (Calderon et al., 2015). At all ages, most of the intra-islet infiltrate was represented by T cells, B cells, and cDCs (Fig. 1, E and F). These were present as early as 4 wk of age and exhibited strong expansion as diabetes progressed. Importantly, examination of islets from 2-wk-old mice showed only the resident macrophages. The presence of the major intra-islet populations was validated by flow cytometry (Fig. 1 G). The pattern of infiltration identified by sc-RNASeq recapitulated the flow cytometry data (Fig. 1 H).

The intra-islet lymphoid population in addition to B cells and CD4/CD8 α/β T cells included a number of minor populations: γ/δ T cells, natural killer (NK), and NKT cells and innate lymphoid cells (ILCs), characterized by expression of cell-associated markers (Fig. 2, A and B). The γ/δ T cells were identified based on the coexpression of genes encoding T cell receptor γ constant region 1, 2, and Cd3e (Fig. 2 C). Additionally, NK cells (positive for Eomes, Ncr1, encoding NKP46, Klra8, encoding Ly-49h, and Klra1) and NKT cells (positive for Cd3e, Klra1; negative for Cd4, Cd8a, Klra8, Eomes) were found (Fig. 2 C). While the CD4 and CD8 T cells showed strong expansion, the minor lymphoid cells expanded to a very minor degree, more noticeable at the late 15-wk stage (Fig. 2 D). The presence of some of these minor populations, such as ILC2, ILC3, γ/δ T cells, and NK cells, was confirmed by flow cytometry (Fig. 2, E and F) and was also reported in other studies (Falcone et al., 2004; Markle et al., 2013; Dalmas et al., 2017).

As shown in Fig. 1 E, a small number of B cells were found in islets at 4 and 8 wk, increasing by week 15. We also identified a low number of plasma cells (highly positive for genes encoding immunoglobulin constant regions of heavy chains and J chain; see Fig. S1, C and D) along with cycling B cells (Fig. 2 G), having up-regulated cell cycle–related transcriptional programs (Fig. 2, H and I). B cells are important in the diabetic autoimmune process (Serreze et al., 1996), most likely by presenting islet antigens (Noorchashm et al., 1999; Silveira et al., 2002; Wan et al., 2020) and recognizing insulin, the major protein of the β cells (Kendall et al., 2007; Henry-Bonami et al., 2013; Felton et al., 2018). In sum, from the very start of the diabetogenic process, islets harbor a range of various infiltrating leukocytes.

We analyzed 2,649 CD4 T cells pooled from the three time points. Leveraging annotated immunologic gene pathways (see Materials and methods) allowed us to dissect them into seven distinct populations: CD4-0, CD4-1,4,5, CD4-2, CD4-3, CD4-6, CD4-7, and CD4-8 (Fig. 3 A and Fig. S2 A). (Populations were referred to by a number with one or two marker genes in parentheses.) These sets differed by expression of a number of markers corresponding to various T cell activation states (Lef1, Ccr7, Ikzf2, Foxp3, Pdcd1, etc.; Fig. 3 B and Fig. S2 B) and were identified based on their unique transcriptional signatures (Table S1), as naive, effector, memory, regulatory, and anergic sets (Fig. 3 C; see also Materials and methods and Fig. S2 C). These various sets were clearly distinguished over the time course (Fig. 3 D).

Although in small numbers, the CD4 T cells at 4 wk were already heterogeneous. Their level increased from the week 8 to week 15, but varied among the various populations (Fig. 3 E, left panel). The salient issue to emphasize is that during the three time points, the islets contained regulatory/anergic cells together with effector/memory ones.

The CD4-0(Lef1) (expressed Lef1) was considered a “naive” set, which may well include nonspecific T cells that enter inflamed, “receptive” islets (Calderon et al., 2011a). The effector cycling set CD4-1,4,5(Ptma) was consistently found at all times (Fig. S2 D and Fig. 3 E).

The regulatory population consisted of two distinct sets: CD4-3(Foxp3) and CD4-6(Lag3, IL10). Both expressed Tnfrsf18 (glucocorticoid-induced tumor necrosis factor receptor [GITR]), Ikzf2, Arl5a (Fig. S2 E). The CD4-3(Foxp3) represented classic Foxp3 positive regulatory T cells (Fig. 3 F), while CD4-6(Lag3, IL10) was negative/low for Foxp3 and coexpressed genes encoding IL-10 and IFN-γ (Fig. 3 B). This latter set corresponds to the type 1 regulatory cells (Tr1; Roncarolo et al., 2006). Tr1 cells produce high levels of an anti-inflammatory cytokine IL-10, as well as IFN-γ, while undergoing activation in an antigen-specific manner (Chihara et al., 2016; Clemente-Casares et al., 2016). To summarize, the two populations potentially controlling the autoimmune process were present in the islets at all stages, even at 4 wk of age.

We also identified an anergic set of CD4 T cells, group CD4-2(Tnfsf8). This set was identified by the expression of genes encoding negative regulators of immune response (Ctla4, Pdcd1, Lag3, Lgals1, Cd200) and also FR4 (Folr4/Izumo1r; Fig. 3 G), markers associated with CD4 T cell anergy (Crespo et al., 2013; Kalekar et al., 2016). This anergic set expanded by the 15-wk time point, highlighting another control established in the islets during progression (Fig. 3 E).

The pathogenic CD4 T cells were represented by the effector/memory population, meaning cells that experienced previous antigenic stimulation. These sets increased at week 15, when β-cell mass starts decreasing (Fig. 3 E, left). This effector/memory population consisted of two transcriptionally distinct groups: CD4-7(Itgb1,Ccr7) and CD4-8(Itgb1, Cxcr6). Both were positive for memory-associated genes Il7r and Klf2 (Fig. 3 B and Fig. S2 E; Kondrack et al., 2003; Weinreich et al., 2009). Both also expressed Itgb1, the gene encoding integrin β1 (CD29), part of the VLA-4 complex (α4β1) that binds VCAM1 molecule on inflamed islet endothelium (Calderon et al., 2011b). However, on the basis of expression of Ccr7 gene and other activation markers, these two subsets were separated into the T effector memory phenotype (Cxcr6⁺, Ly6c1⁺, and Id2⁺) and T central memory phenotype (Ccr7⁺; Sallusto et al., 1999; Fig. 3 H). Separation of the two memory populations and the presence of naive, regulatory, and effector subsets were validated at the protein level based on markers found in the current study (Fig. S2, F and G). In sum, CD4 T cells in islets show marked transcriptional heterogeneity, with mixtures of effector, regulatory, anergic, and naive T cells all at various stages, most likely influencing each other.

Analysis of the 1,873 CD8 T cells pooled from three time points revealed eight populations (CD8-0, CD8-1, CD8-2, CD8-3, CD8-4, CD-5,7, CD8-6, and CD8-8; see Fig. 4 A). These sets were characterized by expression of genes corresponding to different activation states (Fig. 4 B and Table S2). The various patterns were based on the population-specific gene signatures (Fig. S3 A and Table S2). CD8-3(Junb) was the only cluster defined by <15 genes (13) and was not further characterized. The various sets included early effector CD8-5,7(Hmgb2) and CD8-1(Xcl1); and late effector (cytotoxic) CD8-2(Gzma); one naive state, CD8-6(Lef1); and two exhausted states, CD8-8(Tcf7, Tox) and CD8-0(Lag3, Pdcd1) (Fig. 4 C; see also Fig. S3, B and C).

As with the CD4 T cells, the CD8 T cells at 4 wk of age were heterogeneous and, to our surprise, already contained the CD8-2(Gzma) population (Fig. 4, D and E). This set represented the CTL group positive for cytotoxic molecules, such as Gzma, Gzmb, Gzmk, Klrc1, Klrc2, Klre1, and Klrk1 (Fig. 4 F). Similarly to CD4 T cells, the presence of the CTLs at the early time point indicates that the cells already contained a stigma of activation, perhaps acquired via priming in the peripheral lymphoid tissue, most likely in the draining pancreatic node (Gagnerault et al., 2002). Cycling effector CD8-5,7(Hmgb2) and naive CD8-6(Lef1) were present at this stage as well (Fig. 4, D and E).

At the 8- and 15-wk time points, the CD8 T cells expanded 7- and 22-fold, respectively, along with dramatic changes in composition (Fig. 4 E). These changes were characterized by the emergence of the early effector CD8-1(Xcl1) and memory CD8-4 (positive for Il7r, Cd28, Ifng, and Cd40lg; Fig. 4 G and Table S2) and, at the same time, the exhausted cells. The CD8-1(Xcl1) was similar to the effector cycling group CD8-5,7(Hmgb2) in expressing genes linked to the Myc and MTORC1 pathways but lacking the gene expression associated with the cell-cycle phenotype (Fig. S3, D and E).

The exhausted group consisted of two populations: CD8-0(Lag3, Pdcd1) and CD8-8(Tcf7, Tox). Both were positive for markers associated with the exhausted phenotype: Pdcd1 (PD-1), Lag3, and Tox (Fig. 4, B and H; Yao et al., 2019), but were differentiated by the Tcf7 transcript, which was strongly up-regulated in CD8-8(Tcf7, Tox) (Fig. S3 F) and showed correlation with this phenotype (Fig. S3 G). The Tcf7 gene encodes transcription factor TCF1, previously shown to mark exhausted CD8 T cells that potentially give rise to cells with effector functions in response to checkpoint blockade (Im et al., 2016; Utzschneider et al., 2016; Wu et al., 2016; Siddiqui et al., 2019).

Strikingly, the only set that did not increase at 15 wk was the CTL group, CD8-2(Gzma) (Fig. 4 E, left). Moreover, its abundance relative to the total CD8 cohort dramatically dropped (Fig. 4 E, right), indicating that expansion of this population was under strict regulation.

In sum, heterogeneity was also found for CD8 T cells. Surprisingly, CD8 CTLs were already in islets at the 4-wk period but did not increase. At the time of effector reactions (15 wk), they represented a minor population. Note that we found the exhausted population that expresses Tcf7.

We analyzed 2,031 cDCs pooled from 4-, 8-, and 15-wk-old NOD mice and identified two major populations: cDC1 and cDC2 (Fig. 5 A). Both were present at 4 wk and expanded during progression (Fig. 5 B, left). cDC2 was the predominant set at all stages (Fig. 5 B, right). cDC1 expressed genes corresponding to the Batf3-dependent CD103⁺ cDC1: Xcr1, Cd24a, Irf8, and Batf3 (Fig. 5 C; Murphy, 2013; Grajales-Reyes et al., 2015). This set was shown to be indispensable for diabetes initiation in NOD mice (Ferris et al., 2014).

The cDC2 group expressed the Sirpa gene (Fig. 5 C) but was heterogeneous, comprising three populations: cDC2(Ccr7), cDC2(Mgl2), and cDC2(Ltb). High expression of Ccr7 in cDC2(Ccr7) (Fig. 5, D and E) marked the CCR7⁺ cDC2 group with potential migratory activity from the sites of inflammation to the draining lymph node (Dieu et al., 1998; Allavena et al., 2000). The cDC2(Mgl2) set uniquely expressed Ccl17 and Mgl2 (encoding CD301b protein), pointing to the CD301b⁺ cDC2 (Kumamoto et al., 2013), which shares some features with monocytes and macrophages, such as expression of genes encoding Fc receptors (Fig. 5 D). This group did not belong to the macrophage population, as it was negative for Emr1(F4/80), expressed by islet macrophages, and did not express the transcription factor Mafb (Fig. 5 F). The last group, cDC2(Ltb) represented the largest population over the time course. These cells had elevated levels of Havcr2 gene encoding Tim-3, potential negative regulator of cDC activation (Chiba et al., 2012).

In addition to the cell type heterogeneity, we identified a superimposed inflammatory signature (see Materials and methods), manifested by a number of IFN-γ–inducible genes, such as Cxcl9, Ly6a, and Gbp2 (Fig. 5 G). Such a signature was identified in a fraction of both the cDC1 and cDC2 subsets, reflecting the ability of all cDCs to equally respond to the autoimmune milieu of the islets (Fig. 5 H). The proportion of these activated cDCs was similar between cDC1 and cDC2 and increased as diabetes progressed over time (Fig. 5 I). The presence of both cDC1 and cDC2 populations was confirmed by flow cytometry (Fig. 5 J). To summarize, the analysis of DCs indicated the presence of cDC1 and newly characterized cDC2 subsets, both responding to inflammatory cues.

Macrophages are the only resident myeloid cells found under normal conditions in islets of all mouse strains. In steady state, macrophages support the function of endocrine cells (Banaei-Bouchareb et al., 2004; Geutskens et al., 2005; Calderon et al., 2015). In NOD mice, they have an indispensable role in the initiation of the autoimmune process, likely by controlling the initial entrance of T cells: their ablation results in strong protective effects (Carrero et al., 2017). In our previous study, we identified the first inflammation-induced changes in macrophages (Ferris et al., 2017) by 3 wk of age.

Here we expanded our transcriptional analysis, performing whole-population RNASeq on FACS-purified macrophages from 2-, 4-, 8-, and 12-wk-old NOD mice. The rationale was to get deep transcriptional information, often not available from sc-RNASeq. To distinguish the immune response from age-related changes, we complemented our analysis with samples from nondiabetic age-matched control C57BL/6J (B6) mice (Fig. 6 A).

Principal component analysis (PCA) revealed progressive changes in the transcriptomes, matching disease progression (Fig. 6 B). We found 1,216 differentially expressed genes up-regulated in the 12-wk sample compared with the 4-wk sample (Fig. 6 C). Many of the elevated responses were downstream of IFN-γ signaling, IL-2/STAT5, and IL-6/STAT3 axis (Fig. 6 D). The transcripts with elevated expression corresponded to genes involved in antigen presentation and interaction with T cells, such as genes encoding MHC-II complex (H2-Ab1, H2-Aa, H2-Eb1), MHC-I (H2-Q4, H2-Q6, H2-K1, B2m), costimulatory molecules (Cd40, Cd86), inflammatory cytokines (Ccl2, Ccl5, Ccl7, Ccl8, Cxcl9, Cxcl13), negative regulators of immune response (Cd274/PD-L1, Lag3), and transcription factors (Ciita, encodes MHC-II transactivator, Stat1, Stat3, Hif1a, Klf4, Fos, Fosb). Fig. S4 A shows that differentially expressed genes progressed in time from 4 wk of age. The B6 control group exhibited no autoimmune-related changes (Fig. 6 E).

The sc-RNASeq analysis of 5,975 macrophages revealed transcriptionally distinct clusters at the three time points; five main sets were identified (Fig. 7 A). Each had a unique transcriptional signature and functional markers (Fig. 7 B and Table S3), including markers associated with inflammatory activation (Ccl3, Atf3, Cxcl9, Stat1, Cd40), regulatory and scavenging functions (Il1rn, Lgals1, Lgals3, Cd36, Anxa5), and cell-cycle genes (Stmn1, Mki67).

Mac-1(Apoe) showed limited changes in time (Fig. 7, C and D). It was characterized by elevated levels of Apoe and Trem2 but lacked an inflammatory gene signature (Fig. 7 B). This group consisted of three subsets (Fig. S4 B): its complex composition probably reflected stages of macrophage maturation, based on differences in expression of genes encoding MHC-I, MHC-II, and a number of transcription factors (Fig. S4, C–F).

Mac-5(Stmn1) expressed genes linked to the cell-cycle pathway (Fig. 7 E), indicating the self-replicating capacity of islet macrophages. This finding is consistent with our previous report (Calderon et al., 2015), but now shows that such a replicative capacity is present throughout diabetes progression. We consider the Mac-1(Apoe) and Mac-5(Stmn1) sets to represent macrophages involved in islet functions uninfluenced by an autoimmune response.

We identified two groups of inflammatory activated macrophages, Mac-2(Atf3) and Mac-3(Cxcl9) (Fig. 7, A and B). The Mac-2(Atf3) set did not show major changes over the time course (Fig. 7, C and D). It was characterized by strong NF-κB signaling activation, with elevated levels of NF-κB–inducible transcripts, including genes encoding inflammatory cytokines Ccl3, Ccl4, Cxcl2 (macrophage inflammatory proteins 1-α, 1-β, and 2-α), Atf3, Ccrl2, and Dusp1 (Fig. 7, F and G).

The hallmark feature of Mac-3(Cxcl9), aside from an NF-κB signature (Fig. 7 H), was the elevated expression of genes encoding several cytokines, costimulatory molecules, and molecules involved in antigen presentation. Of note was a high expression of the genes encoding IL12/IL23 subunit p40 (Il12b), CD40, MHC-I (H2-Q6, H2-Q4, H2-K1), and MHC-II molecules (H2-Aa, H2-DMa, H2-DMb1, H2-DMb2; see Fig. 7 I and Fig. S4 F). It also had up-regulated expression of transcripts related to proteasomal peptide processing (Psmb8, Psme2b, Psme2, Psme1, Psmb9) and class I antigen presentation (Tap1, Tap2, Tapbp; Fig. S4 F and Table S3). Importantly, this population was the only one that showed continuous strong expansion during T1D progression: approximately four-fold relative increase from 4- to 15-wk sample (Fig. 7, C and D). Its signature matched the inflammatory transcriptional response found throughout diabetes using whole-population RNASeq (Fig. 6, D and E). Based on the kinetics of this population and its strong inflammatory signature, we identified this group as a highly pathogenic cluster promoting autoimmunity.

This Mac-3(Cxcl9) population comprised two distinct subsets, both having strong up-regulation of IFN-γ–inducible genes (Fig. S4, B and F). The two were distinguished because the Mac3(Cxcl9)-2 expressed a number of genes found in monocytes: Ccr2, Ly6c1, and Plac8; it also had decreased levels of C1qa, C1qb, C1qc, Apoe, and Ms4a7 (Fig. S4 G). The other subset, Mac-3(Cxcl9)-3, lacked the monocytic markers. From these differences, we surmised that Mac3(Cxcl9)-2 was monocyte derived, while in contrast, Mac-3(Cxcl9)-3 developed from the resident macrophages. This observation suggests that during disease progression, monocytes enter islets and may contribute to the inflammatory status.

In contrast to the general inflammatory groups, Mac-4(Prdx1) had a gene signature associated with alternative activation of macrophages: lysosomal activity, extracellular matrix protein expression, and iron processing (Fig. 7 E). This set expressed genes encoding annexins (Anxa1, Anxa2, Anxa4, Anxa5; Fig. 7 B), the scavenging receptor Cd36, and galectin3 (Lgals3). This gene signature strongly points to this set as specialized in efferocytosis, apoptotic cell clearance, an activity that may impose negative responses against inflammation (Arandjelovic and Ravichandran, 2015). Indeed, we also found elevated expression of genes with protective and strong anti-inflammatory function such as Prdx1 encoding peroxiredoxin 1, an antioxidant enzyme (Park et al., 2017); the IL1 receptor antagonist (Il1rn), and galectin 1 (Lgals1) and galectin 3 (Lgals3), which have potent anti-inflammatory activity (MacKinnon et al., 2008; Perone et al., 2009; Volarevic et al., 2010). Low in numbers from the beginning of inflammation, Mac-4(Prdx1) was the only set that further decreased in proportion as autoimmunity moved forward (Fig. 7, C and D).

Macrophages and DCs responded to inflammation in the islets, and the transcriptional programs imprinted by the autoimmune responses were similar between them (Fig. 7 J). Both cells strongly up-regulated IFN-γ–inducible inflammatory transcripts (Stat1, Cxcl9, Cxcl10, Cd40, Cd274/PD-L1), but also had a number of genes specifically elevated in each of them (Fig. 7 J). In particular, macrophages up-regulated expression of several genes encoding inflammatory cytokines and MHC-I and MHC-II complexes, suggesting their nonredundant function in T1D.

We examined islets from two variant NOD mice: NOD.Rag1^−/− mice lack mature lymphocytes and consequently are free of diabetes; NOD.Ifngr1^−/− mice lack the gene encoding the type II interferon receptor. IFN-γ is a major macrophage-activating cytokine. In our colony, NOD female mice lacking the type II interferon receptor signaling have reduced disease incidence, along with a diminished frequency of the islet-infiltrating T cells (Carrero et al., 2018). To examine these additional samples, we supplemented the NOD time-course data with the samples from these knockout mice and reclustered the merged sample (Fig. S5 A). The identification of previously described subpopulations was confirmed by the gene markers (Fig. S5 B).

Examining the macrophages from the islets of the NOD.Rag1^−/− mice was instructive. The two main sets were represented by Mac-2(Atf3) and Mac-1(Apoe) (Fig. 8 B). The finding of Mac-2(Atf3) in NOD.Rag1^−/− mice confirms that the basal activation program is independent of an autoimmune process (Ferris et al., 2017). Two other sets at lower levels included the anti-inflammatory group Mac-4(Prdx1) and a very low level of the cell-cycling group Mac-5(Stmn1).

Mac-3(Cxcl9), the second proinflammatory set, resulted from autoimmunity, because it was completely absent in the NOD.Rag1^−/− sample. Furthermore, this population was not found in NOD.Ifngr1^−/− mice (Fig. 8 B), demonstrating its dependence on IFN-γ signaling.

To obtain a deeper insight into the transcriptional path, corresponding to the changes in macrophages, we performed pseudotime analysis using Monocle3 software (Trapnell et al., 2014; Fig. 8 C). This approach allowed us to deduce the transcriptional trajectory of macrophages at the single-cell level as they progressively changed their gene expression profile during the autoimmune process. The 2D uniform manifold approximation and projection (UMAP) plot revealed transitions among the various sets during diabetes development (Fig. 8 C).

Progression along pseudotime started with the cell cycling group Mac-5(Stmn1) (see Materials and methods), and then through the other sets found in nondiabetic mice, ending with the most distant from the entry point group, represented by the autoimmune-specific Mac-3(Cxcl9) (Fig. 8 C). This pattern was also seen by the expression of marker genes (Fig. 8 D).

Furthermore, the pseudotime analysis revealed a stepwise pattern of macrophage inflammatory activation, consisting at the cellular level of two consecutive stages represented by two proinflammatory populations (Fig. 8 C): stage 1 was represented by Mac-2(Atf3) and stage 2 by Mac-3(Cxcl9). Based on pseudotime analysis, Mac-3(Cxcl9) belongs to the terminal stage of activation.

The result of pseudotime alignment was highly consistent with the real-time changes as mice progressed to the terminal stage of inflammation (Fig. 8 E). The macrophages from 4-wk-old NOD mice were mostly represented by Mac-1(Apoe) and Mac-2(Atf3) (corresponding to stage 1 activation) populations, which were also present in NOD.Rag1^−/− (Fig. 8 E). In striking contrast, macrophages from the 15-wk-old NOD mice were mostly concentrated in stage 2 activation. The 8-wk-old sample reflected the intermediate state between the two extremes. Macrophages from NOD.Rag1^−/− and NOD.Ifngr1^−/− contained virtually no cells belonging to the stage 2 activation group (Fig. 8 E).

The pseudotime analysis of macrophage activation revealed distinct groups of genes driving sequential inflammatory states (Fig. 8 F). Cell-cycle genes were expressed at the start and IFN-γ–inducible genes at the end of the trajectory. Moreover, a comparison between populations of macrophages representing the two stages of activation (Fig. 8 G) revealed detailed characteristics of each. Stage 1 activation, represented by Mac-2(Atf3), had an NF-κB activation signature that was not triggered by autoimmune response. The autoimmune-driven stage 2 (represented by Mac-3(Cxcl9)) was characterized by the additional up-regulation of IFN-γ–inducible genes (Cxcl9, Ly6a, Cd40, Ccl5, etc.).

To summarize, the devised trajectory revealed the stage-wise progression pattern of macrophages during autoimmune development and was confirmed by the real-time biological time course. We defined a program of macrophage activation and progression based on a computational pseudotime approach, biological time-course data, and evaluation of two gene knockouts. This analysis defined distinct functional groups of macrophages in pancreatic islets, and also established stages of their activation and subpopulation remodeling as a function of autoimmune diabetes development.

Single-cell analysis of diabetic autoimmunity identified cellular diversity and transcriptional heterogeneity during the various stages of this chronic and persistent autoimmune process. To our surprise, islets at the initial stage of the autoimmune process contained various effector and regulatory CD4 and CD8 T cells, albeit at a low level. This was followed by an apparent control phase, with the expansion of nonpathogenic and regulatory populations, ending with a prediabetic stage. The single-cell analysis adds considerable evidence to support each stage and shows that each is characterized by a complex mixture of cells with marked transcriptional heterogeneity. A single reactive cell is likely to go through various transcriptional changes as it is influenced by the complex metabolic changes in islets throughout diabetes. Some have distinctly different functions in an apparent competition among them.

This study provides insights that complement and expand information previously derived from the examination of human islets. Importantly, as in NOD, myeloid and lymphoid populations were present in islets of recent T1D patients (Damond et al., 2019; Wang et al., 2019), indicating that the immune cell complexity is equally represented. A hallmark in T1D patients and in NOD mice is the simultaneous recruitment of both CD4 and CD8 T cells into the islets, supporting that islet autoreactivity is not driven by a dominant cell type but originates from cooperative activities among different immune cells. But the additional and highly unique information that the NOD model provides are the findings at the very early stage, long before any overt clinical signs. Here the NOD model shows the active autoimmune process very early in development. This 4-wk stage defines the start of a diabetes susceptibility period characterized by a low level of β-cell apoptosis and more demands for insulin secretion, likely resulting from a growth spurt in the mouse (Bonner-Weir, 2000).

Several conclusions were derived from the transcriptional signatures of the T cells. First is that CD4 T cells with effector/memory phenotype were found at the initial 4-wk stage. Some of the entering CD4 T cells already exhibited a stigma of activation, likely caused by their circulation in the peripheral lymphoid organs, as was just shown for T cells to insulin (Wan et al., 2018). The peripheral sensitization depends on the exocytosis of diabetogenic antigens released from islets (Wan et al., 2018). The crucial role of peripheral priming has been shown in various reports (Höglund et al., 1999; Calderon et al., 2014), including those depicting the profound effects of ablation of lymph nodes (Gagnerault et al., 2002; Levisetti et al., 2004). Alternatively, and not exclusively, one can posit that T cells entered islets as naive but rapidly changed their biology upon reaction with the islet microenvironment, one highly rich in bioactive molecules.

A second finding to note is that the mechanisms controlling the autoimmune process were already evident at the very initial stage accompanying the first effector response. The regulatory T cells were represented by two different subsets: Foxp3⁺ regulatory T cells and Il10⁺Ifng⁺ Tr1 cells, which likely restricted the pathogenic T cells. Such a complex control later manifested itself in the accumulation of anergic CD4 and exhausted CD8 T cells. The control of islet-specific CD8 T cells was previously proposed (Green et al., 2003). Indeed, the cytotoxic cells detected at the initial stage were the only CD8 T cell population that did not undergo expansion over the time course. These findings show the presence of the regulatory self-protecting mechanisms in the islet microenvironment starting from the very initial stage of inflammation, and are supported by the observation that ablation of regulatory molecules greatly accelerates disease progression (Lühder et al., 1998; Keir et al., 2006). We note that the exhausted CD8 T cells contained a Tox⁺Tcf7⁺ subpopulation previously described as a stem-like cell sustaining long-term immune response with an ability to undergo rapid proliferation upon checkpoint blockade (Alfei et al., 2019; Im et al., 2016; Siddiqui et al., 2019; Utzschneider et al., 2016; Wu et al., 2016). These may be the CD8 T cells inducing diabetes in PD-1–deficient NOD mice (Keir et al., 2006).

One limitation of this study is that we have no information on T cell specificity. However, two recent reports indicate that the early T cell response in the islets of NOD mice is dominated by T cells recognizing autoreactive insulin B chain 12–20 epitope. One, from Teyton’s group (Gioia et al., 2019), was based on single-cell qPCR of insulin-specific tetramer-positive CD4 T cells. A second, from our group (Wan et al., 2020), identified the peptides in islets bound to MHC-II and the corresponding CD4 T cells. Interestingly, the dominance of insulin-reactive CD4 T cells in the islets was more pronounced in the early to middle stages (6 wk), in contrast with the later prediabetic stage (12 wk). This indicates that at the beginning of the autoimmune process, insulin autoreactive CD4 T cells prevail over T cells with other specificities, making islets more receptive to the nonspecific ones (Calderon et al., 2011a). At later stages, T cells with more diverse specificity were recruited, similar to the broad T cell specificity in islets of T1D individuals (Babon et al., 2016). Also, novel antigens might be released at the later stages, when more and more β-cells experience stress and undergo apoptosis. As for the CD8 T cells, various reports identified T cells directed to peptides from the B chain of insulin (Wong et al., 1999, 2002) as well as to the islet-specific glucose-6-phosphatase catalytic subunit-related protein (Lieberman et al., 2003; Han et al., 2005). Both are likely to be a component of the CD8 T cells found here in islets.

We have paid much attention to the islet macrophages, the only resident myeloid cell found in all islets, including those from nondiabetic mice, and here tested in NOD.Rag.1^−/− mice. Combining the whole-population RNASeq approach with sc-RNASeq enabled us to identify pathogenic transcriptional programs and subsets of macrophages within a two-step program of proinflammatory activation. The activated macrophages contained multiple populations, including the cycling cells and a novel subset, Mac-4(Prdx1). The macrophages of the stage 1 activation were present in the steady-state (in NOD.Rag1^−/−) and comprised the Mac-2(Atf3). The Mac-2(Atf3) was characterized by activation of the NF-κB pathway. Previously, we showed that islet macrophages produced inflammatory cytokines at gene and protein level; in particular, TNF (Ferris et al., 2017). Because TNF receptor signaling is one of the major inducers of NF-κB activation, our data suggest that the stage 1 activation might be driven via autocrine/paracrine TNF signaling. Importantly, NOD mice deficient of TNF receptor 1 (NOD.Tnfr1^−/−) were completely protected from T1D (Chee et al., 2011), suggesting the indispensable role of the initial stage 1 activation for diabetes initiation.

A novel population also present at steady state, Mac-4(Prdx1), was most likely involved in efferocytosis based on the gene expression pattern. A small number of apoptotic β-cells are present in the islets, especially during the neonatal period (Scaglia et al., 1997; Trudeau et al., 2000). This set most likely has an important homeostatic role in the endocrine tissue under noninflammatory conditions. Efferocytosis has been associated with an anti-inflammatory program (Poon et al., 2014).

Stage 2 of macrophage activation was autoimmune induced and up-regulated both IFN-γ and NF-κB signaling pathways. Our analysis also revealed IFN-γ as an important signal, driving terminal pathogenic (stage 2) macrophage activation. Interestingly, knockout of the IFN-γ receptor delays diabetes in female mice and reduces islet leukocyte infiltration levels but does not completely prevent T1D (Carrero et al., 2018). The explanation might be that the lack of type II interferon signaling was partially compensated by type I interferon activation, since double knockout of the genes encoding receptors to both type I and type II interferons has a synergistic effect (Carrero et al., 2018).

The myeloid compartment also included DCs, represented by both cDC1 and cDC2. The cDC1 subset was minor but was important for the initiation of T1D (Ferris et al., 2014). The cDC2 group shared inflammatory responses with cDC1 (Fig. 5 G) but had about a 10× higher frequency (Fig. 5 B), raising an important question about its role. The cDC2 group was heterogeneous and included a CCR7⁺ cDC subset, known to be migratory from sites of inflammation to the draining lymph node (Jang et al., 2006).

The autoimmune activation program was similar between cDCs and macrophages (Fig. 7 J), indicating that they adapt to the same autoimmune microenvironment. However, macrophages exhibited a number of inflammation-associated genes not activated in DCs, indicating a nonredundant activation program and probably specialized functions of macrophages.

Finally, during the late advanced stage, the control mechanisms are superseded, and dysglycemia as a measure of reduced β-cell mass becomes evident. The cellular analysis is giving us some clues on the process that favors progression. It shows three major changes: (1) a relative reduction in the sets of regulatory cells; (2) a striking increase in the number of activated macrophages, the Mac-3(Cxcl9) set; and (3) an increase in the number of cDC2 in the islets. It was clear that while the two CD4 regulatory T sets were relatively high in comparison to the other sets at 4 and 8 wk, by the advanced 15 wk, they had been superseded by diverse populations with potential pathogenicity. At the same time, an exhausted CD8 T cell set became prominent. But not to be ignored is the progressive increase of B cells in the islets of NOD mice. Higher levels of CD20^high B cell infiltration have been shown in islets of patients having accelerated disease (Leete et al., 2016). What determines this profound change from 8 to 15 wk is important to determine: perhaps a broadening of T cell specificities resulting from epitope spreading; a change in the islet bioactive molecules; a slow but progressive increase in the plethora of proinflammatory molecules; or an increase in the antigen presentation capacity as more DCs and macrophages are present. Overall, an intricate balance between the pathogenic and the regulatory elements is a hallmark of the diabetic autoimmune process, which is reflected by the vast degree of cellular and molecular heterogeneity in the target organ from initiation to progression.

Transcriptional data generated in this study was deposited in National Center for Biotechnology Information Gene Expression Omnibus database under accession no. GSE141786.

C57BL/6J mice (B6), NOD/ShiLtJ (NOD), and NOD.129S7(B6)-Rag1^tm1Mom/J (NOD.Rag1^−/−) were originally obtained from the Jackson Laboratory. NOD.Ifngr1^−/− mice were generated as described before (Carrero et al., 2018). All mice were bred and maintained in our pathogen-free animal facility. All mouse experiments were performed in accordance with the Division of Comparative Medicine of Washington University School of Medicine (Association for Assessment and Accreditation of Laboratory Animal Care accreditation no. A3381-01). Only female mice were used in the current study.

Pancreatic islets were isolated as previously described (Calderon et al., 2015) and dispersed into single-cell suspension using nonenzymatic Cell Dissociation Solution (Sigma-Aldrich) for 3 min at 37°C. To block the engagement of Fc-receptors, the cell suspensions were incubated at 4°C for 15 min in PBS (pH 7.4) supplemented with 1% BSA and 50% of FC-block (made in-house). For surface staining, cells were incubated with fluorescently labeled antibodies (1:200 vol/vol; Table S6) supplemented with Fixable Viability Dye eFluor 780 (1:1,000 vol/vol) at 4°C for 20 min. Cells were then washed and analyzed by flow cytometry or subjected to FACS. For intracellular transcription factors, islet cells were stained with surface antibodies for 20 min, fixed, and permeabilized using Foxp3/transcription factor staining buffer set (Thermo Fisher Scientific) following the manufacturer’s instructions. After that, cells were stained with FoxP3 (or RoRγt, GATA3) antibodies at 4°C for 30 min. All flow cytometry was done on a BD FACSAria II (BD Biosciences) and analyzed using FlowJo software (TreeStar).

Islet macrophages were FACS purified from islets of 6–14 NOD or C57BL/6J mice per one replicate, three to four biological replicates per condition. Approximately 100–150 islets were collected per mouse. Macrophages were sorted as Live⁺CD45⁺ B220⁻CD90⁻CD11c⁺F4/80⁺, and approximately 800–1,500 macrophages were collected per biological replicate. Then, mRNA was isolated with RNAqueous-MicroKit (Thermo Fisher Scientific), and a SeqPlex RNA Amplification kit (Sigma-Aldrich) was used for cDNA preparation following the manufacturer’s instructions. Illumina NovaSeq-2500 was used for sequencing. Libraries were prepared and sequenced at the Genome Technology Access Center (Washington University in St. Louis, MO). Raw reads were aligned with the STAR aligner (Dobin et al., 2013; Table S6), counts were generated with htseq-count utility from the HTSeq Python library (Anders et al., 2015), and differential gene expression was performed using DESeq2 R package (Love et al., 2014) with Benjamini–Hochberg adjusted P = 0.05 as a threshold. Heatmaps were generated with GENE-E software. PCA was based on the top 5% genes ranked by variance.

For sc-RNASeq analysis, live CD45⁺, CD31(Pecam1)⁺, and CD45⁻CD90(Thy1)⁺ cells were sorted using islet cells from 4–12 mice per one sample. After sorting, cells were pelleted and resuspended in 10% FBS (HyClone) in PBS (pH 7.4) at 10³ cells/µl and loaded onto the Chromium Controller (10X Genomics; Zheng et al., 2017). Samples were processed using the Chromium Single Cell 3′ Library & Gel Bead Kit (10X Genomics, v2 and v3) following the manufacturer’s protocol. For NOD 4-, 8-, and 15-wk-old samples and NOD.ifngr1^−/−, the v2 Reagent Kit was used. For NOD.Rag1^−/−, the v3 Reagent Kit was used. The libraries were sequenced on Illumina NovaSeq6000. Raw sequencing data were processed using Cell Ranger (10X Genomics; v2.1.1 and v3.0.2). The CellRanger function “count” was used to align raw reads with the STAR aligner (Dobin et al., 2013) against GRCm38 reference, and after counting nonredundant unique molecular identifiers (UMIs), the single-cell digital-expression matrices were obtained. Library preparation and sequencing was done at the Genome Technology Access Center core facility (Washington University in St. Louis). To ensure reproducibility of the results, islet isolation, FACS, library preparation, and sequencing for two of the time-course samples (4- and 8-wk-old NOD) were done twice in two independent experiments.

The output “cellranger” matrices were loaded into the Seurat package (v2.3.4) environment (Satija et al., 2015; Butler et al., 2018; Stuart et al., 2019) using R (v3.5.1). The UMI barcodes were tagged with the sample ID and merged using the “MergeSeurat” function. The UMI count matrix was normalized with a global-scaling normalization method, “LogNormalize.” The unwanted sources of variation (the UMI coverage, number of ribosomal genes per cell) were regressed out using the “ScaleData” function implemented in the Seurat package by passing variables for regression with the “vars.to.regress” argument. A canonical correlation analysis was run to align datasets using the “RunCCA” function, and calculation of the t-SNE projections and clustering (function “FindClusters”) were performed in the aligned subspace. General cell types were annotated based on the expression of canonical markers, driving separation of clusters (Emr1, C1qa, macrophages; Cd79a, Cd79b, B cells; Cd3e, Cd3g, T cells; Kdr, Pecam1, endothelial cells; and Pdgfrb, Pdgfa, mesenchymal cells; see Fig. S1, C and D). Cellular events identified with overlapping markers between different populations (i.e., Cd3e and Emr1, Cd79a and C1qa, etc.) were considered as artifactual aggregates and removed from the analysis.

The cDCs were first clustered using Seurat as described above. As a result, we found that the inflammatory signal was confounded with the signal driving separation of cDC by subtype. To address this issue, we examined the principal components (PrintPCA function) and identified the top genes contributing to the inflammatory-driven separation: Ifitm3, Cd40, Cxcl9, Ifitm2, Gbp4, Zbp1, Stat1, and Isg15. Then, we reanalyzed the cDC sample, reducing this inflammatory-driven signal, and repeated clustering. To extract the inflammatory signal, we did the opposite, repeating clustering after reducing the signals contributing to the subtype separation.

Unsupervised clustering analysis using Pathway Gene Set Overdispersion Analysis (PAGODA2; Fan et al., 2016) allowed de novo identification of the intra-islet CD4 and CD8 T cell subsets based on coordinated differences of gene expression from annotated pathways across cells in our single-cell sample (Fig. S2 A). For this analysis, we leveraged a publicly available database of immunological gene pathways from the Broad Institute (MSigDB, c7; https://software.broadinstitute.org/gsea/msigdb/index.jsp; Subramanian et al., 2005). We performed PCA dimensionality reduction using the top 15 principal components based on the top 3,000 variable genes from the variance-adjusted expression matrix. The k-nearest neighbor was computed with cosine distance (k = 20), and clusters were identified with a “multilevel community” detection algorithm. The t-SNE projections were calculated with perplexity = 30. For CD4 T cells, clusters 1, 4, and 5 were combined into one based on the coordinated enrichment in effector pathways (Fig. 3 A and Fig. S2 A).

We used Seurat software to carry out differential gene expression analyses, construct violin plots, and plot heatmaps (unless indicated otherwise). Differential expression analysis was conducted using a nonparametric Wilcoxon rank sum test. The gene signatures defining subpopulations inside analyzed cells were generated by differential expression analysis with adjusted P < 0.05 and log2 fold-change >0.3. The resulting gene signatures for CD4 T cells, CD8 T cells, cDCs, and macrophages are listed in Table S1, Table S2, Table S3, and Table S4, correspondingly. Heatmaps were plotted using top (limited to 70 per group) differentially expressed genes using Seurat’s “DoHeatmap” function with following RGB color-scheme: col.low = #0038e6, col.mid = #000000, col.high = #FFFF00. The genes that were present in more than one gene signature were plotted only once on the heatmaps. Differential expression scatter plots were generated using GraphPad Prism software. Hypergeometric pathway analysis (Fig. 2 H, Fig. S2 D, Fig. S3, D and E; and Fig. 7, E, F, and H) was done based on differentially expressed genes with adjusted P < 0.05 and logarithm fold change (LFC) >0.3.

Monocle3 R package (Trapnell et al., 2014; Qiu et al., 2017; McInnes et al., 2018,Preprint; Cao et al., 2019) was used for pseudotime inference analysis. The sc-RNASeq dataset was normalized, dimensionality reduction was done with the top 20 principal components based on the genes driving macrophage heterogeneity (Table S4) with “preprocess_cds” function, and the sample-based batch effect was subtracted from the data. The number of principal components used for downstream analysis was determined based on visual examination of the percentage of variance explained by each component plotted with “plot_pc_variance_explained.” Dimensionality of the data was further reduced with UMAP algorithm with the “cosine” metric. The trajectory graph was learned with the “learn_graph” function using “gamma = 1e-3” and “sigma = 1e-3,” minimum distance 0.18; and the number of neighbors to use during k-nearest neighbor graph construction was set to 20. Once the trajectory graph was learned, cells were ordered in pseudotime with “order_cells” function. The pseudotime algorithm does not have a preference for which position along a trajectory would be a start point. For macrophage trajectory reconstruction, we hypothesized that the most naive state of macrophage, the group of self-renewing macrophages Mac-5(Stmn1), would feature the entry point (root) that feeds cells into the general path (Fig. 8 C). The heatmap showing expression pattern in a trajectory-dependent manner was plotted with the top genes differentially expressed by Moran’s I value, all false discovery rate (FDR) <0.01.

For identification of CD4 and CD8 T cell phenotypes, we cross-referenced the gene signatures of found subpopulations against the corresponding reference bulk transcriptional datasets. The reference datasets, which reflect various T cell phenotypes, were retrieved from the National Center for Biotechnology Information genomics data repository with the following Gene Expression Omnibus accession nos.: GSE89555 (Ibitokou et al., 2018), GSE7852 (Feuerer et al., 2009), and GSE113624 (Alonso et al., 2018) for CD4 T cells; GSE41867 (Doering et al., 2012) and GSE83978 (Utzschneider et al., 2016) for CD8 T cells. First, we ran differential expression analysis on these datasets using Limma R package (Ritchie et al., 2015). The specific samples used for each comparison are listed in Table S5. Then, we generated ranked lists of genes reflecting differences between each pair of conditions in the reference datasets. To do so, we used signal-to-noise statistics as a metric and considered only the top 10,000 genes ordered by mean expression. Finally, gene set enrichment analysis (GSEA; Subramanian et al., 2005) revealed a correlation of the sc-RNASeq subpopulations with previously identified cell phenotypes (Fig. S2 C) based on their unique gene signatures (Table S1 and Table S2).

Fig. S1 shows the main populations of non-endocrine islet cells and their gene markers. Fig. S2 and Fig. S3 show the procedure of identifying CD4 and CD8 T cell subsets, respectively. Fig. S4 details transcriptional characteristics of macrophages at bulk and single-cell levels. Fig. S5 depicts clustering analysis of macrophages pooled from NOD wild-type samples, NOD.Rag1^−/− and NOD.Ifngr1^−/−. Table S1, Table S2, Table S3, and Table S4 list the differentially expressed genes among CD4 T cell, CD8 T cell, cDC, and macrophage subsets, respectively. Table S5 shows the design of differential gene expression analysis based on publicly available datasets, which were used as references for CD4 and CD8 T cell annotation. Table S6 lists the software, reagents, and transcriptional datasets used in the study.

We thank Katherine Fredericks and Orion Peterson for their assistance in the maintenance of our mouse colony and islet isolation. We benefited from much discussion with all our laboratory members including Anthony Vomund, Neetu Sirivastava, Mohammed Zaman, and Cheryl Lichti. We thank the Genome Technology Access Center and Flow Cytometry & Fluorescence Activated Cell Sorting Core, Washington University in St. Louis, School of Medicine.

Our work is supported by National Institutes of Health grants AI114551 and NK508177 and by a general support grant from the Kilo Diabetes & Vascular Research Foundation.

Author contributions: P.N. Zakharov was the main contributor to the design of the various experiments and for the single-cell RNAseq analysis. All authors contributed to the interpretation of the data. H. Hu and X. Wan performed some of the cellular experiments. E.R. Unanue, P.N. Zakharov, X. Wan, and H. Hu wrote, revised, and edited the manuscript.