We used CCR7 and CD45RA markers to separate CCR7^hi CD45RA^hi naïve (Tna), CCR7^hi/CCR7^lo CD45RA^lo central and effector memory (Tcm/Tem) and CCR7^lo CD45RA^hi terminal effector memory CD45RA⁺ (Temra) cells from ≥65-year-old individuals ( Table 1 ) and studied their genome-wide gene expression by microarray. The PCA analysis of all gene signals separated the CD8⁺ T cell subsets, from which the Tna cell population was notably different from Tcm/em and Temra populations. ( Figure 1A ). CD8⁺ Tna cells were distinct from CD8⁺ Tcm/em and Temra cells in many differentially expressed genes. Altogether, 449 and 700 genes were upregulated, whereas 1063 and 1230 genes were downregulated in Tna compared to Tcm/em and Temra cells, respectively ( Figure 1B , Supplementary Table S1 ). Tcm/em and Temra subsets were closer in their expression profiles as they differentially expressed only 303 genes (184 of those upregulated in Tcm/em and 119 upregulated in Temra). Both Tcm/em and Temra cells showed upregulation of multiple membrane receptors associated with NK cell activity (KLRG1, KLRF1, KLRB1, KLRC1, KLRC2, KLRD1, NKTR, NKG7, FCR3A, KIR2D family, CD244), and molecules involved in cytotoxicity (PRF1, GZMA, GZMB, GZMH, GZMK, GNLY) ( Figure 1C , Supplementary Table S1 ). Compared to Tcm/em, Temra cells had higher expression of the genes associated with NK cell and cytotoxicity functions (KLRG1, KLRF1, KLRC1, CD244, PRF1, GZMB, GZMH, GNLY). Furthermore, the NCR1 gene, known as NK-cell marker NKp46, was the top upregulated gene in Temra population. The chemokine receptors CCR4 and CCR7, and SESN3, a stress-induced regulator of the AMPK-mTORC1 pathway were downregulated in Temra cells ( Figure 1C ). The gene ontology analysis also showed that Temra cells upregulate genes that are enriched in cytotoxicity-related pathways ( Supplementary Figure S2 ). We confirmed the downregulation of costimulatory molecules CD27 and CD28, and chemokine receptor CCR7 as well as upregulation of XCL1, KLRD1, and GZMB in Temra cells by real-time PCR ( Figure 2 ). The results show the increased NK cell-like and cytotoxicity function in Temra population.

As chronic CMV infection is considered the main driver of Temra cell accumulation and TCR repertoire inflation in old individuals, we studied the proportion of CMV pp65-responsive T cells (Tcmv) from CD8⁺ T cells after their stimulation with an overlapping peptide pool spanning the sequence of the pp65 antigen. We identified Tcmv cells by measuring activation induced markers CD154 and CD137 and correlated their numbers to the Temra counts. The ratio of Tcmv cells ranged from 0.26% to 26.3% (median:1.46%, IQR: 2.88%) of the total CD8⁺ T cell fraction ( Supplementary Figure S3 ); however, we found no correlation between Temra or Tcm/em cell proportions and pp65-responsive cell numbers ( Supplementary Figure S4A ). The T cells specific for the viral protein pp65 (peptide pool) were a mixed population of Tcm/em and Temra cells (75% - mixed, 20% mostly Temra, 5% - mostly Tcm/em; Supplementary Figures S4B, C ). The transcriptome profile of Tcmv cells was closer to Tcm/em and Temra cells and distinct from the Tna population as estimated by differential gene expression number and fold changes ( Figures 3A, B , Supplementary Table S2 ). There were altogether 564 genes that were differentially expressed in resting T cell subsets compared to Tcmv cells ( Figure 3B ). The transcriptome of CD8⁺ Tcmv cells showed strong upregulation of IFNG and TNF genes, activation-specific T cell markers, TNFRSF9, CD25, and TNFRSF4, chemokines (XCL1, XCL2, CCL3, CCL4L2, CCL20, CCL1), CRTAM and transcription factors EGR2 and NR4A3 ( Figure 3A , Supplementary Figure S5 , Supplementary Table S2 ). Expectedly, anti-apoptotic genes BCL2L1 and BCL2A1 were upregulated in Tcmv cells. The downregulated genes were functionally diverse and did not form a common group, underlining the activated status of Tcmv cells. However, most downregulated genes in Tcmv cells were SEPT9, PLAC8 and FCMR ( Supplementary Figure S5 ). Tcmv cells had also lower expression of SESN3 when compared to Tna and Tcm/em populations. The gene ontology analysis showed that CD8⁺ Tcmv cells upregulated the genes associated with cytokine and chemokine pathways ( Figure 3C ). We also compared our CMV-induced gene data to recently published dataset of stimulated and non-stimulated Tcm, Tem and Temra cells (20). After subtracting the activation-induced genes from the list of differentially expressed genes (Tcmv vs Temra, and Tcmv vs Tem/cm), we observed that the top genes that differentiate Tcmv from the two memory cell subtypes overlapped to a great extent ( Supplementary Figure S6 ). Thus, our results confirm that Tcmv cells are transcriptionally more similar to Tcm/em and Temra populations than to Tna.

During the late stages of their differentiation, CD8⁺ T cells undergo epigenetic changes in their chromatin structure (27). We used ATAC-seq to explore the chromatin changes occurring in terminal differentiation from Tcm/em to Temra cells. We identified 63920 peaks that were present in at least 2 individuals. The PCA plot of the ATAC-seq peaks showed more uniform clustering of Tcm/em and higher heterogeneity of Temra samples ( Figure 4A ). We identified 1322 peaks in 985 individual genes, which showed decreased chromatin accessibility in Temra compared to Tcm/em population ( Figure 4B , Supplementary Table S3 ). In contrast, only 127 chromosomal regions, annotated to 94 genes, gained openness in Temra cells, indicating decreasing chromatin accessibility during the CD8⁺ T cell terminal differentiation. Moreover, there were 187 gene regions in Tcm/em and 20 in Temra cells that included more than one open region per gene locus indicating lower chromatin accessibility in Temra cells on single gene levels. In total, 14.4% of the differentially accessible regions were located within 1 kb, and another 6% were approximately 1-2 kb from promoters. Most differentially accessible peaks were distal intergenic (25.5%) or located in introns (41.5%) ( Figure 4C ). The clustering of all 10609 differential peaks with FDR <0.05 separated Tcm/em and Temra populations into two distinct groups ( Figure 4D ).

The top open chromatin regions in Temra cells were the genomic loci of metalloproteases ADAM7 and ADAM28 in chromosome 8 and nuclear transcriptional regulator RARB in chromosome 3 ( Supplementary Table S3 ). ADAM28 has been implicated in the shedding of surface target proteins, such as FASL and CD40L in lymphocytes. In Tcm/em cells the top accessible genes were SEMA5A, PGM1 and CCR4 ( Supplementary Table S3 ). However, we found none of these top genes, except CCR4, to be differentially expressed in the CD8⁺ T cell subsets. The open chromatin regions in Tcm/em cells were associated with several pathways like the movement of cells, intracellular signal transduction, and vesicle-related pathways.

Taken together, in our genome-wide analysis we found the Temra cells to be epigenetically more closed than Tcm/em cells as memory cells had approximately 10-fold more open chromatin regions than Temra cells.

We next tested whether changes in chromatin accessibility directly correlated with gene expression and found 51 genes with overlap in their transcriptome and epigenome changes. The majority of the genes were downregulated in Temra cells, in agreement with decreased chromatin accessibility of regions associated with these genes ( Table 2 , Figures 5A, B , Supplementary Figure S7 ). These included well-known markers of memory T cell subset such as costimulatory molecules CD28, CD27, ICOS and chemokine receptors CXCR5, CCR8, CCR4 ( Figure 5B , Supplementary Figure S7 , Table 2 ). Tcm/em population showed higher expression and chromatin accessibility in several genes encoding transcription factors, such as LEF1, KLF7, BACH2, and AHR ( Table 2 ). In Temra cells, we found increased accessibility in the promoter region of NCR1 and PALLD, the two genes, for which expression was also upregulated in comparison to Tcm/em cells ( Figure 5B , Supplementary Figure S7 , Table 2 ). The upregulation of NCR1 in Temra cells is in line with the development of NK-like phenotype in Temra cells. PALDD is a widely expressed protein needed for the formation of the actin cytoskeleton and supports cell motility. Although the role of the PALLD gene in Temra cells remains unknown, it acts as a biomarker in heart-infiltrating immune cells of patients with ischemic cardiomyopathy (28).

In addition, we compared the donors based on their T cell response to CMV antigens. We divided the individuals into two groups of higher or lower responders (36% and 64% of individuals, respectively) based on their mean response to CMV antigens (2.14%) and searched for chromatin regions that are differentially open in Temra population between those groups. We identified 35 chromatin regions [15 of them had log2 fold change ≥abs (2)], many of them in gene promoters, where accessibility differed between CMV responders ( Supplementary Table S4 ). Most of these regions (34 out of 35) were more open in donors of high CMV responders. For example, 4 of those differentially accessible sites were annotated to the NALCN gene that codes nonselective cation channel (29).

We used a program called diffTF (30) to calculate differential transcription factor (TF) accessibility between CD8⁺ Tcm/em and Temra cell samples. Chromatin peaks in Tcm/em were strongly associated with the cellular activators such as Jun/Fos, NFkappaB, and STAT ( Figure 6A ). In contrast, Temra cells had a higher representation of T-box (TBR1, TBXT, TBX19) and zinc finger (ZNF18, ZNF121, ZNF631) transcription factor motifs. As TBR1, TBXT and TBX19 are not expressed in T cells, the motifs of T-box factors, identified by diffTF likely represent the binding sites of T-bet (TBX21) and EOMES, the two reciprocally acting T-box transcription factors in CD8⁺ T cell differentiation (31). We next searched for top TF (based on their activity) binding sites within the vicinity of 51 genes that overlapped in their transcriptome and epigenome changes. Motifs for TF7L1 (more active in Tcm/em) and TBX19 (more active in Temra) located in the intron of CD28 gene and promoter of CD27 gene, respectively ( Figure 6B ). Moreover, in the promoters of KLF7 and MAL genes, there were binding sites for Eomes and TBX21, respectively ( Supplementary Figure S8 ). KLF7 overexpression transforms T cells into less differentiated state (32), while MAL is involved in T cell signal transduction (33). Together, our results show that Temra cells have higher activity of T-box transcription factors compared to Tcm/em cells.

Although defined by CCR7 and CD45RA, the late differentiated CD8⁺ T cells are known to be a heterogeneous population (24, 34). For a more detailed analysis of terminally differentiated CD8⁺ T cells, we applied single-cell transcriptome profiling in a cohort of 11 individuals to characterize Temra-related subpopulations and their relation to age and CMV infection. We sorted the CD8⁺ Temra population marked by CCR7^lo CD45RA^hi from young (n=4) and old (n=7) individuals with high (n=4) and low (n=7) CMV antibody levels and analyzed their single-cell profiles with 10x Genomics technology. After merging and pre-processing the data we identified 10 clusters based on the gene expression patterns, which we annotated by known marker genes ( Figure 7A ; Supplementary Figure S9 ). Single-cell clustering showed cellular heterogeneity and did not directly correspond to the known populations identified in flow cytometry as the differentiation-related marker genes varied and showed intermediate expression among the clusters. This is expected as in contrast to the expression levels in single cell analysis, the flow cytometry applies the splicing difference of CD45 (RA vs RO) as a marker to identify Temra cells. Nevertheless, four cell clusters (0, 1, 2, and 6), including the three largest subpopulations, expressing CD3, CD8A, CCL5, KLRD1 and NKG7 likely corresponded to Temra cells (Temra 1-4). Temra 1, but also 2 and 4, weakly expressed EOMES. We performed a differential expression analysis to identify genes that distinguished these cell subsets (35) ( Supplementary Table S5 ). The most distinguishing gene in Temra 1 and 2 was CMC1, which was also highly expressed in the transcriptome of Temra population ( Supplementary Table S1 ). Temra 3 was negative for CMC1 but had high expression of granzymes. Clusters 5 (FCGR3A⁺, KLRD1⁺, PRF1⁺) and 9 (KLRC2⁺, KLRC3⁺, NKG7⁺) were assigned to NKT-like populations, as they expressed low levels of CD3, CD8A, CD8B, and TYROBP (DAP12), but also shared granzyme and perforin gene markers with CD8⁺ Temra 3 cells. We stained the clusters with 18 Temra and NK cell markers confirming the higher expression of cytotoxic and lack of IL7R genes in Temra and NKT-like cells ( Supplementary Figure S10 ). In addition to Temra and NKT-like cells, cluster 3 expressed genes SELL, CCR7 and CD28, which we named Tcm/em-like. This population was positive for IL7R that was also highly expressed in MAIT cells, which were represented by a distinct KLRB1⁺ population (cluster 4, SLC4A10⁺). Cluster 7 expressed early activation genes CD69 and NR4A2, representing activated T cells (Tact). Among the sorted cells, we also found gamma-delta T cells (cluster 8, TRDC⁺, TRGC1⁺) and a small population corresponding to IKZF2 Temra precursors (cluster 10, CD27⁺, LEF1⁺, TCF1⁺, IKZF2⁺) ( Figure 7A , Supplementary Figure S9 ). Our analysis of CD3⁺, CD8⁺, CCR7^lo CD45RA^hi flow cytometry sorted cells show a heterogeneity of these cells consisting of several CD8⁺ Temra and NKT subset populations with overlapping marker characteristics as well as related cells such as memory cells, MAIT and gamma-delta T cells that might come along with the sorting gating.

We next compared cell abundance across clusters to determine whether age and CMV altered the subpopulation distributions within the sorted CCR7^lo CD45RA^hi CD8⁺ T cells. For this, we calculated the density of cells from different conditions and overlaid it on the UMAP plot ( Figure 7B–C ). We also compared each cluster separately between the age and CMV groups ( Supplementary Figures S11 , 12 ). However, the studied individuals had substantial variation in their subset proportions and most of the differences were not statistically significant. Among main clusters, we observed altered frequencies of Temra 1, 2 and 3, and MAIT cells in aged individuals. CMC1⁺ Temra 1 and 2 were overrepresented among older whereas Temra 3 population was more enriched in young individuals. We also found statistically significant decrease of IL7R-expressing Tcm/em-like and MAIT cells in CMV^hi individuals, however, this association was not present in comparison between old and young individuals. Our pseudobulk analysis of differentially expressed genes in age and CMV group did not reveal substantial differences, suggesting that the cell subsets and not individual genes are altered among the sorted cell population. The results suggest that Temra populations in old and CMV^hi individuals are altered and shifted towards CMC1⁺ Temra population whereas younger individuals tend to have more Tcm/em-like and MAIT cells.

All study participants ( Table 1 ) were recruited from Tartu University Hospital from Internal medicine or Dermatology clinics and signed the informed consent. The study was conducted according to the Declaration of Helsinki and approved by the Ethics Review Committee of Human Research of the University of Tartu according to permissions no 272/T-12, 275M-17, and 368/M-4. As the study focused on aging, all participants were 65 years or older. Whole blood samples from donors were collected on several occasions during the routine medical appointments. PBMC-s were isolated by Ficoll-Paque density gradient centrifugation method and stored using CTL-Cryo ABC Media Kit (CTL) in a −150°C freezer. All study participants were studied for their CMV antibodies using luciferase immunoprecipitation system (LIPS) assay as reported in (22).

PBMC-s were thawed in X-Vivo 15 cell medium (Lonza). After washing, mix of diluted costimulants [anti-CD28 (BD) and anti-49d (BD)], pure CD40 antibody (for blocking CD40-CD154 binding; Miltenyi Biotec) and also CMV pp65 peptide pool (1 μg of each peptide per ml, PepTivator, Miltenyi Biotec) were added to each sample. As a negative control, cells were incubated only with costimulants and CD40. PBMCs were stimulated for 16-24h in CO2 enriched thermostat at 37°C.

2mM EDTA-PBS solution (to detach clumped cells) was added to the stimulated PBMC-s, followed by incubation at 37°C for 10 minutes. After washing, diluted FcR Blocking Reagent (Miltenyi Biotec) and antibody mix were added to each sample and incubated for 30 minutes at 4°C. Antibodies used for sorting were as follows: CD137 PE (Miltenyi Biotec), CD154 BV421 (BioLegend), CD3 AF488 (BioLegend), CD4 AF700 (BioLegend), CD8-BV605 (BioLegend), CD45RA PE-Cy7 (BioLegend) and CCR7 PE-Dazzle (BioLegend; Supplementary Table S6 ). 7-AAD (Miltenyi Biotec) was added to the samples to exclude nonviable cells. Samples were sorted with MA900 Multi-Application Cell Sorter (Sony Biotechnology; gating strategy shown in Supplementary Figures S3 , S13A ) and were collected at 4°C into collections tubes containing TRIzol reagent (for transcriptome analysis; Thermo Fisher Scientific) or pre-coated with 4% BSA (for epigenome analysis). Publication-grade FACS plots were generated with the help of FCS Expess 7 (De Novo Software).

RNA was isolated with RNeasy Micro kit (Qiagen) and the concentration and quality of RNA was measured with NanoDrop (Thermo Fisher Scientific) and Bioanalyzer (Agilent) using RNA 6000 Pico kit (Agilent), respectively. Samples that had RIN (RNA integrity number) ≥ 6.6 were deemed suitable for further analysis. Reverse transcription, DNA labeling, and hybridization of the labeled target to the microarray were done using GeneChip 3’ IVT Pico Kit (Thermo Fisher Scientific) and analyzed by the GeneChip Clariom S Human HT platform (Thermo Fisher Scientific). Microarray plates were read on the GeneTitan (Thermo Fisher Scientific) instrument and analyzed with the TAC 4.0 software. After the analysis, we used a threshold of at least a two-fold expression change (FDR-adjusted p-value < 0.05) to define the differentially expressed genes among CD8⁺ T cell subsets. Gene set enrichment analysis (GSEA) was performed using gProfiler webtool or Genekitr R package. For each sample, 11 μl of total RNA was used as a template for cDNA synthesis. Gene expression results were validated by qPCR (ViiA7, Thermo Fisher Scientific) using SYBR Green/ROX qPCR Master Mix (Fermentas). Primers used are shown in Supplementary Table S7 . The expression of genes of interest was normalized to the housekeeping gene β2M expression and analyzed using the comparative Ct method (53). Tna cells were used as reference group.

To see if any differentially expressed genes between Tcmv cells and resting T cell subsets (Tna, Tcm/em, Temra) are truly differentially expressed or are so due to stimulation status of Tcmv cells, we conducted an additional bioinformatic analysis. For that, we compared our list of differentially expressed genes to the list of differentially expressed genes from Rose et al. study (20), who used RNA-seq to compare transcriptomes of stimulated and unstimulated CD8⁺ T cell subsets.

Specifically, we compared

Differentially accessible chromatin regions were analyzed by ATAC-seq according to the protocol by Buenrostro et al. (54). Briefly, 50 000 CD8⁺ Temra and Tcm/em cells were isolated by FACS, followed by transposition and DNA purification steps (GeneJET Gel Extraction and DNA Cleanup Micro Kit, Thermo Fisher Scientific). Then the experiment was put on hold and DNA samples were stored in -20°C freezer until all the samples were collected. The correct number of cycles for PCR was calculated with the help of qPCR and additional 5 cycles were added. The quality of libraries was evaluated by Qubit (Thermo Fisher Scientific), and TapeStation (Agilent). All libraries were sequenced on the Illumina NextSeq500 with 75 bp paired-end reads. After sequencing, preprocessing of the data was done using the Galaxy web platform. First, Illumina adapters were trimmed by Cutadapt. In addition, Cutadapt was used to filter out low quality (quality less than 20) and short reads (< 20 bp). Samples were aligned to hg38 with Bowtie2 using –dovetail and –very-sensitive options. Next, reads that mapped to the mitochondrial genome, had low mapping quality, and didn’t pair properly, were removed. The distribution of fragment lengths ( Supplementary Figure S13B ) and mapped reads ( Supplementary Figure S13C ) was acquired from all ATAC-seq libraries and indicated a good quality. Peak calling was done by Genrich using ATAC-seq mode. DiffBind was used for group comparison (log2 fold change ≥ abs (2); FDR-adjusted p-value < 0.05). Peaks only present in at least two samples were included in the analysis. The consensus peakset was generated gathering all peaks from both cell types. Differentially open regions were annotated to genes using ChipSeeker package.

diffTF (30) was used to calculate differential transcription factor accessibility between CD8⁺ Tcm/em and Temra cells. diffTF uses ATAC-seq and transcriptional data such as RNA-seq or microarray data and classifies TFs into transcriptional activators or repressors by quantifying and correlating TF motif accessibility. diffTF was ran in default mode using 1,000 permutations. Normalized microarray data, peakset with 65,669 peaks (representing the consensus peakset of all samples that overlaps at least 10 different samples) together with recommended and pre-calculated transcription factor data (based on HOCOMOCO v11 motifs + PWMScan, see diffTF documentation) were used as input.

scRNAseq was carried out in two different reactions consisting of 5 or 6 samples ( Supplementary Figure S14 ). Samples were combined and barcoded so that every reaction contained one young CMV^hi and CMV^lo sample and two (or one) elderly CMV^hi and CMV^lo samples. For single-cell analyses, previously extracted and frozen PBMCs were thawed using RPMI 1640 media (Corning) supplemented with Penicillin-Streptomycin (Corning) and Fetal Bovine Serum (Gibco). After washing with media, cells were transferred to a cell culture dish for overnight incubation in previously described RPMI 1640+FBS+PS media. The next morning cells were collected and labeled simultaneously with TotalSeq-B Hashtags (BioLegend), CD8 Microbeads (Miltenyi Biotech) and FACS antibodies: CD3 Alexa Fluor 488, CD4 Alexa Fluor 700, CCR7 PE-Dazzle 594 and CD45RA PE-Cy7 (all from BioLegend) in 1% BSA in PBS. CD8⁺ cells were extracted using MS columns and MACS separator (both from Miltenyi Biotech). Cells were washed according to the 10X Genomics protocol CG000149 Rev C. After the washes, Temra cells were sorted using Sony MA900 Cell Sorter (Sony Biotechnology) and counted using Luna FL Automated Cell Counter and Acridine Orange/Propidium Iodide Stain (Logos Biosystems). Approximately 5000 Temra cells from 5-6 individuals with different barcodes were combined into one reaction and loaded onto the Chromium controller and cDNA was generated using Chromium Single Cell 3’ Reagent Kit (v3.1 Chemistry Dual Index) and Chromium Next GEM Chip G Single Cell Kit (both from 10X Genomics) according to manufacturer’s instructions. Gene expression and hashtag libraries were paired-end sequenced on an Illumina NovaSeq 6000 instrument in Novogene.

The raw sequencing data were processed using the 10x Genomics Cell Ranger pipeline (version 6.1.2) with the GRCh38 reference. The scRNA-seq alignment and quantification were carried out on the Tartu University High Performance Computing Center Rocket cluster. The samples were demultiplexed using the HTODemux function from the Seurat package (35). The barcodes having reads for more than one molecular hashtag were designated as doublets and excluded from the downstream processing. Further filtering was applied to remove the cells with less than 1100 and more than 12% of mitochondrial genes. We have observed cell populations that expressed high levels of CD68, which is a monocyte marker. These cells were excluded from further analyses.

The transcriptomes of the retained high-quality cells were further processed according to the Seurat data integration workflow to account for the between-batch technical differences. In brief, the data were normalized using the SCTransform function (version 2). The integration features were selected as the intersection of the top within-batch highly variable genes. Next, the integration was performed as previously described (55). The first 30 principal components were selected for building the UMAP projection and clustering the data. To visualize the cell density in CMV and aging, the data were exported for processing in Scanpy package. Gaussian kernel density estimation was used to calculate the density of cells in a UMAP space (56).

The cell identities were assigned based on a combination of the differentially expressed genes (DE) for each cluster and were identified using the FindAllMarkers function in Seurat (Wilcoxon Rank Sum test). The DE marker genes of each cluster were determined based on the following criteria (1): expressed in more than 10% of the cells within either or both two groups (2); absolute value of the log fold change > 0.25 (3); adjusted p-value < 0.01.