Heterogeneity and clonality of kidney-infiltrating T cells in murine lupus nephritis

Abstract

We previously found that kidney-infiltrating T cells (KITs) in murine lupus nephritis (LN) resembled dysfunctional T cells that infiltrate tumors. This unexpected finding raised the question of how to reconcile the "exhausted" phenotype of KITs with ongoing tissue destruction in LN. To address this, we performed single-cell RNA-Seq and TCR-Seq of KITs in murine lupus models. We found that CD8+ KITs existed first in a transitional state, before clonally expanding and evolving toward exhaustion. On the other hand, CD4+ KITs did not fit into current differentiation paradigms but included both hypoxic and cytotoxic subsets with a pervasive exhaustion signature. Thus, autoimmune nephritis is unlike acute pathogen immunity; rather, the kidney microenvironment suppresses T cells by progressively inducing exhausted states. Our findings suggest that LN, a chronic condition, results from slow evolution of damage caused by dysfunctional T cells and their precursors on the way to exhaustion. These findings have implications for both autoimmunity and tumor immunology.

Keywords: Autoimmune diseases; Autoimmunity; Immunology; T cells.

Figure 1. Low-resolution clustering of T cells from lupus-prone mice.

(A) Schematic of experimental design encompassing tissue source, sorting algorithm, marking with HTOs, number of mice used in each experimental cohort, and a description of downstream analytic techniques. (B) UMAP of T cells from mice in Main-Seq, outlining 14 clusters. (C) HTO-based assignment of cell source was mapped onto the UMAP. MRL/lpr, MRL.Fas^lpr; B6, C57BL/6; UMAP, uniform manifold approximation and projection; TF, transcription factor.

Figure 2. High-resolution clustering of CD4⁺ T cells uncovers unique transcriptional programs in KITs.

(A) UMAP of CD4⁺ T cells from Main-Seq outlining 13 clusters, with the right panel exhibiting assignment of cell source as determined by HTO. (B–D) Gene set enrichment analysis (GSEA) performed in each cell by Wilcoxon’s test (–log₁₀ [P value]) using published reference gene signatures (Supplemental Table 1) and related TF expression were overlaid onto the UMAP to identify CD4 phenotypes. This included (B) Treg gene signature and Foxp3 expression, (C) Th1 signature and Tbx21/Tbet expression, and (D) Th2 signature and Gata3 expression.

Figure 3. TF analysis suggests overarching transcriptional regulation of infiltrating CD4⁺ T cell clusters.

(A) Heatmap representing z-scored regulon activity of top TFs inferred by SCENIC and association with CD4⁺ T cell clusters. (B) Expression of selected TFs overlaid onto the CD4 UMAP as depicted in Figure 2.

Figure 4. CD4⁺ KITs exhibit a progressive transcriptional phenotype from cytotoxicity to hypoxia/dysfunction through pseudotime.

GSEA performed in each cell by Wilcoxon’s test (–log₁₀ [P value]) using published reference genes signature (Supplemental Table 1) overlaid onto the UMAP from Figure 2A. (A) Hypoxia signature, exhaustion signature, and cytotoxic CD4 signature. (B) Dot plots show the distribution of exhaustion score in each CD4⁺ T cell, grouped by cluster number. Dots are colored according to the source of cells they represent. Statistics were calculated by Kruskal-Wallis rank test. (C) Monocle trajectory mapping of CD4⁺ KITs wherein time 0 (dark purple) represents lineage origination with progression to most differentiated (yellow), with cell source mapping. (D) Gene signature mapping for the indicated signatures as defined in A.

Figure 5. Kidney-infiltrating Tregs exhibit features of tissue reprogramming.

(A) High-resolution reclustering of Tregs identified as cluster 6 in Figure 2 (Main-Seq) identifies 3 unique clusters as illustrated by color-coding. (B) HTO-based identification of splenic (B6 and MRL/lpr) and KIT Tregs. (C) Volcano plot shows top significant (FDR < 0.05) DEGs that are upregulated or downregulated in KITs compared with splenic MRL/lpr Tregs, with genes in pink having been previously associated with tissue-resident Tregs (30).

Figure 6. High-resolution clustering of CD8⁺ T cells identifies unique functional phenotypes.

(A) UMAP of CD8⁺ T cells from Main-Seq delineating 9 clusters. (B) UMAP with overlay exhibiting assignment of cell source as determined by HTO. (C) Heatmap shows top significant (FDR < 0.01) DEGs associated with each CD8⁺ cluster and their expression at single-cell level in columns. (D–I) GSEA performed in each cell by Wilcoxon’s test (–log₁₀ [P value]) using published reference gene signatures (Supplemental Table 1), overlaid onto the UMAP to identify CD8 phenotypes. Clusters are outlined as per A. This included gene signatures for (D) resident memory (T_RM), (E) circulating/lymphoid (T_periph), (F) exhaustion (T_EX), (G) effector memory (T_EM), (H) central memory (T_CM), and (I) hypoxia.

Figure 7. TF and lineage progression analysis of CD8⁺ KITs.

(A) Heatmap represents z-scored regulon activity of top TFs inferred by SCENIC (rows) and association with CD8⁺ T cell clusters (columns). (B) Representative TFs were mapped onto CD8 UMAPs; TF selection was based on known regulatory functions or due to identification via SCENIC analysis. Outlines highlight spleen-derived, kidney-infiltrating, and exhausted cells, with dot red color intensity representing log₂ expression. (C) CD8⁺ T cells grouped into 9 distinct clusters and ordered by Slingshot pseudotime trajectory. (D) Slingshot lineage overlay on CD8⁺ T cell UMAP.

Figure 8. CD8⁺ KITs are clonally expanded, with clones and proliferation spanning the exhausted and transitional compartments.

(A) UMAP of all 3 scRNA-Seq cohorts (see Figure 1 for cohorts), integrated using Harmony, followed by projection of individual cohorts onto this “combined UMAP.” (B) Main-Seq–defined CD8⁺ clusters labeled by cluster name (Figure 6), mapped onto the combined UMAP. (C) Combined UMAP with all putative clusters as outlined in B. (D) High-frequency clones (defined as clones representing the top quartile of expressed TCRs) from each cohort are mapped onto the combined UMAP. (E) Exhaustion gene set enrichment calculated using Wilcoxon’s test overlaid onto the combined UMAP. (F) Dot plots represent exhaustion scores for cells grouped based on clonal frequency among MRL/lpr and FcγR2B^–/–.Yaa KITs, with exhaustion scores for B6 naive, early T_EX, and terminal T_EX shown at left for reference (*P = 0.05, **P < 0.01, ****P < 0.0001 as determined by 1-way ANOVA with Tukey’s test for multiple comparisons). Horizontal bars represent medians. (G) Pie charts represent the relative cluster distribution of unique TCR T cells as compared with high-frequency TCRs from MRL/lpr (top) and FcγR2B^–/–.Yaa (bottom) models with the relative distribution of these clones within the putative T cell clusters as defined in B and C. (H) Cell cycle state was assessed over all cells for each cluster using GSEA for genes indicative of G1, G2/M, and S phases. Proliferative potential was analyzed for enrichment of G2/M and S phase genes, comparing individual clusters with the total T cell population. (G and H) Data were analyzed using χ² analysis and corrected for multiple comparisons for 9 comparison groups (*P = 0.05, ***P < 0.005, ****P < 0.0001). Green indicates enrichment of G2/M/S phase genes, and red indicates enrichment of G1 genes.