Replicative history marks transcriptional and functional disparity in the CD8+ T cell memory pool.

Clonal expansion is a core aspect of T cell immunity. However, little is known with respect to the relationship between replicative history and the formation of distinct CD8+ memory T cell subgroups. To address this issue, we developed a genetic-tracing approach, termed the DivisionRecorder, that reports the extent of past proliferation of cell pools in vivo. Using this system to genetically 'record' the replicative history of different CD8+ T cell populations throughout a pathogen-specific immune response, we demonstrate that the central memory T (TCM) cell pool is marked by a higher number of prior divisions than the effector memory T cell pool, owing to the combination of strong proliferative activity during the acute immune response and selective proliferative activity after pathogen clearance. Furthermore, by combining DivisionRecorder analysis with single-cell transcriptomics and functional experiments, we show that replicative history identifies distinct cell pools within the TCM compartment. Specifically, we demonstrate that lowly divided TCM cells display enriched expression of stem-cell-associated genes, exist in a relatively quiescent state, and are superior in eliciting a proliferative recall response upon activation. These data provide the first evidence that a stem-cell-like memory T cell pool that reconstitutes the CD8+ T cell effector pool upon reinfection is marked by prior quiescence.

Introduction

The CD8+ T cell compartment serves to provide protection against intracellular pathogens and also acts as a modifier of cancer growth. Upon antigen encounter, naïve T cells (TN) undergo extensive gene-expression alterations, while entering a highly proliferative state, dividing every 4h to 6h[1,2] in mice. This phase of clonal expansion gives rise to a phenotypically and functionally diverse pool of effector T cells (TEFF) that exceeds its precursor population size by >10,000-fold[3,4]. Unlike TN, these TEFF have the capacity to disseminate to peripheral tissues, and scan for and kill infected or transformed cells. Upon antigen clearance, around 95% of the TEFF pool succumbs to apoptosis, leaving behind a small long-lived pool of memory T cells (TM) that is equipped to provide long-term protection against recurring pathogens.

The central role of proliferation in the T cell response has inspired many to study the relationship between replication and T cell state. While earlier work hinted that memory precursor T cells have undergone limited clonal expansion[5,6], more recent work studying acute T cell responses in human subjects demonstrated that TM, as a whole, are derived from precursor cells that have undergone an extensive number of divisions[7]. Furthermore, prior work has shown that cell cycle speed can differ substantially between phenotypically distinct T cell subsets at different time-points in the T cell response. Specifically, central memory T cells (TCM), a subgroup of memory cells that are endowed with a high level of multipotency, have been documented to undergo homeostatic proliferation after pathogen clearance, while effector memory T cells (TEM) have a low turnover rate[8,9]. In contrast, during the effector phase, a TCM-like state has been linked to lower division speed and reduced clonal burst size compared to their TEM-like and terminally differentiated counterparts[10-13].

The phase-dependent association of proliferative activity within specific cell states, in combination with the reported phenotypic instability of certain T cell subsets[14,15], makes it difficult to deduce the replicative history (i.e., the cumulative number of prior divisions) of different memory T cell populations, and the possible relationship between such replicative history and functional properties. Here, we develop a genetic-tracing approach—termed DivisionRecorder—that allows the measurement of prior division of cell pools over extensive rounds of division, and apply this approach to determine to what extent replicative history identifies distinct memory T cell states and behaviors. In this effort, we focus on three central issues: (1) What are the differences in replicative history between (precursor-)TCM and TEM in the effector and memory phase? (2) Is there heterogeneity in prior division within the TCM pool? (3) If so, does replicative history of cells within the TCM pool predict their capacity to mount a secondary T cell response?

Results

Division-linked genetic labeling of cell pools

The genome contains a large number of hypervariable short tandem nucleotide repeats (STRs) that accumulate intra-allelic length mutations through DNA polymerase slippage during cell division. Such slippage mutations in endogenous STRs have been used to study lineage trees in various organisms and tissues[16,17], and synthetic STRs have previously been employed in a probabilistic labeling approach to define stem cells in the intestinal epithelium and the mammary gland[18,19]. To investigate the replicative history of memory T cells, we engineered a synthetic STR-reporter system to continuously ‘record’ proliferation in cell pools. This genetically encoded system, termed DivisionRecorder, utilizes a synthetic STR domain to achieve a division-linked low-probability acquisition of a fluorescent mark (Fig. 1a). The DivisionRecorder consists of two separate elements: (1) a retroviral-vector encoded module that contains a synthetic STR linked to an out-of-frame CRE recombinase gene; (2) A CRE-activity reporter module that irreversibly induces the expression of a red fluorescent protein (RFP). In its base configuration, all cells that contain the DivisionRecorder only express GFP (hereafter referred to as DRGFP cells). As cells undergo successive divisions, slippage mutations that occur within the synthetic STR yield in-frame variants of the downstream CRE recombinase gene at a fixed, division-dependent, probability (p). The resulting CRE activity induces an irreversible activation of the RFP gene, giving rise to GFP+RFP+ cells (hereafter referred to as DRRFP) that pass this genetically encoded label on to subsequent generations, resulting in a cumulative increase in the DRRFP cell fraction within the DivisionRecorder+ (DR+, i.e., the sum of DRGFP and DRRFP) population as the cell pool expands (Fig. 1b, Supplementary Note 1). Importantly, when p is small (< 0.01) the DivisionRecorder yields a near-linear relationship between the DRRFP fraction and the average number of divisions over dozens of population doublings (Fig. 1c)[20], thereby allowing analysis of replicative history—at the population level—far beyond what can be achieved with classical cell labeling dyes[21] (Fig. 1d).

Fig. 1

DivisionRecorder activation is a proxy for replicative history.

a, Schematic overview of the DivisionRecorder system. b, Cartoon depicting progressive DivisionRecorder activation in a proliferating cell pool. c, Simulation of the minimal ODE model (See Supplementary Note 2 for detailed description and equations), depicting DRRFP acquisition as a function of population doublings for the indicated values of DRRFP acquisition probability (p). d, Maximal number of theoretically recordable population doublings, approximated by calculating the amount of division events required to reach a 99% DRRFP population. Approximate maximums for selected values of p are indicated, colors correspond to legend in panel c. e-f, Percentage of DRRFP cells over time in cultured DivisionRecorder+ (DR+) CRE-activity reporter HEK 293T cells (n=3 replicates per group) in which the CRE recombinase gene was preceded by either a stable nucleotide region (indicated as “no STR”) or a repeat of 24 guanines (indicated as “with STR”). Representative plots (e) and summarizing line graphs (f) are shown. Lines connect experimental replicates g, Percentage of DRRFP cells across population doublings in DR+ CRE-activity reporter HEK 293T cells (n=3 replicates per group) in which the CRE recombinase gene was preceded by either a low stability STR ([G]24) or a high stability STR ([CA]30). Dots indicate individual samples, lines represent fitted linear regression, dotted lines indicate bounds of the 95% confidence interval. h-i, Percentage of DRRFP cells across population doublings in immortalized DR+ mouse embryonic fibroblasts. Representative flow cytometry plots (h) and summarizing graph (i) are shown. Best fits of the minimal ODE model are depicted (100 bootstraps per experimental replicate, Supplementary Note 2). Blue line represents the median of the bootstraps, grey lines represent individual fits, dots indicate experimental measurements (n=3 replicates). p indicates the estimated DRRFP acquisition probability. Depicted experimental data are representative of at least two independent experiments. P values (g) were determined by two-sided ANCOVA test.

To test the utility of the DivisionRecorder, we established a reporter cell-line carrying a lox-STOP-lox-RFP cassette. Following retroviral introduction of the GFP-STR-CRE module, a progressive increase in DRRFP cells was observed over time, whereas no label acquisition was observed when the STR was replaced with a stable DNA sequence (Fig 1e, f). Moreover, the rate at which DRRFP cells accumulated was dependent on the sequence stability of the STR[22,23], underpinning that p is linked to the likelihood of STR slippage (Fig 1g). Similarly, upon introduction of the DivisionRecorder into immortalized embryonic fibroblasts from the Ai9 mouse strain—that carry an endogenous lox-STOP-lox-RFP cassette[24]—a low and predictable DRRFP acquisition was observed, with a [G]33 STR conferring a p of 0.0052 ±0.00074 (Fig. 1h, i), thereby enabling the measurement of replicative history over many cell divisions (in theory >1,500 population doublings, Fig. 1d).

To test whether the DivisionRecorder can be used as a proxy for replicative history in the CD8+ T cell compartment in vivo, we generated Ai9;OT-I mice, in which all T cells recognize the OVA257-264 epitope, thereby allowing examination of T cell pools in the context of equal TCR affinity. Ai9;OT-I T cells were isolated, modified with the DivisionRecorder to obtain DR+ OT-I T cells, transferred into Listeria monocytogenes-OVA (Lm-OVA) infected mice, and the fraction DRRFP cells was measured over time (Fig. 2a). At early time-points post cell transfer (d1-d4), a rapid increase in DRRFP cells was observed (Fig. 2b, c), coinciding with the proliferative burst of the antigen-specific CD8+ T cell pool. To determine whether the observed accumulation of DRRFP cells formed an accurate measure of prior cell division, DR+ OT-I T cells were stained with CellTrace Violet (CTV) prior to cell transfer. Notably, analysis of the fraction DRRFP cells within cell pools with different degrees of CTV dilution revealed a close correlation (Fig. 2d, e, rrm = 0.94), providing direct evidence that in vivo DRRFP acquisition reflects the extent of past division in the CD8+ T cell pool. In conclusion, these data establish that the DivisionRecorder allows the long-term measurement of division history in cell pools in vivo, in a way that is compatible with down-stream methodologies such as single cell sequencing (see below).

Fig. 2

The DivisionRecorder can be applied to study T cell division kinetics in vivo.

a, Overview of experimental setup. b-c, DR+ OT-I T cells were transferred into recipient mice 24 hours post infection with Lm-OVA. Spleen samples were analyzed for the percentage of DRRFP cells at day 1-4 post cell transfer. Representative pseudo-color density plots (b), and boxplots (c) in which the boxes indicate group median and 25th/75th percentiles, whiskers represent min/max, dots represent individual samples (n=8 mice for day 1 and 2; n=7 mice for day 3 and 4). d-e, CTV-stained OT-I T cells were retrovirally transduced with the DivisionRecorder and transferred into recipient mice (n=4) 24 hours post infection with Lm-OVA. 48 hours post-transfer, splenic DR+ OT-I T cells were assessed for CTV dilution (d), and the percentage of DRRFP cells within each division peak was analyzed (e). All depicted data are representative of at least two independent experiments, lines and symbols indicate individual mice or samples. P values were determined by two-sided Kruskal-Wallis test, with Dunn’s multiple comparisons test (c), or two-sided repeated measurement correlation test (e).

CD8+ TCM cells are derived from replicative mature T cells

Having validated the utility of the DivisionRecorder to record T cell division, we next sought to determine the replicative history of the total CD8+ TM pool relative to that of the TEFF pool. Analysis of the size of the DR+ OT-I T cell compartment in blood following Lm-OVA infection showed the characteristic rapid expansion phase, with T cell numbers peaking around day 6, and subsequent contraction into a stable memory pool (Fig. 3a). Notably, DRRFP cells remained detectable following formation of T cell memory, thus allowing analysis of replicative history at late time points after infection (Fig. 3b).

Fig. 3

The multipotent memory T cell pool is formed by replicative ‘mature’ cells.

a-c, Kinetics of DR+ OT-I T cells (a) and the percentage of DRRFP relative to day 4 (c) in response to Lm-OVA, measured in peripheral blood (n=6 mice). Representative flow cytometry plots (b) showing DRRFP and DRGFP frequencies at indicated time points, and line graphs (a, c) depicting kinetics of single mice (grey) and group median (black). d, DRRFP percentages within blood at day 5/6 (TEFF) and day >60 (TM) following LCMV-OVA infection (n=7). e, Representative plots depicting DRRFP frequencies in blood (Bl), spleen (Spl) and liver (Liv). f, Percentage of DRRFP detected in indicated organs of recipient mice at the indicated time points (n=6 mice per time point; response to Lm-OVA). Boxplots indicate group median and 25th/75th percentiles, whiskers represent min/max, dots represent individual samples. g, Moving average of surface marker expression level on splenic DR+ cells plotted against the percentage of DRRFP within each window during memory (day 86; n=6), means are shown in black. DRRFP percentages within each window are corrected for the total percentage of DRRFP detected in that sample. h, Gating strategy (left) and DRRFP percentages (right) of CD27HIKLRG1LO and CD27LOKLRG1HI cells in spleen during effector (d6, top) and memory phase (d86, bottom; n=6) in response to Lm-OVA. i, DRRFP percentages within the CD27HIKLRG1LO and CD27LOKLRG1HI cell populations in blood, comparing effector (day 5/6) and memory (day >60) phases. Data shown for Lm-OVA (top; n=22) and LCMV-OVA (bottom; n=7) infections. Lines connect individual mice. j-k, Ki67 expression by CD27HIKLRG1LO and CD27LOKLRG1HI OT-I cells in blood in response to Lm-OVA. Representative flow cytometry plot (j), and line graphs (k) where solid lines indicate population means, shaded areas indicate 95% confidence interval (n=11 mice). All depicted data are representative of at least two independent experiments, lines and symbols indicate individual mice or samples. P values were determined by two-sided Kruskal-Wallis test with Dunn’s multiple comparisons test (f), or two-sided Wilcoxon’s signed-rank test (d, h, i).

In case TM would primarily be derived from T cells that had undergone limited proliferation upon primary antigen encounter, the fraction of DRRFP cells would be expected to decay during the contraction phase, due to the decline in the number of clonally expanded TEFF (Extended Data Fig. 1, Supplementary Note 3). However, analysis of DRRFP frequencies in blood demonstrated that the fraction of DRRFP cells did not decline, but instead continued to increase during contraction and memory phase (an increase of 2.07% ±0.77% between day 13 and 59, Fig. 3c). This increase in DRRFP frequencies post pathogen-clearance was not restricted to T cell responses induced by Lm-OVA infection, but was also observed upon infection with LCMV-OVA[25] (Fig. 3d), and was not due to anatomical redistribution of cells with distinct division histories, as the fraction of DRRFP cells increased concurrently in peripheral blood and the primary sites of Lm-OVA infection (spleen/liver; Fig. 3e, f). Thus, in line with work by Akondy et al.[7], our results support the notion of a replicative ‘mature’, rather than ‘nascent’, CD8+ TM pool, and extends this observation beyond the peripheral blood compartment to the sites of infection.

Extended Data Fig. 1

Simulation of different scenarios of memory T cell formation.

Simulated data depicting a responding antigen-specific T cell population (blue), comprised of TEFF undergoing clonal expansion and subsequent contraction (red), plus memory precursor T cells (MP, green) that develop into TM. Activated TEFF are modeled to divide rapidly for 6 days (expansion phase), die at a fixed rate throughout the response, and can differentiate into MP cells only during the expansion phase. Cell numbers (top row) and DRRFP percentages (bottom row) are shown for 3 scenarios: (left) TEFF can give rise to MP cells during the entire expansion phase, irrespective of the number of prior divisions, (middle) only TEFF that have gone through at most 24 divisions can give rise to MP cells, or (right) only TEFF that have gone through at most 10 divisions can give rise to MP cells. Note the strong decay in DRRFP percentage that is observed during memory formation in case T cell memory is founded by T cells that have undergone few divisions. See Supplementary Note 3 for detailed description and equations.

It has been well documented that TCM are able to maintain the memory pool through infrequent homeostatic cell division[15,26,27], and recent work has shown that precursor-TCM slow down their replicative cycle early during the expansion phase[10], suggesting limited clonal expansion of these cells during the early phase of the T cell response. However, it is difficult to translate cell-cycle activity at a given time-point into cumulative proliferative history, and we therefore wished to directly test the relationship between cell state (e.g., TCM or TEM) and replicative history during different stages of the T cell response. To this end, the fraction of DRRFP cells within the TM pool was calculated at varying expression levels of proteins associated with either multipotency or terminal differentiation (Fig. 3g). This analysis revealed a positive correlation between replicative history and the expression of the TCM-associated proteins CD27 (rrm = 0.81, P = 6.2·10−14) and CD62L (rrm = 0.62, P = 5.6·10−7)[15,28,29], and a negative relationship between prior division and the expression of the TEM- associated proteins KLRG1 (rrm = −0.83, P = 9.0·10−15) and CX3CR1 (rrm = −0.75, P = 4.5·10−11)[14,15,30]. Likewise, defining multipotent TCM and terminally differentiated TEM subsets by joint expression or absence of CD62L and CD27, respectively, (Extended Data Fig. 2a), and further partitioning based on the expression of KLRG1 or CX3CR1, revealed a positive association between division history and a less differentiated cell state (Extended Data Fig. 2b). Furthermore, the division history of CD27HIKLRG1LO TCM present in lymph nodes equaled that of TCM in the spleen, implying that division history is dictated by cell state rather than anatomical location (Extended Data Fig. 2c).

Extended Data Fig. 2

Evaluation of the division history of T cell subsets throughout a response to Lm-OVA.

a, Gating strategy used to identify indicated TM populations (d86) in spleen samples. b, DRRFP percentages within splenic TM populations (n=6 mice) as identified in panel a. c, DRRFP percentages within the CD27HIKLRG1LO TCM subset in spleen and lymph nodes (LN) and within the CD27LOKLRG1HI TEM subset in spleen. d, Cell surface expression of CX3CR1, CD62L, and CD43 within splenic CD27LOKLRG1HI and CD27HIKLRG1LO populations at the peak of the TEFF phase (day 6 post infection) and in memory phase (day 86 post infection). e, Moving-average of surface marker expression of splenic DR+ OT-I T cells during effector phase (day 6), depicted as in Fig. 3g. f, Boxplots depicting DRRFP percentages within TEFF (day 6 post infection) subsets in spleen (n=6 mice), relative to the total DRRFP percentage. g, Kinetics of DRRFP percentages within CD27LOKLRG1HI (left) and CD27HIKLRG1LO (right) DR+ OT-I T cell populations in blood. Values are relative to the percentage of DRRFP cells detected at the peak of the response (day 6). Grey lines represent individual mice (n = 22), red and blue lines indicate group mean. h, Simulation of the phenotype model (See Supplementary Note 5 for details) illustrating a scenario in which conversion of CD27HIKLRG1LO to CD27LOKLRG1HI cells occur only after the peak of the response at a low rate. Depicted are the overall cell numbers (left), and the percentage DRRFP cells of DR+ OT-I T cells (right) in CD27HIKLRG1LO cells (blue), CD27LOKLRG1HI cells (red) and the total T cell population (green). Note that in this scenario the fraction DRRFP within the terminally differentiated CD27LOKLRG1HI population would increase to almost twice the experimentally observed frequency. All depicted data are representative of at least two independent experiments. Boxplots (c, d, g) represent group median and 25th/75th percentiles, whiskers indicate the interquartile range multiplied by 1.5 (c, d) or min/max (g), dots indicate individual samples. P values were determined by one-way ANOVA followed by Tukey’s HSD post-hoc test (c and d), two-sided Student’s T test (c), two-sided repeated measurement correlation test (h), or two-sided Friedman test (g). All significant (< 0.05) P values are indicated in the plots.

Next, to delineate at which point the divergence in replicative history between T cells with a TCM-like multipotent and TEM-like terminally differentiated phenotype developed, we assessed the link between phenotypic marker expression and DRRFP fractions throughout the T cell response. Notably, replicative history varied minimally across TEFF cell states at the peak of the antigen-specific T cell response (d6 post transfer, Fig. 3h, Extended Data Fig. 2d-f), followed by selective accumulation of DRRFP within the CD27HIKLRG1LO early-TCM pool directly after the peak of the expansion phase (Fig. 3h, i, Extended Data Fig. 1g), due to continued replicative activity of this subset (Fig. 3j, k). The observation that the division history of CD27LOKLRG1HI T cells stays constant post effector phase (Fig. 3i) suggests that, in addition to the previously documented lack of proliferative activity of this cell pool[15,26,31], this terminally differentiated subset also does not receive significant replenishment by the replicative active CD27HIKLRG1LO T cell pool (Extended Data Fig. 1h). The substantial number of divisions that we observe in the CD27HIKLRG1LO cell pool at the peak of the response appears at odds with the proposed limited clonal expansion of precursor-TM. However, these observations may either be reconciled by the reported trans-differentiation between TEFF cell states[14,15,30], or by the fact that a reduced proliferative activity may form a property of only a small part of the memory precursor pool[10,11,32]. In summary, the above data indicate that the high amount of prior division of the TCM pool results from both strong proliferative activity during the effector phase and selective proliferative activity after pathogen clearance.

Replicative history identifies distinct TCM cell states

Increasing evidence suggests that the TCM pool is highly heterogeneous in terms of both gene expression profiles and prior and ongoing replicative behavior[14,15,33], providing an incentive to test for possible associations between division history and transcriptional states within this cell pool. To this end, we carried out single-cell mRNA sequencing (scRNAseq) on DRGFP and DRRFP memory OT-I T cells (75-85 days post Lm-OVA infection; Extended Data Fig. 3). In addition, to test whether DR+ OT-I TM assume the same spectrum of transcriptional states as non-modified T cells, we performed scRNAseq on OT-I TM that were generated through adoptive transfer of a small number (2,000) of naïve OT-I T cells followed by Lm-OVA infection 24 hours later. DR+ OT-I and unmodified OT-I memory T cells were jointly grouped into 23 transcriptionally distinct MetaCells (MCs)[34] that included 4 TEM and 19 TCM MCs based on the expression of a small set of multipotency- and effector-associated genes (Fig. 4a,b). Notably, while memory T cells derived from small numbers of unmodified OT-I T cells showed a proportionally greater contribution to TEM MCs—consistent with the relationship between precursor frequency and TEM formation[35]—DR+ OT-I T cells and unmodified OT-I T cells were equal in their potential to yield the 19 distinct TCM MCs (Extended Data Fig. 4), indicating that the introduction of the DivisionRecorder did not measurably impact the ability of T cells to differentiate into different TCM states.

Extended Data Fig. 3

Single cell mRNA sequencing of DivisionRecorder+ and unmodified memory T cells.

Single cell mRNA sequencing was performed on DivisionRecorder modified and unmodified OT-I memory T cells (Day 75 and 85 post Lm-OVA infection), isolated from spleens (n=7 mice with DR+ memory T cells; n=4 with unmodified memory T cells). Obtained data were aggregated from two independent experiments (Experiment 1: M1-3; Experiment 2: M4-11). All cells were jointly analysed and clustered. a, Cell count per sample. b, Total cell count per MC. c, Sample composition of each MC. d, Relative contribution of DRGFP and DRRFP to the total DR+ pool within each MC.

Fig. 4

Replicative history identifies distinct transcriptional states within the TCM pool.

Single cell transcriptomic profiling of DR+ T cells obtained from spleen in memory phase (Day 75 and 85 post Lm-OVA infection). a, 2D projection of all profiled cells, colors indicate MCs (left), or relative expression of effector- and multipotency-associated genes (right). Gene list in Supplementary Table 1. b, Hierarchical clustering of MCs by their expression of effector- and multipotency-associated genes used in A. MCs are divided into 3 clusters based on Euclidean distance. c, Expression of selected genes by each MC cluster. d, DRRFP/DRGFP ratio within each MC, depicted as waterfall plot (left) and boxplot (right). e-f, Enrichment of gene signatures from MSigDB (C7, collections deposited by Goldrath (GR) and Kaech (KA), Supplementary Table 2) by gene set enrichment analysis comparing ldTCM and hdTCM (e), and enrichment plots (f) of 2 representative gene sets. g, Heatmaps depicting genes involved in immune function that were significantly (P < 0.05) depleted (left) or enriched (right) within ldTCM (See Extended Data Figure 4d, Supplementary Table 3). Selected genes are annotated, complete gene lists in Supplementary Table 4. h, QstemScore of all TCM MCs depicted as waterfall plot (left) and boxplot (right). QstemScore is based on marker genes of quiescent stem cells (Supplementary Table 5)[30], see methods for calculation. Data depicted were accumulated in two independent experiments (3-4 mice per experiment). Boxplots (c, d, h) indicate group median and 25th/75th percentiles, whiskers indicate the interquartile range multiplied by 1.5, dots signify individual MCs. The phenotype clusters TEM, TCM(eff) and TCM(mult) contain 4, 9 and 10 MCs, respectively. P values were determined by two-sided Tukey’s HSD test (c), two-sided Student’s T test with false-discovery rate correction (d, h), the FGSEA algorithm followed by the Benjamini-Hochberg procedure (e), or two-sided Wilcoxon Rank Sum test with Bonferroni correction (g). Significant P values (< 0.05) are indicated in the plots.

Extended Data Fig. 4

TCM transcriptional states are preserved in DR+ OT-I T cells.

Comparison of transcriptional states of splenic memory T cells generated by either DivisionRecorder modified, or unmodified OT-I T cells (Day 75 and 85 post Lm-OVA infection). a-b, Memory OT-I T cells cluster into TCM (blue) and TEM (red). 2D projection colored by subset (a), and violin plots depicting normalized UMI counts of selected genes (b) are shown. c, 2D projection of either DR+ (left) or unmodified (right) memory OT-I T cells. d, Contribution of DR+ and unmodified memory T cells to the TCM and TEM subsets. e, Contribution of DR+ and unmodified OT-I T cells to the 19 MCs that jointly make up the TCM subset. Dots indicate individual mice (n=3 per condition). Note that all TCM states are generated in near-equal proportions by DR+ and unmodified memory T cells. Depicted scRNAseq data was obtained from 6 individual mice, and was aggregated from 2 independent experiments. P values were determined by two-sided Student’s T test followed by Bonferroni correction for multiple testing (d and e). P values < 0.05 are indicated.

Amongst the observed TCM MCs, two transcriptionally distinct subgroups could be identified (Fig. 4b). Specifically, while all TCM showed the expected high expression of Bcl2, Sell and Cd27, and minimal expression of Cx3cr1, Zeb2, Gzma and Prdm1 (Fig. 4c, Extended Data Fig. 5a), a dichotomy was observed in the expression of multipotency-associated (e.g. Myb, Ccr7) and effector-associated (e.g. Tbx21, Lgals1) genes within the TCM pool (denoted as TCM(mult) and TCM(eff), respectively in the figures; Fig. 4b, Extended Data Fig. 5a). Next, we assessed the relation between transcriptional state and replicative history within the memory T cell pool. In line with the flow cytometry data, the replicative history of TCM—as a whole—exceeded that of TEM, thereby validating the scRNAseq approach. Strikingly, TCM enriched for effector genes had overall higher DRRFP/DRGFP ratios compared to TCM enriched for multipotency genes, demonstrating that stemness-related transcriptomic features are inversely associated with division history within the TCM pool (Fig. 4d). Correspondingly, comparison of the three TCM MCs with the highest and lowest level of prior division (hdTCM and ldTCM, respectively) revealed that ldTCM were marked by the expression of key multipotency-associated genes, including Tcf7, Sell, Myb and Eomes, and several survival factors (Gimap and Birc family members, Extended Data Fig. 5b, c). Moreover, one ldTCM MC was highly enriched for transcripts involved in inhibitory function (Lag3, Cd160, Tox), suggesting a possible analogy with the inhibitory signaling-dependent TCM- precursor subset identified by Johnnidis et al.[33] (Extended Data Fig. 5c). In contrast, hdTCM commonly expressed genes related to terminal differentiation, such as Lgals1 and S100 family members, and showed increased transcript levels for cytotoxicity-associated genes (Nkg7, Ctsw; Extended Data Fig. 5b, c). This link between replicative history and a multipotency versus effector-associated gene expression signature within the TCM pool was further validated by differential gene expression analysis and gene set enrichment analysis (Fig. 4e-g, Extended Data Fig. 5d). In line with this association, ex vivo antigen stimulation of TCM harvested from Lm-OVA memory mice showed that TCM that had undergone more prior divisions were more likely to degranulate and less likely to produce IL-2, as compared to their less divided TCM counterparts (Extended Data Fig. 5e, f).

Extended Data Fig. 5

Replicative history identifies distinct transcriptional states within the TCM pool.

Single cell transcriptomic profiling of DR+ T cells obtained from spleen in memory phase (Day 75 and 85 post Lm-OVA infection). a, Log2 enrichment of selected genes in each MC cluster. Boxplots indicate group median and 25th/75th percentiles, whiskers indicate the interquartile range multiplied by 1.5, dots signify individual MCs. The phenotype clusters TEM, TCM(eff) and TCM(mult) contain 4, 9 and 10 MCs, respectively. For definition of TCM(eff) and TCM(mult), see Fig. 4B. b, Top and bottom marker genes of ldTCM (Top, MC2, 11, 14) and hdTCM (Bottom, MC6, 8, 18), see Fig. 4D for ldTCM and hdTCM definitions. c, Heatmaps depicting z-score transformed enrichment values of genes related to cell survival (left), cytotoxicity and effector function (middle), inhibitory markers (top-right), and transcription factors involved in T cell multipotency (bottom-right). Expression is depicted for the 3 ldTCM and 3 hdTCM MCs. d, Volcano plot depicting differentially expressed genes in ldTCM versus hdTCM. Significantly (adjusted P-value < 0.05) differentially expressed genes are depicted in red. Selected genes are highlighted. e, Cytokine release of CD27HIKLRG1LO DR+ T cells (isolated from spleen at day >60 post infection) 4 hours post ex vivo stimulation. Percentage DRRFP cells within cytokine producers (+) and non-producers (-), relative to the average DRRFP percentage within each sample, is depicted. Lines connect individual ex vivo stimulated samples samples (n=12), obtained from 3 mice. f, Ex vivo degranulation of CD27HIKLRG1LO DR+T cells (isolated from spleen at day >60 post infection) 4 hours post ex vivo stimulation. Percentage DRRFP cells within the CD107a/b positive (+) or negative (-) cell populations is depicted. Lines connect individual samples ex vivo stimulated samples (n=17), obtained from 5 mice. g, Enrichment of gene signatures from MSigDB (Hallmark) by gene set enrichment analysis comparing ldTCM and hdTCM. Data depicted was accumulated in two independent experiments (3-4 mice per experiment). P values were determined by Tukey’s HSD test (a), Wilcoxon Rank Sum test with Bonferroni correction (d), two-sided Wilcoxon signed-rank test (e, f), the FGSEA algorithm followed by the Benjamini-Hochberg procedure (g). P values < 0.05 are indicated.

The observed divergence in replicative history between distinct TCM states potentially reflects the selective quiescence of a subset of TCM with a less differentiated state. Of note, ldTCM showed reduced expression of Myc targets and genes involved in cell metabolism (Extended Data Fig. 5g), suggesting that these cells exist in a transcriptionally-enforced replicative quiescent state. To test for such a transcriptional state, we scored the expression of a core gene set of quiescent stem cells from various tissues[36] (hereafter referred as QstemScore). Notably, TCM that showed increased expression of multipotency-associated genes were marked by a higher QstemScore than TCM with increased expression of effector-associated genes. (Fig. 4h). Moreover, variation in QstemScore could also be detected in gp33-specific P14 TCM from an external data-set[37], and those P14 TCM that prominently expressed this gene set transcriptionally resembled the multipotency-signatureHI, effector-signatureLO OT-I ldTCM described here (Extended Data Fig. 6). Together, these data suggest a link between TCM quiescence and the expression of multipotency-associated genes, driving the divergence in replicative history between distinct TCM states.

Extended Data Fig. 6

gp33-specific P14 TCM with increased expression of genes associated with replicative quiescence resemble OT-I ldTCM.

Re-analysis of scRNAseq profiled splenic of P14 memory T cells, published in Kurd et al. (Kurd et al., Science Immunology, 2020). a-b, 2D projection of P14 memory T cells 90 days post LCMV infection, colors indicate individual MCs (a), or the relative expression of effector- and multipotency-associated genes (b). Gene list in Supplementary Table 1. c, P14 memory T cells cluster into TCM (blue) and TEM (red). 2D projection colored by subset (top), and violin plots depicting normalized UMI counts of selected genes (bottom) are shown. d, QstemScore of all TCM MCs in the Kurd et al. dataset. e, Pearson correlations between the Kurd et al. P14 TCM MCs that score high (MC1, 3) or low (MC6, 7) for QstemScore, and all OT-I TCM MCs described here. Data are depicted as waterfall plots, asterisks indicate significant correlations. TCM(eff), TCM(mult), ldTCM and hdTCM MCs are defined in Figure. 4. P values were determined by two-sided Pearson correlation test followed by Bonferroni correction (e). P values < 0.05 are indicated in the plots.

To directly test whether replicative behavior in the TCM pool is associated with a multipotency- associated state and relates to the functional capacity of TCM to re-expand upon secondary activation, we established a DivisionRecorder-independent, CTV-based serial-transfer approach (Fig. 5a). Naïve OT-I and GFP;OT-I T cells were transferred into primary recipients that were subsequently exposed to Lm-OVA infection. At day 30 post-infection, early memory T cells were harvested, CTV labeled and transferred into infection-matched secondary recipients. 75 days later, CTVHI (div0-2) and CTVLO (div5+) TCM were isolated, and the resulting TCM populations were then profiled by scRNAseq, or transferred at a 1:1 ratio into tertiary recipients that were subsequently challenged with Lm-OVA. Strikingly, comparison of quiescent (div0-2) cells and proliferative (div5+) TCM by gene set enrichment analysis revealed a clear negative association between quiescence and an effector-like transcriptional state, while quiescence was positively associated with multipotency-associated gene expression (Fig. 5b, c, Extended Data Fig. 7a). Likewise, inspection of MCs (Extended Data Fig. 7b-e) that were enriched in the div0-2 cells, showed a prominent expression of multipotency-associated genes (Myb, Tcf7, Id3), whereas those enriched in div5+ cells showed increased expression of effector-associated genes (Id2, S00a4, Lgals1) (Fig. 5d, e). Furthermore, comparison of the expansion potential of div0-2 and div5+ TCM demonstrated that quiescent TCM were superior in generating offspring upon renewed infection (Fig. 5f, g), further demonstrating that replicative heterogeneity in the TCM pool is both linked to transcriptional state and functionality.

Fig. 5

Replicative history is linked to recall-potential within the TCM pool.

a, Experimental setup. Primary recipient mice received 5*105 naïve OT-I and 5*105 naïve GFP;OT-I T cells. 30 days after Lm-OVA challenge, CD8 T cells were enriched, labelled with CTV and transferred into infection-matched secondary recipient mice (1 primary recipient per secondary recipient). At d105 post infection, splenic CD27+KLRG1- memory T cells that had either divided 0-2 or 5+ times and were either GFP+ or GFP- were isolated by FACS. b, Enrichment of gene signatures from MsigDB (C7, collections deposited by Goldrath (GR) and Kaech (KA), Supplementary Table 2) between Div0-2 and Div5+ cells. Top and bottom 5 pathways are depicted. c, Enrichment plots of representative pathways detected in by gene set enrichment analysis. d, Ratio of normalized counts between Div0-2 and Div5+ cells within each MC separately calculated for GFP+ and GFP- populations. Bars indicate averages, dots indicate ratios of either GFP+ or GFP- OT-I T cells. Red dotted lines indicate a fold change of 2. e, Waterfall plots depicting top and bottom 6 marker genes for selected MCs, filtered for genes involved in immune function (Supplementary Table 3). f, Flow cytometry plots depicting pre-transfer mixes of Div0-2 and Div5+ TCM. g, 8,000-12,000 total memory T cells as described in f were transferred into infection-naïve mice, following Lm-OVA challenge 24 hours later. Ratios between Div0-2 and Div5+ derived cells was determined from peripheral blood samples at indicated days post infection. Lines connect populations from individual mice (Experiment 1 n = 3; Experiment 2 n = 5). Depicted scRNAseq data was collected from 4 mice, data describing recall potential was obtained from 8 mice. P values were determined by the FGSEA algorithm followed by the Benjamini-Hochberg procedure (e).

Extended Data Fig. 7

Single cell mRNA sequencing analysis of highly divided and less divided splenic TCM.

a, Volcano plot depicting differentially expressed genes in Div0-2 versus Div5+ TCM. Significantly differentially expressed genes (Adjusted P < 0.05) are depicted in red. Selected immune-related genes are highlighted. b, Cell count per MC. c, Number of sequenced cells per sample included in the analysis. d, Sample composition of each MC. e, 2D projection, colors indicate different MCs.Depicted scRNAseq data was collected from 4 individual mice. P values were determined by Wilcoxon Rank Sum test with Bonferroni correction (a).

Re-expansion potential of TCM is linked to prior division

Having observed a link between prior division and recall potential in adoptive transfer experiments, we set out to verify this relationship without disruption of the TM niche, through re-challenge of recipient mice carrying DR+ memory OT-I T cells. In case the capacity for renewed expansion would primarily be restricted to replicative quiescent TCM cells, the fraction of DRRFP cells should show an initial decay upon reinfection—due to the increased preponderance of offspring derived from this previously quiescent population—followed by a gradual recovery throughout the contraction phase, as a result of novel division-dependent label acquisition. Notably, analysis of the fraction of DRRFP T cells in blood revealed a steep decline during the first days post-secondary infection, followed by a gradual recovery during secondary memory formation (Fig. 6a, Supplementary Note 4). This transient reduction in the DRRFP fraction was observed in multiple anatomical compartments (blood, spleen, liver), occurred independent of cell phenotype, and was also observed in LCMV-OVA induced TM pools responding to secondary challenge (Fig. 6b-d). Of note, DRRFP cell accumulation during the secondary contraction phase occurred at a comparable rate as during the primary response (Fig. 6e), yielding a secondary TM pool that—despite extensive renewed clonal expansion—had undergone a similar number of divisions as the initial memory pool (Fig. 6f, median fold difference = 1.03). Thus, the replicative histories of the TEFF and TM pools of the secondary T cell response mimic those of the primary T cell response, supporting the notion that the secondary expansion wave is mounted by a group of TCM that has undergone limited prior division. Furthermore, this low-division TCM pool is able to repeatedly reconstitute the effector T cell pool, as the same decrease in the fraction of DRRFP cells was observed upon tertiary infection of mice (Fig. 6g).

Fig. 6

The secondary TEFF pool is predominantly generated by previously quiescent memory T cells.

a, Kinetics of the percentage of DRRFP cells in blood upon secondary Lm-OVA infection. Values are relative to the DRRFP percentage within the respective memory pools (n=6 mice), black line represents group mean. b-c, DRRFP percentages in indicated organs (b), or within splenic CD27LOKLRG1HI and CD27HIKLRG1LO populations (c) at indicated time points (n=6 mice per time point) post-secondary infection. Boxplots indicate group median and 25th/75th percentiles, whiskers represent min/max, dots represent individual samples. d, DRRFP percentages in blood at memory (day >60) and at the peak of the secondary response (day 4/5 post-recall). Memory pools were generated with LCMV-OVA, recall infection was performed with Lm-OVA. e, DRRFP acquisition in blood following primary and secondary infection. Values are relative to DRRFP percentage at the peak of the primary or secondary response. Lines represent group medians (n=6 mice per group), greyed areas represent 95% confidence intervals. f, DRRFP percentages in blood during effector and memory phases of the primary and secondary responses. Lines connect data of individual mice (n=6). g, DRRFP percentages in blood (n=5 mice) upon tertiary infection. Mice were challenged twice with Lm-OVA with a >60 days interval, and subsequently infected with LCMV-OVA >60 days post-secondary infection. Depicted data are representative of at least 2 independent experiments. P values were determined by two-sided Kruskal-Wallis test with Dunn’s multiple comparisons test (b, c), two-sided Wilcoxon signed-rank test (d, f), or repeated-measures one-way ANOVA followed by Dunnett correction (g).

To determine whether the observed data are consistent with re-expansion being driven by a memory T cell subset that becomes quiescent early in the immune response, we simulated T cell responses in which a fraction of TCM precursors acquires replicative quiescence during the primary T cell response (see Supplementary Note 5, Extended Data Fig. 8a). Specifically, T cell responses were simulated that yielded quiescent T cells at a frequency of either ~0.1% or ~1% of the TEFF pool, resulting in TM pools in which quiescent TCM accounted for ~3 and ~25 percent of the memory population (Fig. 7a). Modeling of DRRFP labeling rates during recall responses in which the potential to re-expand was either abruptly lost as a function of the number of prior divisions (fun 1 and 2), or was lost more gradually across division history (fun 3), demonstrated that the transient drop in DRRFP fractions is only consistent with models in which the capacity to re-expand is restricted to cells that have undergone limited clonal expansion (Fig. 6a, b). Furthermore, the stringency of this relation is strongly dependent on the relative size of the quiescent TCM pool (Fig. 6b).

Extended Data Fig. 8

Modelled T cell responses are consistent with the presence of a replication-competent quiescent TCM population.

a, Cartoon of the phenotype model depicting phenotypes, the considered interactions among them and the parameters associated with the interactions. Arrows indicate various events occurring during the response, such as cell division (denoted with λ), differentiation to a different phenotype (denoted with δ), cell death during contraction (denoted with μ), and recruitment toward the secondary response during recall infection (denoted with r). Subscripts indicate the phenotype of the cell that the parameter is affecting. Full list of parameters can be found in Supplementary Note 5. b-d, Best fit of the modelled T cell response to the experimental measurements depicting either cell numbers (top plot in each panel), or DRRFP percentages (bottom plot in each panel). The total number of quiescent T cells generated was either capped at 1% (b) or 0.1% (c, d) of the TEFF pool. Lines depict the modeled populations; Dots indicate the experimental measurements obtained from peripheral blood (b, d) or spleen (c). See Supplementary Note 5 for more details and calculations. Experimental data points are representative of at least two independent experiments, dots indicate individual mice (n=6 mice per time point).

Fig. 7

Modelled T cell responses are consistent with the presence of a replication-competent quiescent TCM population.

a, Division history of TCM and TEM pools generated by modelled T cell responses (see Supplementary Note 5) during which a high (capped at 1% of the TEFF pool size) or low (capped at 0.1% of the TEFF pool size) fraction of T cells acquire quiescence during the effector phase (top). 3 re-expansion functions were used to restrict which fraction of TCM with a given number of prior divisions will re-expand during recall (bottom). For reference, the division history of TCM is shown as a shaded area. b, Modelled DRRFP percentages within the CD27LOKLRG1HI and CD27HIKLRG1LO populations during secondary responses, with each re-expansion function applied to a memory pool containing either a high or low number of quiescent TCM. Black dots indicate experimental measurements. j, Best fit of the modelled T cell response (number of quiescent T cells capped to 1% of TEFF) experimental data obtained from spleen, depicting either cell numbers (left) or DRRFP percentages (right). See Supplementary Note 5 for details. Lines indicate the modeled populations; dots indicate experimental measurements.

Taken together, our data establish that replicative state is not homogeneously distributed within the TCM pool and is linked to distinct transcriptional and functional properties. Specifically, our observations are consistent with a dichotomy in the TCM pool, in which a self-renewing TCM population maintains the TM pool but marginally contributes to secondary expansion, and a replication-competent quiescent TCM population is required to form the TEFF pool that arises upon renewed infection (Fig. 6c, Extended Data Fig. 8, 9).

Extended Data Fig. 9

Model describing replicative behaviors in the CD8 +memory T cell pool.

Upon infection, antigen-specific CD8+ T cells activate and rapidly expand (phase 1, p1). Following pathogen clearance (p2), a subset of memory T cells continues to divide, resulting in a progressive increase in the replicative history of the overall T cell memory pool (dotted line). Within this population, three separate behaviors of transcriptionally disparate memory T cell pools can be distinguished.

Top) Terminally differentiated TEM cells that cease division after the inflammation phase (p1) and that are marked by high transcription of effector- and minimal expression of multipotency-associated genes ([E], [M]). Upon reactivation, these cells exert rapid effector functions, but lack the potential to re-expand.

Middle) A subgroup of TCM that continues to proliferate in the memory phase, exhibits diminished levels of multipotency-associated transcripts, and that abundantly expresses effector-associated genes. Although the functionality of these cells upon reinfection requires further study, their heightened expression of effector-associated genes suggests that these cells exert cytotoxic activity upon reinfection. The contribution of these cells to the secondary TEFF pool is limited.

Bottom) A subgroup of TCM cells that shows low expression of effector-associated genes but increased expression of multipotency-associated genes, and that exists in a near-quiescent state after the inflammation phase. Upon renewed infection, this cell pool is primarily responsible for the generation of a new wave of secondary TEFF. Based on their transcriptional profile, these cells are expected to have limited immediate cytotoxic functions.

Discussion

Here, we report the development and application of the DivisionRecorder to dissect the replicative history of cell pools in vivo. We show that this approach allows longitudinal examination of division history, and how it may be combined with technologies such as flow cytometry and scRNAseq to couple cell state to division history. In the application presented here, the DivisionRecorder requires viral transduction to introduce one of its modules. While this did not significantly disrupt cell behavior in our study, development of a fully germline encoded DivisionRecorder system will be attractive, for instance, to follow replicative behavior of cell pools that are not amenable to adoptive transfer.

Using the DivisionRecorder, we demonstrate that, as a whole, the multipotent CD8+ T cell pool has undergone substantial proliferation at the peak of the expansion phase, and continues to proliferate following pathogen clearance, resulting in a cumulative replicative age of the TCM pool that exceeds that of the TEFF and TEM pool. Previous work has shown that a fraction of CD62LHI precursor-TM divide at a lower rate than terminally differentiated effector subsets[10,11,32]. In line with this, we observed a lower fraction of Ki67HI cells within the multipotent effector pool compared to the terminally differentiated pool, early post infection. At the same time, our data indicate that this difference does not result in a reduced cumulative number of past divisions within the entire CD62LHI TEFF pool. Conceivably, these findings may be reconciled by the ability of highly proliferative CD62LLO TEFF to phenotypically convert to a less differentiated CD62LHI state[14,15,30]. Alternatively, the precursor-TCM pool may harbor a heterogeneity in replicative history that is not revealed by the phenotypic markers used.

In line with the latter possibility, by combining the DivisionRecorder with scRNAseq we reveal that, while the TCM pool has undergone substantial prior division as a whole, replicative history is heterogeneous within this pool and is associated with specific transcriptional states. First, our data demonstrate the presence of TCM that bear transcriptional similarities to TEM cells but, in contrast to TEM, remain highly proliferative in the absence of inflammation (Extended Data Fig. 9). Second, we identify a population of quiescent TCM that expresses reduced levels of effector-associated genes, and high levels of pro-survival genes and genes associated with quiescent stem cells[36]. Several recent studies have reported the early emergence of TCF-1HI and CD62LHI effector cells that develop into memory T cells exhibiting stemness features[38,39]. Moreover, Johnnidis et al.[33] propose early expression of inhibitory receptors as a mechanism preserving hallmark memory features. Although these early T cell subsets bear similarities to the quiescent TCM observed here, further investigations into the developmental origin of distinct TCM states are necessary to better understand the lineage relationships between the TCM states described here, and those present during the early phases of the T cell response.

A hallmark of immunological memory is the ability to efficiently generate a new wave of TEFF upon renewed infection. Our data demonstrate that this ability is predominantly confined to a subgroup of replicative nascent TM cells. The combined observations of a less differentiated quiescent TCM population, and the reconstitution of the secondary and tertiary TEFF pool by the output of these nascent progenitors, make a compelling argument for the presence of a bona fide stem cell population within the TM pool. A growing body of work has examined a stem cell-like memory T cell (TSCM) population[40,41], generally using cell phenotype to enrich and study these cells ex vivo. Using a function-driven, phenotype-agnostic, approach that does not require removal of cells from their niche, we observe a cell behavior that fits the profile of stem cell-like memory T cells in situ.

In high turnover tissues, such as the bone marrow[42,43], the intestinal epithelium[44,45] and skin epidermis[46,47], two distinct behaviors of multipotent progenitor cells have been described: Actively dividing cells that promote normal tissue homeostasis, and quiescent cells that have been documented to break their dormancy upon tissue injury and exhibit profound re-population capacity[42,45,48,49]. We propose that the two TCM behaviors we describe provide the T cell compartment with the same capacity for renewal. Thus, the T cell pool can be viewed as an autonomous tissue that abides by organizing principles akin to those of the hematopoietic system and solid organs.

Methods

DivisionRecorder vector generation

In order to prevent expression of Cre recombinase during bacterial cloning, a synthetic intron—containing a splice donor, a branch site, a pyridine rich region, and a splice acceptor— was inserted into the Cre gene through three-fragment isothermal assembly. To prevent low level Cre translation occurring from alternative start sites, two ATG codons (position 78 and 84) were replaced by TGT codons. Finally, the Cre start codon was replaced by an EcoRI-spacer-XhoI site, to facilitate subsequent introduction of synthetic STRs. To generate the DivisionRecorder vector, two lox511 sites were introduced into the multiple cloning site of the pMX retroviral vector. Subsequently, an eGFP gene and the modified Cre recombinase gene were introduced directly upstream and downstream of the 5’ lox511 site, respectively. Finally, a P2A element was inserted directly in between the eGFP gene and the 5’ Lox511 site. Together, this resulted in a cassette comprising from 5’ to 3’: Kozak, an eGFP gene, a P2A site, a lox511 site, an EcoRI restriction site, spacer, an XhoI restriction site, a Cre recombinase gene, and a lox511 site. In its base configuration, Cre recombinase is out of frame. Synthetic STR domains were ordered as oligonucleotides (Invitrogen) and subsequently dimerized. STR dimers were inserted via the EcoRI and XhoI sites. Full sequences of all oligonucleotides are supplied in Supplementary Table 6.

Cre-activity reporter vector generation

LoxP sites were introduced into the multiple cloning site of the pCDH-CMVp-MCS-PGK-BlastR vector. In addition, a Katushka open reading frame was introduced, resulting in a vector containing from 5’ to 3’; The CMV promoter, a floxed scrambled open reading frame, a Katushka open reading frame, the PGK promoter, and a blasticidin resistance gene.

Establishment of cell lines

The Cre-activity reporter cell line used in Figure 1 was generated by retroviral transduction of HEK 293T cells (ATCC) with the Cre-activity reporter plasmid and subsequent Blasticidin selection (2 μg/ml, InvivoGen). Transduced cells were seeded at 1% confluency, and resulting single cell-derived colonies were transferred to individual wells. Clones were then examined for efficiency of induction of Katushka expression upon transfection with Cre recombinase, and the best-performing clone was selected. Cre-activity reporter cells were cultured in IMDM (Gibco) supplemented with 8% fetal calf serum (FCS, Sigma), 100 U/ml penicillin (Gibco), 100 μg/ml streptomycin (Gibco) and 2 mM Glutamax (Gibco). A mouse embryonic fibroblast (MEF) cell line from the Ai9 mouse strain was generated by modification of E14.5 embryonic fibroblasts with a retroviral vector encoding short-hairpin RNA directed against the p53 mRNA. Resultant cells were cultured in IMDM supplemented with 8% FCS, 100 U/ml penicillin, 100 μg/ml streptomycin and 2 mM Glutamax.

Mice

C57BL/6J-Ly5.1, OT-I, UBC-GFP and Ai9 mice were obtained from Jackson Laboratories, and strains were maintained in the animal department of The Netherlands Cancer Institute (NKI). Ai9 and OT-I, and UBC-GFP and OT-I mice were crossed to obtain the Ai9;OT-I and GFP;OT-I strains, respectively. Between 5-10 mice, both male and female, of the age of 6 to 15 weeks were used for each experiment. All animal experiments were approved by the Animal Welfare Committee of the NKI, in accordance with national guidelines.

Generation of DivisionRecorder+ OT-I T cells

Platinum-E cells (Cell Biolabs Inc) cultured in IMDM supplemented with 8% FCS, 100 U/ml penicillin, 100 μg/ml streptomycin, and 2 mM Glutamax were transfected with the DivisionRecorder vector using FuGene™ 6 (Promega). Retroviral supernatant was harvested 48h after transfection and stored at -80°C. Spleens from Ai9;OT-I mice were harvested and mashed through a 70 μm strainer (Falcon) into a single cell suspension and resulting splenocytes were subsequently treated with NH4Cl to remove erythrocytes. Subsequently, splenocytes were cultured in T cell medium (RPMI (Gibco Life Technologies) with 8% FCS, 100 U/ml penicillin, 100 μg/ml streptomycin, Glutamax, 10mM HEPES (pH 7.4), MEM Non-Essential Amino Acids (Gibco), 1mM Sodium pyruvate (Gibco), 50 μM 2-mercaptoethanol, supplemented with 1 ng/ml recombinant murine IL-7 (PeproTech) and 2 μg/mL ConcanavalinA (Merck)). After 48h, splenocytes were re-seeded on RetroNectin (Takara) coated plates in T cell medium supplemented with 60 IU/mL human IL-2 and DivisionRecorder virus, and were centrifuged for 90min at 400g to allow spinfection. Virus concentration was chosen such that a transduction efficiency of approximately 10-15% was achieved, in order to minimize the occurrence of multiple retroviral integrations (Supplementary Note 6). Cells were harvested 24h later and a small aliquot was stained with anti-CD8-PercpCy5.5, anti-Vb5-PeCy7, anti-CD45.2-AF700 and DAPI to determine the fraction viable OT-I T cells (DAPI-CD8+Vb5+CD45.2+) by flow cytometry (Fortessa, BD Bioscience), which generally was around ~80%. CD8+Vb5+CD45.2+ cells that expressed GFP were considered as DivisionRecorder+ OT-I cells. Within the initial population of DivisionRecorder+ OT-I cells, the fraction of cells that already showed reporter activation (as inferred by tdTomato expression) 24h after transduction was consistently between 0.4 and 0.8%. Activated splenocytes were prepared for adoptive transfer (see below).

Infection, adoptive transfer and cell recovery

C57BL/6J-Ly5.1 mice were infected with 5,000-10,000 CFU of a recombinant Listeria monocytogenes strain that expresses ovalbumin or with 5,000 PFU artLCMV-OVA[25]. Approximately 24h later, infected mice received 5,000-40,000 DivisionRecorder+ OT-I T cells through intravenous tail vein injection. To analyze OT-I T cell responses in peripheral blood over time, 25-50 μL blood samples were obtained from the tail vein at the indicated time points, and were treated with NH4Cl supplemented with 0.2 mg/ml grade-II DNaseI (Roche) to remove erythrocytes (see Methods, Flow Cytometry). To obtain spleen and liver samples, mice were sacrificed, organs were harvested, and single cell suspensions were prepared by means of mashing through a 100μM or 70μm strainer (Falcon), respectively. Subsequently, erythrocytes were removed by treatment with NH4Cl. To purify leukocytes from single cell suspensions of liver tissue, cell suspensions were separated over a 37.5% Percoll (Sigma) density gradient. Obtained blood, spleen and liver samples were further processed for flow cytometric analysis, scRNA-sequencing or functional in vitro assays, as indicated. Samples were monitored for the occurrence of retroviral silencing; which was not observed in any of the examined samples (Supplementary Note 7)

Validation of DivisionRecorder functionality

To assess the ability of the DivisionRecorder to faithfully report on the replicative history of T cell populations using dilution of cell dyes as a reference, as described in Fig. 2d-e, we employed an experimental approach that was optimized to obtain sufficient DRRFP events within the limited number of cell divisions that can be followed using cell dyes such as CTV (i.e., by transferring a high number of cells modified at a high transduction efficiency). Conclusions from this experiment are restricted to the validation of the functionality of the DivisionRecorder in dividing CD8+ T cells. Splenic CD8+ T cells were isolated using the Mouse CD8 T Lymphocyte Enrichment Set (BD Biosciences) and were subsequently stained with CellTrace™ Violet (Thermofisher). Next, cells were activated for 16h in T cell medium supplemented with 0.05 μg/mL SIINFEKL peptide and 60 IU/mL IL-2. Following this activation step, cells were seeded onto RetroNectin® (Takara Bio) coated plates and were transduced with DivisionRecorder virus by spinfection for 4h in the presence of IL-2 and SIINFEKL peptide. Analysis of CellTrace™ Violet signal by flow cytometry indicated that the cells had not undergone a full cell division post labeling. Subsequently, 6×106 OT-I T cells were transferred into Lm-OVA infected recipients. Spleens were harvested 48h after adoptive transfer, processed into single cell suspensions and prepared for flow cytometric analysis. In order to accurately determine the fraction of DRRFP cells per division during the initial stages of the proliferative burst when cumulative switching rate is still low, analysis of a large number of DivisionRecorder+ OT-I T cells events is required. For this reason, a transduction efficiency of ~60% was chosen in these experiments, instead of the 10-15% transduction efficiency used in other experiments. Note that a high transduction efficiency will result in the more frequent occurrence of cells that carry multiple retroviral integrations. The presence of cells with multiple integrations will result in a higher, yet stable, DRRFP acquisition rate, as compared to the experimental set-up used in the remainder of the study.

Ex vivo analysis of degranulation and cytokine secretion potential of memory T cells

Spleens were harvested from recipient mice at >60 days post-infection, and CD8 T cells were isolated using the Mouse CD8 T Lymphocyte Enrichment Set (BD Biosciences). Following isolation, T cells were plated at 106 cells per well in 96-well U bottom plates in T cell medium supplemented with 0.05 μg/mL SIINFEKL peptide to selectively activate OVA-specific T cells. Following a 4hr incubation, capacity of indicated T cell populations to either produce the indicated cytokines or to degranulate was assessed. To allow analysis of cytokine production, Brefeldin A (GolgiPlug™, BD Biosciences) was added 30 minutes after initiation of T cell stimulation. To allow analysis of degranulation, T cell medium was supplemented with anti-CD107a and anti- CD107b antibodies at the initiation of T cell stimulation, and Brefeldin A (GolgiPlug™, BD Biosciences) and Monensin (GolgiStop™, BD Biosciences) were added 30 minutes after initiation of T cell stimulation. At the end of the T cell stimulation period, cells were stained for KLRG1 and CD27 and prepared for flow cytometric analysis (see below).

Flow cytometric analysis

Cells were taken up in PBS (Invitrogen) supplemented with 0.5% bovine serum albumin (BSA, Fisher Scientific), and stained with antibodies directed against the indicated cell surface proteins (1:200 dilution), for 30min on ice. To allow detection of intracellular cytokine production, cells were fixed and permeabilized with CytoFix/CytoPerm™ (BD Biosciences) according to the manufacturer’s protocol and subsequently stained using antibodies against IL-2, TNFα and IFNγ. To detect intranuclear Ki-67 expression, the Foxp3/Transcription factor Staining buffer set (eBioscience) was used. See Supplementary Table 7 for list of antibodies used in the study. All samples were acquired on a BD LSR Fortessa™ (BD Bioscience); DRGFP and DRRFP cells were identified as CD8+Vβ5+CD45.2+GFP+tdTomato- and CD8+Vβ5+CD45.2+GFP+tdTomato+, respectively. Flow cytometry data analysis was performed using FlowJo V10. An example of the used gating strategy is depicted in Extended Data Fig. 10.

Extended Data Fig. 10

Gating strategy.

General gating applied to flow cytometry data presented in the study. Single lymphocytes were first selected using morphology gates, and were subsequently gated on CD8+ T cells and transferred OT-I T cells (Vβ5+CD45.2+). Next, DRRFP and RFGFP could be directly selected, or first separated by phenotype depending on the analysis. The data presented here was analyzed from blood of a recipient of DR+ cells, and was acquired 6 days post infection with Lm-OVA. Phenotype gates other than those shown here are defined in their respective figures.

For the moving average analysis depicted in Fig. 3g and Extended Data Fig. 2e, CD8+Vβ5+CD45.2+GFP+ events were exported and further processed using the R package FlowCore[50]. In brief, outlier events (i.e., antibody aggregates/cell doublets) were removed, fluorescence intensities of each of the cell surface proteins were normalized using an inverse hyperbolic sine transformation and subsequently scaled between 0 and 1. To obtain the depicted moving averages, the fraction of DRRFP cells was calculated within windows that each contained 10% of total cells, starting with the 10% of cells with the lowest expression levels for the indicated marker, and with subsequent windows moving up by steps of 2.5%.

Single cell RNA sequencing and data analysis of DivisionRecorder modified cells

The scRNAseq dataset of DivisionRecorder modified and unmodified OT-I memory T cells was obtained in two independent experiments, comprising 11 mice in total (See Extended Data Fig. 3). Experiment 1 included 3 mice containing DR+ memory T cells (mouse 1-3), which were processed in a single batch. Experiment 2 included 4 mice containing DR+ memory T cells (mouse 4-7) and 4 mice containing memory T cells derived from naïve OT-I T cells (unmodified, mouse 8-11), which were processed in two separate batches (batch 1: mouse 4-5 and mouse 8-9, batch 2: mouse 6-7 and mouse 10-11).

Spleens of DivisionRecorder+ OT-I T cell recipient mice (n=7) or naïve OT-I T cell recipient mice (n=4) were harvested >65 days post-infection. Splenocytes were stained with fluorochrome-conjugated antibodies directed against CD8, CD45.2 and Vβ5 (See Supplementary Table 7), to allow purification of transferred cells by FACS using the BD FACSAria™ Fusion Flow Cytometer (BD Biosciences). DR+ cells were subsequently FACS purified based on their expression of RFP and GFP. Following the isolation of DRGFP and DRRFP memory T cells by FACS (FACSAria Fusion, BD Biosciences), obtained cell populations were barcode-labeled with distinct anti-mouse TotalSeq™Hashtag antibodies (TotalSeq-A0301-0306, Biolegend), and pooled, with an equal number of cells from each mouse to form the total pool of cells for scRNA-sequencing. If the amount of sorted DRRFP cells from a particular sample was limited, it was pooled together with another DRRFP sample to reduce cell loss during cell hashing (as indicated in Extended Data Fig. 3). Single-cell RNA isolation and library preparation was performed according to the manufacturer’s protocol of the 10X Genomics Chromium™ Single Cell 3’ kit, and the cDNA library was sequenced on a NextSeq™ 550 Sequencing System (Illumina). Cumulative data tallied to a total of ~15,000 cells. Feature-barcode matrices were generated using the Cell Ranger software of the 10X Genomics Chromium™ pipeline. Cells that could be ascribed to multiple samples or to no sample (inferred from the detection of multiple or no Hash tags), cells with a transcript (UMI) count lower than 1,500 and cells with a mitochondrial-gene fraction higher than 0.12 were excluded from downstream analysis. Next, cells were further filtered based on gene counts, setting upper and lower thresholds separately for each sample-batch to control for differences in sequencing depth (gene-count-thresholds: Experiment 1 [1,200-3,000], experiment 2 batch 1 [800 -2,500], experiment 2 batch 2 [1,000-3,000]). Subsequent analysis of the remaining 11,767 cells was performed using the Seurat[51] and MetaCell[34] R packages.

To examine enrichment or depletion of DRRFP cells within the different MetaCells, cell counts were first normalized across hashtags. Data obtained from the different mice were subsequently aggregated and used to calculate the ratio of DRRFP versus DRGFP cells in each MetaCell. The immune signature gene list used in several analyses was composed of gene clusters involved, or proposed to be involved in T cell function. The full gene list is described in Supplementary Table 3.

Differential gene-expression testing was performed using the FindMarkers function (Wilcoxon Rank Sum test) implemented in Seurat, comparing all ldTCM to all hdTCM. Significantly differentially expressed genes (P < 0.05) were subsequently used for gene-set enrichment analysis using the R package fgsea[52], testing for enriched gene-sets from the C7 immunologic or the H Hallmark gene-sets from Molecular Signatures Database (only including sets that consisted of >10 genes). Results from this analysis were filtered for collections deposited by Kaech and Goldrath (Supplementary Table 2), focusing on relevant CD8+ T cell biology.

To calculate the QstemScore, the log2 enrichment values of genes that were positively or negatively associated with stem cell quiescence (Supplementary Table 5) were first summed within each MetaCell resulting in a positive and a negative score. QstemScore was then obtained by subtracting the negative-score from the positive-score.

Re-analysis of LCMV specific memory T cell scRNAseq dataset

Single cell transcriptomes from P14 memory T cells (harvested from spleen at day 90 post infection) were obtained from the Gene Expression Omnibus (accession GSE131847, sample GSM3822202). All single cells from this dataset were clustered applying the MetaCell algorithm. Next, TCM MetaCells were determined based on the expression levels of core effector- and multipotency-related genes (Supplementary Table 1). QstemScores were then calculated for each of the TCM MetaCells, and the 2 highest and 2 lowest scoring MetaCells were selected. Pearson correlations were subsequently calculated between each of these 4 TCM MetaCells and all of the TCM MetaCells from the OT-I dataset described here.

CTV-based serial transfer experiment and analysis

Spleens from OT-I and GFP;OT-I mice were harvested and CD8+ T cells were isolated using the Mouse CD8 T Lymphocyte Enrichment Set (BD Biosciences) according to the manufacturer’s protocol. The obtained cells were mixed in a 1:1 ratio and transferred to 4 primary recipient C57BL/6J-Ly5.1 mice (1.5x106 T cells per recipient), and 24 hours later recipients were infected with 5,000-10,000 CFU Lm-OVA. 30 days post-infection, spleens and lymph nodes were harvested and CD8+ T cells were enriched using the Mouse CD8 T Lymphocyte Enrichment Set (BD Biosciences), replacing the supplied antibody-cocktail with a mixture of anti-mouse CD19, CD20 and CD4 biotinylated antibodies (used 1:200 each, See Supplementary Table 7 for information on antibody clones). The enriched cell pool was subsequently stained with CellTrace™ Violet (Thermofisher) and re-transferred into 4 infection-matched secondary C57BL/6J-Ly5.1 recipients. 74 days after secondary transfer (104 days post-infection) spleens and lymph nodes were harvested from the secondary recipients and stained with anti-mouse KLRG1-PE, CD27-APC, and CD45.2-AF700 (See Supplementary Table 7 for information on antibody clones). Next, stained cell-pools were first enriched for transferred cells (i.e., CD45.2+) through FACS using the BD FACSAria™ Fusion Flow Cytometer (BD Biosciences), and subsequently sorted again to obtain 4 populations of TCM based on both GFP expression and CTV dilution: KLRG1-CD27+GFP+Division0-2, KLRG1-CD27+GFP+Division5+, KLRG1-CD27+GFP-Division0-2, KLRG1-CD27+GFP-Division5+. These cell pools were then further processed for tertiary transfer or single-cell RNA sequencing.

For tertiary transfer, GFP-Division0-2 cells were mixed 1:1 with the GFP+Division5+ cells (experiment 1), or GFP+Division0-2 cells were mixed 1:1 with GFP-Division5+ cells (experiment 2), thereby controlling for potential confounding effects of the donor strain. Next, 10,000 cells of either obtained cell pool were transferred in naive tertiary recipient C57BL/6J-Ly5.1 mice (3 mice for experiment 1, 4 mice for experiment 2). 24 hours later recipients were infected with 10,000 CFU Lm-OVA and the ratio of GFP+ over GFP- cells within the transferred population (Ly5.2+) in blood was monitored by flow-cytometry over time.

For scRNAseq analysis, cell pools obtained by cell-sorting were barcode-labeled with distinct anti-mouse TotalSeq™ Hashtag antibodies (TotalSeq-A0301-0304, Biolegend), and subsequently pooled. Single-cell mRNA isolation and library preparation was performed according to the manufacturer’s protocol of the 10X Genomics Chromium™ Single Cell 3’ kit, and the cDNA library was sequenced on a NextSeq™ 550 Sequencing System (Illumina). Feature-barcode matrices were generated using the Cell Ranger software of the 10X Genomics Chromium™ pipeline, resulting in 13,064 single-cell transcriptomes. Cells that could be ascribed to multiple samples or to no sample (inferred from the detection of multiple or no Hashtags), cells with a transcript (UMI) count lower than 2,000 and cells with a mitochondrial-gene fraction higher than 0.12 were excluded from downstream analysis. Finally, cells with a gene-count of >2,800 were additionally excluded from further analysis. Subsequent analysis of the remaining 9,702 cells was performed using the Seurat[51] and MetaCell[34] R packages.

Differential gene-expression testing was performed using the FindMarkers function (Wilcoxon Rank Sum test) implemented in Seurat, comparing all CTVHI (division0-2) cells to all CTVLO (division5+) cells. Significantly differentially expressed genes (P < 0.05) were subsequently used for gene-set enrichment analysis using the R package fgsea[52], testing for enriched gene-sets from the C7 immunologic gene-sets (only including sets that consisted of >10 genes). Results from this analysis were filtered for collections deposited by Kaech and Goldrath (Supplementary Table 2), focusing on relevant CD8+ T cell biology.

For the MetaCell-based analysis, the number of cells within each hashtag-MetaCell combination was counted, and subsequently normalized to 1,000 cells within each hashtag. The ratios of CTVHI over CTVLO was then calculated separately for the GFP;OT-I and OT-I derived cells.

Statistical analysis

Flow cytometric data was acquired using BDFACSDiva (v8.0) software. Flow cytometric data was analyzed using Flowjo (v10.4.2), R (v6.3.1, ‘Action of the Toes’), and FLowCore (v1.52.1). Single cell RNA sequencing data was analyzed using R (v 6.3.1), Seurat (v3.1.1), and MetaCell (v0.3.41). Data was visualized using Graphpad (V8.4.1, Prism software) and GGplot (v3.2.1). No statistical methods were used to pre-determine sample sizes, and sample sizes were chosen based on those reported in previous publications[13,53]. Data distribution was assumed to be normal but this was not formally tested. Mice were stratified according to age and sex where appropriate. Data collection and analysis were not performed blind to the conditions of the experiments. No data points were excluded from the analyses.

Materials

All commercially available reagents are listed in Supplementary Table 8.