Single-cell RNA sequencing identifies inflammatory tissue T cells in eosinophilic esophagitis

Abstract

T cell heterogeneity is highly relevant to allergic disorders. We resolved the heterogeneity of human tissue CD3+ T cells during allergic inflammation, focusing on a tissue-specific allergic disease, eosinophilic esophagitis (EoE). We investigated 1088 single T cells derived from patients with a spectrum of disease activity. Eight disparate tissue T cell subtypes (designated T1-T8) were identified, with T7 and T8 enriched in the diseased tissue. The phenotypes of T7 and T8 resemble putative Treg (FOXP3+) and effector Th2-like (GATA3+) cells, respectively. Prodigious levels of IL-5 and IL-13 were confined to HPGDS+ CRTH2+IL-17RB+FFAR3+CD4+ T8 effector Th2 cells. EoE severity closely paralleled a lipid/fatty acid-induced activation node highlighted by the expression of the short-chain fatty acid receptor FFAR3. Ligands for FFAR3 induced Th2 cytokine production from human and murine T cells, including in an in vivo allergy model. Therefore, we elucidated the defining characteristics of tissue-residing CD3+ T cells in EoE, a specific enrichment of CD4+ Treg and effector Th2 cells, confinement of type 2 cytokine production to the CD4+ effector population, a highly likely role for FFAR3 in amplifying local Th2 responses in EoE, and a resource to further dissect tissue lymphocytes and allergic responses.

Keywords: Allergy; Immunology.

Figure 1. The schematic overall strategy to study tissue lymphocytes.

(A) Flow chart of the overall single-cell platform used to study human tissue lymphocytes, which consists of 3 basic modules: flow cytometry (FACS), bulk CD3⁺ RNA-seq, and C1 Fluidigm–based single-cell RNA-seq. (B) CD3⁺ lymphocytes isolated from esophageal biopsy tissue and autologous blood were analyzed by FACS and substratified into CD4⁺ and CD8⁺ populations.

Figure 2. FACS analysis of CD3⁺ T cells in blood and tissue.

(A) The number of CD4⁺, CD8⁺, and pan CD3⁺ T cells isolated from a single biopsy were enumerated from 4 different cohorts: normal (N, n = 9), complete remission (CR, n = 6), partial remission (PR, n = 12), and active inflammation (A, n = 13). Two-way ANOVA, EoE factor: P < 0.001; CD4/CD8 factor: P < 0.001; interaction, nonsignificant. (B) With the same 4 disease cohorts, the CD4/CD8 ratio was quantified from autologous blood and tissue samples. Two-way ANOVA, EoE factor nonsignificant; tissue/blood factor: P < 0.001; interaction, nonsignificant. (C) The percentages of CD25⁺ events were quantified by normalizing to gated CD4⁺ or CD8⁺ events. CD25⁺/CD4⁺ 2-way ANOVA, EoE factor: P < 0.001; tissue/blood factor: P < 0.001; interaction, P < 0.01. Bonferroni test: P < 0.05, normal controls (N) versus PR; P < 0.001, N versus A. (D) The percentage of CD25⁺CD4⁺ cells was correlated with tissue eosinophilia, but not autologous blood eosinophilia, with the P value and Spearman r shown. (E) Tissue and autologous blood CD3⁺CD4⁺ events (blue) plotted in the context of multiple proteins of interest, with CD25⁺ events labelled in red. Two-way ANOVA was used for multi-comparisons in C, D and E. Box and whisker plots are shown indicating mean and minimum/maximum values.

Figure 3. Bulk CD3⁺ RNA-seq analysis of tissue lymphocytes.

(A) CD3⁺ T cells were sorted from autologous blood and tissue from control (normal [N], n = 3) individuals and patients with EoE with active (active [A], n = 4) or inactive (complete remission [R], n = 4) EoE for bulk RNA-seq analysis. The bidirectional expression profile difference between blood and tissue T cells are displayed in the heat diagram, with each column representing a patient sample. A set of 331 genes reflected tissue-specific functions (paired moderated t test; FDR-adjusted P < 0.05, fold change >10). (B) Gene expression (150 genes) tracking with disease activity, 147 of which were significantly upregulated and 3 of which were downregulated in active EoE tissue lymphocytes compared with healthy controls (moderated t test, FDR-adjusted P < 0.05, fold change >10). (C) Among the 147 genes in B, the expression trends of 6 Th2-related genes across tissue and disease states are shown (blue: autologous blood CD3⁺; red: tissue CD3⁺). (D) The top 10 genes (fold change [FC] = 323–1449, NL tissue CD3⁺ vs. EoE tissue CD3⁺) from the volcano plot are displayed in a heat diagram. The expression pattern among the 6-comparison group is summarized in the adjacent dot plots showing 3 distinct trends labelled by roman numerals I–III. The fold change (FC, blue) and the FDR-adjusted P value (FDR-P, orange) were exhibited in a radar map in the context of the 10 top genes.

Figure 4. scRNA-seq identification of 8 lymphocyte subpopulations.

(A) Left panel: with t-SNE dimensionality reduction, the T1–T8 clusters were 2D plotted with each tissue T cell subclass color coded as indicated. Right panel: a 2D plot was color coded in the context of patients with and without EoE (non-EoE group includes NL and EoE remission). (B) With the large rectangle representing all single T cells isolated from NL, EoE remission, and active EoE subjects, the cellularity presence of T1–T8 subclasses was proportionally graphed within the disease status bin that was set to the equal area as 100%. (C) To emphasize the expression profiles of T7 and T8, a heatmap was selectively generated focusing on the T7-specific Treg genes and the T8-specific Th2-functionality genes, as well as those genes common to both clusters. Four patterns were identified on the basis of differential expression. Patterns I, II, III, and IV were enriched with distinct T clusters—I with T7+T8, II with T7 only, III with ubiquitous T clusters but T8 enriched, and IV with T8 only. (D) A volcano plot comparing top differentially expressed genes between T7 and T8, with key significant genes labeled. Colorized entities represent passing the filter of FDR-adjusted P < 0.05 and fold change >16 with red-blue gradient on FDR-adjusted P values. (E) A series of radar plots based on relevant functional gene sets was generated on logarithm scale to depict the expression characteristics of T1–T8 restricted to Th2 cell markers and cytokines, Treg properties, migration molecule expression related to chemotaxis and adhesion, and T cytotoxic cells properties, respectively. T7 (blue line) and T8 (red line) are shown in bold to emphasize the 2 disease-associated clusters.

Figure 5. Single-cell Th2 cytokine mRNA analyses.

To characterize the Th2 cytokine expression at mRNA and protein levels, relevant gene expression of 1088 tissue T cells at the single-cell level was analyzed along with independent FACS analyses. (A) Scrutinizing the human 5q31-33 Th2 cytokine cluster from IL5 to KIF3A, the RPKM scatter plots of all 5 adjacent genes in the context of disease activity (N, normal; R, remission; A, active EoE) are shown with each dot representing a single cell. The expression pattern for Th1 cytokine INFG (interferon γ) was plotted as a function of disease activity (mean ± SEM, *P < 0.05, **P < 0.01, ***P < 0.001; 1-way ANOVA with Bonferroni’s multiple comparison test). (B) Boolean logical analysis of Th2 cytokine–producing pattern of all Th2-capable cells (n = 72, center red circle represents 7 triple-cytokine–positive cells). (C) Three-way quantification of cytokine load by CD3⁺ T cells from allergic inflamed tissues encompassing the 1088 T cells. (D) Pie-chart break down of tissue T cell cytokine abundancy on the basis of single-cell RPKM enumeration.

Figure 6. Single-cell Th2 cytokine protein analyses.

(A) Top row: IL-4, IL-5, IL-13, and IFN-γ protein was assessed by FACS in tissue T cells from 29 subjects across 3 disease states. Major Th1/Th2 cytokine abundancy quantified (as % CD4⁺) in the context of disease activities. N, normal; R, remission; A, active EoE (mean ± SEM, *P < 0.05, **P < 0.01, ***P < 0.001, 1-way ANOVA with Bonferroni’s multiple comparison test). Bottom row: linear correlation between Th2 cytokine level and tissue eosinophilia, blood versus tissue. P values and Spearman r shown next to each correlation line. (B) All CD3⁺ tissue lymphocytes from EoE were pooled for a SPADE (spanning-tree progression analysis of density-normalized events) analysis, identifying a dichotomy of tissue lymphocytes (TL) driven by CD8-CD4. The developmental trees for major EoE tissue cytokines were juxtaposed to representative FACS plots to reveal their levels, T cell subtype source, and developmental patterns of each cytokine, with the heat color indicating expression levels and dot size indicating the cellular abundancy.

Figure 7. TCR clonotype analyses.

(A) EoE TCR repertoire diversity (as represented by Shannon entropy normalized to the total number of cells) is shown as a function of disease activities (normal control [N], remission [R], active EoE [A]) (mean ± SEM, 1-way ANOVA with Bonferroni’s multiple comparison test). (B) Stacked bar charts representing stratified diversity for each TCR repertoire, in which the clonal occurrence of tissue TCR repertoire for all patients across 3 different disease activities were enumerated: unique (n = 1, cyan), duplicated (n = 2, red), and clonal (n > 2, blue). (C) The TCRB CDR3 lengths from sequence samples from NL, remission, and active EoE individuals were analyzed (statistical legend same as A). (D) The TCR repertoires of cells classified as T7 or T8 are distinct from T1–T6 as measured by the Morisita-Horn dissimilarity index. Additionally, the T7 and T8 repertoires are distinct from each other in that they have no TCR clonal overlap. The Morisita-Horn dissimilarity index measures how dissimilar any 2 populations are based on overlapping species and the species’ respective hierarchies. (E) Principal component analysis (PCA) of TRBV-TRBJ usage shows that remission samples largely segregate from patients with active allergic inflammation and NL control samples based on TRBV-TRBJ usage along dimension 1 (Dim1) and dimension 3 (black outlined shape, centroid). (F) On the basis of the PCA in E, a superimposed projection of all detected TRBV-TRBJ pairs indicating their disease state enrichment (remission sector: upper 2 quadrants; active and NL sectors: lower quadrants), shown by individual pair’s Eigen vectors.

Figure 8. FFAR3 induction and Th2-enhancing effects of butyrate.

(A) A magnified view in the upper panel emphasizing the IL5-correlating genes specifically present in the T8 cluster, with FFAR3 tracking with IL5 expression shown in the lower panel (orange, IL5 producers; blue, FFAR3⁺ T cells). (B) Upper panel: the single-cell expression patterns of FFAR3 across 3 disease activity states (N, normal; R, remission; A, active EoE), with each data point representing 1 of the 1088 tissue T cells. Lower panel: with 3 days of Th2 differentiation (anti–CD3-CD28 [activated] with or without IL-4), human blood CD4⁺ T cells isolated from patients with EoE upregulated FFAR3 transcript induced by αCD3-CD28 (activation) plus IL-4 but not αCD3-CD28 (activation) alone. (C) IL5 transcript production by human Jurkat cells following short-chain fatty acid exposures (C2, C3, and C4 each at 10 mM for 24 hours). (D) Jurkat cells express FFAR3 on the cell surface by immunofluorescence and FACS. Original magnification, ×400. (E) In a murine asthma model, eGFP–IL-4 reporter mice were aspergillus allergen–challenged (ASP-challenged) intranasally with and without C4 (1 mg coadministration). IL-4–eGFP expression in CD4⁺ cells was analyzed by FACS. (F and G) Lung tissue CD4⁺ lymphocytes were assayed for Th1 and Th2 cytokine production by FACS. (H) The bronchoalveolar fluid (BALF) cells were quantified for major leukocyte populations. EOS, eosinophils; PMN, polymorphonuclear neutrophils; LYM, lymphocytes; MΦ, macrophages. All scatter plots are presented as mean ± SEM, and all experiments were repeated at least 3 times (*P < 0.05, **P < 0.01, ***P < 0.001; 2-tailed t test).