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 cell (TCM) pool is marked by a higher number of prior divisions than the effector memory T cell pool, due 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 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 (T EFF ) that exceeds its precursor population size by >10,000-fold 3,4 . Unlike TN, these T EFF have the capacity to disseminate to peripheral tissues, and scan for and kill infected or transformed cells. Upon antigen clearance, around 95% of the T EFF pool succumbs to apoptosis, leaving behind a small long-lived pool of memory T cells (T M ) 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 T M , 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 (T CM ), 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 (T EM ) have a low turnover rate 8,9 . In contrast, during the effector phase, a T CM -like state has been linked to lower division speed and reduced clonal burst size compared to their T EM -like and terminally differentiated counterparts [10][11][12][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-)T CM and T EM in the effector and memory phase? (2) Is there heterogeneity in prior division within the T CM pool? (3) If so, does replicative history of cells within the T CM pool predict their capacity to mount a secondary T cell response?

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 DR GFP 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, divisiondependent, probability (p). The resulting CRE activity induces an irreversible activation of the RFP gene, giving rise to GFP + RFP + cells (hereafter referred to as DR RFP ) that pass this genetically encoded label on to subsequent generations, resulting in a cumulative increase in the DR RFP cell fraction within the DivisionRecorder + (DR + , i.e., the sum of DR GFP and DR RFP ) 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 DR RFP 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).
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 DR RFP 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 DR RFP 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 DR RFP 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 DR RFP cells was measured over time (Fig. 2a). At early time-points post cell transfer (d1-d4), a rapid increase in DR RFP 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 DR RFP 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 DR RFP 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 DR RFP 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).

CD8 + T CM 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 T EFF 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, DR RFP cells remained detectable following formation of T cell memory, thus allowing analysis of replicative history at late time points after infection (Fig. 3b).
In case T M would primarily be derived from T cells that had undergone limited proliferation upon primary antigen encounter, the fraction of DR RFP cells would be expected to decay during the contraction phase, due to the decline in the number of clonally expanded T EFF (Extended Data Fig. 1, Supplementary Note 3). However, analysis of DR RFP frequencies in blood demonstrated that the fraction of DR RFP 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 DR RFP 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 DR RFP 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 + T M pool, and extends this observation beyond the peripheral blood compartment to the sites of infection. 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., T CM or T EM ) and replicative history during different stages of the T cell response. To this end, the fraction of DR RFP cells within the T M 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 T CM -associated proteins CD27 (r rm = 0.81, P = 6.2·10 −14 ) and CD62L (r rm = 0.62, P = 5.6·10 −7 ) 15, 28,29 , and a negative relationship between prior division and the expression of the T EM -associated proteins KLRG1 (r rm = −0.83, P = 9.0·10 −15 ) and CX 3 CR1 (r rm = −0.75, P = 4.5·10 −11 ) 14,15,30 . Likewise, defining multipotent T CM and terminally differentiated T EM 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 CX 3 CR1, revealed a positive association between division history and a less differentiated cell state (Extended Data Fig. 2b). Furthermore, the division history of CD27 HI KLRG1 LO T CM present in lymph nodes equaled that of T CM in the spleen, implying that division history is dictated by cell state rather than anatomical location (Extended Data Fig. 2c).
Next, to delineate at which point the divergence in replicative history between T cells with a T CM -like multipotent and T EM -like terminally differentiated phenotype developed, we assessed the link between phenotypic marker expression and DR RFP fractions throughout the T cell response. Notably, replicative history varied minimally across T EFF 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 DR RFP within the CD27 HI KLRG1 LO early-T CM 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 CD27 LO KLRG1 HI 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 CD27 HI KLRG1 LO T cell pool (Extended Data Fig.  1h). The substantial number of divisions that we observe in the CD27 HI KLRG1 LO 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 transdifferentiation between T EFF 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 T CM pool results from both strong proliferative activity during the effector phase and selective proliferative activity after pathogen clearance.

Replicative history identifies distinct T CM cell states
Increasing evidence suggests that the T CM 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 DR GFP and DR RFP memory OT-I T cells (75-85 days post Lm-OVA infection; Extended Data Fig. 3). In addition, to test whether DR + OT-I T M assume the same spectrum of transcriptional states as non-modified T cells, we performed scRNAseq on OT-I T M 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 T EM and 19 T CM 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 T EM MCs-consistent with the relationship between precursor frequency and T EM formation 35 -DR + OT-I T cells and unmodified OT-I T cells were equal in their potential to yield the 19 distinct T CM 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 T CM states.
Amongst the observed T CM MCs, two transcriptionally distinct subgroups could be identified (Fig. 4b). Specifically, while all T CM 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 multipotencyassociated (e.g. Myb, Ccr7) and effector-associated (e.g. Tbx21, Lgals1) genes within the T CM pool (denoted as T CM (mult) and T CM (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 T CM -as a whole-exceeded that of T EM , thereby validating the scRNAseq approach. Strikingly, T CM enriched for effector genes had overall higher DR RFP /DR GFP ratios compared to T CM enriched for multipotency genes, demonstrating that stemness-related transcriptomic features are inversely associated with division history within the T CM pool (Fig. 4d). Correspondingly, comparison of the three T CM MCs with the highest and lowest level of prior division (hdT CM and ldT CM , respectively) revealed that ldT CM 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 ldT CM MC was highly enriched for transcripts involved in inhibitory function (Lag3, Cd160, Tox), suggesting a possible analogy with the inhibitory signaling-dependent T CM -precursor subset identified by Johnnidis et al. 33 (Extended Data Fig. 5c). In contrast, hdT CM 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 T CM 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 T CM harvested from Lm-OVA memory mice showed that T CM that had undergone more prior divisions were more likely to degranulate and less likely to produce IL-2, as compared to their less divided T CM counterparts (Extended Data Fig. 5e, f).
The observed divergence in replicative history between distinct T CM states potentially reflects the selective quiescence of a subset of T CM with a less differentiated state. Of note, ldT CM 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, T CM that showed increased expression of multipotency-associated genes were marked by a higher QstemScore than T CM with increased expression of effectorassociated genes. (Fig. 4h). Moreover, variation in QstemScore could also be detected in gp33-specific P14 T CM from an external data-set 37 , and those P14 T CM that prominently expressed this gene set transcriptionally resembled the multipotency-signature HI , effectorsignature LO OT-I ldT CM described here (Extended Data Fig. 6). Together, these data suggest a link between T CM quiescence and the expression of multipotency-associated genes, driving the divergence in replicative history between distinct T CM states.
To directly test whether replicative behavior in the T CM pool is associated with a multipotency-associated state and relates to the functional capacity of T CM 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, CTV HI (div0-2) and CTV LO (div5+) T CM were isolated, and the resulting T CM 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+) T CM 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+ T CM demonstrated that quiescent T CM were superior in generating offspring upon renewed infection ( Fig. 5f, g), further demonstrating that replicative heterogeneity in the T CM pool is both linked to transcriptional state and functionality.

Re-expansion potential of T CM 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 T M 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 T CM cells, the fraction of DR RFP 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 DR RFP 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 DR RFP fraction was observed in multiple anatomical compartments (blood, spleen, liver), occurred independent of cell phenotype, and was also observed in LCMV-OVA induced T M pools responding to secondary challenge ( Fig. 6b-d). Of note, DR RFP cell accumulation during the secondary contraction phase occurred at a comparable rate as during the primary response (Fig. 6e), yielding a secondary T M 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 T EFF and T M 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 T CM that has undergone limited prior division. Furthermore, this low-division T CM pool is able to repeatedly reconstitute the effector T cell pool, as the same decrease in the fraction of DR RFP cells was observed upon tertiary infection of mice (Fig. 6g).
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 T CM 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 T EFF pool, resulting in T M pools in which quiescent T CM accounted for ~3 and ~25 percent of the memory population (Fig. 7a). Modeling of DR RFP 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 DR RFP 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 T CM pool (Fig. 6b).
Taken together, our data establish that replicative state is not homogeneously distributed within the T CM pool and is linked to distinct transcriptional and functional properties. Specifically, our observations are consistent with a dichotomy in the T CM pool, in which a self-renewing T CM population maintains the T M pool but marginally contributes to secondary expansion, and a replication-competent quiescent T CM population is required to form the T EFF pool that arises upon renewed infection (Fig. 6c, Extended Data Fig. 8, 9).

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 T CM pool that exceeds that of the T EFF and T EM pool. Previous work has shown that a fraction of CD62L HI 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 Ki67 HI 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 CD62L HI T EFF pool. Conceivably, these findings may be reconciled by the ability of highly proliferative CD62L LO T EFF to phenotypically convert to a less differentiated CD62L HI state 14,15,30 . Alternatively, the precursor-T CM 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 T CM 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 T CM that bear transcriptional similarities to T EM cells but, in contrast to T EM , remain highly proliferative in the absence of inflammation (Extended Data Fig. 9). Second, we identify a population of quiescent T CM 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-1 HI and CD62L HI 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 T CM observed here, further investigations into the developmental origin of distinct T CM states are necessary to better understand the lineage relationships between the T CM 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 T EFF upon renewed infection. Our data demonstrate that this ability is predominantly confined to a subgroup of replicative nascent T M cells. The combined observations of a less differentiated quiescent T CM population, and the reconstitution of the secondary and tertiary T EFF pool by the output of these nascent progenitors, make a compelling argument for the presence of a bona fide stem cell population within the T M pool. A growing body of work has examined a stem cell-like memory T cell (T SCM ) 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 repopulation capacity 42,45,48,49 . We propose that the two T CM 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.

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.
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 NH 4 Cl 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 NH 4 Cl. 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 DR RFP 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×10 6 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 DR RFP 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, DR RFP 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 10 6 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); DR GFP and DR RFP cells were identified as CD8 + Vβ5 + CD45.2 + GFP + tdTomatoand 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.
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 DR RFP 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 DR GFP and DR RFP 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 DR RFP cells from a particular sample was limited, it was pooled together with another DR RFP sample to reduce cell loss during cell hashing (as indicated in Extended Data Fig. 3). Singlecell 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 samplebatch 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 DR RFP 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 DR RFP versus DR GFP 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 ldT CM to all hdT CM . 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, T CM 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 T CM MetaCells, and the 2 highest and 2 lowest scoring MetaCells were selected. Pearson correlations were subsequently calculated between each of these 4 T CM MetaCells and all of the T CM 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.5x10 6 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 antibodycocktail 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 T CM 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 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 CTV HI (division0-2) cells to all CTV LO (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 CTV HI over CTV LO 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.

Extended Data
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 T EFF undergoing clonal expansion and subsequent contraction (red), plus memory precursor T cells (MP, green) that develop into T M . Activated T EFF 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 DR RFP percentages (bottom row) are shown for 3 scenarios: (left) T EFF can give rise to MP cells during the entire expansion phase, irrespective of the number of prior divisions, (middle) only T EFF that have gone through at most 24 divisions can give rise to MP cells, or (right) only T EFF that have gone through at most 10 divisions can give rise to MP cells. Note the strong decay in DR RFP 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. Fig. 2. Evaluation of the division history of T cell subsets throughout a response to Lm-OVA. a, Gating strategy used to identify indicated T M populations (d86) in spleen samples. b, DR RFP percentages within splenic T M populations (n=6 mice) as identified in panel a. c, DR RFP percentages within the CD27 HI KLRG1 LO T CM subset in spleen and lymph nodes (LN) and within the CD27 LO KLRG1 HI T EM subset in spleen. d, Cell surface expression of CX3CR1, CD62L, and CD43 within splenic CD27 LO KLRG1 HI and CD27 HI KLRG1 LO populations at the peak of the T EFF 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 DR RFP percentages within T EFF (day 6 post infection) subsets in spleen (n=6 mice), relative to the total DR RFP percentage. g, Kinetics of DR RFP percentages within CD27 LO KLRG1 HI (left) and CD27 HI KLRG1 LO (right) DR + OT-I T cell populations in blood. Values are relative to the percentage of DR RFP 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 CD27 HI KLRG1 LO to CD27 LO KLRG1 HI cells occur only after the peak of the response at a low rate. Depicted are the overall cell numbers (left), and the percentage DR RFP cells of DR + OT-I T cells (right) in CD27 HI KLRG1 LO cells (blue), CD27 LO KLRG1 HI cells (red) and the total T cell population (green). Note that in this scenario the fraction DR RFP within the terminally differentiated CD27 LO KLRG1 HI 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 25 th /75 th 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. . Note that all T CM 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.   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 T EFF pool is limited. Bottom) A subgroup of T CM 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 T EFF . Based on their transcriptional profile, these cells are expected to have limited immediate cytotoxic functions.

Extended Data
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, DR RFP and RF GFP 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.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.

Data Availability
Transcriptomic data presented in the manuscript have been deposited to the Gene Expression Omnibus (GEO), and can be accessed under the GEO accessions GSE169154 and GSE184947. The gp33-specific P14 T cell scRNAseq dataset was retrieved from GEO (accession GSE131847, sample GSM3822202). All statistical source data of the figures presented in the present study are provided with this paper. Indicated gene sets used in gene set enrichment analyses were retrieved from the Molecular Signatures Database (MSigDB) at http://www.gsea-msigdb.org/gsea/msigdb. Any additional data supporting the findings of this study are available from the corresponding authors upon request.

Code availability
R scripts that were used to produce the main and extended data figures in the manuscript are available from GitHub (https://github.com/kasbress/DivisionRecorder_analysis).       Table 3). Selected genes are annotated, complete gene lists in Supplementary  Table 4. h, QstemScore of all T CM 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  Table 3). f, Flow cytometry plots depicting pre-transfer mixes of Div0-2 and Div5+ T CM . 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).  test (b, c), two-sided Wilcoxon signed-rank test (d, f), or repeated-measures one-way ANOVA followed by Dunnett correction (g).