Epigenomics of human CD8 T cell differentiation and aging^*

Genome-wide mapping of chromatin accessibility in CD8 T cell subsets

Using the assay for transposase-accessible chromatin (ATAC-seq) (16, 17), we generated genome-wide maps of chromatin accessibility in three human CD8 T cell subsets, purified by sorting peripheral blood mononuclear cells (PBMC): naïve (CD28^hiCD45RA^hiCD62L^hi), central memory (CD45RA^loCD62L^hi), and effector memory (CD45RA^loCD62L^lo; Figures 1A and S1). The PBMC were isolated from ten healthy individuals, comprising five donors under age 35 (“young;” 2 males and 3 females) and five over age 60 (“old;” 3 males and 2 females). We excluded terminally differentiated CD28^loCD45RA^hi effector T cells, which mostly contain oligoclonally expanded herpes-virus-specific cells, to focus on age-related chromatin changes rather than changes induced by chronic viral stimulation. To support the functional relevance of the chromatin openness patterns we identified, we also performed RNA-seq in each of the T cell subsets.

Sequencing read density at phenotype-defining genes confirmed successful transposition in each subset, with accessible chromatin at the CD8A locus and minimal accessibility at the CD4 locus (Fig. 1B). IFNG and GZMB, which are known to be more rapidly inducible in memory T cells after activation, were more open in central (CM) and effector memory (EM) than naïve cells.

Across all three subsets, we identified a total of 64,089 distinct chromatin accessibility peaks. Using predictions of chromatin states in naïve CD8 T cells from the Roadmap Epigenomics Project (18), we observed that the peaks were concentrated at promoters and enhancers, and depleted at heterochromatic regions (Fig. 1C). Hierarchical clustering of the peaks segregated the three T cell subsets (Fig. 1D), with the two memory populations clustering more closely with each other than with the naïve population. To mitigate the impact of donor heterogeneity, we excluded sex chromosomes, which have been shown to be a major site of inter-individual variation in chromatin accessibility (19).

To uncover the chromatin accessibility landscapes underlying naïve and memory cell identity in the young, we identified loci of differential accessibility between each pair of subsets using DESeq2 (20). We observed greater numbers of differentially open loci between naïve and memory subsets than between the two memory subsets (Fig. 2A); comparing young naive to young EM yielded 18,489 differentially accessible peaks (28.8% of those tested), while comparing young CM to young EM resulted in only 914 (1.4%). To determine which molecular pathways were subject to regulation, we performed k-means clustering of the differentially open peaks (Fig. 2B) and used GREAT (21) to calculate gene ontology (GO) enrichment for each cluster, separately for peaks for which openness was positively versus negatively correlated with expression of the nearest gene. We found that each cluster was enriched for GO terms reflecting the underlying T cell biology (Figures 2B and S2A). The accessibility of the peaks most open in CM and closed in naive (cluster 1; red) was positively correlated with adaptive immune response pathways, particularly those regulating leukocyte cell-to-cell adhesion. Accessibility of the peaks most open in EM and closed in naive (cluster 2; blue) was correlated with cytokine production, concordant with the increased effector cytokine production of EM cells. Moreover, peaks in this cluster were associated with the up-regulation of genes that repress cell activation, and the down-regulation of genes that facilitate activation, suggesting the induction of negative regulatory pathways in effector cells. This pattern held for genes regulating proliferation, consistent with the loss of proliferative potential following differentiation. Accessibility of the peaks more open in naive (cluster 3; green) was associated with the up-regulation of T cell proliferation, T cell receptor signaling, and V(D)J recombination, representing features of T cell biology unique to naïve cells and their progenitors. Thus, changes in chromatin accessibility across CD8 T cell differentiation were associated with the regulatory pathways that define the T cell subsets.

We next tested whether changes in chromatin openness were directly correlated with changes in gene expression. We found that peaks more open after differentiation also neighbored genes more expressed in the differentiated cells. Peaks that were closed after differentiation likewise neighbored genes less expressed after differentiation (Fig. 2C). Furthermore, when the openness changes were large, neighboring genes were more likely to change in expression (GSEA p = 0.03), and the expression changes became more extreme (Fig. S2B). This finding held both for canonical naïve and effector genes (correlation between accessibility and expression changes = 0.71, p = 1.1 * 10⁻⁶; Fig. 2D) and overall (p < 2.2 * 10⁻¹⁶, chi-squared test; Fig. S2C). These observations indicate that chromatin accessibility changes across CD8 T cell differentiation are associated with the regulation of functionally relevant pathways and genes.

Maintenance of naïve and memory T cell states is associated with distinct transcription factor families

TF binding can instigate chromatin remodeling and drive differentiation (22, 23). To identify TFs putatively associated with differentiation in T cells, we calculated the enrichment of TF binding motifs in the loci that changed in accessibility after differentiation. At peaks more open in EM than naive (Fig. 3A), we found motifs for several members of the bZIP, T-box, and RUNT TF families among the most significantly enriched (Table S1). Our top hits encompassed known regulators of T cell differentiation and development, such as Eomes and TBX21 (T-bet) (FDR < 10⁻⁶⁶; (24)). The TF family with the greatest representation was bZIP, which includes Jun-AP1 and BATF, both critical to T cell differentiation (25). These TFs were also more expressed in EM than naive. The only bZIP TF more expressed in naive cells than effector cells was Bach2, which is known to suppress effector-related genes to maintain a naive state (26); this TF also shares strong motif similarity with other bZIP TFs. Some TFs known to be overexpressed in EM, such as PRDM1 (27, 28), showed smaller, yet significant, binding motif enrichment. The motif enrichment in peaks more accessible in EM than in CM was modest (Table S1), as expected due to their similarity in accessibility patterns. We did observe minor enrichment for the binding motifs of T-bet and Eomes, which are known to be more expressed in effector than central memory populations (29). Overall, the finding of known regulators confirmed that our inference of TF binding is concordant with current models of human CD8 T cell differentiation.

Peaks that closed following differentiation were enriched for motifs recognized by the ETS and zinc finger (Zf) families of TFs (Fig. 3B), and most ETS and Zf TFs were more expressed before differentiation. ETS TFs have been shown to be required for normal T cell development, survival, and activation (30, 31), but they have not previously been reported to be central to the epigenetic maintenance of the naïve T cell state. Taken together, the motif enrichment, gene expression levels, and gene expression changes nominate BATF and Eomes as key factors in the epigenetic changes in CD8 T cell differentiation, and ETS1 and Sp1 in maintaining the chromatin accessibility landscape of naïve cells.

We hypothesized that greater TF expression plus greater chromatin openness at the TFs’ binding motifs would lead to greater expression of the TFs’ target genes. We evaluated this hypothesis using published TF target sets from knockdown experiments in human-derived cell lines (25, 32). We found that the majority of BATF targets increased in expression following differentiation (p = 0.001, young CM versus naive, and p = 0.006, young EM versus naive, Wilcoxon rank-sum tests; Fig. 3C). The trend was similar for Eomes targets (p = 0.001, young CM versus naive, and p = 0.01, young EM versus naive). Expression of BATF and Eomes targets did not differ uniformly between the two memory populations (p = 0.43 and 0.65, respectively), consistent with the high degree of epigenetic similarity between the memory populations.

In the reciprocal analysis, we evaluated ETS1 and Sp1 target genes (33–35). We found that the majority of ETS1 target genes decreased in expression following differentiation (p = 0.001, young naive versus CM, and p = 0.06, young naive versus EM). Targets of Sp1 decreased in expression on average (Fig. 3C), but these changes did not achieve statistical significance, possibly because available lists of curated Sp1 targets are less comprehensive than those for the other TFs. Overall, our analyses directly implicate ETS1 and Sp1 as maintaining the naïve state in CD8 T cells, and BATF and Eomes as contributing to the epigenetic changes induced by CD8 T cell differentiation. While this analysis suggests that these TFs are key regulators of the gene expression patterns that define T cell subset identities and functions, other family members sharing their motifs may also serve as important regulators of T cell identity or as redundancy in the transcription network circuitry. The statistical significance of the enrichment of all TF binding motifs is included in Table S1.

Accessibility changes with age in CD8 T cells are linked to TF families involved in memory differentiation

We next investigated whether the TF binding and gene expression changes associated with differentiation were preserved across age cohorts. We found that differentiation in old was similar to differentiation in young, with a strong overlap in the differential peaks and motifs identified (Fig. S3).

To examine age-related dysfunction in the maintenance of T cell states, we compared the chromatin landscapes between young and old individuals within the T cell subsets. We identified 2,864 differentially open peaks (4.5% of all peaks) across these three comparisons (Fig. 4A). Interestingly, the majority of the differentially accessible peaks associated with aging were within the naive and CM subsets, with comparatively few changes observed in the EM.

As verification of the specificity of the peak set we found to be differentially accessible with age, we performed principal component analysis (PCA) on those peaks and evaluated the sample clustering. The samples were effectively separated by age (Figure S4A), but not by donor or sex (Figure S4B), indicating that our results were not meaningfully affected by sample heterogeneity. To assess more directly whether interpersonal variability impacted our estimates of age-associated differences, we quantified their stability by randomly dropping out all possible pairs of samples (one from each age cohort), then rerunning our differential analyses. We found that no single donor drove our results (Figures S4C and S4D). We also compared our findings to those attained by using alternative normalization strategies designed to correct for variability across samples in overall levels of accessibility (36), or in GC and peak-length biases (37). These analyses confirmed that our estimates were highly robust (correlations of fold-change estimates > 0.97, p < 2.2 * 10⁻¹⁶; Figures S4E and S4F).

Using PCA to obtain a broad-based perspective on T cell aging and differentiation, we found that the first PC was sufficient to delimit the six age-specific T cell subsets (Fig. 4B). We observed that the old naive samples were closer to the CM populations than were the young naive, and we noted a similar trend in old CM relative to the EM, indicating that the epigenetic landscape of old T cells is globally biased toward differentiation. Similarly, applying k-means clustering to the peaks changing in accessibility with differentiation, we found that the old cells, and particularly the old naive, attained more effector-like patterns of chromatin structure (Fig. S5A). Furthermore, genes more expressed with age were enriched for GO categories closely aligned with those enriched in genes more expressed following differentiation (p = 4.2 * 10⁻¹³, Fisher’s exact test; Fig. S5B).

We found a similar relationship between aging and differentiation in the transcriptome for the naive subset (Fig. 4C). The first PC distinguished the three differentiation stages, and the transcriptome of the old naive was closer to the CM than was the young naive. This differentiation-oriented shift in naive cell gene expression with age was confirmed by RNA-seq of an independent validation set of young and old individuals. In contrast to the ATAC-seq results, no age-associated shift was seen in the CM subset.

When we performed k-means clustering of peaks changing in accessibility with age, we observed that age-associated closures mirrored differentiation-associated closures in the young, with peaks more open in the young naive than young CM and EM uniformly more open in young than old (top two clusters, Fig. 4D). Likewise, peaks more open in EM than naive (bottom cluster, Fig. 4D) increased in openness with age in all three subsets, indicating a widespread progression toward effector-like chromatin accessibility patterns with age. Furthermore, these peaks bordered genes with effector-like GO signatures, associated with pathways related to the regulation of immune response. GO results for the other clusters are shown in Fig. S5C.

The trend of changes in aging resembling those in differentiation was also present in the gene expression levels, but less apparent than in the chromatin accessibility patterns. Applying k-means clustering, we found that genes more expressed in the young naive than young CM and EM declined with age (green cluster, Fig. 4E), and that genes more expressed in EM than naive increased with age (dark blue cluster, Fig. 4E). However, these effects were less striking than in the chromatin accessibility results, and GO analysis of this latter cluster only yielded terms that were of low enrichment and not obviously related to effector function. The observation that the link between aging and differentiation was more evident in chromatin openness than gene expression patterns was unsurprising, as naive and memory T cell populations have been shown to differ significantly in their ability to induce gene expression upon activation, but only minimally in their transcription profiles while in the resting state (12).

To examine whether the drivers of T cell differentiation also contribute to the immune aging of CD8 T cells, we identified the TF binding motifs enriched in loci that changed in accessibility with age. We found high concordance between the enrichment in age-associated and differentiation-associated peaks (Fig. 5A), as was expected because the peak sets overlapped considerably (Fig. S5D). Specifically, the motifs enriched in regions that gained accessibility with age were dominated by bZIPs, notably BATF; in the loci that lost accessibility with age, motifs for ETS and Zf TFs were enriched, with steady-state expression levels of the TFs again nominating ETS1 and Sp1 as having key roles (Fig. 5B).

The age-associated enrichment for BATF binding motifs was of particular interest because BATF is an essential regulator of T cell differentiation (25), and its overexpression has been shown to lead to T cell exhaustion (25, 38). We found that BATF target genes were significantly more expressed in old compared to young in every T cell subset (naive, CM, EM: p = 2.48 * 10⁻⁷, 6.02 * 10⁻⁵, and 0.0028, respectively, Wilcoxon rank-sum tests; Fig. 5C).

Other classical differentiation-associated TFs, including Eomes, TBX21, RUNX1, and PRDM1, appeared to be specific to differentiation and did not exhibit aging-associated binding motif enrichment. Concordant with this finding, we determined that age-driven expression changes were significantly greater for BATF target genes than for Eomes target genes (Fig. S5E). Additionally, canonical effector genes such as IFNG and GZMB did not increase in accessibility with age (Fig. 1B), suggesting that the mechanisms driving aging only partially resemble those underlying antigen-induced differentiation. Overall, our findings indicate that changes in chromatin accessibility patterns, TF binding motif utilization, and gene expression all support the model that naïve and central memory CD8 T cells become more epigenetically primed for memory and effector cell function with age, while losing naive- and central memory-associated features.

Loss of promoter openness during aging associated with decline in NRF1 activity

Regions of greater chromatin accessibility in young tended to be highly open promoter regions, while enhancers constituted a greater percentage of peaks more open in old (Fig. 6A, Fig. S6A). Compared to all regions in our peak set, the ratio of promoters to enhancers was significantly greater for peaks more open in young (p < 10⁻⁹¹, binomial test) and significantly lesser for peaks more open in old (p < 10⁻⁶⁰, binomial test; Fig. 6B, upper pane), a disparity which was not attributable to technical biases or systematic imbalances in data quality (Figures S6B–E). None of the comparisons between T cell subsets showed a shift from the overall promoter-to-enhancer ratio as great as that of the aging comparisons (Fig. 6B, lower pane). This discrepancy suggests that loss of chromatin openness at promoters is a unique feature of immune aging.

Because we did not observe a similar shift in the chromatin openness changes associated with differentiation, we hypothesized that a molecular driver of the chromatin signatures of aging would be a TF with a binding motif highly enriched in loci closing with age, but not highly enriched in loci closing with differentiation. We found that NRF1 best matched this pattern (Fig. S7A). Among the top twenty most significantly enriched TF binding motifs from the aging and differentiation comparisons, sites for NRF TFs were the most skewed toward promoter regions relative to enhancers (Fig. 7A; NRF in yellow). Enrichment for NRF1 motifs in loci closing with age was maintained even after controlling for the overall promoter-to-enhancer shift, both for the naive (Fig. S7B) and the CM (Fig. S7C).

Footprinting the NRF1 motif, we found greater chromatin openness in young naive compared to old naive around NRF1 sites in loci that close with age (Fig. 7B). Using NucleoATAC (39), we determined that aggregate occupancy of nucleosomes at NRF1 sites was 30% greater in the old than in the young (Fig. S7D) – the greatest shift among the TFs with the top 10 most significantly enriched motifs in the naive aging comparison. These results indicated that NRF1 binding activity is a key feature of the epigenetic landscape of young naïve CD8 T cells. To further investigate this inference, we performed NRF1 ChIP-seq in naive CD8 T cells from young donors. Peaks more open in young than old naive were significantly associated with NRF1 binding (p < 2.2 * 10⁻¹⁶, Fisher’s exact test; Fig 7C). This observation confirmed the specificity of NRF1 binding in maintaining chromatin openness at loci characteristic of youth in the naive subset.

NRF1 is a critical regulator of mitochondrial biogenesis (40), specifically mitochondrial respiratory chain complex (MRC) genes (41, 42). MRC promoters closed with age in naïve CD8 T cells at a greater rate than did other loci (p = 6.32 * 10⁻⁵, Fisher’s exact test), and MRC genes decreased in expression in naïve cells with age, both in our original cohort (p = 1.3 * 10⁻⁴, Wilcoxon rank-sum test; Fig. 7D) and an independent validation set (p = 3.7 * 10⁻⁹, Wilcoxon rank-sum test; Fig 7E). Furthermore, NRF1 binding was significantly enriched both at MRC promoters in general (p = 1.5 * 10⁻¹², Fisher’s exact test; Fig 7F), and at those that exhibited a loss in accessibility in the naive cells with age (p = .012, Fisher’s exact test).

NRF1 gene expression decreased slightly, though not significantly, in the naive cells with age (Fig. 5A). To test whether reducing NRF1 expression induces a corresponding loss in NRF1 binding, we knocked down NRF1 in young naive cells. We detected only a minimal decrease in occupancy at NRF1 binding sites (Fig. 7G), despite a substantial decrease in NRF1 protein levels (Fig. 7H). Accordingly, the knockdown induced a decrease in MRC gene expression (p = .002, Wilcoxon rank-sum test; Fig 7I), but one that was more modest than the age-associated difference.

We investigated the functional impact of the drop in MRC gene expression in naive cells with age by comparing the rates of mitochondrial respiration in young and old naïve CD8 T cells. Both basal and maximal oxygen consumption rates (OCR) were decreased in the old naive relative to the young (p = 0.022 and 0.021, respectively, F-tests; Fig. 7J). Interestingly, mitochondrial mass was increased in the old compared to the young (p = 0.029, Wilcoxon rank-sum test; Fig. 7K), resulting in an even larger difference in OCR when normalized to mitochondrial mass. Overall, our findings suggest that the failure of NRF1 to maintain promoter openness with age may compromise the metabolic state of old naive CD8 T cells, resulting in increased cell loss and an inability to mount a vigorous, energy-demanding clonal expansion in immune response.

Study design

This study was designed to characterize the epigenetic landscapes and molecular drivers associated with human CD8 T cell aging, as well as those associated with maintaining the naive and memory states.

Our experimental design utilized repeated measures, with donors as the blocking variable: for every donor, each of the naive, central memory, and effector memory T cell subsets were subjected to ATAC-seq. Having multiple data points for each donor allowed us to statistically control for interpersonal variability in chromatin openness.

Cells were collected from five young (aged 20–35 years) and five elderly adults (aged 60–85 years). All individuals were healthy, regular platelet donors without a history of autoimmune disease, diabetes mellitus, or cancer.

ATAC-seq analysis and peak annotation

We performed ATAC-seq on purified T cell subsets using standard protocols (see Supplemental Methods). We aligned the data to hg19 using Bowtie 2, and called a total of 64,089 peaks with MACS2 (50). Gene and TSS annotations were drawn from RefSeq and GENCODE. ChromHMM calls for naive CD8 T cells were downloaded from the Roadmap Epigenomics Project.

Differential openness and TF binding motif enrichment calls

To assess differential openness, we applied DESeq2, using a model with a main effect for age, an interaction term for age with subset, and an interaction term for age with donors (blocking on donors in each age cohort). Peaks were considered to be differentially open if the corrected p-value was less than 0.1.

For highly significant peaks (corrected p-value < 0.05), we performed k-means clustering using the median z-score for each T cell subset. The number of clusters was determined with the gap statistic. Gene ontology analysis was performed using GREAT (21).

Transcription factor binding motif enrichment was calculated with HOMER (51). TFs of interest were chosen as those that had both enriched motifs and significant changes in expression. TF target sets were drawn from published results on genome-wide expression studies based on knockdown experiments in human-derived cell lines (25, 32–35).

RNA-seq

RNA-seq libraries were constructed according to standard protocols (see Supplemental Methods). The data were aligned with TopHat2 and differential expression was tested using DESeq2. Gene ontology analysis was performed with topGO.

When running PCA, we used the counts from the 1,000 genes with the highest overall variance in the original cohort. The validation cohort was then projected onto these dimensions. To enable direct comparisons between the two cohorts, within each cohort the median PC 1 value of the young naive was subtracted from each sample’s coordinate.

ChIP-seq

ChIP-seq was performed on naïve T cells using two NRF1-specific antibodies, according to standard protocols (see Supplemental Methods). Reads were aligned to the genome using Bowtie 2 and a total of 3,375 peaks were called with MACS2.

Seahorse Assay

The Seahorse assay was performed according to the manufacturer’s instructions (see Supplemental Methods). Statistical significance was assessed by analyzing each time-point as separate repeated measures ANOVA experiments.