Single-Cell Profiling Reveals Divergent, Globally Patterned Immune Responses in Murine Skin Inflammation.

Inflammatory response heterogeneity has impeded high-resolution dissection of diverse immune cell populations during activation. We characterize mouse cutaneous immune cells by single-cell RNA sequencing, after inducing inflammation using imiquimod and oxazolone dermatitis models. We identify 13 CD45+ subpopulations, which broadly represent most functionally characterized immune cell types. Oxazolone pervasively upregulates Jak2/Stat3 expression across T cells and antigen-presenting cells (APCs). Oxazolone also induces Il4/Il13 expression in newly infiltrating basophils, and Il4ra and Ccl24, most prominently in APCs. In contrast, imiquimod broadly upregulates Il17/Il22 and Ccl4/Ccl5. A comparative analysis of single-cell inflammatory transcriptional responses reveals that APC response to oxazolone is tightly restricted by cell identity, whereas imiquimod enforces shared programs on multiple APC populations in parallel. These global molecular patterns not only contrast immune responses on a systems level but also suggest that the mechanisms of new sources of inflammation can eventually be deduced by comparison to known signatures.

Introduction

The term “inflammation” serves as an abstraction that summarizes the activity of dozens of types of immune cells, whose global reprogramming by exogenous antigens counteracts pathologic intrusions. Until recently, the heterogeneity of the immune response has been addressed by applying flow cytometry and immunohistochemistry to interrogate specific subpopulations such as innate lymphoid cells (Sonnenberg and Artis, 2015), regulatory B cells (Rosser and Mauri, 2015), and regulatory T cells (Malhotra et al., 2018). Although these studies assigned functional repertoires to cellular subtypes, their scope limited systematic comparisons between cell types and the discovery of rare, functionally active cell subpopulations. The recent emergence of single-cell sequencing approaches has enabled the characterization of multiple, divergent immune cell populations in parallel (Villani et al., 2017).

Active inflammation typifies the skin, both in host defense and pathogenic processes that cause rashes. Studies using single-cell RNA-sequencing (scRNA-seq) have recently refined the contributions of T cells (Hughes et al., 2019), myeloid cell subsets (Jordão et al., 2019), keratinocytes (Cheng et al., 2018), and fibroblasts (Croft et al., 2019) to cutaneous inflammation. Because the skin is persistently, naturally exposed to exogenous provocation, it offers a tractable means to recapitulate physiologic inflammation. Indeed, mouse skin has been used to model common human rashes, using topical and systemic administration of inflammation-inducing agents. For example, application of oxazolone to mouse skin activates cytotoxic CD8+ T cell-mediated delayed hypersensitivity, in which mast cells and neutrophils play a critical role in priming the inflammatory reaction (Dudeck et al., 2011; Vocanson et al., 2009; Weber et al., 2015). Oxazolone-mediated inflammation also appears to involve various dendritic and T helper (Th) cell types, the roles of which remain incompletely understood (Vocanson et al., 2009). In contrast, imiquimod, an agonist of Toll-like receptor 7, activates the IL-23/IL-17/IL-22-axis and interferon signaling pathways, modeling an established pathogenic mechanism for psoriasis vulgaris. However, imiquimod also likely activates additional pathways, reportedly producing molecular and histopathologic changes suggestive of systemic lupus erythematosus (Liu et al., 2018) and contact dermatitis (Garzorz-Stark et al., 2018).

We sought to understand how single-cell analysis could enhance our understanding of the molecular effects of topical provocations such as oxazolone and imiquimod on the cutaneous immune system. To what extent are inflammatory transcriptional programs shared across related CD45+ subpopulations, rather than being compartmentalized by cell identity? We were especially interested in exploring whether this patterned distribution of such transcriptional responses could be quantitatively described. An in-depth understanding of these molecular signatures might be used to help dissect the causes of pathogenic inflammation in which classical immunopathology or conventional flow cytometry may be limiting, ultimately allowing us to better classify complex skin pathologies in clinical settings such as immunological senescence (Dulken et al., 2019).

Results

Unbiased Classification of Single-Cell Transcriptomes Readily Recapitulates Known, Cutaneous Immune Cell Populations

To explore and compare inflammatory responses, we applied two distinct topical inflammatory agents, oxazolone and imiquimod, in parallel with their vehicle controls, to the ear skin of C57BL/6J mice. Each of the four treatments was applied over the course of a week, and each treatment condition was performed on three different mice (biological replicates, see Transparent Methods for details). We refer to these triplicate samples hereafter as OXA-C (oxazolone control), OXA (oxazolone), IMQ-C (imiquimod control), and IMQ (imiquimod). For each replicate, two ∼0.5 × 1.0-cm samples of treated ear skin were harvested 24 h after the last challenge. Histological analysis revealed an acanthotic epidermis with parakeratosis and spongiosis and a dermal inflammatory infiltrate with oxazolone treatment, whereas imiquimod-treated skin displayed an acanthotic epidermis with mounds of parakeratosis, underlying attenuation of the granular layer, and a dermal inflammatory infiltrate (Figure S1A). We enzymatically digested these skin samples, flow sorted for individual CD45+ cells, and then performed scRNA-seq using droplet-based microfluidics (10X Genomics). The number of total single cells initially isolated from each treatment condition was OXA-C: 20,120; OXA: 21,316; IMQ-C: 10,247; and IMQ: 7,632 (Table S1).

Using Seurat (Stuart et al., 2019), data were log-normalized and canonical correlation analysis was performed using 3,000 input variable genes to identify integration anchors among the four different treatment datasets (see Transparent Methods for details). Following integration, dimensional reduction and unsupervised Louvain modularity-based clustering was performed on 59,315 total cells to yield 21 initial cell clusters (see Methods). Visualizing these subpopulations using t-distributed stochastic nearest neighbor embedding (t-SNE, Figure S1B) confirmed their distinct identities. To define these populations, we utilized the Seurat FindAllMarkers and FindConservedMarkers functions (Stuart et al., 2019), which combine p values using meta-analysis methods to identify differentially expressed genes between clusters (conserved across treatment conditions for the latter function; Table S2). After thus identifying transcripts enriched in and typifying each of our 21 initial cell classes (p.adj <0.05), we manually assigned class identity based on comparison to well-established marker genes (Table S3). We excluded small numbers of apparent contaminant keratinocytes expressing high levels of Krt15, Col17a1, Krt2, Krt5, and Krt17; fibroblasts expressing Col1a2, Col1a1, Col3a1, Col6a1, Col6a2, and Col6a3; and endothelial cells expressing Selp, Sele, Cdh5, Cd34, and Cd93 (Hughes et al., 2019) (Figure S1B and Table S2). We also observed, as have others, that actively dividing cells cluster independently, regardless of cell identity, because of the large number of shared mitotic transcripts (Hsiao et al., 2019). Dividing cells were represented in our data by three groupings with elevated expression of Pclaf, Top2a, Mrc1, Mki67, Birc5, and numerous histone genes. Finally, we also excluded two other small CD45+ groupings. The first was enriched in the lncRNAs Gm42418 and Gm26917 associated with ribosomal RNA contamination. The second was characterized by only a few transcripts (Neb, Ppp1r16b, Rora, Odc1, Uhrf2, Fnbp1), for which a clear immune cell class could not be defined (Figure S1B and Table S2).

These exclusions produced 13 distinct immune cell subpopulations, each showing robust representation across both our oxazolone and imiquimod induction and control treatments (Figures 1A and 1B Tables S3 and S5). Six cell classes harbored markers, such as MHC II transcripts, suggestive of antigen-presenting cells (APCs). One subpopulation appeared consistent with macrophages (“Mac”), based on its expression of Adgre1, Itgam, Fcgr1, and Cd68 (Yu et al., 2016). This population also expresses the so-called tissue-remodeling macrophage markers (similar to a class also referred to as M2 macrophages or alternatively activated macrophages) such as Cd163, Mrc1, Folr2, and intermediate levels of Arg1 and Chil3 (Rőszer, 2015) (Figures 1C, 1D, and S2 and Tables S2 and S3). The second cluster appears to represent a combination of monocytes and monocyte-derived macrophages (designated as “M/MdM” in this report). This population expresses not only macrophage markers (Adgre1, Fcgr1, Itgam) but also monocyte markers (e.g., Lyz2, Ly6c2, Plac8, and Cd14) (Figures 1C and S2 and Tables S2 and S3). A smaller population of cells with elevated Cd207, Cd24a, and Epcam expression defines Langerhans cells (referred to hereafter as “LC”; Figures 1C, 1D, and 2A and Tables S2 and S3).

Figure 1

Single-Cell Profiling Demarcates Key Immune Cell Populations in Murine Skin

(A) t-SNE map shows 13 immune cell classes conserved across treatment conditions, delineated by Louvain clustering. Each dot represents one of 44,130 profiled cells.

(B) Relative proportions for each immune cell population for each paired treatment condition and its control. Top panel shows imiquimod (IMQ, blue), and imiquimod control (IMQ-C, green). Bottom panel shows oxazolone (OXA, red) and oxazolone control (OXA-C, turquoise).

(C) Immune cell population marker transcript expression levels (x axis) for the 13 immune cell populations (y axis). Size of dots represents the fraction of cells expressing a particular marker, and color intensity indicates mean normalized scaled expression levels.

(D) Violin plots show normalized transcript expression distribution on a per cluster basis for selected immune cell population marker genes that distinguish major populations. Cd163, Mrc1, and Folr2 for macrophages; Clec9a for cDC1; Sirpa in macrophages, monocytes, and cDC2; Cd207 for LC; and Cd4, Cd8a, and Foxp3 for Thet. Each dot represents gene expression of a single cell, and the kernel density plot shape represents the overall distribution.

See also Figure S1 and Tables S1, S2 and S3.

Figure 2

Cluster-Specific Single-Cell Transcript and Protein Expression Validates Immune Cell Population Assignments

(A) Cluster-specific heatmaps show normalized RNA expression count values of selected marker genes on the x axis and single cells on the y axis. Cells are ordered on y axis by treatment condition (IMQ-C, IMQ, OXA-C, OXA).

(B) Projection of protein epitope (CITE-seq) expression “gated” immune cell populations from OXA-C/OXA datasets onto transcript-based t-SNE plot from Figure 1A. Depicted are CD3−I-A/I-E+ APC cells, CD3−I-A/I-E+CD24+CD326+ Langerhans cells, CD3+TCRγ/δ+CD90.2+ DETCs, and CD3−CD117+ mast cells.

See also Figure S2 and Tables S2 and S3.

The three remaining APC populations display increased expression of H2-Ab1, H2-Eb1, and Flt3 and lack the monocyte/macrophage marker Fcgr1 and are consistent with dendritic cells. To better understand these three populations, we identified differentially expressed genes for each, in comparison with the aforementioned macrophage and Langerhans cell classes. One subpopulation shows elevated Xcr1, Irf8, and Clec9a and absence of Irf4, Itgam, and Sirpa, identifying them as likely conventional type 1 dendritic cells (referred to hereafter as “cDC1”) (Kashem et al., 2017) (Figures 1C, 1D, and S2 and Tables S2 and S3). A second population showed upregulated transcripts consistent with conventional type 2 dendritic cells (cDC2s), including Mgl2, Sirpa, Itgam, and Irf4, and absence of cDC1 markers (Irf8, Batf3, and Cd207) (Figures 1C, 1D, and S2 and Tables S2 and S3). The third population is characterized by elevated migration and activation markers Fscn1, Cacnb3, Ccr7, Cd40, Tmem123, and Cd274, demarcating a group of migratory/mature dendritic cells (referred to hereafter as “mDC”) (Bros et al., 2011; Ma and Clark, 2009; Riol-Blanco et al., 2005; Takekoshi et al., 2010; Versteven et al., 2018; Yamakita et al., 2011; Figures 1C and S2 and Tables S2 and S3).

Of the remaining seven non-APC CD45+ cell classes, three harbored elevated Cd3e/g expression suggestive of T cells. The first of these classes was characterized by high relative expression of Cd3e/g, Trdc, Tcrg-C1, Thy1, Nkg7, Fcer1g, and Il2rb, and lack of Cd4/Cd8a expression, consistent with dendritic epidermal T cells (DETCs), a population of embryonically derived, tissue-resident γδ T cell receptor-expressing cells that function in cutaneous immune surveillance (O'Brien and Born, 2015; Turchinovich and Pennington, 2011) (Figures 1C, 1D, and 2A and Tables S2 and S3). We identified a second cluster as dermal γδ T cells (dγδT) based on strong Rora expression and intermediate expression of Cd3e/g, Trdc, and Tcrg-C1/2/4, as well as expression of either Scart1 (Cd163l1; indicative of Vγ6 cells) or Scart2 (5830411N06Ri; indicative of Vγ4 cells) and lack of Cd4 and Cd8a expression (Tan et al., 2019) (Figures 1C, 1D, and 2A and Tables S2 and S3).

The remaining T cell cluster was consistent with ⍺β T cells displaying broad expression of Cd3d, Trac, Trbc1, Cd4, and Foxp3, suggesting a heterogeneous mix of T regulatory and conventional Cd4 cells. We thus denoted this category as “heterogeneous T cells,” or Thet (Figures 1C and 1D, and S2 and Tables S2 and S3). To better define this mixed Thet population, we executed a transcriptional gating strategy based on classic T cell markers (see Methods) and were able to identify Foxp3+ Tregs, conventional Cd4Foxp3- T cells (Tconv), Cd8+ T cells, and double-negative T cells (Figure S4; Table S5).

We identified another cluster as type 2 innate lymphoid cells (ILC2) based on high expression of the transcription factors Gata3, Rora, Thy1, as well as Il7r, Il5, and Il13 coupled with very low levels of Cd3d/e/g, Eomes, Rorc, Tbx21, and Cd4 (Figures 1C, 1D, and S2 and Tables S2 and S3). This group was also negative for natural killer (NK) and NK T cell markers, e.g., Klra7, Klra8, and Klrb1 (Haug et al., 2019; Juvet and Zhang, 2012; Dutton et al., 2018; Vivier et al., 2018) (Figure S2). An NK cluster with increased relative expression of NK markers, e.g., Gzma, Klra8, Klra7, Klrb1c (NK1.1), Eomes, Ncr1 was identified (Figures 1C, 1D, and S2 and Tables S2 and S3). Finally, a population of granulocytes likely composed of both mast cells and basophils (hereafter referred to as “M/B”) was defined by Gata2, Ms4a2, Kit, Mcpt4, Itgam, and Mcpt8 expression (Figures 1C and S2 and Tables S2 and S3; Li et al., 2015), whereas neutrophils, termed “Neu”, were typified by Cd14, S100a8, S100a9, Csf3r, Cebpd, Slc11a, and Spi1 (Park et al., 2018).

Epitope and Biaxial Transcript Analysis Validate and Refine Transcriptional Classification

To compare transcript-based immune cell classifications with more traditional surface epitope-based approaches (e.g., flow cytometry or CyTOF), we utilized Cellular Indexing of Transcriptomes and Epitopes by Sequencing (CITE-seq) on the OXA-C and OXA samples to quantitate single cell-level protein expression utilizing a sequencing-based output (Stoeckius et al., 2017). We bioinformatically “gated” populations based on CITE-seq expression to define protein epitope-based cell populations (Zhang et al., 2020), as would be done with flow cytometry, and overlaid these on our transcript-based clusters. Consistent with our APC transcript based-assignments (Mac, M/MdM, cDC1, cDC2, mDC, and LC), we find strong overlap with CD3−I-A/I-E+ CITE-seq gated cells (Figure 2B). CD3−I-A/I-E+CD24+CD326+ cells localize to the transcript-defined LC population. CD3+TCRγ/δ+CD90.2+ cells show ample representation within the RNA-based DETC cluster but not dγδT cells (which show lower γ/δ expression than DETCs). The transcript-based mast cell population harbors CD3−CD117+-expressing cells (Figure 2B).

Oxazolone or Imiquimod Treatment Induce APC-Dominant Shifts into Mouse Skin

Although the OXA and IMQ datasets contained the same basic cell types, their composition by type showed substantial differences. To evaluate the relative abundance of 13 immune cell types in two conditions, we examined the percentage of each subpopulation in each of the OXA-C, OXA, IMQ-C, and IMQ datasets (Figure 3A and Table S5). However, because a subpopulation can only be assessed as a proportion of total cells in a given profile, it is not possible to formally ascertain if changes in population abundance represent absolute increased numbers of a cell type or relative decreases in other major cell types. Therefore, the largest shifts are more likely to represent new, infiltrating cells during treatment. Imiquimod and oxazolone both lead to marked increase in relative percentages of macrophage (102% and 273% for OXA and IMQ, respectively), M/MdM (340% and 476%), mDC (133% and 147%), NK (2,087% and 44%), and neutrophil (1,643% and 555%) cells. Population shift decreases after oxazolone treatment include dermal γδ T cells (−87%), ILC2 (−90%), and LC (−65%), as well as decreases for DETC for both oxazolone (−92%) and imiquimod (−68%) (Figure 3A and Table S5).

Figure 3

Oxazolone and Imiquimod Divergently Reprogram Immune Cells

(A) Stacked bar plots showing relative percentages of each immune cell population for each of the treatment conditions (IMQ-C, IMQ, OXA-C, and OXA).

(B) Unsupervised hierarchical clustering heatmap for treatment-induced differentially expressed genes on a per cluster basis. Columns depict average logFC of oxazolone- or imiquimod-treated cells (versus control) for each cluster. Rows depict selected treatment-specific DEGs for OXA versus OXA-C and IMQ versus IMQ-C.

(C–F) (C) Cluster-specific scatterplots showing OXA or IMQ treatment-induced differentially expressed genes for macrophages, (D) M/MdM cells, (E) dγδT cells, and (F) M/B populations (avg_logFC for OXA DEGs relative to OXA-C on y axis and IMQ DEGs relative to IMQ-C on x axis). Pseudocolored dots represent significant DEGs (p.adj <0.05) for OXA (red), IMQ (blue), or both treatments (purple).

See also Figures S3 and S4 and Tables S4 and S5.

Oxazolone and Imiquimod Divergently Reprogram CD45+ Cell Transcriptomes toward JAK2 or Interferon Responses

Our initial classification of CD45+ single-cell transcriptomes into 13 subpopulations conflates some immune cell populations that are differentiated through classical serial gating strategies, particularly for T cell classes. However, we hypothesized that transcriptional commonalities and differences in these initial classes would enable an informative, high-level survey of gene expression differences caused by oxazolone and imiquimod. We compared transcripts differentially expressed (avg_logFC >0.5, p.adj <0.05) in either treatment condition in at least one cluster (1,088 transcripts in total, Figures 3B and S3A). This approach immediately accentuated striking treatment- and cluster-specific transcriptional shifts. Globally, IL17/22 and IFN-γ signatures emerge in multiple subpopulations after imiquimod treatment. In contrast, oxazolone generates a more complex picture, with upregulation of type 1 and type 2 cytokines, and IFN-responsive genes in selected cell classes, consistent with previous reports (Liu et al., 2019; Pantelyushin et al., 2012; Thomson et al., 1993; Wenzel et al., 2005) (Figures 3B, S3A, and S3B and Table S4).

The widespread IMQ-specific induction of IL17a, Il17f, and Il22 transcripts is most robust in dγδT cells (Figures 3B, 3E, and S3A and Table S4). We also observe pervasive upregulation of interferon-responsive genes, including a gene cluster composed of Ifitm3, Isg15, Zbp1, Rtp4, and Isg20 as well as transcripts not classically known to be interferon responsive (e.g., Ly6c2, Ms4a4c, and Pld4). This interferon signature is most prominent in imiquimod-treated APC populations (such as Mac, M/MdM, and cDC2) and seen only faintly in the corresponding oxazolone-treated APC populations (Figures 3B, S3A, and S3B and Table S4). Imiquimod also upregulates chemokines associated with type 1 inflammatory responses (Ccl4 and Ccl5) most prominently in M/B cells and APC populations (Figures 3B, 3D, 3F, S3A, and S3B). Oxazolone, in contrast, primarily elevates only Ccl4, in M/B cells and neutrophils.

Oxazolone widely elevates Jak2 and Stat3 expression, as well as Il4ra, effects absent after imiquimod treatment (Figures 3B and S3A and Table S4). Oxazolone-treated cells also show marked upregulation of type 2 cytokine transcripts (Il4, Il6, Il13) and CCR1 chemokine ligands (Ccl6 and Ccl9), most prominently in M/B cells (Figures 3B and 3F and Table S4). We also observe robust upregulation in OXA-treated APC populations of genes such as the Ccr1 and Ccr5 chemokine receptors, the Ccl8 and Ccl24 beta-chemokines, as well as matrix metalloproteinases (Mmp12/14/19) and Fabp5 (Figures 3B–3D and S3 and Table S4). Except for Ccr5, these shifts are nearly absent from the IMQ dataset. The neutrophil chemoattractant Cxcl3 is also elevated in oxazolone-treated neutrophils (Figures 3B and S3 and Table S4).

To better characterize the effects of imiquimod and oxazolone on the Thet cell population, we subsetted these cells based on a transcriptional gating approach into the four T cell types described above: Foxp3+ Tregs, conventional Cd4+Foxp3- T helper cells (Tconv), Cd8a/b1+ T cells, and Cd4-Cd8a/b1- double-negative T cells (DNT) (Figure S4 and Tables S4 and S5; Methods) (Gao et al., 2011). We then assessed imiquimod or oxazolone-induced RNA expression changes for the 265 most variable genes, as well as selected type 1, type 2, type 17 cytokine, and interferon response-related genes (Figure S4). We find key treatment-induced differentially expressed genes previously identified in the Thet population in all four Thet subpopulations, i.e., IMQ-specific upregulation of Il17a, Il17f, Il22, and OXA-specific upregulation of Jak2 and Stat3.

Oxazolone induced a marked increase in the number of Tregs (Table S5), which was accompanied by upregulation of Ets1 and Ikzf2, transcription factors known to both regulate Foxp3 expression and mediate the suppressive activity of this cell type (Kim et al., 2015; Mouly et al., 2010) (Figure S4, Table S4). In oxazolone-treated cells, there is a trend toward T cell activation gene upregulation in Cd8+ cells, e.g., the co-inhibitory receptors Pdcd1 (Pd-1), Lag3 and Havcr2 (Tim-3), granzyme molecules Gzma and Gzmb, and Ccr5 (Okoye et al., 2017; Viola and Luster, 2008) (Figure S4 and Table S4). We also observe OXA-induced upregulation of Havcr2 in Cd4+T cells, comporting with its known expression in CD4+ Th1 cells, where it negatively regulates IFNG production (Das et al., 2017; Hastings et al., 2009) (Figure S4 and Table S4).

Imiquimod treatment induced the upregulation of MHC I and interferon response genes across all T cell types, although only Ifi27l2a was significant across all T cell subpopulations (Figure S4 and Table S4). Interestingly, the cathepsin Ctss, which is known to be expressed in barrier skin cells and is mechanistically implicated in both atopic dermatitis and psoriasis, was strongly upregulated in all imiquimod-treated T cell types (Figure S4 and Table S4; Ainscough et al., 2017; Kim et al., 2012; Schönefuss et al., 2010). Consistent with our transcriptional data, flow cytometry analysis revealed a trend toward increased IL17A and IL22 expression in Tconv, CD8, and dγδT cells, as well as increased CCL5 expression in CD8 and Treg cells with imiquimod treatment (Figure S5).

Basophils Infiltrate Skin as Primary Producers of IL4 and IL13 following Oxazolone but Not Imiquimod Treatment

We observed Il4 and Il13 induction mainly restricted to the M/B cluster, which expresses known granulocyte markers (Figure 3 and Table S4). To more deeply understand this phenomenon, we subjected these cells to κ-means clustering, which revealed two distinct classes (Figure 4A). Although M/B granulocyte cells globally express Gata2 (Figure 4B), the putative mast cell cluster was clearly demarcated by Mcpt4 and Kit (the latter on both the transcript and epitope levels), whereas the presumptive basophils uniquely expressed Mcpt8 and Itgam (Akahoshi et al., 2011; Siracusa et al., 2013) (Figure 4B). Surprisingly, basophils were found almost exclusively in the OXA dataset, and not in OXA-C, IMQ-C, or IMQ (Figure 4B).

Figure 4

Infiltrating Basophils Produce Il4 and Il13 after Oxazolone but Not Imiquimod Treatment of Murine Skin

(A) Heatmap displaying count values for the top 100 most variable genes (y axis) defined by pseudotime (x axis) for M/B population cells. Two nonredundant and non-covariate k-means-defined subpopulations are denoted by the color bar at the top of the panel (mast cells, red, left; basophils, blue, right). Rows represent genes and columns represent single cells of the M/B cluster.

(B) Feature plots showing mast and basophil marker expression projected onto the t-SNE M/B cluster from Figure 1A. Each row depicts transcript expression of a different marker gene for mast cells (c-Kit RNA and protein [CITE-seq] and Mcpt4) and basophils (Itgam, Mcpt8). Each column shows M/B population t-SNE plots representing only the cells from the treatment condition listed at the top.

(C) Principal-component analysis (PCA) of mast cells (orange/red, left) and basophils (blue/green, right) demonstrating these subpopulations do not share a differentiation trajectory.

(D) Differential gene expression volcano plot between mast cells and basophils. x axis depicts the average logFC for basophils relative to mast cells. y axis depicts –log10(p.adj). Significant genes are defined as |avg_logFC |≥1 and p.adj ≤0.05 [–log10(p.adj) ≥1.3]. Blue dots represent basophil, and red dots represent mast cell upregulated genes.

See also Table S4.

The restriction of basophils to this dataset suggests that they were uniquely induced to infiltrate skin by oxazolone (Figure 4B). Given the transcriptional relatedness of granulocytes, we first evaluated the alternative hypothesis that our putative basophils might instead represent reactive mast cells. Utilizing global lineage structure and inference of pseudotime variables (Street et al., 2018) (Slingshot), we were unable to discern an evolutionary trajectory connecting these subpopulations (Figure 4C). We next analyzed gene expression differences between mast cells (found in both OXA and IMQ) and basophils (found only in OXA). Basophils infiltrating after oxazolone treatment appear not only specifically enriched in key marker transcripts such as Mcpt8 but also appear to serve as the exclusive source of new Il4, Il6, and Il13 (Figure 4D).

Oxazolone Provokes a More Compartmentalized Pattern of Global Transcriptional Response, Relative to CD45+ Cell Identity, Than Imiquimod Stimulation

As noted above, oxazolone activates sharply distinct transcriptional programs in different immune cell types. Imiquimod, in contrast, appears to activate related patterns across multiple populations. One of the opportunities conferred by single-cell data is the ability to quantitate such differing global molecular responses, allowing their comparison and classification by algorithms. Such approaches could eventually enable us to deduce the mechanism of new sources of inflammation, based solely on analyzing single-cell patterns of immune response and placing them in context of known agents.

As an initial test of this hypothesis, we first systematically identified global, single-cell, transcriptional response patterns produced by oxazolone and imiquimod. We next examined whether these patterns showed simple compartmentalization between our identified immune cell types, or whether they showed more complex distribution throughout our CD45+ subpopulations. In other words, we sought to ask whether the molecular responses to a given agent are wholly defined by immune cell type, or are they more heterogeneously distributed across cell types? We aimed to develop unbiased, quantitative means to make such an assessment.

To technically execute this approach, for each cell in the oxazolone and imiquimod-treated samples, we initially identified one cell in its respective control group with the most similar, harmonized transcriptomic pattern, hypothesizing that this control cell could be used as a proxy for a pre-treatment version of the same cell type (Figure 5A). We only considered control cells that were in the top 50 nearest neighbors. We then calculated gene expression differences over the 1,000 most variable genes between thousands of such pairs of cells, generating a vector of pairwise single-cell gene expression comparisons between a treatment and its control:DiVNCE(Tx) =

Figure 5

DiVNCE Analysis Yields High Degree of Compartmentalization by Cell Identity for Oxazolone but not Imiquimod Transcriptional Response in CD45+ Cells

(A) Schematic representation by which DiVNCE identifies the most similar, “neighboring” control cell for each treatment cell and then catalogs their transcriptional differences.

(B) UMAP representation showing partitioning of oxazolone-treated cells based on 12 DiVNCE profiles (six highlighted for comparison, upper left panel), with same clusters colored by cell identities as depicted in Figure 1A (upper right). Lower left panel shows imiquimod-treated cells partitioned on 11 DiVNCE profiles (six highlighted for comparison) and same clusters colored by cell identity (lower right). Mac clusters 1–5 represent five DiVNCE patterns finely subdividing macrophages in the Mac classification, highlighting greater compartmentalization. In contrast, in the IMQ dataset, interferon-responsive DiVNCE signatures (IRS 1–5) are each shown by multiple APC types, showing mixing and lower relative compartmentalization. Also displayed for validation are unbiased detection of a basophil Mcpt8/Il13 DiVINCE signature in the OXA M/B granulocyte cell identity cluster and an IL17 signature in a subset of IMQ dermal γδ T cells (Figure S6 and Table S6).

(C) Log2-fold elevations from DiVNCE comparisons (represented by circle size) for inflammatory transcripts (y axis) for MacS S and IRS signatures (x axis), all p.adj <1 × 10−3.

See also Figure S5 and Table S6.

where Tx represents a treated cell, C(Tx) represents the nearest control cell, and sctx() is a vector of scaled expression values for 1,000 genes in a given cell. We refer to this library of pairwise comparisons as “Difference Vs. Nearest Control cEll,” or “DiVNCE” profiles. The DiVNCE profiles were clustered, producing 11 discrete response signatures for the OXA dataset and 13 for IMQ (Table S6). We next compared DiVNCE patterns for OXA and IMQ against our initial 13 CD45+ cell classifications. This unsupervised approach re-identifies a basophil Mcpt8/Il13 DiVNCE signature in the OXA M/B granulocyte cluster and an Il17 signature in a subset of IMQ dermal γδ T cells, validating our approach (Figures 5B and S6 and Table S6).

In Figure 5B, the upper left UMAP shows OXA CD45+ cells clustered based on 12 DiVNCE signatures; the upper right panel shows the same cell organization colored by our original 13 CD45+ cell identities. Similarly, the lower left shows the IMQ dataset clustered based on 13 DiVNCE patterns, then colored on the lower right panel by cell identity. It is readily apparent from this representation that DiVNCE signatures more precisely delineate cell identities for OXA than for IMQ. For example, the OXA Mac subpopulation is subdivided into five single-cell DiVNCE subpopulations (MacS1–5), each with distinct inflammatory patterns (Figure 5C). Even in the Mac DiVNCE cluster that includes other APCs groups such as mDC and M/MdM (MacS 5), these other APCs are spatially distinct (Figure 5B), suggesting that they express the MacS 5 pattern differently. Overall, the DiVNCE clusters recapitulate cell identities in the OXA data, in cases at higher resolution. In contrast, APCs in the IMQ data are mixed in five DiVNCE clusters based on interferon response signatures (IRS 1–5) (Figures 5B and 5C). We generated a quantitative measure for positive compartmentalization by calculating the adjusted rand index (ARI) between the DiVNCE cluster assignment and the CD45+ cell type assignments. The ARI for OXA is 0.272, 95% confidence interval (CI) = (0.268–0.275), compared with 0.128, 95% CI = (0.124–0.131) for IMQ, reflecting significantly enhanced compartmentalization of global expression changes in the oxazolone response.

Discussion

We present here one of the first single-cell comparisons of inflammatory responses to distinct provocations in mouse skin. While confirming well-established molecular features of the response to oxazolone and imiquimod, our global CD45+ cell capture also allows us to directly examine the molecular behavior of multiple cell populations, some of which have not been interrogated directly in standard reports. In each of our datasets, we readily identify major APC and T cell classes, many of which we further refine with successive rounds of multiparametric segregation. To validate our transcript-based classifications, we also utilized sequencing-based protein epitope analysis (CITE-seq), which showed strong correlation for transcript- and protein epitope-defined APC, Langerhans cell, DETC, and mast cell populations. This approach validates most functionally important immune cell subpopulations, including DETC and cDC2, as well as Langerhans and mast cells.

These single cell profiles offer unprecedented resolution in dissecting divergent cutaneous inflammatory immune response. Continuous oxazolone application is classically associated with immediate type hypersensitivity followed by a shift toward production of Th-2-associated cytokines (Man et al., 2008; Wang et al., 2000). Oxazolone has been used to model both allergic contact dermatitis and atopic dermatitis, although transcriptomic studies suggest that IL23 injection more closely resembles the latter (Ewald et al., 2017). Our data was striking for elevation of Jak2 and Stat3 expression in most CD45+ cell populations, suggesting that a global activation of this pathway is central to some instances of allergic contact dermatitis. Successful JAK inhibition of the oxazolone response has been observed previously, primarily using JAK3-selective agents (Fujii and Sengoku, 2013; Mahajan et al., 2015), leaving open the possibility that JAK2-selective agents may be even more effective in counteracting a subset of delayed-type hypersensitivity.

While granulocytes have long been known to be generally capable of Il4 and Il13 production (Gessner et al., 2005), our analyses demonstrate the power of scRNA-seq to precisely differentiate cell identity and demarcate transcript expression as a function of that cell identity. Our single-cell analyses detect newly arriving basophils after oxazolone treatment, which then appear to serve as the sole detectable source of elevated Il4 and Il13 production. The restriction of this induction to infiltrating basophils is nearly complete, echoing type 2 results in human patients with atopic dermatitis (Mashiko et al., 2017). Thus, agents inhibiting basophil migration or activity may be especially useful in ameliorating some type 2 responses. Although mast cells are known to produce these type 2 cytokines under some conditions, they do not demonstrate this role in our study. Imiquimod provokes no comparable Il4 or Il13 production. Future investigation will determine whether other inducers of delayed-type hypersensitivity produce a similar profile, or whether the landscape of type 2 cytokines production is instead antigen-specific.

Topical provocation with imiquimod is known to produce Il17, mainly from IL23-induced Th17 cells (Zheng et al., 2007). More recently, skin γδT cells, mainly Vγ4+ γδT cells, as well as Rorγt+ ILCs cells, have been established to produce Il17a, Il17f, and Il22 (Pantelyushin et al., 2012). Consistent with these reports, we detect dγδT cells as the main producer of these cytokines (Figure 3A), but there is also significant Il17/Il22 production by Vγ6+ γδT cells (Figure 2A). Along with type 17 cytokine upregulation in macrophage, M/MdM, and Thet cell subpopulations, our data suggest that imiquimod impacts target tissue through a broader set of immune cells than previously appreciated.

We also observed the pan-upregulation of Ccl4 and Ccl5 in most imiquimod-treated CD45+ subpopulations, a less-described phenomenon. A recent study indicates that these cytokines may also be involved intimately in the imiquimod response, quenching inflammation through Treg recruitment (Oka et al., 2017). However, Ccl5 and Ccr5 (encoding the receptor for both CCL4 and CCL5 ligands) are upregulated in human psoriasis and downregulated after treatment, supporting a role in pathogenicity. The discrepancy may derive from the differences between humans and mice (Mack et al., 2001), as well as complex ligand-receptor interactions between CCL4/CCL5 and CCR1/CCR3. Further work will be required to deconvolute these roles.

We describe and demonstrate a novel method, DiVNCE, that quantitatively describes discrete expression programs generated by stimuli and allows us to assess their distribution among immune cell populations. We show that whereas oxazolone produces a highly cell type-specific pattern, imiquimod enforces a more heterogeneous response in T cells, consisting of multiple, distinct, interferon-related programs (Figure 5B). Interestingly, APCs show greater plasticity in response to the agents in this study than T cell subpopulations. We describe and apply the ARI metric, which quantitates this greater degree of molecular response compartmentalization for oxazolone at statistical significance. Even more nuanced pattern-based matching approaches may allow the immunological consequences of new chemical agents to be rapidly assessed by comparison against a library of existing compounds. Further investigation, based on profiling additional types of immune stimulation, will be required to validate and extend our findings.

Limitations of the Study

Transcriptional similarity-based single-cell computational classification may lead to conflation of similar immune cell types. In particular, this was noted in the heterozygous T cell and mast/basophil populations, which required further computational subdivision to identify cell populations that better correlate with classic flow cytometry-defined populations. Furthermore, the inherent biases of droplet-based scRNA-seq and utilization of RNA expression for immune cell type classification will lead to variation from protein epitope-based flow cytometry-defined populations.

Large treatment-induced increases in APC populations and their known anatomic site-specific diminishment in mouse cutaneous ear led to a relative paucity of certain T cell subtypes (e.g., CD4, CD8, Treg). Consequently, statistical power for analyses with these populations was limited.

Resource Availability

Lead Contact

Further information and requests for resources should be directed to and will be fulfilled by the Lead Contact, Jeffrey Cheng (jeffrey.cheng@ucsf.edu).

Materials Availability

This study did not generate new unique reagents.

Data and Code Availability

Single-cell RNA-seq data have been deposited in the NIH Gene Expression Omnibus (GEO) database, under accession number GSE149121. All statistical analysis and plotting of scRNA-seq and cell surface protein data were performed using Rstudio software (Version 1.2.5033). All analysis scripts are available upon request to corresponding authors.

Methods

All methods can be found in the accompanying Transparent Methods supplemental file.