Two subsets of regulatory CD8+ T cells with differential transcriptome revealed by single cell analysis

Abstract

Although CD8⁺ regulatory T cells (Tregs) were described in the 1970's, they remain poorly defined compared to CD4⁺ Tregs. Their phenotypic heterogeneity and lack of consensus markers have hindered mechanistic studies and slowed clinical development despite their therapeutic potential. In this study, we performed single-cell RNA sequencing coupled with CITE-seq and TCR-seq on peripheral blood CD8⁺ T cells from four healthy donors, including the CD45RC marker to distinguish pro-inflammatory from pro-regulatory subsets. We analyzed ∼7000 freshly isolated, non-stimulated CD8⁺ T cells and identified two distincts CD8⁺ Tregs subsets, defined by HELIOS or TNFR2 expression, with unique transcriptional and surface marker profiles. Functional assays revealed potent suppressive capacity of the TNFR2⁺CD29^lowCD45RC^low/- subset. These findings were independently validated using a publicly available single-cell dataset from four additional individuals. This work provides the most comprehensive profiling of human CD8⁺ Tregs to date and supports their translation to clinical application.

Keywords: Immunology; Transcriptomics.

Graphical abstract

Figure 1

Single cell RNA-seq coupled to CITE-seq reveals CD8⁺ T cell subpopulations (A) Schematic representation of 5′ single-cell RNA-sequencing workflow (illustrated with Biorender). (B) Identification of CD8⁺ T cell subpopulations by clustering based on gene expression using uniform manifold approximation and projection (UMAP) (resolution = 1). Ten distinct clusters were defined identifying naive (N), central memory (CM), effector memory (EM) types 1 (EM1), 2 (EM2), 3 (EM3), and Pe2 (TEMRA) CD8⁺ T cells. Each dot corresponds to a single cell. (C) Heatmap showing the expression of genes defining cell cluster’s identities. Columns represents individual cells grouped by clusters and rows correspond to specific genes. Expression levels were scaled per gene with blue color indicating lower expression and red indicating higher expression. (D) Proportion of cells per cluster among CD8⁺ T cells. Data represent concatenated single-cell data from all four donors. (E) Distribution of CD45RA protein marker across clusters using UMAP. (F) Frequency of gene expression across CD8⁺ T cells clusters. Each dot corresponds to an individual donor. Each color corresponds to a specific cluster. Mean +/− standard deviation (SD).

Figure 2

Characterization of CD8⁺ regulatory T cell clusters (A) Distribution of CD45RC protein marker across clusters using UMAP. (B) Volcano plot (min.pct = 0.25, logfc.threshold = 0.2) comparing gene expression between CD45RC^high (cluster 0, 5 and 7) vs. CD45RC^low/- (clusters 1–4, 6, 8, 9) CD8⁺ T cells. Thresholds were set to p value adjusted <0.05 and fold change >1.5. Red dots indicate genes meeting both criteria and blue dot indicates genes with significant adjusted p values. Some differentially expressed genes are labeled. (C) Violin plots displaying the expression of regulatory T cell-associated genes across the clusters. Violin plots depict gene expression distributions across the clusters, with each color corresponding to a specific cell cluster and each dot representing a single cell. The central line indicates the median, the box the interquartile range (IQR), and the whiskers extend to the minimum and maximum values. (D) Clonotype distribution across clusters with red indicating clonotypes found in 10%–100% of cells within a cluster, light blue in 1%–10% and dark blue in 0%–1%. (E) UMAP showing the affiliation of the CD8⁺ T cells to individual healthy volunteers. Each dot represents a single cell, and colors indicate the corresponding healthy volunteer, highlighting the distribution and affiliation of cells across the dataset. (F) Feature plot showing KIR2DL3 gene expression in CD8⁺ T cells. Gene expression is represented on a color scale from gray (low expression) to blue (high expression) with each dot corresponding to a single cell.

Figure 3

Identification of markers specific to HELIOS⁺CD8⁺ T regulatory cells (A) Volcano plot showing differentially expressed genes between cell cluster 8 vs. all other clusters. The thresholds have been set at adjusted p value <0.05 and fold change >1.5. Red dots indicate genes that exceeded both thresholds, while blue dots correspond to genes with a significant adjusted p value but not fold change. (B) Pathway analysis of differentially expressed genes in HELIOS⁺CD8⁺ Tregs compared to other clusters. Pathways are ranked by significance using an adjusted p value threshold of <0.05. Bars represent the degree of enrichment, with the x axis showing the fold enrichment and the y axis listing the biological processes or pathways. The key genes involved in each pathway are labeled on the right of the bar. (C) Violin plots displaying the expression of in HELIOS⁺CD8⁺ Tregs-associated genes across the clusters. Violin plots depict gene expression distributions across the clusters, with each color corresponding to a specific cell cluster and each dot representing a single cell. The central line indicates the median, the box the IQR, and the whiskers extend to the minimum and maximum values.

Figure 4

Identification of markers specific to TNFR2⁺CD8⁺ T regulatory cells (A) Volcano plot showing differentially expressed genes between cell cluster 2 vs. all other clusters. The thresholds have been set at adjusted p value <0.05 and fold change >1.5. Red dots indicate genes that exceeded both thresholds, while blue dots correspond to genes with a significant adjusted p value but not fold change. (B) Pathway analysis of differentially expressed genes in TNFR2⁺CD8⁺ Tregs compared to other clusters. Pathways are ranked by significance using an adjusted p value threshold of <0.05. Bars represent the degree of enrichment, with the x axis showing the fold enrichment and the y axis listing the biological processes or pathways. The key genes involved in each pathway are labeled on the right of the bar. (C) Violin plots displaying the expression of in TNFR2⁺CD8⁺ Tregs-associated genes across the clusters. Violin plots depict gene expression distributions across the clusters, with each color corresponding to a specific cell cluster and each dot representing a single cell. The central line indicates the median, the box the IQR, and the whiskers extend to the minimum and maximum values.

Figure 5

TCR diversity of CD8⁺ regulatory T cell clusters (A) Shannon diversity index (left) and Inverse Simpson index (right) scores for each CD8⁺ T cell cluster, representing TCR clonotype diversity. Each dot corresponds to a single cluster, and error bars represent the mean ± SEM for all groups. The Shannon index reflects the evenness and richness of TCR clonotypes within each cluster, with higher values indicating greater diversity. The inverse Simpson index measures the probability that two randomly selected clonotypes belong to different clonotypes, with higher values indicating less dominance by a single or few clonotypes. (B) Chord plot showing the extent of TCR clonotypes sharing between cell clusters for all healthy volunteers. Each colored segment represents a distinct cell cluster, and the connecting bands indicate shared TCR clonotypes between clusters. The thickness of the bands corresponds to the number of shared clonotypes, with thicker bands representing higher levels of clonotype overlap. (C) Relative frequency of variable domain of the alpha (TRAV, upper) and beta (TRBV, lower) TCR chains in cluster 2 and 8. Results are expressed as mean ± SEM (n = 4). (D) Kmer analysis of the CDR3 region identifies enriched motifs within TCR sequences across cluster 2 and 8. The heatmap displays the frequency of specific amino acid kmer motifs in cluster 2 and 8. The relative frequency of each kmer is represented by color intensity, with lighter shades indicating higher enrichment.

Figure 6

TNFR2⁺CD29^lowCD8⁺CD45RC^low/- Tregs exhibit enhanced suppressive function highlighting TNFR2’s role (A) UMAP projection of CD8⁺CD45RC^low/- T cells showing gene and protein expression of TNFRSF1B (CD120b) and ITGB1 (CD29). Gene expression is scaled from gray to blue for genes, while protein expression is scaled from gray to green. Each dot corresponds to a single cell. (B) Left: Flow cytometry gating strategy to identify CD45RC^low/- and CD45RC^high CD8⁺ T cells and their subsets based on TNFR2 and CD29 surface expression. Right: Frequency of the four subsets (TNFR2⁺CD29^low, TNFR2⁺CD29^high, TNFR2⁻CD29^low, TNFR2⁺CD29^high) in CD45RC^low/- and CD45RC^high CD8⁺ T cells. Results are expressed as mean +/− SD. Mann-Whitney two-tailed test. ∗, p < 0.05. (C) Normalized MFI of TNFR2 and CD29 in FOXP3⁺ and FOXP3^- CD45RC^low/−CD8⁺ T cell subsets, relative to isotype control. Results are expressed as mean +/− SEM. Mann-Whitney two-tailed test. ∗, p < 0.05. (D) Left: suppressive assay workflow. Responder T cells (CD4⁺CD25⁻ T cells) were labeled with CFSE and co-cultured in presence of irradiated allogenic Antigen Presenting Cells (APCs) with or without subsets of CD8⁺ Tregs for 5 days at 37°C with 5% of CO₂. Proliferation was measured on CFSE dilution. Right: bar plot showing the percentage of proliferation of responder T cells with different subsets of CD8⁺ Tregs (ratio Teff:Treg 1:1) normalized to proliferation in absence of Tregs. Results are expressed as mean ± SEM (n = 6). Mann-Whitney two-tailed test. ∗, p < 0.05; ∗∗, p < 0.01. (E) Left: percentage of proliferation of responder T cells co-cultured with CD8⁺CD45RC^low/-TNFR2⁺ T cells in the presence of saturating 10 μg/mL blocking anti-TNFR2 or isotype control Ab (Left) or in the presence of TNC-scTNF(143N/145R) STAR2 (Right), in a range of Teff:Treg ratios, normalized to proliferation without Tregs. Two-way row matched Anova test. ∗p < 0.05. Results are shown as mean ± SEM (n = 4–5).

Figure 7

Identification of transcriptionally distinct CD8⁺ Treg subsets in an independent scRNA-seq dataset (A) Left, UMAP plot showing clustering of CD8⁺ T cells from a public dataset (10× Genomics) comprising ∼60,000 cells per individual and 4 individuals. Nine transcriptionally distinct clusters were identified. Right: proportion of cells per cluster among CD8⁺ T cells. Data represent concatenated single-cell data from all four donors. (B) Frequency of cells positive for the expression of canonical genes used to define naive (SELL, CCR7, CD27, LEF1), effector (FGFBP2, GZMB, GNLY, KLRB1), and memory (CCR7, CD27, FAS) T cell subsets. Each dot represents an individual donor. Results are expressed as mean +/- SD. (C) Violin plots displaying gene expression levels of selected markers across clusters. Colors correspond to cell clusters, and each dot represents a single cell. The central line indicates the median, the box the IQR, and the whiskers extend to the minimum and maximum values. (D) Volcano plots showing the transcriptional profiles of cluster 4 (TNFR2⁺) and cluster 5 (HELIOS⁺).

Figure 8

TCR repertoire analysis of CD8⁺ Treg subsets in an independent public dataset (A) Shannon diversity index (left) and inverse Simpson index (right) scores for each CD8⁺ T cell cluster from the public dataset, representing TCR clonotype diversity. Each dot corresponds to a single cluster, and error bars represent the mean ± SEM for all groups. (B) Relative frequency of variable domain usage for the alpha (TRAV, upper) and beta (TRBV, lower) TCR chains in TNFR2⁺ (cluster 4) and HELIOS⁺ (cluster 5) CD8⁺ Treg subsets. Results are expressed as mean ± SEM. (C) Heatmap displaying the frequency of specific amino acid k-mer motifs in CDR3 region in TNFR2⁺ (cluster 4) and HELIOS⁺ (cluster 5) CD8⁺ Treg subsets. The relative frequency is represented by color intensity.