Fig 2: Expansion of CD161−CD8+TEMRAcells in SLE. (a) Percentages of CD161−CD8+ TEMRA and CD161+CD8+ TEMRA cell subclusters in HCs and SLE samples. (b) Two-dimensional plots showing the changes of CD161 expression over pseudotime in HCs and SLE samples. P value is calculated by Spearman correlation. (c) CD161 expression in CD161−CD8+ TEMRA and CD161+CD8+ TEMRA cell subcluster (left), and relative abundance of these two subsets in HCs and SLE subjects from public dataset (GSE162577) (right). P value is calculated by Wilcoxon rank-sum test. (d) CD161 expression of PBMCs from a new case–control cohort (HCs, n = 30; SLEs, n = 25). P value is from Unpaired t-test (e and f) Flow cytometry (left) and quantification (right) of the frequency of CD161−CD8+ TEMRA cells (e) cytoplastic levels granzyme B of CD161−CD8+ TEMRA (f) from a validation cohort (HCs, n = 23; SLEs, n = 33). P value is from Unpaired t-test. (g) Representative histogram showing granzyme B expression in HC and SLE samples. (h) Comparison of granzyme B expression in CD161−CD8+ TEMRA and CD161+CD8+ TEMRA cells from SLE samples (n = 33). P value is from Paired t-test. (i) Boxplot showing SLEDAI scores of SLE patients with high levels of CD161 and low levels of CD161 from public dataset (GSE121239). P value is calculated by Wilcoxon rank-sum test (j) Sankey plot representing the relative abundance of CD161−CD8+ TEMRA cells and CD161+CD8+ TEMRA cells from case–control cohorts of pSS (left, GSE177278) and MS (right, GSE193770).

Fig 3: CD161−CD8+TEMRAcells contribute to SLE by LIGHT signaling. (a) The number (top) and strength (bottom) of interaction among all cells in HCs and SLE samples. (b) Heatmap of differential interactions between HCs and SLE samples in cell–cell communication network. The top bar indicates the sum of incoming signaling and right bar indicates the sum of outgoing signaling. Red indicates increased signaling and blue indicates decreased signaling in SLE. (c) Number of significant ligand-receptor pairs between CD161−CD8+ TEMRA cells (outgoing) and other cell subclusters (incoming) in HCs (left) and SLEs (right). The relative number of ligand-receptor pairs is represented by the edge width. (d) Comparison of the significant outgoing signaling from CD161−CD8+ TEMRA cells between HC and SLE. Empty space means the communication probability is zero. P-values are computed from one-sided permutation test. ∗∗∗P < 0.001. (e) qPCR analysis of TNFSF14 expression in primary CD8+ T cells electroporated with empty vector, DTHD1-WT or DTHD-mut plasmid cocultured with THP1 cells (n = 3). P values are determined by unpaired t-test. (f) Receiver operating curve for out-of-sample prediction of case–control state by a logistic regression model trained on DEGs in CD161−CD8+ TEMRA cells. DEGs include IFI27, IFI44L, RSAD2, IFI44, FAM118A, LGALS9, MX1, EPSTI1, USP18, OAS3, LAIR2, IFIT1, XAF1. Inset depicts the changes of DEGs in the public transcriptome profile and CD161−CD8+ TEMRA cells. (g) Graphic abstract showing the expansion of CD161−CD8+ TEMRA in patients with SLE and DTHD1 downregulation promotes MYD88-mediated expansion and cytotoxicity of this pathogenic CD161−CD8+ TEMRA subset in SLE.

Fig 4: DTHD1 deficiency increased the cytotoxicity of CD161−CD8+TEMRAcells. (a) Heatmap showing the correlation of co-expression modules with different cell subsets. (b) Functional annotations of DTHD1-related genes. (c) Quantification (left) and flow cytometry (right) of the apoptosis percentage of P815 cells after cocultured with primary CD8+ T cells electroporated with si-NC or si-DTHD1 siRNAs for 5 h at the ratio of 1:5. P value is determined by Unpaired t-test, P = 0.0124. (d) Flow cytometric analysis of mean fluorescence intensity (MFI) of Granzyme B of CD161−CD8+ T cells electroporated with si-NC or si-DTHD1 siRNAs after coculture with P815 cell. P value is calculated by Unpaired t-test, P = 0.0012. (e) Structural prediction of DTHD1-MYD88 interaction by ZDOCK. Grey, DTHD1; blue, MYD88; red, mutated regions in DTHD1. (f) Western blotting analysis of HA in HEK293T cells co-transfected with myc-MYD88 and HA-empty vector or HA-DTHD1 (WT) or HA-DTHD1 (Mutant) expressing plasmids after immunoprecipitated with anti-Myc beads. (g) Boxplot indicating the average expression of MYD88 in CD161−CD8+ TEMRA cells in HCs and SLE samples from our dataset. P value is from Wilcoxon rank-sum test. (h) Luciferase activity analysis of lysates of HEK293T cells co-transfected luciferase reporter plasmid for NF-κB, pRL-TK-renilla-luciferase plasmid, MYD88 plasmid and WT-DTHD1 plasmid or mutant DTHD1 plasmid (n = 6). P values are determined by Unpaired t-test. (i) Quantification (left) and flow cytometry (right) of the apoptosis percentage of P815 cells after cocultured with primary CD8+ T cells pre-treated with DMSO or TAK-242 for 5 h at the ratio of 1:5. P value is determined by Unpaired t-test, P = 0.0024. One representative experiment of three is shown (f). ∗P < 0.05, ∗∗P < 0.01, ∗∗∗P < 0.001.

Fig 5: Analysis of CD8+T cell subsets transition trajectory in SLE. (a) Pseudotemporal analysis of CD8+ T cell subsets. Trajectory of CD8+ T Cell subclusters is inferred using monocle3 and subclusters are marked by colors (left). Pseudotime-ordered variables are inferred using monocle2 (right). Dotted lines and arrows indicate inferred differentiation trajectory and direction. (b) Correlation between pseudo-time and cytotoxicity of CD8+ T cell subsets. The solid line represents loess fitting of the relationship between cytotoxicity scores and Monocle components. P values are calculated by Spearman correlation. (c) Distribution of CD8+ subpopulations during the transition, along with the pseudo-time (upper). Heatmap showing the co-expression modules with the highest average expression in each CD8+ T cell subcluster (lower). (d) Volcano plots (left) and GO/KEGG enrichment analysis (right) of differentially expressed genes between CD8_GZMB and CD8_GZMH (left). Red represents upregulation in CD8_GZMB, and blue represents upregulation in CD8_GZMH. (e) Boxplot indicating the average expression of exhaustion and activation gene signatures in CD8_GZMB and CD8_GZMH. The P values are from a Wilcoxon test. Loess, locally weighted regression.