Finding a TRAIL to allergy control
House dust mites (HDMs) are a major source of aeroallergens that trigger human allergic responses including asthma. To characterize the heterogeneity of human CD4+ T helper (TH) cell responses to HDM antigens, Seumois et al. stimulated blood T cells with HDM peptides and isolated the freshly activated T cells for single-cell RNA sequencing. A TH subset expressing an interferon response gene signature (designated THIFNR) and the cell death–inducing TRAIL ligand was more frequent in individuals without HDM allergy. This TRAIL+ THIFNR subset may function to restrain allergic TH2 responses. Individuals with both HDM allergy and asthma were enriched for interleukin-9–producing TH2 cells. These findings provide a basis for designing more precise approaches to thwart pathogenic TH responses contributing to allergic asthma.
Abstract
CD4+ T helper (TH) cells and regulatory T (Treg) cells that respond to common allergens play an important role in driving and dampening airway inflammation in patients with asthma. Until recently, direct, unbiased molecular analysis of allergen-reactive TH and Treg cells has not been possible. To better understand the diversity of these T cell subsets in allergy and asthma, we analyzed the single-cell transcriptome of ~50,000 house dust mite (HDM) allergen–reactive TH cells and Treg cells from asthmatics with HDM allergy and from three control groups: asthmatics without HDM allergy and nonasthmatics with and without HDM allergy. Our analyses show that HDM allergen–reactive TH and Treg cells are highly heterogeneous and certain subsets are quantitatively and qualitatively different in individuals with HDM-reactive asthma. The number of interleukin-9 (IL-9)–expressing HDM-reactive TH cells is greater in asthmatics with HDM allergy compared with nonasthmatics with HDM allergy, and this IL-9–expressing TH subset displays enhanced pathogenic properties. More HDM-reactive TH and Treg cells expressing the interferon response signature (THIFNR and TregIFNR) are present in asthmatics without HDM allergy compared with those with HDM allergy. In cells from these subsets (THIFNR and TregIFNR), expression of TNFSF10 was enriched; its product, tumor necrosis factor–related apoptosis-inducing ligand, dampens activation of TH cells. These findings suggest that the THIFNR and TregIFNR subsets may dampen allergic responses, which may help explain why only some people develop TH2 responses to nearly ubiquitous allergens.
INTRODUCTION
Asthma is characterized by aberrant type 2 immune responses to common inhaled aeroallergens such as house dust mite (HDM), grass pollen, animal dander, and mold (1–6), leading to “asthma attacks” in sensitized asthmatic individuals in response to inhalation of such allergens (7). The hallmarks of asthma, namely, airway narrowing and sputum eosinophilia, have been shown to result from the specific activation of [major histocompatibility complex (MHC)] class II–restricted CD4+ T helper (TH) cells by challenging asthmatics with synthetic allergen-derived peptides (8–12). Further evidence of the centrality of TH cells in asthma pathology is that their depletion reduces allergic airway inflammation in animal models (13) and that inhibition of TH cell–derived type 2 cytokines [interleukin-5 (IL-5), IL-13, and IL-4] is clinically beneficial in patients with asthma (14–17). However, despite the central role of allergen-reactive TH cells and their products in driving airway inflammation, the full spectrum and function of TH cell subsets that respond to common allergens has yet to be defined. Similarly, although an imbalance between regulatory T (Treg) cells and TH cell responses to allergens is associated with the development of allergy and asthma (18–23), the heterogeneity of allergen-reactive Treg cells remains unstudied.
Previous studies of allergen-reactive T cells have characterized their phenotype on the basis of the expression of cell surface markers or canonical cytokines (24–26). Because of their relative rarity, analyses of these cells usually require in vitro expansion, which can alter their molecular properties, thus limiting the value of unbiased transcriptomic studies (27–29). Furthermore, transcriptomic studies performed at the whole population level fail to capture cellular heterogeneity and also lack the resolution to detect biological differences associated with asthma or allergy (30). A recent single-cell analysis of TH cells in mouse models of allergic airway inflammation revealed substantial heterogeneity and also identified TH subsets that had not been previously described (20).
Characterizing the various subsets of TH and Treg cells in asthmatic individuals and comparing their frequency and properties with those in individuals without asthma is ideally achieved at single-cell resolution. Single-cell transcriptomic analysis can help define the molecular properties of allergen-reactive TH cells associated with pathology and assess whether these features are the result of an expansion of a preexisting population of cells or the result of their aberrant differentiation in response to environmental signals (31, 32). To address the latter issue, the subsets of allergen-reactive TH cells must also be defined in individuals without asthma and allergy. Such allergen-reactive TH cells are present even in nonallergic individuals (33–36), although it is not known why or how these cells fail to cause overt allergic responses.
To address these questions in a hypothesis-free manner, we performed single-cell transcriptomic analysis of TH and Treg cells that react to HDM allergen. HDM is one of the most common and ubiquitous allergens, and sensitization is associated with both the onset of allergic asthma and its severity (37–40). The relatively high abundance of HDM-reactive T cells in the blood makes it possible to isolate sufficient number of cells for high-throughput single-cell transcriptomic analysis. Here, we report on the single-cell transcriptomes of >50,000 HDM-reactive T cells from allergic asthmatic individuals and relevant control groups. Our analysis revealed multiple distinct subsets of TH and Treg cells that are either preferentially expanded or depleted in asthmatic individuals with and without HDM allergy, defined the pathogenic properties of TH subsets associated with allergic asthma, and uncovered a unique HDM-reactive TH subset that is expanded specifically in individuals without HDM allergy.
RESULTS
Bulk RNA sequencing analysis of HDM allergen–reactive T cells does not identify asthma-specific features
To comprehensively characterize the molecular properties of allergen-reactive TH and Treg cells from patients with asthma, we isolated pure populations of HDM-reactive memory TH and Treg cells ex vivo (see Materials and Methods and fig. S1) from asthmatics with HDM allergy (n = 6) and performed both bulk and single-cell RNA sequencing (RNA-seq) (Fig. 1A). To distinguish the molecular features that are specific to asthma as opposed to HDM allergy, we performed similar assays in HDM-reactive TH and Treg cells isolated from HDM-allergic individuals without asthma (n = 6). Because allergen-reactive T cells are present even in nonallergic individuals (33–36), we also isolated HDM-reactive TH and Treg cells from asthmatic (n = 6) and healthy participants (n = 6) without HDM allergy to uncover features that may contribute to the lack of HDM allergy, i.e., immunoglobulin E (IgE) reactivity. In total, we performed 95 bulk RNA-seq and ~50,000 single-cell RNA-seq assays on T cells from a total of 24 individuals (Fig. 1A and table S1).
(A) Schematic representation of the study summarizing the participant groups, activation assay, and sorting and sequencing strategies. (B) Dot plot showing the percentage of CD4 memory HDM-reactive T cells that were CD154+ (left) or CD154− CD137+ (right) for each participant group. Horizontal line mean; error bar, SE. (C) tSNE plot of bulk RNA-seq samples: 6 HDM− T cells (from asthma-allergic patients), 24 HDM+ TH (from all patients), and 18 HDM+ Treg (from all patients). (D) Heat map of row-wise z-score–normalized expression for 724 genes differentially expressed between the three groups of cells in (C). Each column represents data from one individual. DEG, differentially expressed genes. (E) GSEA of grouped bulk RNA-seq datasets presented in (C). q, false discovery rate; NES, normalized enrichment score; RES, relative enrichment score; a list of genes is provided in table S4. (F) tSNE plot of bulk RNA-seq datasets for HDM− T samples (n = 6) and TH samples colored by disease group, where each dot represents one RNA-seq data from one individual (n = 12 per group). (G) Scatter plot shows log2 fold change in expression of significantly differentially expressed genes between asthma (n = 24) versus nonasthma (n = 24) (x axis) and allergic (n = 24) versus nonallergic (n = 24) (y axis) TH samples. (H) Scatter plots show coexpression of the indicated canonical TH1 and TH2 cytokines in TH samples coded by disease group.
HDM-reactive TH cells (0.2 to 3% of all memory TH cells) and Treg cells (1 to 5% of all memory Treg cells) were detected in all four participant groups, including the HDM-allergic and nonallergic individuals (Fig. 1B). Bulk transcriptome analysis showed that HDM-reactive TH and Treg cells clustered separately from one another and from HDM-nonreactive cells (HDM− T cells) (Fig. 1C). A total of 724 transcripts were differentially expressed between both HDM-reactive, activated T cell populations, TH, and Treg (after stimulation with HDM peptide/MHC complex from antigen-presenting cells) and HDM− T cells (not stimulated by HDM allergen–derived peptides) (adjusted P < 0.01, log2 fold change >2; Fig. 1D and table S2). As expected, these differentially expressed transcripts were highly enriched for genes in the T cell receptor (TCR) signaling pathway (Fig. 1E, top). Allergen-activated HDM-reactive TH cells expressed greater amounts of several transcripts encoding cytokines (IL-2, IL-13, IL-5, IL-4, IL-9, IL-31, IL-17F, IL-22, TNF, IFN-g, and GM-CSF) and chemokines (CCL20 and CXCL10) linked to effector functions (Fig. 1E, middle). HDM-reactive Treg cells expressed higher levels of genes linked to Treg function, such as IL2RA, FOXP3, CTLA4, IKZF2, and TNFRSF8, when compared with HDM− T cells [Fig. 1E (bottom), fig. S2, and table S2].
Clustering analysis of HDM-reactive TH cells by disease group showed separation based on HDM allergy status rather than asthma phenotype (Fig. 1F). For example, in HDM-reactive TH cells from HDM-allergic individuals, expression of canonical TH2 cytokines was increased compared with those from HDM-nonallergic individuals (Fig. 1G), whereas no significant differences were observed between the HDM-reactive TH cells from asthmatic versus nonasthmatic individuals with HDM allergy (Fig. 1G). The heterogeneity observed within the HDM-reactive TH population, reflected in the coexpression of transcripts encoding canonical TH1, TH2, and TH17 cytokines (Fig. 1H), is likely to have limited the resolution of bulk transcriptome data to distinguish asthma-specific features.
Single-cell RNA-seq analysis reveals heterogeneity among HDM-reactive TH cells
Single cells from all six individuals in each disease group were pooled for droplet-based single-cell RNA-seq (10x Genomics platform), and genotype-based deconvolution was used to obtain individual-specific single-cell transcriptomes and to exclude potential cell doublets (see Materials and Methods and fig. S3). Our cell isolation strategy, based on the CD154 activation marker, primarily enriches for HDM-reactive TH cells (41–43). Analogous to flow cytometry–based approaches, single-cell transcriptome analysis allows discrimination of activated (true positives) from nonactivated (false positives) TH cells. On the basis of a TH activation signature, derived by comparing HDM-reactive TH and HDM-nonreactive (HDM− T cells) single cells (Fig. 2A and fig. S3), we eliminated potential false-positive cells from the HDM-reactive TH cell population (Fig. 2A and fig. S3).
(A) Top: Heat map of row-wise z-score–normalized expression for 110 genes used to establish the activation score (Materials and Methods). Rows are ordered by hierarchical clustering. Each column represents a single-cell RNA-seq data, ordered by activation score. Right: Examples of genes included in the TH activation score. Bottom: The histogram shows the density function of HDM− T (n = 3075; gray) and HDM+ TH (n = 31,105; red) cell activation scores (Materials and Methods). The red line indicates the threshold of selection (activation score of −0.27, 5% of HDM−). (B) tSNE visualization of Seurat clustering analysis of approximately 28,313 single HDM+ TH cell transcriptomes obtained from all 24 individuals. IFNR, interferon response; ACT, biologically uncharacterized activated T cells (three groups). Top right: The pie chart shows the cell number proportion of each cluster. (C) Heat map of row-wise z-score–normalized mean expression of cluster-specific differentially expressed genes. Columns represent the average expression for each cluster, ordered on the basis of biological relevance. Right: Lists of biologically relevant example genes for each cluster. Equal numbers of cells were sampled from each cluster. (D) Violin plots show log2(CPM + 1) normalized expression in each cluster (three THACT clusters merged) for 24 cluster-specific signature genes (six per cell type). The color scale represents the fraction of cells within each cluster expressing the given gene, excluding cells with no expression. (E) GSEA between each cluster and other clusters of single-cell transcriptome datasets presented in (B). The list of genes is provided in table S4. a.u., arbitrary units.
Analysis of the single-cell transcriptomes of HDM-reactive TH cells (nondoublet and activation signature positive) revealed seven clusters (Materials and Methods and fig. S4) present at varying frequency among individuals, highlighting the importance of studying cells from multiple individuals (Fig. 2B and fig. S5). To understand the molecular properties unique to each cluster, we performed multiple pairwise single-cell differential gene expression analyses (Materials and Methods and table S3). Several hundred genes (n = 687) were especially highly expressed by each cluster, allowing classification into specific TH subsets (Fig. 2C). Cells in cluster 1 were highly enriched for transcripts encoding canonical type 2 cytokine genes (IL5, IL13, and IL4), the TH2 master transcription factor GATA3, and receptors (IL1RL1 and IL17RB) for the TH2-polarizing cytokines IL-33 and IL-25, indicating that this cluster represented TH2 cells (Fig. 2D). The TH2 subset only represented ~6.3% of the HDM-reactive TH cell population (Fig. 2B). Cluster 2 was enriched for TH1 phenotype- and function-related genes such as IFNG, CXCR3, and PRF1 (Fig. 2, C and D) (44–47). In addition, we found that the expression of genes encoding the chemokines X-C motif chemokine ligand 1 (XCL1) and XCL2 was correlating with expression of IFNG (fig. S6), suggesting a potentially important role of these chemokines in the function of TH1 cells (46–49). Cluster 3 was enriched for TH17 phenotype- and function-related genes such as IL17A, IL17F, CCR6, IL22, CTSH, and CCL20 (50–52). The characteristics of cell clusters 1 to 3 were confirmed by gene set enrichment analysis (GSEA) using curated lists of signature genes (Fig. 2E and table S4). Cells in cluster 4 were very highly enriched for type I and type II interferon response genes (IFI6, MX1, ISG20, OAS1, IFIT1, and IFI44L) (53, 54), indicating that they represent a previously uncharacterized TH subset, which we termed TH subset expressing the interferon response signature (THIFNR) (Fig. 2, C and D). GSEA analysis confirmed enrichment of interferon response genes in this cluster (Fig. 2E and table S4). Cells in clusters 5 to 7 were enriched for genes linked with cell activation; cluster 5 (THACT1) was enriched in genes linked to ribosomal proteins and RNA translation (RPLx and RPSx; see table S3), cluster 6 (THACT2) was enriched with genes linked to endocytosis and membrane trafficking (ARL6IP5, ARPC5, and BIN1; see table S3), and cluster 7 (THACT3) was enriched in genes linked to chromosome maintenance (NPM1 and NHP2) and apoptosis (GADD45B, NFKB1, ATF4, and PMAIP1; see table S3). Overall, our single-cell transcriptome analysis uncovered substantial heterogeneity among HDM-reactive TH cells.
Proportions of HDM-reactive TH subsets differ between HDM-allergic and nonallergic individuals
We next asked if the proportions of any of the HDM-reactive subsets varied between individuals with or without HDM allergy or asthma. As expected, the TH2 cluster (cluster 1) was present only in the HDM-allergic groups (Fig. 3, A and B, and fig. S7), consistent with the central role of TH2 cells and type 2 cytokines in IgE class switching and allergy and asthma pathogenesis (55–57). On the other hand, the TH1 cluster, although observed in all participant groups, was present at greater proportions in participants without HDM allergy, consistent with the reciprocal role of TH1 cells in dampening TH2 differentiation. Several other clusters, including the TH17 cluster, showed no significant differences in their proportions among disease and control groups (Fig. 3, A and B, and fig. S7).
(A) tSNE visualization of Seurat clustering analysis shown in Fig. 2B, using equal cell numbers for each disease group (n = 3720), obtained from all six individuals in each group. Cells are colored according to cluster as in Fig. 2B. (B) Pie chart illustrating the relative proportion of cells from each disease group in the four biologically relevant clusters. (C) Scatter plot shows the log2 fold change of expression of THIFNR signature genes between asthmatic (n = 12) versus nonasthmatic (n = 12) (x axis) and allergic (n = 12) versus nonallergic (n = 12) (y axis) individuals among cells in the THIFNR cluster. Equal numbers of cells were sampled per group. Dotted lines indicate the threshold value of fold change for gene filtering. (D) Violin plots show log2(CPM + 1) normalized expression of CXCL10 and TNFSF10 in each TH cluster. Cells with no expression were excluded. (E) Scatter plots show coexpression of TNFSF10 with the products of the THINFR signature genes IFI6 and ISG15 by THINFR cells (left) or HDM− T cells (right). Each dot represents one cell. (F) Contour plots show the expression of CD69 versus TRAIL in memory CD4+ T cells before (left) and after 6 hours (center) of TCR stimulation with anti-CD3 and anti-CD28. Both plots are representative of five independent experiments. Numbers indicate the percentage of cells in each quadrant. Right: Quantification of each of the six experiments; bars represent mean and SE. (G) Diagram of the TRAIL functional assay. rhTRAIL, human recombinant TRAIL. (H) Left: Contour plots show the expression of the cell activation markers CD154, CD69, and CD137 in memory CD4+ cells after 6 hours of stimulation in the presence or absence of TRAIL. Data shown are from a representative experiment. Right: Quantification of each of the six experiments; bars represent mean and SE. **P < 0.01 and ***P < 0.001.
Individuals without HDM allergy (both the asthmatic and healthy control groups), despite displaying a substantially broad TH response to HDM allergen, failed to generate TH2 cells that respond to HDM ex vivo. The large majority of HDM-reactive TH cells expressing the interferon response signature (THIFNR, cluster 4) were observed in individuals without HDM allergy (Fig. 3, A to C). This negative association raised the possibility that the THIFNR subset plays a role in dampening TH2 responses to allergens. Cells in the THIFNR subset expressed the highest levels of CXCL10 and TNFSF10 (Fig. 3D). CXCL10 encodes CXCL10, a chemokine that recruits TH cells expressing the chemokine receptor CXCR3, which mainly comprises TH1 cells (58, 59). Thus, CXCL10 expression by THIFNR cells is likely to promote selective recruitment of TH1 cells (60). TNFSF10 encodes tumor necrosis factor–related apoptosis-inducing ligand (TRAIL), which can drive apoptosis in cells expressing its receptors (TRAIL-R) (61, 62). More recently, both surface-bound and soluble TRAIL have been shown to dampen TCR signaling by inhibiting the phosphorylation of downstream kinases (63–66). Given that activated TH cells express TRAIL-R (65, 66), TRAIL produced by the HDM-reactive THIFNR cells may play an important role in blocking TH cell responses to HDM in vivo. A small fraction of cells expressing THIFNR signature genes (IFI6 and ISG15) and TNFSF10 were observed even in resting cells not reactive to HDM, suggesting persistence of this population within peripheral blood mononuclear cells (PBMCs) (Fig. 3E). We confirmed that after TCR stimulation TRAIL was expressed by population of TH cells, which is likely to be enriched for the THIFNR subset (Fig. 3F). In published single-cell datasets (67), we found that CD4+ T cells in human lungs expressed interferon response signature genes and TRAIL, which indicated that the THIFNR subset is also present in the human lung tissue (fig. S8). We next experimentally tested TRAIL’s function and found that recombinant TRAIL inhibited TCR-dependent activation of TH cells ex vivo, as measured by the surface expression of the activation markers CD154, CD69, and CD137 (4-1BB) (Fig. 3, G and H) (68). These findings support a potential regulatory role of HDM-reactive THIFNR cells in dampening allergic responses.
A subset of HDM-reactive Treg cells express the interferon response signature
We also investigated whether HDM-reactive Treg cells differed between HDM-allergic and nonallergic individuals. As shown previously, we confirmed that the proportion of HDM-reactive Treg cells was not related to HDM allergy status (Fig. 1B) (43). Furthermore, transcriptomic analysis of bulk populations of HDM-reactive Treg cells revealed no major disease-related differences (fig. S9). To determine whether specific subsets of HDM-reactive Treg cells varied with disease state, we performed single-cell transcriptomic analysis of ~10,000 HDM-reactive Treg cells across the four participant groups, which separated this population into three distinct clusters (Fig. 4, A and B, and table S3). The proportion of cells in cluster 2 was greater in asthmatic individuals without HDM allergy compared with HDM-allergic asthmatics, suggesting preferential expansion of this subset in asthmatic individuals without HDM allergy (Fig. 4, C and D, and fig. S10). GSEA analysis of transcripts enriched in this cluster (n = 248) revealed significant enrichment of interferon response genes (Fig. 4E). The features of this Treg cluster were similar to those of the THIFNR cluster, which was also present at higher proportions in individuals without HDM allergy; for example, interferon-responsive Treg cells (TregIFNR) also expressed higher levels of transcripts encoding for TRAIL (Fig. 4F). Overall, our findings indicate preferential expansion of HDM-reactive Treg and TH cells, expressing the interferon response signature in asthmatic individuals without HDM allergy, and that expression of TRAIL by these subsets is likely to play an important role in dampening TH2 responses to HDM allergens, although further studies in animal models would be required to confirm this hypothesis.
(A) tSNE visualization of Seurat clustering analysis of the transcriptomes of 10,526 single HDM-activated Treg cells obtained from all 24 individuals. (B) Heat map showing hierarchical clustering of row-wise z-score–normalized mean expression of cluster-specific differentially expressed genes (n = 1559). Columns represent each Treg cluster. Right: Lists of biologically relevant examples genes for the TregIFNR cluster. Equal numbers of cells were sampled per disease group. (C) tSNE visualization of Seurat clustering analysis shown in Fig. 4A, using equal cell numbers for each participant group (n = 2180) obtained from all six individuals. (D) Pie charts illustrate the relative proportion of cells from each participant group within each of the three Treg clusters. (E) GSEA between each TregIFNR and the two other Treg clusters. The list of genes is provided in table S4. (F) Violin plots show log2(CPM + 1) normalized expression of TNFSF10 in each Treg cluster. Cells with no expression are excluded (Materials and Methods).
HDM-reactive TH2 cells are enriched for transcripts linked to enhanced functionality
Given the important role of TH2 cells in the pathogenesis of allergy and asthma, we analyzed the genes enriched in the TH2 cluster to gain insights into their functional properties. Gene coexpression analysis is a powerful method to find genes that are likely to play important roles in the differentiation or function of a given cell type (69, 70). Hierarchical gene clustering (Fig. 3A) and weighted gene coexpression network analysis (WGCNA) (Materials and Methods and Fig. 3B) (71) of the 214 “TH2 enriched” transcripts, defined from single-cell transcriptome analysis (Fig. 2, A and B), revealed five modules of highly coexpressed genes (Fig. 5A). Among these five modules, two modules (green and purple) contained genes linked to cellular metabolism, protein trafficking, active transcription, and oxidative phosphorylation (EIF3J, EIF5B, CALM3, FKBP1A, PTPN11, ATP13A3, PSMD13, UBE2S, and DUSP4), indicating increased metabolic and transcriptional activity in TH2 cells; a third (blue) module contained genes encoding for important transcription factors linked to TH2 cell differentiation such as GATA binding protein 3 (GATA3), interferon regulatory factor 4, and special A-T sequence-binding protein 1 (SATB1) (72–74).
(A) Hierarchical clustering of Spearman correlation coefficient matrix for SAVER-imputed expression values of the 214 genes uniquely up-regulated in the TH2 cluster. Values are clustered with complete linkage. The dotted red line indicates the Euclidean distance (Dist) threshold value used to define the five modules of coexpressed genes. Right: List of example genes for each module (modules 1 and 2 merged). (B) Gephi visualization of WGCNA for genes coexpressed in module 3 (top) and module 4 (bottom) from (A). Nodes correspond to a given gene and are sized on the basis of the number of edges (connections); edges thickness correlates to strength degree of correlation. (C) tSNE visualization of Seurat clustering analysis of single TH2 cluster cell transcriptomes (n = 1751) obtained from 12 allergic individuals regardless of asthma status [red, TH2 cluster 1 (n = 1440); purple, TH2 cluster 2 (n = 311)]. Red and purple circling lines represent limits of each TH2 clusters 1 and 2, respectively. (D) Heat map showing row-wise z-score–normalized mean expression of genes shown in (A) between both TH2 subclusters (columns). (E) Violin plots show log2(CPM + 1) normalized expression for genes biologically relevant between both TH2 subclusters. Color code represents the fraction of cells expressing the given gene in each TH2 subcluster; cells with no expression are not included. (F) GSEA between TH2 subclusters. Plots illustrate significative enrichment of module genes shown in (A) between both TH2 subclusters. The list of genes is provided in table S4.
The module including the canonical type 2 cytokine genes (IL5 and IL13) [red in Fig. 5 (A and B)] likely includes other genes that play an important role in driving the effector functions of TH2 cells. One of the most highly coexpressed transcripts encodes for the effector cytokine IL-9, which has recently been shown to be produced by a subset of TH2 cells that expressed peroxisome proliferator–activated receptor–γ (PPAR-γ) after transforming growth factor–β (TGF-β) signaling (75). We found that transcripts encoding for PPAR-γ and the TGF-β receptor 3 (TGFBR3) were also enriched in TH2 cells and highly coexpressed with IL9 (Fig. 5, A and B), suggesting that the IL-9 differentiation pathway is active in TH2 cells. Last, transcripts linked to cytotoxicity function (GZMB and RAB27A) (76–78) and differentiation of cytotoxic TH cells (ZEB2 and RUNX3) (46, 77, 79) were also highly coexpressed with IL5 and IL13, indicating that HDM-reactive cells may include cytotoxic TH2 cells (Fig. 5, A and B). Cytotoxic TH cells are known to contribute to antiviral immunity (77, 80) and autoimmunity (81), and our findings bring up the possibility that they may also play a role in asthma pathogenesis.
The gene for another canonical type 2 cytokine, IL-4, also linked to the function of T follicular helper cells (TFH) and IgE class switching, was present in a fourth module (yellow). Important molecules encoded by genes in this module include IL-31, a member of the IL-6 family of cytokines that is produced by activated TH2 cells and leads to itching in skin inflammation (Fig. 5, A and B) (82, 83), and IL-3, which is linked to hematopoietic progenitor proliferation and recruitment (84–86). We recently showed that IL-3 plays an important role in the activation and survival of eosinophils (87). Other genes in this module encode for products such as inducible T-cell costimulation protein (ICOS) and IL-21, which is linked to B cell help and immunoglobulin isotype class switching (Fig. 5, A and B), suggesting that this module was enriched for genes linked to TFH cell function. The presence of gene modules with distinct coexpression patterns indicated potential heterogeneity in the TH2 population. To address this issue, we reclustered cells only from the TH2 population; this analysis revealed two distinct subclusters (Fig. 5C), each highly enriched for genes in modules 3 and 4 (yellow and red) (Fig. 5, D to F).
Overall, these results show that cells in the TH2 cluster were enriched for the expression of transcripts encoding for several costimulatory and inhibitory receptors as well as transcription factors and molecules known to promote T cell survival. The expression of the first class of molecules, including CD28, ICOS, BTLA, CTLA-4, PD-1, and HVEM receptor (LIGHT and TNFSF14) (88), suggests that these molecules could be targeted to dampen the proinflammatory function of TH2 cells in asthma. Prosurvival factors included several in the nuclear factor κB signaling pathway, including NFKBID, NIPAI, MAP3K8, FOSL2, NEDD9 (89), ZEB2, BCL2A1 (90), BIRC3 (91), DUSP4/MKP-2 (92, 93), and CFLAR (cFLIP/CASPER) (94). Together, these expression patterns suggest that these cells are endowed with properties that allow them to exert sustained and strong type 2 inflammatory responses in asthma.
IL-9–expressing HDM-reactive TH2 cells are increased in asthma
We next sought to identify potential asthma-specific changes in HDM-reactive TH2 cells from individuals with HDM allergy. Single-cell differential gene expression analysis of HDM-reactive TH2 cells from HDM-allergic asthmatics versus nonasthmatics revealed that among the TH2-enriched transcripts, IL9 was the most up-regulated gene in asthmatic individuals (Fig. 6A). Furthermore, as shown in Figure 5 (D and E), subclustering of the TH2 subset showed that IL9-expressing cells were highly enriched in the larger TH2 subcluster (Fig. 6B). The relative proportion of cells in the IL9-enriched TH2 subset (TH2 cluster 1) was greater in asthmatics compared with (75) nonasthmatic individuals (Fig. 6C). Thus, enrichment of IL9 expression in TH2 cells appears to be associated in the development of asthma.
(A) Volcano plot shows statistical significance (−log10 adjusted P value) in function of the log2-fold difference in expression for filtered genes (see Materials and Methods), when comparing expression between TH2 cells from asthma allergic (n = 6) versus nonasthma allergic (n = 6) individuals. Dots are colored accordingly to the average of expression (log2) and sized on the basis of the fraction of cells expressing the given gene, both derived from the group in which the gene is up-regulated. Equal numbers of cells where sampled in each group (n = 661). Gray dotted lines represent the threshold value for fold change [vertical, log2(|fold change|) > 0.5-fold] and significance [horizontal, −log10(adjusted P value) > 2]. FDR, false discovery rate. (B) tSNE visualization of TH2 cluster cell transcriptomes shown in Fig. 5C in which each dot is a cell colored according to expression for IL9 (gray, none). Outlines represent TH2 subcluster limits. (C) Box-and-whisker plot shows percentage of TH2 cells in each subcluster in asthma allergic (n = 6) and nonasthma allergic (n = 5) individuals. Center line, median value; box, quartiles; whisker lines, 10th and 90th percentiles. ***P < 0.01. (D) Volcano plot similar to (A) comparing expression between IL9-positive cells (n = 444) and IL9-negative cells (n = 444) in the asthma allergic TH2 cluster 1.
To determine the properties of IL9-expressing HDM-reactive TH2 cells in asthmatic individuals, we compared the single-cell gene expression profiles of IL9-expressing and nonexpressing cells contained within the IL9-enriched TH2 cluster 1. We were surprised to find that expression of several transcripts encoding products linked to pathogenicity and survival of TH2 cells was increased in IL9-expressing cells (Fig. 6D and table S3). These included transcripts encoding canonical TH2 cytokine IL-5, cytotoxicity molecules [granzyme B, zinc finger E-box binding homeobox 2 (ZEB2), and EF-hand domain family member D2 (EFHD2)], TH2 polarizing and survival-related signaling receptor (IL-33R) (95, 96), and CD109, a membrane-anchored molecule described as negative regulator of TGF-β signaling (97, 98) but also as a coactivator of the Janus kinase/signal transducers and activators of transcription 3 signaling pathway (98–100) that is important for TH2 cell development (Fig. 6D, fig. S11, and table S3) (98, 101). Overall, these findings suggest that IL9-expressing HDM-reactive TH2 cells displayed greater pathogenic properties that could play an important role in driving asthma pathogenesis.
DISCUSSION
In this study, we present large-scale single-cell transcriptome analysis of allergen-reactive TH and Treg cells (n ~ 50,000) from individuals with asthma and/or allergy and healthy controls. Our characterization of HDM-reactive TH cells identified substantial heterogeneity as well as quantitative and qualitative changes related to allergy and asthma and revealed a unique subset of TH cells with a strong interferon response signature.
Our single-cell study addresses some of the important unanswered questions in the field of allergy and asthma research. Most fundamentally, exposure to common allergens, such as HDM, is nearly ubiquitous, and TH responses to these allergens are seen in both allergic and nonallergic individuals (33–35), yet why do only some people develop TH2 responses to allergens? By comparing HDM-reactive T cells from asthmatics with and without HDM allergy (TH2 responses versus no TH2 responses), we found that in individuals without HDM allergy (but sensitized to other allergens; see table S1), a subset of TH and Treg cells expressing the interferon response signature was expanded. These cells expressed higher levels of TRAIL, a molecule that can inhibit TCR signaling, activation of TH cells, and inflammation in model systems (102). Therefore, we hypothesize that these TRAIL-expressing HDM-reactive T cells could play an important role in dampening TH2 inflammation in allergy and asthma. A recent single-cell transcriptomic study identified TH cells expressing the interferon response signature in the lung tissue of mice sensitized and challenged with HDM (20). This finding implies that THIFNR cells can be generated and sustained in vivo and that HDM sensitization of mice may be an appropriate system to test the role of these cells in allergic inflammation. Future studies are also required to determine the molecular mechanisms and signals that drive the differentiation, maintenance, and persistence of THIFNR cells, ideally through the analysis of their molecular and chromatin landscape.
Another important question is why only some patients with allergy develop asthma. Do TH cells from allergic patients with asthma display distinct molecular features from those of allergic patients without asthma? By comparing HDM-reactive TH cells from allergic individuals with and without asthma, we defined a subset of IL9-expressing TH2 cells that are enriched in asthmatic individuals. We show that IL9-expressing HDM-reactive cells display several features likely to enhance their pathogenicity and persistence, which may contribute to asthma pathogenesis. In the context of peanut allergy, IL-9 was shown to best differentiate children with peanut allergy from children with peanut sensitization, who tolerate peanut, suggesting a potentially important role in food allergy (103). Studies blocking IL-9 activity in animal models of asthma have indicated that this may be a promising treatment approach (104, 105), but the relative importance of IL-9–producing T cells has not been fully explored. A recent report in cancer showed that IL-9–producing murine TH cells are more cytolytic, hyperproliferative, and less exhausted (106, 107); these properties conferred potent antitumor activity for these cells when tested in adoptive transfer experiments. Studies in mouse models of allergic airway inflammation are required to demonstrate the relative pathogenicity and persistence of IL-9–producing TH2 cells in vivo. Moreover, this cell population should also be characterized in humans with severe asthma, including those that reside in the airways. In summary, our single-cell transcriptomic study of HDM allergen–specific T cells has identified TH subsets that may contribute to the pathogenesis of allergy and asthma.
MATERIALS AND METHODS
Study design
The goal of this study was to use bulk and single-cell RNA-seq assay to capture the transcriptome of HDM allergen–reactive CD4+ T cell memory subsets from PBMCs of six asthmatic individuals allergic to HDM, six asthmatic individuals nonallergic to HDM, six nonasthmatic individuals allergic to HDM, and six nonasthmatic individuals nonallergic to HDM. To isolate HDM-reactive CD4+ cells, PBMCs were stimulated with HDM and CD4+ memory cells were sorted on the basis of CD154 and CD137 expression: CD154+ (HDM-reactive TH), CD154− CD137+ (HDM-reactive Treg), and CD154− CD137− (HDM-nonreactive T cells). For bulk RNA-seq, we collected 200 cells in triplicates, and for single-cell RNA-seq assay, we collected between 1500 and 2000 cells per cell type per patient (see table S1).
Participant recruitment, ethical approval, and characteristics
Recruitment of participants included in this study followed Institutional Review Board (La Jolla Institute for Immunology, La Jolla, CA) approval, and study participants gave written informed consent. Twelve nonsmoking individuals with mild asthma treated only with inhaled bronchodilators (mild asthma); six individuals with allergic rhinitis but no asthma, meeting established diagnostic criteria; and six healthy nonatopic individuals were studied (Fig. 1A and table S1). Participants with asthma underwent pulmonary function tests and/or methacholine challenge to establish diagnosis (bronchodilator response of ≥12%, or ≥ 200 ml, and/or methacholine challenge with a provocative concentration causing a drop of the forced expiratory volume in 1 s of 20%, 8 mg/ml). All participants were classified as allergic to HDM on the basis of skin test reactivity to HDM allergens [Dermatophagoides pteronissinus (Der p) and Dermatophagoides farinae (Der f); table S1].
Sample processing
We used PBMCs obtained from blood samples by density gradient centrifugation according to the manufacturer’s instructions and cryopreserved in liquid nitrogen.
Antigen-reactive T cell enrichment assay
The HDM peptide pool was generated as described and contained 75 peptides at a total concentration of 20 mg/ml (188.7 ng/ml for each peptide) (108, 109). The assay to isolate HDM-reactive T cells based on HDM peptide pool stimulation, magnetic activated cell sorting (MACS)–based enrichment, and fluorescence-activated cell sorting (FACS) of CD154+ memory CD4+ T cells from PBMCs was adapted from the work of Bacher et al. (43) and is outlined in fig. S1. For each donor, PBMC cryovials were thawed, washed, and plated overnight in six-well culture plates at a concentration of 10 × 106 cells/ml in 2 ml of serum-free TexMACS medium (Miltenyi Biotec) (5% CO2, 37°C). In the presence of a blocking CD40 antibody (1 μg/ml; Miltenyi Biotec), cells were then stimulated by the addition of HDM peptide pool (1 μg/ml) for 6 hours. Subsequently, cells were stained with fluorescence-labeled antibodies and a biotinylated CD154 antibody (clone 5C8; Miltenyi Biotec). Anti-biotin microbeads (Miltenyi Biotec) were added to allow MACS-based enrichment of CD154+ cells using MS columns (Miltenyi Biotec). Cells (5%) were kept as control sample (“input”) and used for FACS of HDM− T cells and analysis of cell frequencies before enrichment. Positively selected cells (CD154+) were eluted from the column and used for FACS of CD154+ memory CD4+ T cells. The flow-through from the MACS column was collected, stained with a biotinylated CD137 antibody (clone REA765; Miltenyi Biotec) and anti-biotin MicroBeads, and applied to a second MS column. Positively selected cells (CD137+) were used for FACS of CD137+ cells. All cell populations were FACS-sorted using a FACSAria II (Becton Dickinson); the gating strategy is shown in fig. S1. All flow cytometry data were analyzed using FlowJo software (version 10).
Cell isolation for bulk and single-cell RNA-seq assay
For bulk assays, cells of interest were directly collected by sorting 200 cells into 0.2-ml polymerase chain reaction (PCR) tubes (low retention; Axygen) containing 8 μl of ice-cold lysis buffer [Triton X-100 (0.1%; Sigma-Aldrich)] containing ribonuclease (RNase) inhibitor (1:100; Takara). Once collected, tubes were vortexed for 10 s, spun for 1 min at 3000g, and stored at −80°C. For single-cell RNA-seq assays (10x Genomics), 1000 to 2000 HDM-reactive T cells per individual were collected by sorting in low-retention and sterile ice-cold 1.5-ml collection tubes containing 500 μl of phosphate-buffered saline (PBS):fetal bovine serum (FBS) (1:1 v/v) with RNase inhibitor (1:100). HDM-reactive T cells from six individuals in each of the four groups (asthma with HDM allergy, asthma without HDM allergy, HDM allergy without asthma, and healthy without HDM allergy) were collected in the same tube. Collection tubes with ~9000 to 12,000 sorted cells per study group were inverted a few times; ice-cold PBS was added to reach a volume of 1400 μl, and they were centrifuged for 5 min at a speed of 600g at 4°C. Supernatant was cautiously removed, leaving 5 to 10 μl of volume. Pellets were then resuspended with 25 μl of 10x Genomics resuspension buffer [0.22-μm filtered ice-cold PBS supplemented with ultrapure bovine serum albumin (0.04%; Sigma-Aldrich)]. A total of 33 μl of cell suspension was transferred to an 8-PCR-tube strip for downstream steps as per the manufacturer’s instructions (10x Genomics, San Francisco).
Bulk RNA library preparation for sequencing
For full-length bulk transcriptome analyses, we used the Smart-seq2 protocol (110), adapted for samples with small cell numbers (111, 112). We followed the protocol as described previously (111, 112) with following modifications: (i) The preamplification PCR cycle for T cells was set at 22 cycles. (ii) To eliminate any traces of primer-dimers, the PCR-preamplified cDNA product was purified using 0.8× AMPure XP beads (Beckman Coulter) before using the DNA for sequencing library preparation. Preamplified cDNA (1 ng) was used to generate barcoded Illumina sequencing libraries (Nextera XT library preparation kit, Illumina) in 40-μl reaction volume. Samples failing any quality control step [DNA quality assessed by capillary electrophoresis (Fragment Analyzer, Advanced Analytical) and quantity (PicoGreen quantification assay, Thermo Fisher Scientific)] were eliminated from further downstream steps. Libraries were then pooled at equal molar concentration and sequenced using the HiSeq 2500 Illumina platform to obtain 50–base pair (bp) single-end reads (HiSeq SBS Kit v4; Illumina). In total, 1.7 billion uniquely mapped reads were generated with a median ± SD of 17.8 ± 3.8 million uniquely mapped reads per sample.
10x Genomics single-cell RNA library preparation for sequencing
Samples were processed using 10x Genomics 3′TAG v2 chemistry as per the manufacturer’s recommendations; eleven cycles were used for cDNA amplification and library preparation, respectively (77). Barcoded RNA was collected and processed following the manufacturer’s recommendations. After quantification, equal molar concentration of each libraries was pooled and sequenced using the HiSeq 2500 Illumina sequencing platform to obtain 26- and 100-bp paired-end reads using the following read length: read 1, 26 cycles; read 2, 100 cycles; and i7 index, 8 cycles.
Bulk RNA-seq analysis
Bulk RNA-seq data were mapped against the hg19 reference using TopHat (113) [v1.2.1 (--library-type fr-unstranded --no-coverage-search) with FastQC (v0.11.2), Bowtie (v1.1.2) (114), and SAMtools 0.1.18.0] (115), and we used htseq-count -m union -s no -t exon -i gene_name (part of the HTSeq framework, version 0.7.1) (71, 116). Cutadapt (v1.3) was used to remove adapters (117). Values throughout are displayed as log2 TPM (transcripts per million) counts; a value of 1 was added before log transformation (pseudocount). We performed principal component (PC) analysis and clustering analysis using t-distributed stochastic neighbor embedding (tSNE) dimensional reduction algorithm (based on three PC using the top 200 most variable genes) (118). To identify genes expressed differentially between groups, we performed negative binomial tests for paired comparisons by using DESeq2 (1.16.1) (119) with default parameters. We considered genes to be expressed differentially by any comparison when the DESeq2 analysis resulted in a Benjamini-Hochberg–adjusted P value of at most 0.01 and a log2 fold change of at least 2. GSEA were performed as previously described (112, 120) using the Qlucore visualization software (version 3.5) (121). Gene lists used for GSEA analysis are shown in table S4.
Single-cell RNA-seq analysis
Analysis of 3′ single-cell transcriptomes using the 10x Genomics platform
Raw data were processed as previously described (70, 77), merging multiple sequencing runs using count function from Cell Ranger (table S3) and then aggregating multiple cell types with cell ranger aggr (v3.0.2). The merged data were transferred to the R statistical environment for analysis using the package Seurat (v3.0.2) (122).
Doublet cell filtering
Barcoded single-cell RNA-seq was demultiplexed patient-wise using Demuxlet (123) with the following parameters: alpha = 0, 0.5 and --geno-error = 0.05. Cells called as doublet by Demuxlet were removed from downstream analyses (fig. S3A and table S3). Identities were inferred by analyzing VCF files from the genotyping analysis containing the corresponding individual for each particular library. Each cell was assigned a donor ID or marked as a doublet and then incorporated to the annotation table. We did not observe major changes in singlets/doublets proportion between the different 10x Genomics libraries (reflecting cell type and participant groups), suggesting optimal processing of cells during 10x GEM (Gel Bead-In Emulsions) generation and downstream steps (fig. S3A). All downstream analyses were performed using singlet cells.
Activation score and cell filtering
To filter out cells with low level of activation or no activation by the HDM peptide pool, i.e., HDM nonreactive cells, we performed pairwise single-cell differential gene expression analysis using a model-based analysis of single-cell transcriptomics (MAST) algorithm (124) between HDM-reactive TH cells (CD154+ cells, n = 3075, random sampled) and HDM− T cells (CD154neg cells, n = 3075). We defined a gene set (n = 110 genes; called HDM+ TH activation genes) that captured the transcripts up-regulated in the CD154+ versus CD154neg cells (table S3) using the following parameters (fig. S3B): Benjamini-Hochberg–adjusted P ≤ 0.05, log2 fold change ≥ 2, log2 mean of expression ≥ 0.75 counts per million (CPM), and percentage of expressing cells (>0 CPM) in HDM+ TH cells > 37.5% (fig. S3B). We scored each cell using AddModuleScore from Seurat (122). Briefly, the module score is calculated by binning the genes by the average expression level; then, the average expression of each gene is subtracted by the aggregated expression of the control gene sets (100) randomly selected per bin. Last, on the basis of the distribution of cells based on their activation score (Fig. 2A and fig. S3C), we applied a threshold for defining activated cells. The proportions of cells expressing important canonical genes such as IL4, IL5, IFNG, and IL17A pre- and postactivation filtering indicated that cell eliminated because of low activation score did not up-regulate transcripts for these genes. Similarly, an independent HDM+ Treg activation score was also calculated using a similar approach to analyze HDM+ Treg single-cell datasets (table S3 and fig. S3, right). In total, 3505 cells with low activation score were eliminated from downstream analysis.
Transcriptome-based clustering analysis
The merged data were transferred to the R (v3.5) statistical environment for analysis using the package Seurat (v3.0.2) (122). Only cells expressing more than 200 genes and with a total mitochondrial gene expression less than 5%, and genes expressed in at least three cells were included in the analysis. The data were then log-normalized per cell and a list of the most variable genes with a mean expression > 0.1 (unique molecule identifier, UMI) and explaining 30% of the cumulative standardized variance given by the FindVariableFeatures function were used for clustering analysis (fig. S4). We performed clustering analysis using distinct lists of most variable genes for HDM+ TH, HDM+ Treg and HDM+ TH2 clusters (table S3). In regard to TH2 subclustering, because of the limited number of cells, we considered most variable genes that were expressed by more than 10% of the cells and with a standardized variance greater than 2. We also limited the selection of the most variable genes to 10% of the cumulative standardized variance (fig. S4, right). Normalized single-cell transcriptomic data were then further scaled by the number of UMI-detected and percentage of mitochondrial reads. We then performed PC analysis with RunPCA algorithm (122) using the determined most variable gene lists. We followed Seurat procedure to determine the number of PCs to select for downstream analyses, using the SD of PCs. We applied FindNeighbors and FindClusters functions from Seurat with default settings (table S3 and fig. S4) to identify clusters. All clusters had more than 50 cells and none were excluded from the downstream analysis. Cluster specific markers were obtained by the FindAllMarkers function with default parameters, test.use = MAST (124). Further visualizations of exported normalized data such has “violin” plots were generated using the Seurat package and custom R scripts. Our violin plots show Seurat normalized expression for a particular gene [log2(CPM + 1 pseudo-count)] only for the cells expressing the gene of interest. Violin shape represents the distribution of cell expressing transcript of interest (based on a Gaussian Kernel density estimation model) and is colored according to the percentage of cell expressing the transcript of interest.
Single-cell differential gene expression analysis
Pairwise single-cell differential gene expression analysis was conducted after conversion of data to count per million base pairs (CPM + 1) using MAST algorithm (q < 0.05, v1.2.1) (R package) (124). We used equal number of cells from each participant group and random sampling performed when necessary. A gene was considered differentially expressed when Benjamini-Hochberg–adjusted P value was <0.05 and a log2 fold change was more than 0.25. For cluster-specific signatures, a gene was considered significantly different (unique to a group), only if the gene was enriched in every pairwise comparison for a particular cluster with other clusters.
Single-cell coexpression analysis and WGCNA
To perform coexpression analysis, given the high levels of genes drop-out associated with single-cell analysis, we performed a data imputation using a single-cell analysis via expression recovery (SAVER) imputation algorithm (125). Briefly, SAVER analysis was implemented on the Cell Ranger UMI matrix output for HDM+ TH using the function SAVER (v1.1.1) with pred.genes.only = TRUE. Then, we calculated Spearman correlations coefficients using the cor function and determined the cluster-modules through hclust on Euclidean distances and cutree functions k = 5 according to the within groups sums of squares elbow (similar to TH single-cell clustering analysis, fig. S4). We also performed weighted correlation analysis using WGCNA algorithm (v1.66) (71) using the function TOMsimilarityfromExpr, power = 3, and exportNetworkToCytoscape, weighted = TRUE, threshold = 50th quantile of the topological overlap matrix. Network plots were generated by Gephi (0.9.2) using Fruchterman-Reingold and Noverlap layouts (126). The size and color of the nodes were defined according to the degree, while the edge width and color were scaled according to the weight value.
Genotyping
For each patient, genomic DNA was isolated from PBMCs using the DNeasy Blood and Tissue Kit (Qiagen) and used for genotyping using the Infinium Multi-Ethnic Global-8 Kit (Illumina) following the manufacturer’s instructions. Raw data from the genotyping analysis, data quality assessment, and single-nucleotide polymorphism identification were performed as previously described (127).
Stimulation of memory CD4+ T cells with human recombinant TRAIL
Memory CD4+ T cells were isolated from PBMCs using the Memory CD4+ T Cell Isolation Kit (Miltenyi Biotec) and cultured in Iscove’s modified Dulbecco’s medium (IMDM; Invitrogen) supplemented with 5% (v/v) heat-inactivated FBS and 2% (v/v) human AB serum (Cellgro). The memory CD4+ T cells were stimulated with precoated human recombinant TRAIL (5 μg/ml), anti-CD3 antibodies (2.5 μg/ml), and soluble anti-CD28 antibodies (1 μg/ml) in the presence of IL-7 (5 ng/ml; Miltenyi Biotec). The expression of surface markers (CD69, CD154, and CD137) was analyzed by flow cytometry after 6 hours.
Expression of TRAIL on memory CD4+ T cells
Total PBMCs were thawed, washed, and plated overnight in serum-free TexMACS medium (Miltenyi Biotec)/complete IMDM (as described above). In the presence of a blocking CD40 antibody (1 μg/ml in culture; clone HB14; Miltenyi Biotec), cells were then left untreated or stimulated by the addition of the control reagent CytoStim (1:500 dilution of stock; Miltenyi Biotec). The expression of surface markers (CD69, CD154, and TRAIL) was analyzed by FACS after 6 hours.
Statistical analysis
We used nonparametric Kruskal-Wallis one-way analysis of variance (ANOVA) test to compare unpaired data for more than two conditions and Kolmogorov-Smirnov test when comparing two groups of data. We used paired Student’s t test for time course flow cytometry analysis. We used GraphPad Prism 7.0. All source datasets and statistical tests used are detailed in table S5.
SUPPLEMENTARY MATERIALS
immunology.sciencemag.org/cgi/content/full/5/48/eaba6087/DC1
Fig. S1. Experimental design and gating strategy to isolate HDM-reactive memory TH and Treg cells.
Fig. S2. Expression of differentially expressed genes in HDM-reactive cells.
Fig. S3. Single-cell filtering of doublets and low HDM-reactive T cells.
Fig. S4. Single-cell clustering analysis using Seurat.
Fig. S5. Distribution of cell frequency for the seven HDM+ TH clusters for the 24 individuals.
Fig. S6. Coexpression of TH1- or TH17-specific signature genes.
Fig. S7. Donor and disease groups cell distribution in each of the HDM+ TH clusters.
Fig. S8. Coexpression of THIFNR-specific signature genes in publicly available dataset.
Fig. S9. Treg disease-related differences.
Fig. S10. Proportions of HDM-reactive Treg subsets among disease groups.
Fig. S11. TH2 single-cell analysis.
Table S1. Study participant details (Excel spreadsheet).
Table S2. Differential gene expression analysis in bulk populations (Excel spreadsheet).
Table S3. Single-cell analysis (Excel spreadsheet).
Table S4. Gene signature lists (Excel spreadsheet).
Table S5. Raw data file with statistics (Excel spreadsheet).
This is an article distributed under the terms of the Science Journals Default License.
REFERENCES AND NOTES
- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works
