Defining molecular circuits of CD8⁺ T cell responses in tissues during latent viral infection

Naïve CD8⁺ T cells possess an extraordinary ability to differentiate into a diverse array of effector and memory CD8⁺ T cell subtypes upon encountering pathogens (Kaech and Cui, 2012). However, their capacity to do so has been primarily studied in the context of acute and/or chronic infections. For instance, upon the resolution of an acute infection, CD8⁺ T cells differentiate into various subsets of memory cells that can provide long-lasting immunity and can rapidly respond upon secondary exposure to the same pathogen. In contrast, infections that result in chronic persistence of the pathogen, such as those caused by HIV or hepatitis virus in humans or lymphocytic choriomeningitis virus (clone 13) in mice, among others, can result in a dysfunctional CD8⁺ T cell state with curtailed effector capabilities ineffective in combating the infection (Fenwick et al., 2019; Iannacone and Guidotti, 2022; Wherry and Kurachi, 2015; Zehn et al., 2022).

Certain viruses, however, can evade the immune system and sporadically reactivate to cause a life-long latent infection that must be controlled by an intact immune system throughout the affected individual’s life. Failure to continuously keep these latent viruses in check due to perturbations of the immune system can be detrimental to the overall fitness of the host. Among viruses that can cause life-long latent infections, cytomegalovirus (CMV) is one of the most common and best characterized, infecting >60–100% of individuals, depending on geographical location (Griffiths et al., 2015). Interestingly, both in human and in mice, CMV infection results in the generation of distinct epitope-specific CD8⁺ T cell responses that can give rise to either conventional memory or to unconventional inflationary memory responses (Munks et al., 2006; O’Hara et al., 2012; Zangger and Oxenius, 2022). During the latter, certain CMV epitope-specific CD8⁺ T cells are maintained at an elevated level in an “effector-like” state in the circulation and peripheral organs, a phenomenon that is thought to be driven by sporadic reactivation of the virus and may contribute to the surveillance of viral latency (Snyder et al., 2008).

Whether conventional memory and inflationary memory CD8⁺ T cell responses utilize shared and/or distinct molecular programs is not well understood. Moreover, how these two distinct CD8⁺ T cell responses are regulated in various tissues during latent infection also warrants further investigation. In this study, we sought to understand the dynamics and molecular relationships between the differentiation of conventional and inflationary memory CD8⁺ T cell responses following mouse CMV (MCMV) infection. To achieve this goal, we generated a new MCMV-specific CD8⁺ T cell mouse transgenic (Tg) line that allows us to track distinct epitope-specific CD8⁺ T cell responses and to directly compare the differentiation trajectories between conventional memory and inflationary memory CD8⁺ T cells within the same host experiencing the same infection across tissues and time points.

To understand the dynamics of distinct CD8⁺ T cell responses to MCMV infection, we first profiled endogenous CD8⁺ T cells for MCMV-specific M45 and M38 epitopes using peptide-MHC class I tetramers (M45^H2-Db+ and M38^H2-Kb+), which exhibit conventional and inflationary memory responses, respectively. Consistent with previous observations, the immunodominant M45-specific CD8⁺ T cell response peaked after 7 days postinfection (DPI) and contracted thereafter, whereas the inflationary M38-specific CD8⁺ T cell response peaked after 14 DPI and its proportion remained high in the blood up to 150 DPI (Fig. 1 A) (Munks et al., 2006).

Phenotypic profiling of M45- and M38-specific CD8⁺ T cells by flow cytometry during effector phase (7 DPI) revealed that M45^H2-Db+ CD8⁺ T cells showed a slight but significant increase in effector-memory precursor and a marked decrease in terminal effector (TE) proportions, with no notable difference in central-memory precursor cells, compared with M38^H2-Kb+ CD8⁺ T cells, based on the expression of CD62L and CD27 (gated on CD44⁺) (Fig. 1 B). Similarly, the expression of memory-associated markers, such as CD127 (or IL-7Ra) and CXCR3, was slightly elevated in M45^H2-Db+ CD8⁺ T cells, while short-lived TE markers KLRG1 and CX₃CR1 were enriched in M38^H2-Kb+ CD8⁺ T cells despite comparable expression of the effector molecule granzyme A and the transcription factors (TFs) TCF1 and EOMES in both CD8⁺ T cell populations (Fig. 1 C, upper panel). Moreover, while a slight preferential bias toward a “memory” phenotype was already observed on 7 DPI in both populations, stark differences were observed at 35 DPI, where M38^H2-Kb+ CD8⁺ T cells still largely showed a TE phenotype while M45^H2-Db+ CD8⁺ T cells exhibited characteristics of both central- and effector-memory cells (Fig. 1, B and C).

Although the inflationary M38^H2-Kb+ CD8⁺ T cell response largely resembled effector cells by flow cytometry, principal component analysis (PCA) of bulk RNA sequencing (RNA-seq) data suggested that splenic effector and inflationary memory M38^H2-Kb+ CD8⁺ T cells are transcriptionally distinct (Fig. S1 A). Comparative pairwise differential RNA-seq analysis between naïve, effector and memory M45^H2-Db+ and M38^H2-Kb+ CD8⁺ T cells from spleen revealed five distinct clusters (Fig. 1 D and Fig. S1 B). Cluster 1 represents genes shared by naive and memory cells, but downregulated in effector cells, including CD127 (encoded by Il7ra) and the anti-apoptotic protein BCL2 (encoded by Bcl2). Cluster 2 and Cluster 3 represent genes gradually downregulated or upregulated as naïve CD8⁺ T cells differentiate to effector and memory cells, respectively. Genes in these clusters include the TFs Stat1 and Ikzf1 in Cluster 2 and Bhlhe40 and Id2 in Cluster 3, all of which have been implicated in the differentiation of effector and memory T cells (Cannarile et al., 2006; Hoshino et al., 2022; Li et al., 2019; O’Brien et al., 2014; Quigley et al., 2008; Salmon et al., 2022; Yang et al., 2011). In addition, the effector cytokine Ifng and Il2rb (which encodes a subunit of the IL-15 receptor required for memory cell maintenance [Becker et al., 2002; Schluns et al., 2002]) are also upregulated in Cluster 3.

Cluster 4 represents genes upregulated during the naïve to effector transition, peaked in effector and decreased in memory, but maintained at an elevated level in M38-specific inflationary response compared with M45-specific memory response (Fig. 1 D). These include the effector molecule Gzmb, the negative regulator of T cell activation Ctla4 (Krummel and Allison, 1995; Leach et al., 1996), cell proliferation–associated marker Mki67, and a TCR signaling component, Plcg2 (Fu et al., 2012). Furthermore, genes downregulated during this naïve-to-effector transition and remain downregulated in inflationary M38-specific response but re-expressed by memory M45-specific CD8⁺ T cells belonging to Cluster 5, such as TFs associated with stemness, namely Tcf7 (Zhou et al., 2010) and Foxo1 (Kim et al., 2013), and the secondary lymphoid organ homing molecules Ccr7 and Sell, consistent with the observed phenotypes of these cells (Fig. 1 D). Finally, in total there are 436 differentially expressed genes (DEGs) (|log₂FC| > 0.5, false discovery rate [FDR] < 0.05) between M38- and M45-specific CD8⁺ T cells on 35 DPI, and these DEGs are enriched for effector and memory signatures, respectively (Fig. 1, E and F). Altogether, different CMV-specific CD8⁺ T cell subsets show similar phenotypic and transcriptional programs during the effector phase that eventually diverge to adopt distinct transcriptional fates upon memory acquisition.

To understand whether the underlying transcriptomic differences between conventional and inflationary memory CD8⁺ T cell responses could be explained by differences in their epigenetic program, we subjected naïve (CD62L⁺ CD44⁻ and 0 DPI), effector (CD44⁺ and 7 DPI), and memory (CD44⁺ and 35 DPI) M45^H2-Db+ and M38^H2-Kb+ CD8⁺ T cells to assay for transposable-accessible chromatin with sequencing (ATAC-seq). MCMV-specific CD8⁺ T cells undergo substantial and dynamic chromatin accessibility changes, and the chromatin accessibility profiles between the two distinct epitopes at the effector stage (7 DPI) is almost identical (Fig. 2 A). Moreover, conventional memory M45^H2-Db+ CD8⁺ T cells clustered away from effector M45- and M38-specific cells, as well as from inflationary memory M38-specific cells (Fig. S1 A). Strikingly, however, although RNA-seq suggested that effector and memory M38^H2-Kb+ CD8⁺ T cells are transcriptionally distinct, PCA of ATAC-seq data suggests that they are epigenetically quite similar (Fig. S1 C). Inflationary M38^H2-Kb+ CD8⁺ T cells (35 DPI) clustered together with both effector CD8⁺ T cell populations from 7 DPI (Fig. S1 C), indicating a decoupling between epigenetic and transcriptomic regulations in the inflationary memory but not conventional memory response, and that the epigenetic landscape of inflationary CD8⁺ T cells is more poised for an effector response despite being transcriptionally distinct from effector cells. For instance, two intergenic regions upstream of an effector chemokine Ccl3 maintain their accessibility in inflationary M38-specific CD8⁺ T cells as these cells transition from effector cells, but not during the transition of effector to memory M45-specific CD8⁺ T cells, although the overall decrease in Ccl3 transcripts can be observed in these two responses (Fig. S1 D).

Given the potential for distinct epigenetic and transcriptomic remodeling during these transition processes, we further dissected these datasets by focusing on an epitope-specific pairwise analysis of differentially accessible chromatin regions (DACRs) between M45^H2-Db+ and M38^H2-Kb+ CD8⁺ T cells during their transition from naïve to effector and to memory cells. A tabulation of DACRs during the transition of each of these epitope-specific responses suggests that M45- and M38-specific CD8⁺ T cells experience almost identical epigenetic changes during naïve-to-effector transition (Fig. 2 B), with regions experiencing greatest numerical and magnitude changes (|log₂FC| > 1, adjusted P value [P_adj] < 0.05) located within the intergenic and promoter regions. However, as M45- and M38-specific cells transition from effector cells into either conventional or inflationary memory cells, respectively, they exhibit a distinct pattern of epigenetic changes, where M38-specific CD8⁺ T cells experience a large degree of promoter “opening” during this transition and exhibit more modest changes in overall accessibility; in comparison, M45-specific CD8⁺ T cells show higher fold changes in intergenic and intronic regions (Fig. 2 B). As an example, the Tcf7 locus shows high accessibility in several intronic regions in the conventional memory M45^H2-Db+ CD8⁺ T cells compared with inflationary memory M38^H2-Kb+ CD8⁺ T cells, whereas the Gzma locus is more poised in inflationary memory M38-specific CD8⁺ T cells (Fig. 2 C), which corresponds to the transcript and protein levels between these two populations (Fig. S1 E).

De novo motif analysis of DACRs (P_adj < 0.05) between conventional memory M45^H2-Db+ versus inflationary memory M38^H2-Kb+ CD8⁺ T cells identified TCF/LEF (HMG, high mobility group) as the top enriched motif in conventional memory compared with inflationary memory cells, consistent with the accessibility of regions associated with the Tcf7 locus and the high expression of TCF1 in conventional memory M45^H2-Db+ CD8⁺ T cells (Fig. 2, C and D; and Fig. 1 C). Additionally, other motifs belonging to the IRF family and the RUNT family were also enriched in memory M45-specific cells (Fig. S1 F). In contrast, inflationary memory M38^H2-Kb+ CD8⁺ T cells were enriched in zinc finger motifs, such as those belonging to the Kruppel-like factor (KLF) family members known to be important in controlling various T cell functions, including migration and differentiation (Fig. 2 D) (Hart et al., 2012; Nah and Seong, 2022). We also found the presence of motifs belonging to the AP-1 and T-BOX family members in memory M38-specific CD8⁺ T cells, which may contribute to their “effector”-like epigenetic state (Fig. S1 F).

To further identify the strongest signature that distinguishes conventional and inflationary memory responses, we examined the differences between M45^H2-Db+ and M38^H2-Kb+ CD8⁺ T cells on 35 DPI based jointly on differences in their gene expression and associated chromatin regions (Fig. 2 E). This pairwise analysis revealed that most of the DEGs (63%, 276 of 436 genes) between conventional and inflationary subsets map to changes in chromatin accessibility in 766 regions, indicating a coordinated transcriptomic and epigenetic changes that distinguish these two responses, whereas other DEGs are likely already present in a “poised” chromatin state. Among the top signatures that demarcate coordinated changes between gene expression and chromatin states, we found TFs known to regulate memory formation (e.g., Foxo1 and Tcf7 [Kim et al., 2013; Zhou et al., 2010]), markers associated with memory cells (e.g., Ccr7, Cxcr3, Sell, and Il7r), as well as various co-stimulatory molecules (e.g., Icos, Cd28, Cd27, and Btla) were enriched in memory M45^H2-Db+ CD8⁺ T cells. In contrast, inflationary memory M38^H2-Kb+ CD8⁺ T cells were marked by loci associated with regulators of effector cells, namely the TF Prdm1 (Kallies et al., 2009), and genes important for effector functions of CD8⁺ T cells (e.g., Gzma, Gzmb, and Gzmk), as well as other markers associated with cellular proliferation (e.g., Cdc20b and Mki67), which is in line with their effector-like state. We also identified a handful of genes that have been associated with CD8⁺ T cells but have not been previously associated with inflationary responses, including Matk, Fgl2, and Birc5. Collectively, both conventional and inflationary memory responses are clearly distinguished by their divergent epigenetic trajectories during their transition from effector to memory cells, potentially driven by distinct transcriptional programs. Furthermore, these conventional versus inflationary memory responses exhibit hallmarks of memory versus effector-like transcriptomic and epigenetic signatures, respectively.

To dissect clonal MCMV epitope-specific CD8⁺ T cell responses in greater detail, we generated a TCR Tg mouse line that specifically recognizes the endogenous immunodominant MCMV epitope derived from viral protein M45 (Fig. S2, A–D; see Materials and methods). This Tg model allows us to track monoclonal M45-specific CD8⁺ T responses following adoptive transfer (hereafter referred to as M45^Tg cells). The co-adoptive transfer of M45^Tg and the M38-specific “Maxi” Tg cells (Torti et al., 2011) (hereafter called M38^Tg) recapitulates the kinetics and the phenotypic characteristics of endogenous CD8⁺ T cell responses of these two distinct epitopes (Fig. 1, A–C; and Fig. 3 A). These TCR^Tg models also mimic the contraction and predisposition of a memory phenotype in the endogenous M45-specific CD8⁺ T cells, as well as the preference of endogenous M38-specific CD8⁺ T cells to inflate and patrol the circulation and peripheral tissues (Fig. 3, B and C). Moreover, both M45^Tg and M38^Tg cells possess the ability to infiltrate tissues (as assessed by i.v. labeling) during the effector response (7 DPI), and both can establish tissue residency in various peripheral tissues, albeit to different degrees in the salivary glands where CMV takes up latency (Fig. 3, D and E). Thus, M45^Tg and M38^Tg cells resemble endogenous CD8⁺ T cell responses against MCMV infection, and this co-adoptive transfer model allows us to longitudinally monitor these distinct monoclonal epitope-specific CD8⁺ T cells within the same host.

Although we observed that the molecular programs governing the transition of conventional and inflationary memory CD8⁺ T cell responses from effectors cells diverge in the spleen (as assessed by bulk RNA- and ATAC-seq), it remained unclear how these two distinct CD8⁺ T cell responses are regulated within tissues. Moreover, both in humans and in mice, CMV is thought to primarily reside and persist in the salivary gland, and therefore the ability of conventional and inflationary CD8⁺ T cells to survey and control local virus reactivation may be influenced by tissue-specific signals and controlled by distinct molecular programs (Christo et al., 2021; Reddehase et al., 2008; Smith et al., 2015; Thom et al., 2015).

To dissect the molecular regulation underlying these distinct responses across tissues and time points, we performed multiome (paired single-cell RNA-seq [scRNA-seq] and ATAC-seq) analysis (Mimitou et al., 2021) to simultaneously measure the transcripts and chromatin accessibility of M45^Tg and M38^Tg cells isolated from spleens, livers, and salivary glands of MCMV-infected mice on 7 and 35 DPI, as prototypes for effector and conventional memory/inflationary memory responses, respectively. We also included naïve (CD62L⁺ CD44⁻) TCR^Tg cells from the spleen of M45^Tg and M38^Tg mice as day 0 after infection (Fig. S3 A). Dimensionality reduction by uniform manifold approximation and projection (UMAP) and unsupervised Leiden clustering of scRNA-seq (k_n = 30) identified five different clusters that primarily delineate time point– and tissue-specific responses (Fig. 4 A). Specifically, cluster 1 represents naïve splenic M45^Tg and M38^Tg cells, and clusters 0 and 3 primarily represent splenic and liver effector (7 DPI) and memory cells (35 DPI), respectively. Moreover, salivary gland–specific responses from 7 DPI and 35 DPI dominate cluster 2 and 4, suggesting splenic and liver CD8⁺ T cells may be more closely related to each other than to salivary gland–specific responses (Fig. 4 A and Fig. S3 B).

To further understand the relationships among these cell populations, we performed partition-based graph abstraction (PAGA) analysis (Wolf et al., 2019), using defined cell populations based on combination of time point, tissue origin, and Tg status variables as nodes. This analysis reveals a complex landscape of connectivity that is primarily shaped by tissue origin (Fig. 4 B and Fig. S3 C). Splenic naïve M45^Tg and M38^Tg CD8⁺ T cells exhibit robust connectivity to one another and to splenic M45^Tg cells from 35 DPI, indicating the relatedness of naïve CD8⁺ T cells to conventional memory CD8⁺ T cells. Notably, splenic M45^Tg cells from 35 DPI are also related to splenic M38^Tg cells from 35 DPI, hinting at some relatedness between conventional memory and inflationary memory responses in the spleen compared with other organs. Moreover, effector cells from spleen and liver, irrespective of their TCR^Tg identity, form a hyperconnected network with each other and with splenic M38^Tg cells from 35 DPI (Fig. 4 B), highlighting a shared effector program across these responses and tissues and indicating a close relationship between effector and inflationary memory cells. This spleen–liver effector network, however, does not extend as robustly to effector cells from salivary glands (Fig. 4 B), suggesting an early tissue-specific transcriptional divergence. Moreover, this analysis also reveals that splenic M45^Tg memory cells from 35 DPI act as a “central hub,” bridging connections across various T cell subsets, including effector and memory cells from spleen, liver, and salivary gland. Notably, cells from salivary gland, regardless of their TCR^Tg status and time points, also form their own network that is not as robustly interconnected with the nodes from other tissues, indicating a unique and dominant program tailored to adapt to this specific tissue environment.

To investigate the defining programs that drive the diversity of this transcriptional network, we performed a global pairwise differential expression analysis between each node. This analysis captures three primary gene expression patterns: one enriched in effector cells (n = 1,435 genes), another found in naïve/memory cells (n = 2,013 genes), and a third that specifically marks salivary gland–specific responses (n = 517 genes) (Fig. S3 D). It also highlights shared and distinct programs between the two TCR^Tg cells across different tissues and time points (Fig. S3 E). To dissect which genes are the primary drivers of these distinct programs, we focused our analysis on the top 10 most DEGs from pairwise comparisons between each node and further filtered this list into top 100 DEGs across all pairwise comparisons to perform unsupervised hierarchical clustering. Consistent with the global analysis, this refined gene list also delineates three main archetypes (Fig. 4, C and D). Archetype 1 highlights genes highly expressed by naïve populations and splenic memory M45^Tg cells, delineating naïve/central memory programs and driven by TF genes such as Satb1, Bach2, Tcf7, and Lef1, as well as surface molecules including Sell (encodes for CD62L or L-selectin), Il7r, and Slamf6 (Fig. 4, C and D). Meanwhile, genes in Archetype 2 are intimately linked to effector cells and TE differentiation. In this cluster are the chemokine receptors Cx3cr1 (Gerlach et al., 2016) and Ccr2 (Terwey et al., 2005), the cytotoxic molecules Gzma and Gzmb, TFs driving effector differentiation (e.g., Tbx21 [Intlekofer et al., 2005], Id2 [Cannarile et al., 2006; Yang et al., 2011], and Zeb2 [Guan et al., 2018; Omilusik et al., 2015]), genes involved in cytokine signaling (e.g., Il18r1 and Jak1), integrin subunits Itga4 and Itga1, and a sodium bicarbonate transporter Slc4a7 that is highly expressed across tissues but has not been implicated in effector differentiation of CD8⁺ T cells (Fig. 4, C and D). These genes particularly define effector M45^Tg and M38^Tg cells from the spleen and liver, as well as splenic M38^Tg cells from 35 DPI (Fig. 4 C), highlighting a shared signature of effector and inflationary memory cells, as previously implicated by phenotypic characterization and PAGA analyses. Lastly, Archetype 3 is strongly represented in the salivary gland for both TCR^Tg lines, irrespective of time point, and from the liver at 35 DPI (Fig. 4 C). Key genes in this cluster have been implicated in establishing and maintaining tissue residency (e.g., the activation and retention marker Cd69 [Shiow et al., 2006; Walsh et al., 2019]), the co-stimulatory molecule important for residency Icos (Peng et al., 2022), several AP-1 TFs that have been implicated in tissue-resident memory formation (e.g., Junb, Fos, and Fosb [Buquicchio et al., 2024]), a regulator of TGF-β signaling Smad7 (Nakao et al., 2000), as well as a few regulators of NFκB pathway (e.g., Rel and Nfkb1), whose contributions to tissue-residency program remain to be examined (Fig. 4, C and D). Hence, scRNA-seq captures a spectrum of CD8⁺ T cell states, primarily defined by three transcriptional programs or “archetypes,” that may be shaped by signals from the tissue microenvironment during MCMV infection.

Using bulk ATAC-seq of endogenous responses, we observed that endogenous inflationary memory CD8⁺ T cells are epigenetically distinct from conventional memory CD8⁺ T cells, with the epigenetic differences being more pronounced than the transcriptomic differences (Figs. 1 and 2). Additionally, scRNA-seq captured various spectrums of these distinct responses across tissues (Fig. 4), highlighting cross-tissue complexity and suggesting a more intricate epigenetic regulation governing these processes. Thus, we focused on single-cell ATAC-seq (scATAC-seq) to examine these processes further.

Analysis of scATAC-seq identified a greater number of clusters compared with scRNA-seq, with eight clusters as opposed to the five clusters identified by scRNA-seq, consistent with the ability of ATAC-seq to capture differences in these distinct CD8⁺ T cell responses more clearly (Fig. 5 A). Of note, scATAC-seq captured distinct epigenetic states even between naïve M45^Tg and M38^Tg cells (Cluster 5 and 6, respectively), likely reflecting the differences in the “naïve-ness” of these cells stemming from their thymic selection. Importantly, these fixed TCRs recognize different epitopes and likely have differing affinities, and thus may undergo distinct positive selection, a phenomenon that has been shown to contribute to heterogeneity of naive T cell pool (Rogers et al., 2021). Therefore, it is entirely possible that thymic selection induces differences in the epigenetic landscape between the different TCR^Tg T cells at steady state. Nonetheless, any differences in the naïve state became less pronounced during infection, as splenic and liver effector cells from these two TCR^Tg lines occupied Cluster 3, 4, and 8, with M38^Tg cells more abundant in Cluster 3 and 4, and M45^Tg cells in Cluster 8, likely attributable to the differences in conventional and inflationary responses during the effector phase, as observed in the endogenous response as well (Fig. S4 A, Fig. 1, and Fig. 3 B). Genes that mark these three clusters are associated with an effector phenotype, including Klrg1, Gzma, Cx3cr1, and S1pr5, among others (Fig. S4 C). Moreover, while Cluster 3 is predominantly composed of M38^Tg cells, it also includes a balanced contribution from splenic and liver cells across time points, highlighting that epigenetic heterogeneity exists in the inflationary response (Fig. S4 A). TCR^Tg cells isolated from the salivary gland, however, are minimally represented in these clusters, but instead found primarily in Clusters 1 and 2, which also consist of cells derived from 35 and 7 DPI, respectively. Importantly, Cluster 1 also contains a handful of cells isolated from the liver. These two clusters are enriched in genes associated with tissue-homing, retention, and homeostasis that have been implicated as hallmarks of tissue-resident memory cells, including Cdh1 (Best et al., 2013; Buquicchio et al., 2024) and P2rx7 (Borges da Silva et al., 2018; Borges da Silva et al., 2020; Stark et al., 2018) (Fig. S4 C). In addition, Cluster 7 is mostly comprised of cells from the spleen and liver at 35 DPI, likely representing circulating memory cells from both TCR^Tg lines. Altogether, these data mirror our scRNA-seq findings, albeit with increased granularity, indicating that salivary gland–specific responses significantly diverge from splenic and liver responses.

Moreover, genes associated with the three transcriptional archetypes identified by scRNA-seq also faithfully occupy their expected clusters in the scATAC-seq data (Fig. 5 B). For instance, Ccr7, a naïve/central memory gene, is enriched in naïve and memory cells (Cluster 5, 6, and 7), while Gzma, indicative of effector signature, scores highly in predominantly effector clusters (2, 3, 4, and 8) (Fig. 5 B). Moreover, consistent with NFκB signaling genes identified by scRNA-seq (e.g., Rel and Nfkb1), genes involved in this pathway, namely Tnfaip3 and Tnfrsf1b (encoding for TNFR2), mark salivary gland–specific responses (Clusters 1 and 2) (Fig. 5, B and C). The precise roles of these NFκB-related genes in regulating this tissue-specific response or residency program, however, remain to be explored. Hierarchical clustering of differential gene activity scores (n = 2,796 features) and differential peaks (n = 3,636 features) (log₂FC > 0.5, FDR < 0.05) across nodes (defined by combined DPI, TCR^Tg type, and tissue origin) also supported this finding as that M45^Tg and M38^Tg cells isolated from the salivary gland at 7 and 35 DPI cluster separately from their counterparts derived from the spleen and liver (Fig. 5 C and Fig. S4 B).

Consistent with the idea that the phenotypic diversity of these populations is primarily driven by these three transcriptional archetypes, open chromatin regions surrounding key regulatory TFs like Tcf7 (naïve/central memory), Tbx21 (effector), and Smad7 (salivary gland/residency) exhibited dynamic and co-accessible regulation that aligns with the phenotypic characteristics of the cell population in which these genes are most accessible (Fig. 5 D). Thus, the diversity of conventional and inflationary CD8⁺ T cell responses across various tissues is also mirrored by their chromatin landscape, and tissue-specific signals may drive the convergence of conventional and inflationary responses in the salivary gland.

The differentiation of distinct subsets of CD8⁺ T cells requires a stepwise and coordinated network of transcriptional and epigenetic circuitry along their trajectories. Although many factors can modulate this differentiation process, the tissue microenvironment can play a significant role by enforcing tissue-adaptation programs that may or may not influence the core differentiation program required for the generation of effector and memory cells.

To investigate the dynamic differentiation trajectories of conventional memory M45^Tg and inflationary memory M38^Tg cells across tissues, we performed pseudotime trajectory analyses on each TCR^Tg population within distinct tissues. Naïve cells from each TCR^Tg line were selected as a starting point given that specification of memory subsets may have already occurred during early priming events. Our analysis revealed that M45^Tg and M38^Tg cells in the spleen and liver follow distinct trajectories, diverging into different paths as they progress toward the 35 DPI endpoint (Fig. 6 A). However, in the salivary gland, M45^Tg and M38^Tg cells share a similar trajectory, suggesting a convergence of their differentiation pathways within this tissue. By breaking down the pseudotime trajectory into distinct clusters, we also observed that conventional M45^Tg and inflationary M38^Tg cells from the spleen and liver occupy unique clusters at midpoint and terminal stages of the trajectories. Specifically, naïve splenic M45^Tg cells transited to Cluster 3 before predominantly moving through Cluster 8 and 7 toward the terminal stage. In contrast, splenic M38^Tg cells primarily occupy Cluster 3 from midpoint through near-terminal stages of the trajectory (Fig. S4 D).

These differences in trajectories are reflected by the changes in gene score activity for the canonical effector marker Cx3cr1 and the memory marker Il7r. In M45^Tg cells, the Cx3cr1 gene score activity returns to basal levels in the spleen and liver, whereas it remains elevated in M38^Tg cells in these tissues. Conversely, the Il7r gene score activity rebounds in M45^Tg cells but gradually decreases in M38^Tg cells (Fig. S4 E). Of note, while cells differentiating in the liver showed a trajectory similar to those in the spleen, a subset of cells at the terminal end of the liver trajectory also appeared in Cluster 1. This cluster, along with Cluster 2, primarily represent the midpoint and terminal stages of conventional and inflationary responses in the salivary gland, and both of which are enriched for a “tissue-residency” signature (Fig. S4 D). This is consistent with the ability of these distinct CD8⁺ T cell responses to form tissue-resident memory cells in the salivary gland and other peripheral tissues (Smith et al., 2015; Thom et al., 2015; Welten et al., 2021), as illustrated by the gradual increase of the gene score activity of Fosb in both tissues (Fig. S4 E). Altogether, our data suggest that the continuum of molecular reprogramming that governs conventional and inflationary memory responses within CD8⁺ T cells is evident at the epigenetic level and are distinct across tissues.

To elucidate the transcriptional programs involved in the differentiation of these CD8⁺ T cell responses across tissues during MCMV infection, we examined the accessibility of regions containing TF motifs along the pseudotime trajectory. We first confirmed that regulatory TFs for specific transcriptional archetypes are enriched in the expected clusters. For example, a TCF motif, crucial for the naïve/central memory program, is enriched in naïve/memory cells, whereas T-BOX and SMAD motifs are present in clusters associated with effector responses and salivary gland/tissue-resident responses, respectively (Fig. 6 B). Across tissues and TCR^Tg populations, the differentiation toward conventional and inflationary memory cells is associated with a loss of TCF/LEF motifs compared with their naïve starting points, although these motifs re-appear at terminal trajectory stages in the spleen and liver, but not in the salivary gland (Fig. 6 C).

Furthermore, motifs from the T-BOX family are prominently enriched at the midpoint of the trajectory in all tissues, likely highlighting how T-BET and EOMES govern the effector response (Intlekofer et al., 2005). The enrichment of T-BOX motifs, however, becomes less pronounced in the spleen at terminal trajectory stages, despite being maintained in the liver and salivary gland. Remarkably, the presence of AP-1 motifs in various tissues mirror those of T-BOX motifs, suggesting a potential coordination between these two TF families in regulating effector and memory programs across tissues (Fig. 6 C). Moreover, consistent with bulk ATAC-seq data (Fig. 2 D), the differentiation of inflationary memory M38^Tg cells in the spleen and liver is characterized by an abundance of various KLF motifs (Fig. 6 C), whereas conventional memory M45^Tg cells are distinctly enriched for RUNX motifs (Fig. 6 C and Fig. S1 F), highlighting the potential differences in transcriptional regulation between the response of inflationary and conventional CD8⁺ T cells, respectively, and corroborating our bulk ATAC-seq analyses of the endogenous CD8⁺ T cell responses. In addition to the potential divergent regulation of these two responses by a distinct set of TFs, the NFκB motif exhibits a distinct tissue-specific pattern between the spleen and liver, particularly in M45^Tg cells. In the spleen, NFκB motifs are relatively enriched at the start of the trajectory but only moderately enriched at terminal stage, whereas in the liver, NFκB motifs show strong enrichment only at the terminal trajectory stage, underscoring the potential of this TF family in regulating tissue-specific responses and programs.

Unlike in the spleen and liver, however, SMAD motifs are present early and persist throughout the trajectory in both TCR^Tg populations from the salivary gland, alongside other TFs observed across tissues (Fig. 6 C). This observation is consistent with the requirement of TGF-β signaling for driving the formation of CD103⁺CD69⁺ tissue-resident memory cells in the salivary gland and other barrier tissues (Christo et al., 2021; Crowl et al., 2022; Hirai et al., 2021; Mackay et al., 2013; Zhang and Bevan, 2013). Additionally, the NFκB motif is also relatively abundant in salivary gland, in line with our observation of the elevated expression of several NFκB signaling components (e.g., Tnfrsf1b, Tnfaip3, and Nfkb1) in this tissue (Fig. 4, C and D; Fig. 5, B and C; and Fig. 6 C). Consistent with this analysis, we also observed that both TCR^Tg populations showed elevated levels of phospho-SMAD2/3, two TFs downstream of TGF-β, and the NFκB subunit p65 in the salivary gland (Fig. 6 D). Hence, the differentiation of both conventional and inflationary memory responses is associated with the presence of distinct TFs and their motifs along the differentiation trajectories, and these differentiating CD8⁺ T cells possibly integrate tissue-specific signals that may influence their differentiation in specific tissues.

To better understand the coordinated molecular regulation that drives these two distinct CD8⁺ T cell responses during MCMV infection, we integrated the scRNA- and scATAC-seq datasets. Dimensionality reduction by UMAP on the gene integration matrix revealed six clusters that highlight tissue-specific and time point–specific responses, largely preserving the structure observed in the individual scRNA-seq and scATAC-seq analysis (Fig. 7 A). Notably, the time point–specific response in the salivary gland, regardless of the TCR^Tg status, converged into a single cluster (Cluster 3), which reinforces our observation of a shared regulatory program in this tissue and suggests that conventional and inflationary CD8⁺ T cell responses may adopt a unified differentiation path in the salivary gland where MCMV resides in latency (Fig. 7 A). However, despite adopting a unified program in the salivary gland, the inflationary versus conventional memory potentials of these salivary gland–specific cells are still primarily dictated by their respective viral-derived epitopes (Fig. S4 F). Furthermore, peak-to-gene link analysis also provided insights into the coordinated regulatory program that defines the three transcriptional archetypes across tissues (correlation cutoff = 0.45, FDR < 0.05): genes that define the naïve/central memory signatures (e.g., Tcf7, Ccr7, Bach2, Il7r, Klf2, Sell, Lef1, Bcl2, Themis, and Satb1), the tissue residency signature (e.g., Icos, Smad7, Runx3, Tnfaip3, Nfkbia, Jun, Jund, Fos, Cd69, and Itgae), and the effector signature (e.g., Il18r1, Tbx21, Gzma, Zeb2, Klrg1, Cx3cr1, and S1pr5), all regulated in a coordinated manner (Fig. 7, B and C).

Using this integrated data, we also identified several putative TFs whose motifs positively correlate with their gene expression and the accessibility of regions containing these motifs (Fig. 7 D). Among those are known regulators of CD8⁺ T cell differentiation, such as Batf, Eomes, Fos, Klf2, Tcf7, and Lef1, alongside TF family members, including NFκB (Nfkb1, Rela, and Relb), SMAD (Smad2), FOXO (Foxo1 and Foxo3), and ZEB (Zeb1). Additionally, we uncovered potential novel regulators of the antiviral CD8⁺ T cell response, such as Smarcc1, Jund, and Creb1, among others. Of note, expression of these putative TF regulators varied between TCR^Tg lines and exhibited highly dynamic patterns along the pseudotime trajectory across tissues. For instance, Batf is highly expressed at the endpoint of trajectory in all tissues, except in splenic M45^Tg cells (Fig. S5 A). In contrast, Klf2, a known regulator for CD8⁺ T cell migration that is downregulated during CD8⁺ T cell transition to residency (Skon et al., 2013), remained at high levels in the spleen and liver but showed minimal expression in the salivary gland (Fig. S5 A).

Motif analysis offers insights into the potential for TFs to bind accessible regions, but it does not predict how the binding of TFs to regulatory elements might modulate the expression of the associated genes. To better examine how these TFs might regulate various CD8⁺ T cell differentiation programs across tissues, we assessed the presence of TF motifs in peaks linked to DEGs based on peak-to-gene link analysis (Granja et al., 2021). We then summarized the list of genes associated with the putative binding of each TF motif into a module score and plotted the score along the pseudotime trajectory (Fig. 7 E). In addition, we also defined phenotypic signatures for effector, memory, and residency, based on the shared and distinct DEGs identified by our scRNA-seq (Fig. S3 F; see Materials and methods). By comparing these three archetypes with TF module signatures along the pseudotime, we aimed to capture the dynamics of transcriptional regulation conferred by these TFs as they relate to the phenotypic changes of CD8⁺ T cells during differentiation, under the assumption that genes co-regulated by these TFs, whether directly (in cis) or indirectly (in trans), would exhibit similar patterns of gene expression.

Consistent with their observed phenotypic signature, the effector module score of splenic M45^Tg cells increased and then declined, whereas splenic M38^Tg cells exhibited a prolonged effector program along the pseudotime compared with M45^Tg cells (Fig. S5 B, upper panel). Furthermore, the memory signature gradually increased toward the end of the trajectory in the spleen, whereas our residency signature captured the early upregulation of the residency program observed in the salivary gland, which was subsequently maintained until the terminal end of the trajectory for both TCR^Tg populations (Fig. S5 B, middle and bottom panels). We next correlated these phenotypic signatures with the modules of identified TFs along the pseudotime trajectory to infer the dynamics of transcriptional coregulation, where positive or negative correlation suggest co-regulation or inverse regulation, respectively.

Our analysis revealed potential mechanisms by which these identified TFs may modulate the observed phenotypic signatures. For example, the effector archetype highly correlated with the CREB1 module, and the pattern of CREB1 module closely aligned with the effector signature throughout the differentiation trajectory across tissues and between conventional and inflationary responses, suggesting the role of CREB1 as a positive coregulator of the effector program (Fig. 7 E and Fig. S5 C). Meanwhile, the memory signature gradually increased toward the terminal cell states along the trajectory of both M45^Tg and M38^Tg cells in all tissues. Unlike with CREB1, the LEF1 and TCF1 modules strongly correlated with the memory signature of M45^Tg cells across tissues but only moderately correlated with memory signature of splenic and liver M38^Tg cells (Fig. 7 E and Fig. S5 C). This suggests that the memory signature of the inflationary response in the salivary gland may be regulated by other TFs and/or that LEF/TCF modules confer distinct transcriptional modulation in these tissues, potentially even inversely regulating the memory program in the salivary gland during the inflationary response. Finally, both M45^Tg and M38^Tg cells upregulate the residency program in the salivary gland relatively quickly, with delayed kinetics in the spleen and liver (Fig. 7 E). Our analysis suggests that RELA may act as a positive modulator of the residency program in the salivary gland for both conventional and inflationary responses, but not in the spleen or liver (Fig. 7 E). In contrast, another NFκB family member, RELB, may appear to function as a universal negative regulator of all phenotypic signatures for these two distinct CD8⁺ T cell responses (Fig. S5, C and D). Additionally, other TFs, including those belonging to AP-1 family, are predicted to exhibit various patterns of transcriptional modulations (Fig. S5, C and D). Collectively, our analysis potentially reveals a complex transcriptional regulation that underscores how different members of the same TF family may regulate distinct transcriptional programs in a context- and tissue-specific manner.

Our extensive understanding of CD8⁺ T cell differentiation has primarily been gleaned from settings of acute or chronic infections often involving a single epitope-specific response. However, this process is relatively less understood in the context of latent/reactivating infections, where T cells responding to different epitopes can behave in very distinct patterns (Munks et al., 2006). In this study, we investigated epitope-specific responses during MCMV infection, a model for latent HCMV infection, an infection that is highly prevalent in human populations (Reddehase et al., 2008). We found that both conventional and inflationary memory responses, which are hallmarks of CD8⁺ T cell responses during latent infection, undergo significant and distinct transcriptomic and epigenetic remodeling as they transition from effector to memory cells in a tissue-specific manner at both population and single cell levels. Moreover, these responses can be defined by three molecular archetypes: effector, memory, and residency. Notably, our results suggest that in certain tissues, such as the salivary gland, conventional and inflationary memory programs may converge.

In this study, we also introduce a new TCR Tg mouse line that allows for further dissection of the two major distinct CD8⁺ T cell responses during latent viral infection. Modeling conventional memory CD8⁺ T cells responses with TCR Tgs in MCMV has been challenging because widely used model antigens (e.g., ovalbumin that can be recognized by OT-I T cells or lymphocytic choriomeningitis virus–derived GP33 epitope by P14 T cells) induce an inflationary response similar to the M38^Tg. Moreover, infections by other recombinant MCMV strains do not truly mimic the endogenous T cell responses (Dekhtiarenko et al., 2016). The contracting M45-specific T cell response may be unique, perhaps due to viral evasion mechanisms that restrict presentation of the M45-derived immunodominant epitope during latency (Dekhtiarenko et al., 2016; Grassmann et al., 2020; Holtappels et al., 2004). Thus, our new mouse model adds to existing tools for studying the conventional versus inflationary CD8⁺ T cell responses against MCMV infection.

Furthermore, our study provides insight into how tissue microenvironments may influence these distinct CD8⁺ T cell responses. Immune cells must adapt to their environment to function effectively, a molecular adaptation that involves co-opting transcriptional and epigenetic programs based on the specific tissue microenvironment they infiltrate. In the context of latent infection, where distinct CD8⁺ T cell responses have been primarily characterized in circulation or the spleen, we observed that tissue-specific signals may dictate the convergence of these CD8⁺ T cell responses in certain tissues or divergence in other tissues. Specifically in the salivary gland, a convergence of responses occurs early in the effector stage (7 DPI) and exhibits features of TGF-β signaling (via SMADs), known for its role in establishing tissue-resident memory T cells in various tissues such as the skin, intestine, liver, and adipose tissue, among others (Hirai et al., 2021; Mackay et al., 2013; Zhang and Bevan, 2013). TGF-β signaling has also been previously linked to global modulation of CD8⁺ T cell program that affect their effector functions (Arumugam et al., 2015; Taber et al., 2023), likely mediated by the balance of signals regulated by Smad2, Smad3, and Smad7 (Nakao et al., 2000), which are found in the salivary gland. How this signaling pathway influences the control of viral persistence/latency in the salivary gland by CMV-specific CD8⁺ T cells, however, remains to be explored.

In addition to TGF-β signaling, our findings hint that NFκB signaling may be associated with tissue-specific responses and residency programs. NFκB signaling is crucial for T cell homeostasis and early activation (Silva et al., 2014; Teixeiro et al., 2009), but its role in regulating tissue-specific immunity requires further investigation, although there are studies that have suggested that NFκB signals suppress tissue-resident memory T cell formation in the lung following influenza infection (Pritzl et al., 2023). Many signals, including TCR, co-stimulation, TNFR-mediated signaling, as well as cytokine signaling (e.g., IL-1 and IL-18), can induce NFκB transcriptional programs. Our data suggest that several receptors involved in this signaling pathway, such as Tnfrsf1b and the co-stimulatory molecule Cd28, are elevated in the salivary gland and may play a role in regulating salivary gland–specific immunity and/or residency program. Given the persistence of CMV in the salivary gland, it is also possible that sporadic TCR “tickling” may contribute to this enrichment of the NFκB program. Additionally, the salivary gland–specific CD8⁺ T cell response is also characterized by factors downstream of TCR signaling, such as AP-1 factors and NR4A TF family members, which are known to be part of the core tissue-residency program yet may be unrelated to TCR signaling, as these signatures are present even in the absence of antigen (Crowl et al., 2022).

By integrating scRNA-seq and scATAC-seq datasets, we also highlight the potential different modes in which TFs, even within the same family, might modulate phenotypic programs along the differentiation trajectory of conventional and inflationary CD8⁺ T cell responses in various tissues. Our multiome findings suggest that CREB1 may potentially modulate the effector program throughout the differentiation of these responses. In support of this hypothesis, a previous study implicated its role in regulating the cytotoxic T cell program in a Notch-dependent manner (Maekawa et al., 2008). Additionally, RELA, a member of NFκB family, likely appears to modulate the residency program of both conventional and inflationary responses in the salivary gland, but not in other tissues. This hypothesis is supported by observations that certain TFs, such as HIC1 (Crowl et al., 2022), may impart tissue residency programs only in specific tissues, whereas others, such as RUNX3 (Milner et al., 2017), HOBIT, or BLIMP1 (Mackay et al., 2016), and potentially RELB, may act as universal regulators of tissue residency. Thus, our study underscores the complexity of shared and distinct transcriptional modulations by various TFs that may depend on cell state and tissue-specific microenvironments.

In summary, our study expands our understanding of the molecular remodeling underlying distinct immune memory responses. Infections by latent pathogens, such as CMV, EBV, varicella zoster virus, HSV, and many others, pose significant health threats, especially for individuals who are immunodeficient, including newborns and transplant patients. Given the life-long nature of these infections and the need for constant immunosurveillance, perturbations in the immune system from treatments, including chemotherapy or radiotherapy can lead to severe diseases due to viral reactivation, as exemplified by cases of HCMV reactivation in recipients of bone marrow transplantation (Griffiths and Reeves, 2021; Stern et al., 2019). Thus, insights into how epitope-specific T cell responses are regulated and how they may differentially keep latent viruses at bay in various tissues may improve vaccination and therapeutic strategies against infection.

All mice used in this study were housed and bred under specific pathogen–free conditions with food and water in 12-h light–dark cycles at 72°F with 30–70% humidity at Memorial Sloan Kettering Cancer Center (MSKCC) and handled in accordance with the guidelines of the Institutional Animal Care and Use Committee at MSKCC. The following mouse strains were used in this study: C57BL/6 (CD45.2) (000664; The Jackson Laboratory), C57BL/6 CD45.1 (CD45.1) (Mercier et al., 2016), C57BL/6 Thy1.1 (Thy1.1 or CD90.1) (000406; The Jackson Laboratory), C57BL/6 Thy1.1/1.2 (Thy1.1/1.2 or CD90.1/.2), M45_985–993 Tg (M45^Tg), and M38_316–323 Tg (M38^Tg) (Torti et al., 2011). M38^Tg Maxi mouse strain was a gift from A. Oxenius (ETH Zürich, Zürich, Switzerland). Experiments were conducted using 8–12-wk-old mice in accordance with approved institutional protocols.

MCMV (Smith strain) was serially passaged through BALB/c hosts three times, and salivary gland viral stocks were prepared with a homogenizer to dissociate the salivary glands of infected mice 3 wk after infection.

Biotinylated monomers D^b-M45_983-99 (HGIRNASFI) and K^b-M38_316–323 (SSPPMFRV) were provided by the National Institutes of Health (NIH) Tetramer Core Facility and were conjugated to streptavidin fluorochrome conjugates and tetramers referred to as M45^H2-Db or M38^H2-Kb. Briefly, 10 μl of 1 mg/ml monomers were mixed with 1.2 μg of fluorophore-conjugated streptavidin (catalog no. 405226 and 405206; BioLegend) and incubated for 15 min and repeated for a total of five times on ice to generate tetramers used for staining.

C57BL/6 mice were infected with 1.1 × 10⁴ PFU of MCMV i.p.. 6 mo following infection, splenic CD8⁺ T cells specific for M45_985–993^H2-Db (NIH Tetramer Core Facility) were single-cell sorted using a BD AriaII and sequenced by iRepertoire, Inc. Paired TCR sequences were completed using IMGT database and cloned into a pMSCV-IRES-mCherry vector (a gift from M. Li) and transduced into the 58TCRαβ^KO CD8a-overexpressing cell line (a gift from M. Li). A pMSCV-IRES-mCherry (pMI-mCh) carrying the OT-I construct was used as a control. Briefly, Platinum-E packaging cells were transfected via calcium phosphate precipitation with pMI-mCh-M45 or pMI-mCh-OTI. The supernatant of transfected Platinum-E cells was collected 48 and 72 h after transfection and purified from the remaining cells by centrifugation at 1,500 rpm at 4°C for 7 min. Non-tissue culture-treated 48-well plates were coated with Retronectin (Takara) at 10 μg/ml in PBS overnight at 4°C. Retronectin was discarded, and 500 μl virus supernatant was centrifuged at 3,000 × g at 32°C for 2 h. Upon centrifugation, 50,000 CD8α-overexpressing 58TCRαβ^KO cells were added in 100 μl of DMEM and centrifuged at 800 × g 32°C for 90 min. After 48 h, transduction efficacy was assessed using flow cytometry. Transduced cells were then stained with either OVA_257–264^H2-Kb (SIINFEKL) tetramer or M45^H2-Db tetramer provided by the NIH Tetramer Core Facility. Upon identification of TCR clones that underwent productive pairing and were bound by the M45^H2-Db tetramer, we cloned this construct into the hCD2-LCR vector (a gift from M. Li), and the linearized hCD2-M45_985–993-LCR construct was injected into the blastocyst of C57BL/6 mice by Mouse Genetics Core Facility at MSKCC. Pups carrying the transgene were screened by PCR, and expression of TCR measured by M45^H2-Db tetramer staining.

For adoptive transfer studies, 30,000 and 10,000 of congenically marked M45^Tg and M38^Tg cells were co-transferred into congenically distinct wild-type recipients 1 day prior to MCMV infection. Recipient mice were infected with MCMV by i.p. injection of 1.1 × 10⁴ PFU in 0.5 ml of ice-cold PBS.

Mice were sacrificed in CO₂ chamber according to the Institutional Animal Care and Use Committee guideline. Upon sacrifice, blood was collected using cardiac puncture, lymph nodes and spleens were dissociated, and filtered through 100-μm cell strainer. To isolate lymphocytes from livers, tissues were physically mashed through 100-μm cell strainer and subjected to 40% Percoll gradient before ACK lysis. To isolate lymphocytes from lungs and salivary glands, tissues were dissociated with GentleMACS OctoDissociator with in complete RPMI (10% FBS, 1× L-glutamine, 1× sodium pyruvate, 1× BME, 1× MEM-NAA, and 25 mM HEPES) + 1 mg/ml Collagenase D (catalog no. 11088866001; Roche) using “LIDK_37C” protocol. For i.v. labeling, 3 μg of anti-CD45.2 or anti-CD8β (catalog no. 109822 or 126618; BioLegend) in 1× PBS was injected retro-orbitally per mouse 3 min before sacrifice.

Flow cytometry and cell sorting were done on Cytek Aurora (Cytek Bioscience), and BD AriaII. Endogenous splenic naïve CD8⁺ T cells were defined as CD3ε/TCRβ⁺ CD8α⁺ CD62L⁺ CD44⁻, while M45- or M38-specific CD8⁺ T cells were defined as CD3ε/TCRβ⁺ CD8α⁺ CD44⁺ and gated on M45^H2-Db+ or M38^H2-Kb+, respectively. For intracellular staining, cells were surface stained and fixed using eBioscience Fix/Perm buffer as per the manufacturer’s instruction. For phosphostain of pSMAD2/3 (catalog no. 562586; BD Biosciences), cells were surface stained and fixed in 2% pre-warmed PFA at 37°C for 10 min. Cells were then incubated with prechilled BD Phosflow Perm Buffer III (catalog no. 558050; BD Biosciences) for 30 min at 4°C and subjected to intracellular phospho staining for 50 min at room temperature. For preparation of Tg cells for single-cell multiome sequencing, congenically marked naïve TCR^Tg cells (M45^Tg and M38^Tg) were transferred into congenically distinct naïve recipient mice prior to infection with MCMV (see Materials and methods), as illustrated in Fig. S3 A. On the day of cell isolation, each mouse was injected with 50 µg of anti-ARTC2 nanobody (i.p.) 30 min before sacrifice (catalog no. 149803; BioLegend, Clone S+16a) and processed as described above (see Tissue processing and preparation of single-cell suspension). i.v. labeling was not performed prior to sorting due to its limited ability to distinguish circulating from bona fide tissue-resident cells in the liver, while i.v. labeling showed minimal signal in the salivary gland, comparable with the lymph nodes, which remain unlabeled (data not shown). Upon cell sorting into 1× PBS + 0.04% BSA, TCR^Tg cells were stained with mouse TotalSeq-A Hash Tag Oligonucleotides antibodies (tag number A0301–A0310; BioLegend), oligo-tagged anti-rat/mouse CD90.1 (catalog no. 202547, Clone OX-7; BioLegend), and anti-mouse CD90.2 (catalog no. 105345, Clone 30-H12; BioLegend). Hash-tagged cells were then pooled at approximately equal ratios for multiome sequencing.

5 × 10⁴ cells were sorted for RNA-seq as described previously (see above). For RNA-seq, RNA was isolated from sorted cell populations using Trizol (cat. no. 15596018; Invitrogen) two–three replicates each, and total RNA was amplified using SMART-seq V4 Ultra Low Input RNA kit (Clontech). Afterward, 10 ng of amplified cDNA was used to prepare Illumina HiSeq libraries with KAPA DNA library preparation chemistry (KAPA Biosystems) using eight cycles of PCR. Samples were barcoded and run on HiSeq2500 in a 50/50-bp paired-end run, using the TruSeq SBS Kit v3 (Illumina). For ATAC-seq, 5 × 10⁴ cells were isolated and ATAC-seq was performed as described previously (Buenrostro et al., 2013). Briefly, fresh cells were washed with ice-cold PBS and lysed. Transposition reaction occurred at 42°C for 45 min. DNA was purified using MinElute PCR purification kit (catalog no. 28006; Qiagen) and amplified for five cycles, with additional PCR cycles were added as determined by real-time PCR. AMPureXP Beads (catalog no. BDB564144; Beckman Coulter) at 1.5× ratio was used to clean final product, and libraries were sequenced on HiSeq2500 in a 50/50-bp paired-end run, using TruSeq SBS Kit v3 (Illumina).

For all bulk sequenced data, paired-end reads were trimmed for adaptors and low-quality reads using Trimmomatic (v0.36). Bulk RNA-seq transcripts were quantified using the quasi-mapping–based mode of Salmon (v0.13.1) (Patro et al., 2017) following trimming, based on the mm10 UCSC KnownGene model, correcting for potential GC bias. Transcripts were summarized at gene level using tximport (v1.10.1), and differential analysis was performed with DESeq2 (v1.34.0) (Love et al., 2014). PCA was conducted using normalized log₂ values of the top 500 genes calculated by DESeq2. To generate RNA-seq heatmap, mean-centered normalized log₂ values (z-scores) of all DEGs (P_adj < 0.05 and average TPM > 3) from each pairwise comparison were averaged across replicates and plotted using the ComplexHeatmap package (v2.14.0) with k-means clustering set to 5. Effector and memory signatures were obtained from the mouse version of MSigDB C7 ImmuneSigDB (systematic name M3041 and M3044, respectively). Triwise package (v0.99.5) was used to plot these signatures using the normalized log₂ values of selected genes as input.

For bulk ATAC-seq dataset in this study, peak calling, atlas generation, and annotation were described in detail previously (Lau et al., 2018). PCA of the ATAC peak atlas was performed using normalized log₂ values of the top 500 peaks, similar to bulk RNA-seq. DACR analysis was performed using DESeq2, with peaks considered differentially accessible if they had a P_adj < 0.05. To generate bigwig files, the bamCoverage function from deepTools (v3.5.6) was used, normalizing each bam file using the size factor calculated by DESeq2. For peak-centered heatmaps, size factor normalized bigwig files were averaged across replicates using the bigwigAverage function from deepTools, and a binned matrix was generated using computeMatrix function, spanning ± 500 bp from the peak center and with a bin size of 50. Selected DACR peaks (|log₂FC| > 1.5, P_adj < 0.05) were plotted using ComplexHeatmap package. Raw and processed data are available under Gene Expression Omnibus (GEO) accession number GSE106139.

Motif analyses were performed using HOMER’s findMotifsGenome.pl on all DACRs (P_adj < 0.05) between M45^H2-Db+ and M38^H2-Kb+ CD8⁺ T cells from the day 35 comparison, with parameters “-size given -len 6,8,10,12 -mset vertebrates -mask” (Heinz et al., 2010). Only the top motifs with a motif score >0.85 and P_adj < 0.05 were considered enriched. For meta coverage plots, normalized ATAC-seq reads were counted at a bin size of 50, spanning ± 1,000 bp from the motif center, with median values plotted for each bin.

Single-cell multiome ATAC and gene expression profiling was performed with 10X Genomics Next GEM Single Cell Multiome Reagent Kit A (catalog no. 1000282) and ATAC Kit A (catalog no. 1000280). We followed Chromium Next GEM Single Cell Multiome ATAC + Gene Expression Reagent Kits User Guide, with a modification in single-cell suspension preparation, which was carried out as previously described (Mimitou et al., 2021) (DOGMA-seq protocol with digitonin buffer). Briefly, cells from different conditions were multiplexed using hashtag oligonucleotides (HTO). The pooled samples were then treated with a digitonin lysis buffer (20 mM Tris-HCl, pH 7.4, 150 mM NaCl, 3 mM MgCl₂, 0.01% digitonin, and 2 U/μl of RNase inhibitor) for 5 min on ice, followed by adding 1 ml of chilled DIG wash buffer (20 mM Tris-HCl, pH 7.4, 150 mM NaCl, 3 mM MgCl₂, and 1 U/μl of RNAse inhibitor) and gentle inversion before centrifugation at 500 × g for 5 min at 4°C. The supernatant was discarded, and the cells were resuspended in DIG wash buffer, followed by counting using trypan blue and DAPI on a Countess II FL Automated Cell Counter. Cells were then processed as described by 10X User Guide. Final libraries were sequenced on Illumina NovaSeq S4 platform (R1–28 cycles, i7–8 cycles, and R2–90 cycles).

Multiome reads were mapped to the mm10 reference genome using 10X Genomics cellranger-arc count pipeline (v2.0.0). Downstream processing of CITE-seq and scRNA-seq was done using scanpy (v1.9.8) pipeline. Viable cells were identified based on library size and complexity, whereas cells with <40% of transcripts derived from mitochondria were excluded from analysis. After mitochondrial and doublet clean up (as identified by HTO), CITE-seq count matrix (CD90.1 versus CD90.2) was analyzed to identify correct TCR^Tg status from the multiome data. Briefly, raw count matrix was normalized using center log ratio normalization method adapted from Seurat (Satija et al., 2015). Upon center log ratio normalization, Euclidean distance was used to construct nearest neighbor graph, and Leiden clustering was performed for further filtering irrelevant/contaminating barcode sequences. Raw and processed data are available under GEO accession number GSE291432.

After filtering cells based on viability, mitochondrial, and doublets, the raw gene expression count matrix was normalized by median library size, followed by log transformation. Euclidean distance was used to construct a nearest neighbor graph, and clustering was performed using the Leiden algorithm (k_n = 30). PAGA analysis was performed as described (https://github.com/theislab/paga) (Wolf et al., 2019). Differential gene expression analysis between pseudobulk samples from each condition was performed using MAST (Finak et al., 2015) to identify DEGs (FDR < 0.05 & |coef| > 0). For gene count visualization by UMAP, MAGIC was used to impute transcript counts (van Dijk et al., 2018). To generate DEGs heatmap, unique genes from pairwise MAST analyses between each condition was visualized. The top 10 most DEGs (based on FDR rank) from each comparison were shortlisted, yielding the 100 genes. The average z-score from each condition was used for plotting with ComplexHeatmap.

Processed fragment files were analyzed with ArchR (v1.0.2) using default setting unless stated otherwise (https://www.archrproject.com/index.html) (Granja et al., 2021). Briefly, cells were filtered based on quality metrics with minFrags > 3000, TSS enrichment > 10, and exclusion of mitochondrial chromosomes. Only singlets identified from HTO analysis were considered for downstream analyses. Dimensionally reduction of scATAC-seq data was performed using iterative latent semantic indexing and projected using UMAP. Gene score matrix was generated using the default settings of the addGeneScoreMatrix function, summing up 100-kb upstream and downstream of each gene, and imputed with MAGIC using ArchR’s addImputeWeights function. After clustering, reproducible peaks were called using MACS2 (v2.2.7.1) on pseudoreplicates from each cluster. Cluster markers were identified using the getMarkerFeatures function with a Wilcoxon test, applying a cutoff of FDR < 0.01 and log₂FC > 0.5. The CisBP TF motif set was added using addMotifAnnotation, and TF motif deviations were computed with chromVar using the addDeviationsMatrix function. Peak co-accessibility was calculated using the addCoAccessibility function. Pseudotime trajectory analysis was performed using the addTrajectory function, with splenic naïve TCR^Tg cells specified as the starting point of each path.

scRNA-seq and scATAC-seq were integrated using ArchR’s integration method with slight modifications. Briefly, the gene expression matrix was imported, filtered based on quality control from scRNA-seq analysis, and added to the ArchR object using the addGeneExpressionMatrix function. After dimensionality reduction of the gene expression data, scRNA and scATAC modalities were reduced into a single dimension using the addCombinedDims function. Clustering was performed on this reduced dimension with a resolution of 0.4, and the gene integration matrix was generated using the addGeneIntegrationMatrix function. Peak-to-gene linkage analysis was conducted with a correlation cutoff of 0.45 and a resolution of 1,000 using the addPeak2GeneLinks function. To generate the integrated heatmap in Fig. 6 B, the plotPeak2GeneHeatmap function was used with a correlation cutoff of 0.45 and FDR < 0.05. For identification of candidate TF regulators, we correlated gene integration matrix with the motif matrix and considered only candidate TF regulators with a correlation > 0.3 and FDR < 0.05.

After identifying candidate TF regulators, the peak matrix was imported as SummarizedExperiment (v1.28.0) object and GC bias was corrected. To scan the associated motifs of candidate TFs, the motifmatchr (v1.20.0) package was used based on a previously computed mouse motif position weight matrix with the matchMotifs function, applying a P value cutoff of <5e⁻⁵. Peaks containing motifs for candidate TFs were then overlapped with peaks from the peak-to-gene linkage analysis. The genes associated with the overlapping peak-to-gene regions were considered part of a TF module, which was then filtered based on whether those genes exhibit differential expression changes in pairwise comparisons of all conditions using MAST. The remaining list of genes was used to calculate the module score of the candidate TF using the addModuleScore function based on the gene integration matrix, and the results were plotted along the trajectory. Moreover, to define the transcriptional archetypes in this analysis, we referred to differential expression analysis of scRNA-seq by MAST, as illustrated in Fig. S3 F. For the “effector” archetype, we identified the intersection of shared upregulated genes in effector cells (7 DPI) over naïve from both M45^Tg and M38^Tg across tissues. The “memory” archetype was defined by the shared upregulated genes in memory cells (35 DPI) from M45^Tg cells over naïve across tissues and subtracting it with the effector archetype signature to account for overlapping genes. For the “residency” signature, we defined shared upregulated memory genes from all tissues and salivary gland–specific upregulated memory genes from M45^Tg and M38^Tg and subtracted this signature with the shared upregulated effector genes from all tissues and salivary gland–specific upregulated effector genes. The resulting gene list was used to describe the residency archetype, and these phenotypic archetypes were summarized with addModuleScore function.

Fig. S1 relates to Figs. 1 and 2. Fig. S2 presents a schematic of M45Tg generation along with representative phenotyping of M45Tg cells at steady state. Fig. S3 contains the supporting materials related to Fig. 4, while Figs. S4 and S5 contain the analyses related to Figs. 5, 6, and 7.

All published data used to support the findings in this manuscript have been deposited to the National Center for Biotechnology Information (NCBI)’s GEO with accession number GSE106139. All raw and processed data generated in this manuscript have been deposited to NCBI’s GEO under accession number GSE291432.

We thank all members of the Sun laboratory for their feedback on this manuscript. We also thank M. Li, C. Chou, J. Giacalone, S. James, Y. Furuta (MSKCC, New York, NY, USA), A. Hill (Oregon Health & Science University, Portland, OR, USA), and C. Snyder (Thomas Jefferson University, Philadelphia, PA, USA) for providing reagents and helpful discussions. We acknowledge the use of the Mouse Genetics Core Facility and the Integrated Genomics Operation Core, funded by the National Cancer Institute Cancer Center Support Grant (P30CA008748), Cycle for Survival, and the Marie-Josée and Henry R. Karvis Center for Molecular Oncology.

J.C. Sun was supported by the Ludwig Center for Cancer Immunotherapy, the American Cancer Society, the Burroughs Wellcome Fund, and the NIH (AI100874, AI130043, AI155558, and P30CA008748).

Author contributions: E.K. Santosa: conceptualization, data curation, formal analysis, investigation, methodology, validation, visualization, and writing—original draft, review, and editing. J.M. Zhang: investigation and writing—review and editing. J.C. Sauter: investigation. M.E. Lee: investigation. B.D. Ng: investigation and methodology. S.V. Stulz: formal analysis and investigation. M. Takizawa: investigation. S. Grassmann: investigation. O.-E. Weizman: investigation. N.M. Adams: investigation. R. Chaligné: supervision. A. Oxenius: resources and writing—original draft, review, and editing. G. Gasteiger: supervision. C.M. Lau: data curation, formal analysis, methodology, resources, and software. J.C. Sun: conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, project administration, resources, software, supervision, validation, visualization, and writing—original draft, review, and editing.