Shared and distinct biological circuits in effector, memory and exhausted CD8+ T cells revealed by temporal single-cell transcriptomics and epigenetics

Abstract

Naïve CD8⁺ T cells can differentiate into effector (T_eff), memory (T_mem) or exhausted (T_ex) T cells. These developmental pathways are associated with distinct transcriptional and epigenetic changes that endow cells with different functional capacities and therefore therapeutic potential. The molecular circuitry underlying these developmental trajectories and the extent of heterogeneity within T_eff, T_mem and T_ex populations remain poorly understood. Here, we used the lymphocytic choriomeningitis virus model of acute-resolving and chronic infection to address these gaps by applying longitudinal single-cell RNA-sequencing (scRNA-seq) and single-cell assay for transposase-accessible chromatin sequencing (scATAC-seq) analyses. These analyses uncovered new subsets, including a subpopulation of T_ex cells expressing natural killer cell-associated genes that is dependent on the transcription factor Zeb2, as well as multiple distinct TCF-1⁺ stem/progenitor-like subsets in acute and chronic infection. These data also revealed insights into the reshaping of T_ex subsets following programmed death 1 (PD-1) pathway blockade and identified a key role for the cell stress regulator, Btg1, in establishing the T_ex population. Finally, these results highlighted how the same biological circuits such as cytotoxicity or stem/progenitor pathways can be used by CD8⁺ T cell subsets with highly divergent underlying chromatin landscapes generated during different infections.

Extended Data Fig. 1 |. Flow cytometry gating schemes.

a) Sort strategy of scRNA-seq/scATAC-seq depicted in Fig. 1a,b. b) Sort strategy of scATAC-seq depicted in Fig. 2i. c) Gating strategy for Fig. 3j. d) Gating strategy for Extended Data Fig. 5b. e) Gating strategy for Fig. 3m. f) Gating strategy for Fig. 5b-f. g) Gating strategy for Fig 8 g-m.

Extended Data Fig. 2 |. UMAP analysis of scRNA-seq and scATAC-seq by infection and timepoint.

UMAP from (a) scRNA-seq and (b) scATAC-seq colored by infection and timepoint as indicated.

Extended Data Fig. 3 |. Effector and memory clusters defined by scRNA-seq and scATAC-seq identify shared and non-overlapping cell subsets.

Percentage of cells from Arm infection by timepoint as indicated in (a) scRNA-seq clusters and (b) scATAC-seq clusters. c) scATAC-seq UMAP (left) and scRNA-seq UMAP (right) colored with d30 Arm cells.

Extended Data Fig. 4 |. ZEB1 motif is enriched in non-CTL clusters.

scATAC-seq UMAP of cells from Arm infection colored by ZEB1 motif enrichment. The location of CTL and non-CTL clusters is indicated.

Extended Data Fig. 5 |. CD8⁺ TIL from human melanoma post-PD1 blockade express NK receptors.

a) Sample schematic. b) Representative flow cytometry plots of four patients. Cells are first gated as live single non-naïve (not CD45RA⁺CD27⁺) CD8⁺ T cells. (Extended Data Fig.1d) c) Enumeration of subsets gated in (b). Two-sided paired Student’s t-test. n = 11 patients.

Extended Data Fig. 6 |. scATAC-seq defined clusters Eff-like I and Eff-like II are distinguished by DACRs at gene loci related to migration.

a) Barplot representing the number of DACRs between scATAC-seq clusters Eff-like I and Eff-like II. b) Number of Eff-like II DACRs per gene loci. Genes of interest annotated.

Extended Data Fig. 7 |. Zeb1 is critical for persistence of exhausted CD8⁺ T cells.

a) Experimental schematic for testing the role of Zeb1 in Cl13 infection. b) Frequency of Zeb1 KD versus control (Ctrl) over time in the spleen in Cl13 infection. Data are presented as mean values +/− standard deviation. Enumeration of T_ex subsets gated as in Fig. 3j as percent of parent (c) and total number (d). (b-d) P values calculated with two-sided paired Student’s t-test with Benjamini– Hochberg correction. n = 5 d8 Cl13, 5 d15 Cl13, 5 d30 Cl13, 5 d8 Arm, 5 d15 Arm, 5 d30 Arm mice. Data representative of 2 independent experiments.

Extended Data Fig. 8 |. Identification of Tcf7-expressing progenitor/stem-like CD8⁺ T cell subsets.

a) Gene expression from scRNA-seq of all scRNA-seq defined clusters. b) Motif enrichment from scATAC-seq of all scATAC-seq defined clusters.

Extended Data Fig. 9 |. Btg1 expression is associated with return to quiescence after proliferation.

a) Gene expression of Btg1 compared to cell cycle phase scores in Cl13. b) Correlation of Btg1 with all other expressed genes in within G2M+ cells as indicated. c) Gene ontology of genes positively or negatively correlated Btg1 performed with performed with metascape.org which uses hypergeometric test and Benjamini-Hochberg p-value correction algorithm.

Extended Data Fig. 10 |. Retroviral-mediated knock down of Btg1.

a) Experimental schematic. b) qPCR results of shRNA-mediated knockdown of Btg1. Bar represents mean, points represent independent experiments.

Fig. 1 |. Single-cell transcriptional and accessible chromatin landscape of memory and exhausted CD8⁺ T development.

a, Experimental strategy to capture CD8⁺ T cell differentiation in acute resolving and chronic viral infections. Microfluidic image provided by 10x Genomics. b, Detailed experimental schematic (Extended Data Fig. 1a). c,d, UMAP from scRNA-seq and scATAC-seq colored by infection and time point (c) or by cluster (d). e,f, Enumeration and proportion of cells per cluster as indicated for scRNA-seq (e) or scATAC-seq (f). g, scATAC-seq coverage and tile plots. Sample-specific ACRs are indicated with black boxes below tile plot. Previously identified Pdcd1 enhancer indicated in red. h, Gene expression from scRNA-seq of genes represented in g.

Fig. 2 |. Acute-resolving infection generates two branches of effector and memory CD8⁺ T cells distinguished by epigenetic cytolytic potential.

a, scRNA-seq UMAP; cells from Arm infection are colored by cluster or time point (inset). b, Expression of T cell genes by cluster. c, Number (top) and percentage (bottom) of cells from Arm infection per cluster filled by time point. d, scATAC-seq UMAP; cells from Arm infection are colored by cluster or time point (inset). e, Average gene activity per scATAC-seq cluster. f, scATAC-seq coverage and tile plots. DACRs of CTL versus non-CTL clusters indicated on the bottom. g, Average TF motif enrichment per scATAC-seq cluster of differentially enriched motifs comparing CTL and non-CTL scATAC-seq clusters. h, Differential gene activity comparing CTL and non-CTL scATAC-seq clusters. Gene loci of interest indicated. Calculation performed with two-sided Seurat FindMarkers LR test using Bonferroni correction. i, Experimental schematic of long-term Arm infection experiment (Extended Data Fig. 1b). j, scATAC-seq UMAP of cells from experiment in i colored by time point (top) or cluster (bottom). k, Enrichment score of cluster-specific ACRs from d30 Arm Mem-CTL and Mem scATAC-seq clusters. Two-sided Wilcoxon test of EM (871 cells) or EM-CTL (547 cells) versus the rest (5,741 or 6,065 cells). l, Average gene activity per scATAC-seq cluster with number of cells per cluster indicated on top, filled by time point. m, Enrichment score of gene activity from gene sets derived from TRM or circulating memory cells. Two-sided Wilcoxon test of TRM (419 cells) versus the rest (6,193 cells). n, Number of DACRs between CM d60 and CM d200 clusters. DACRs were calculated with Signac FindAllMarkers two-sided likelihood-ratio (LR) test using Bonferroni correction. o, scATAC-seq UMAP of cells from experiment in i colored by ZEB1 motif enrichment. p, Data summary schematic. In box plots, the median is indicated by the center line; box limits represent upper and lower quartiles; and whiskers extend to 1.5 times the interquartile range.

Fig. 3 |. Exhausted CD8⁺ T cells are transcriptionally heterogeneous and include a distinct subset characterized by expression of natural killer cell receptors.

a, scRNA-seq UMAP; cells from Cl13 infection are colored by cluster or time point (inset). b, Average gene expression per scRNA-seq cluster with proportion of cells per time point in each cluster represented below. c, Phylogenetic tree of scRNA-seq clusters with proportion of cells per time point. Correspondence of clusters with previous nomenclature: α¹¹, β¹⁰, χ¹². d, DEG analysis between Eff and Eff-like clusters. e, Gene Ontology analysis of DEGs in d performed with Metascape, which uses a hypergeometric test and Benjamini–Hochberg P-value correction algorithm. f, Cell cycle S.Score for each cluster. The number of cells in each cluster is available in Supplementary Table 7. g, Predicted cluster identity of proliferating cells shown as the number of cells per cluster and colored by time point (Methods). h, DEG analysis between Exh-Int and Exh-KLR clusters. i, DEG analysis between Exh-Term and Exh-KLR clusters. j, Flow cytometry gating strategy to identify T_ex clusters. Cells were gated as live single CD8⁺ P14 cells (Extended Data Fig. 1c). k, Enumeration of T_ex clusters gated in j. Each point represents a mouse. l, Representative flow cytometry plots gated on Exh-KLR cells as in j from Cl13 infection at d15 and d30. Mean percentage per quadrant is indicated. m, Representative flow cytometry plots from Arm infection at d15 or d30 gated on live singlet CD8⁺ P14 cells (top) or KLRC1⁺KLRD1⁺ P14 cells (bottom) as indicated. Mean percentage per quadrant is indicated. n, DEG analysis between CTL cluster from Arm infection and Exh-KLR cluster from Cl13 infection. o, DEG analysis between Mem-CTL cluster from Arm infection and Exh-KLR cluster from Cl13 infection. In d, h, I, n and o, DEGs were calculated with Seurat FindMarkers two-sided Wilcoxon test using Bonferroni correction. In j–m, n = 5 d8 Cl13, n = 5 d15 Cl13, n = 15 d30 Cl13, n = 5 d15 Arm and n = 5 d30 Arm mice. Data are representative of two independent experiments. In box plots, the median is indicated by the center line; box limits represent upper and lower quartiles; and whiskers extend to 1.5 times the interquartile range.

Fig. 4 |. The accessible chromatin landscape distinguishes fewer exhausted T cell epigenetic cell fates under wider transcriptional diversity.

a, scATAC-seq UMAP; cells from Cl13 infection are colored by cluster or time point (inset). b, Average enrichment score per scATAC-seq cluster of gene sets from scRNA-seq cluster DEGs, using gene activity. c, Phylogenetic tree of scATAC-seq clusters with proportion of cells per time point. d, Percentage of cells from Cl13 infection by time point as indicated in scRNA-seq clusters (top) and scATAC-seq clusters (bottom). e, Average accessibility of DACRs per scATAC-seq cluster with proportion of cells per time point in each cluster represented below. f, Number of DACRs per gene loci for each scATAC-seq cluster. DACRs were calculated with Signac FindAllMarkers two-sided LR test using Bonferroni correction. g, Average TF motif enrichment per cluster for differentially enriched TF motifs. h, Top, ZEB1 motif enrichment. Bottom, Zeb1, Zeb2 and Tox average gene activity per scATAC-seq cluster.

Fig. 5 |. Zeb2 promotes differentiation of epigenetically distinct cytotoxic CD8⁺ T cell subsets in chronic and acute-resolving viral infection.

a, Experimental schematic for testing the role of Zeb2 in Cl13 versus Arm infection. b, Frequency of Zeb2 KD versus control (Ctrl) over time in the spleen in Cl13 infection. Data are presented as mean values ± s.d. c, Enumeration of subsets from Cl13 infection as gated in Fig. 3j. d, Representative flow cytometry plots from d30 Cl13 infection as indicated. e, Frequency of Zeb2 KD versus Ctrl over time in the spleen in Arm infection. Data are presented as mean values ± s.d. f, Enumeration of subsets from Arm infection. g, Representative flow cytometry plots from d30 Arm as indicated. In b–g, n = 5 d8 Cl13, n = 7 d15 Cl13, n = 8 d30 Cl13, n = 5 d8 Arm, n = 5 d15 Arm and n = 5 d30 Arm mice. Data are representative of three independent experiments. Cells are gated as live single CD8⁺ P14 cells KD or Ctrl (Extended Data Fig. 1f). In b, d, e and g, P values were calculated using two-sided paired Student’s t-test with Benjamini–Hochberg correction. h, Venn diagram of overlapping DACRs in scATAC-seq CTL, Mem-CTL and Exh-KLR clusters. i, TF motif enrichment in DACRs comparing scATAC-seq Exh-KLR and CTL (top) or Mem-CTL (bottom) clusters with total number of DACRs represented as bar plot below. TF motif enrichment was calculated with Signac FindMotifs, which uses a hypergeometric test and Benjamini–Hochberg correction. RNP, ribonucleoprotein; sgRNA, single-guide RNA.

Fig. 6 |. PD-1 pathway blockade alters exhausted T cell subset dynamics within the preexisting population structure.

a, Experimental schematic. b, Frequency of blood P14 cells determined by flow cytometry. Data are presented as mean values ± s.d. Each dot is a mouse. d30 with and without αPD-L1 were compared with a two-sided Student’s t-test. n = 10 d8, n = 20 d15, n = 20 d30, n = 15 d30 + αPD-L1 mice. Data are from one scATAC-seq experiment. c, scATAC-seq UMAPs colored by infection, time point and treatment as indicated. d, Difference in the number of cells in each scATAC-seq T_ex cluster with and without αPD-L1 using Fisher’s exact test. e, scATAC-seq UMAP of Cl13 d30 cells with or without αPD-L1, colored by cluster (left) or pseudotime calculated by Monocle (right). f, scATAC-seq UMAP of Cl13 d30 cells with and without αPD-L1 colored by density. g, Pseudotime within each scATAC-seq cluster comparing cells with or without αPD-L1 using two-sided Wilcoxon test. The number of cells in each cluster is available in Supplementary Table 7. h, The schematic shows a summary of results. In box plots, the median is indicated by the center line; box limits represent upper and lower quartiles; and whiskers extend to 1.5 times the interquartile range.

Fig. 7 |. Acute-resolving and chronic infections generate Tcf7-expressing progenitors with divergent accessible chromatin profiles.

a,d, Phylogenetic tree of scRNA-seq (a) or scATAC-seq (d) clusters; bar plots represent the proportion of cells per infection and time point. b,e, UMAP from scRNA-seq (b) or ATAC-seq (e) of Naïve, MP, Mem, Exh-Pre and Exh-Prog subsets, colored by subset. c, Average gene expression per cluster of all pairwise DEGs. f, Average accessibility per cluster of all pairwise DACRs. g, Venn diagram of MP, Exh-Pre and Exh-Prog DEGs. h, Gene Ontology of DEGs in g performed with Metascape, which uses a hypergeometric test and Benjamini–Hochberg P-value correction algorithm. i, DEG analysis between MP and Exh-Pre clusters. Calculation performed with Seurat FindMarkers two-sided Wilcoxon test using Bonferroni correction. j, Venn diagram of DACRs calculated between scATAC-seq naïve and MP or Naïve and Exh-Pre clusters. k,m, Coverage and tile plots from scATAC-seq of the Satb1 locus (k) or Lef1 locus (m); DACRs closing in both MP and Exh-Pre (green) or only in Exh-Pre (blue), compared to Naïve, are indicated on the bottom. l,n, Gene expression of Satb1 (l) or Lef1 (n). o, TF motif enrichment of DACRs comparing MP and Exh-Pre; the total number of DACRs are represented as bar plot below. Enrichment was calculated with Signac FindMotifs, which uses a hypergeometric test and Benjamini–Hochberg correction.

Fig. 8 |. Transition from Exh-Pre to Exh-Prog uncovers Btg1 as a new regulator of exhausted T cell differentiation.

a, Enumerated DEGs from pairwise analysis. DEGs were calculated with Seurat FindMarkers two-sided Wilcoxon test using Bonferroni correction. b, Gene Ontology of DEGs from a performed with Metascape, which uses a hypergeometric test and Benjamini–Hochberg P-value correction algorithm. c, Representative flow cytometry plot (top) and enumeration of MFI of HPG signal (bottom) from in vitro translation assay. Cells were gated as described in Fig. 3j. n = 5 d8 Cl13, n = 5 d15 Cl13, n = 6 d30 Cl13, n = 5 d8 Arm, n = 4 d30 Arm and n = 6 naive mice. Data are representative of two independent experiments. P values were calculated with two-sided Student’s t-test using Benjamini–Hochberg correction. d, Enumerated DACRs from pairwise analysis as indicated. DACRs were calculated with Signac FindAllMarkers two-sided LR test using Bonferroni correction. e, Number of DACRs in d15 Exh-Prog per gene loci. Select genes overlapping with d15 Exh-Prog DEGs are annotated. f, Experimental schematic. g, Flow cytometry plot of input P14 cell mixture containing Btg1 KD and Ctrl. h, Representative flow cytometry plot from d8 p.i. Cells were gated on RV⁺ (GFP) live CD8⁺ P14 T cells (Extended Data Fig. 1g). Mean percentage is indicated. i, Total P14 cells with Btg1 KD or Ctrl RV. P values were calculated with two-sided paired Student’s t-test. j, Representative flow cytometry plot of Ki67 staining; histograms were colored by shRNA target. k, Total Ki67⁺ cells per shRNA target as indicated. P values were calculated with two-sided paired Student’s t-test. l, Representative flow cytometry plots of T_ex subsets. Cells gated on RV⁺ live P14 T cells. Mean percentage as indicated. m, Total RV⁺ live P14 T cells per each subset. Mean fold change is indicated. P values were calculated with two-sided paired Student’s t-test with Benjamini–Hochberg correction. In h–m, n = 6 mice. Each point represents a mouse. Data are representative of three independent experiments.