Dynamic chromatin regulatory landscape of human CAR T cell exhaustion.

Dysfunction in T cells limits the efficacy of cancer immunotherapy. We profiled the epigenome, transcriptome, and enhancer connectome of exhaustion-prone GD2-targeting HA-28z chimeric antigen receptor (CAR) T cells and control CD19-targeting CAR T cells, which present less exhaustion-inducing tonic signaling, at multiple points during their ex vivo expansion. We found widespread, dynamic changes in chromatin accessibility and three-dimensional (3D) chromosome conformation preceding changes in gene expression, notably at loci proximal to exhaustion-associated genes such as PDCD1, CTLA4, and HAVCR2, and increased DNA motif access for AP-1 family transcription factors, which are known to promote exhaustion. Although T cell exhaustion has been studied in detail in mice, we find that the regulatory networks of T cell exhaustion differ between species and involve distinct loci of accessible chromatin and cis-regulated target genes in human CAR T cell exhaustion. Deletion of exhaustion-specific candidate enhancers of PDCD1 suppress the expression of PD-1 in an in vitro model of T cell dysfunction and in HA-28z CAR T cells, suggesting enhancer editing as a path forward in improving cancer immunotherapy.

Chimeric antigen receptor (CAR) T cells have proven to be an effective cancer immunotherapy strategy against B cell malignancies (1–6), but T cell exhaustion has been identified as a major roadblock in development of universally effective CAR T cells (6–11). First noted in models of chronic viral infection, T cell exhaustion occurs in the context of persistent antigen stimulation and is associated with limited T cell proliferation, cytokine production, and cell killing in response to further stimulation. T cell exhaustion is now recognized as a distinct epigenetic state, where fully exhausted T cells exhibit chromatin changes and cannot be reversed by checkpoint blockade but can be remodeled by cessation of cell signaling (12–14). Exhaustion in CAR T cells can arise via high tumor burdens or tonic CAR signaling during ex vivo expansion, which severely limits therapeutic efficacy (8, 15).

Much of the current understanding of the genetic regulation underpinning T cell function (16–18), generally, and T cell exhaustion (19–25), specifically, comes from mouse models of the immune system and chronic infection, while genomic studies in dysfunctional human T cells are only very recently emerging (13, 19–21, 26–32). While we previously identified the dysregulation of AP-1 transcription factors (TFs) as the driver of exhaustion in CAR T cells (3), the dynamics of the chromatin state during the progression of T cell exhaustion and the genetic regulatory program underlying the progression remain unclear. Here, we present the chromatin landscape of exhaustion-prone GD2-targeting HA-28z and nonexhausted CD19-targeting CD19-28z CAR T cells throughout the development of the exhaustion phenotype.

Here we identify and differentiate the gene regulatory networks involved in the early development and later maintenance of the exhausted T cell phenotype, across various T cell subtypes found in CAR T cell cultures. As the genetic networks underlying the T cell response to antigen and T cell exhaustion are largely independent (23), we can dissect the uniquely exhaustion-associated network and uncover coregulated modules of accessible chromatin loci, transcription factors, and dysregulated genes. The genomic features we identify may be used to inform targeted genomic studies using CRISPR-Cas9 to restrict expression of exhaustion markers such that the ex vivo maturation of CAR T cells does not induce T cell exhaustion prior to transfusion.

Results

HA-28z CAR T Cells Become Exhausted during In Vitro Culture and Exhibit a Unique Chromatin Signature.

We have previously reported that human HA-28z CAR T cells, targeting the GD2 surface marker, exhibit tonic signaling and lead to an exhausted T cell phenotype following 10 d of transduction, in contrast to the nonexhausted phenotype of human CD19-28z CAR T cells, targeting the CD19 surface marker (3, 33). To elucidate the dynamics of the exhaustion-associated chromatin state, we profiled the chromatin accessibility landscape and transcriptome of human CAR T cells transduced with either the exhaustion-inducing HA-28z CAR or the nonexhausting CD19-28z CAR (Fig. 1 and ).

Fig. 1.

HA-28z CAR T cells exhibit a differentially accessible chromatin profile during the development of exhaustion. Donor T cells were collected and transfected with either the HA-28z or CD19-28z CAR construct. (A) Schematic view of the two CAR constructs. (B) Experimental schematic outlining the time course of stimulation, transfection, and cell collection for ATAC-seq. (C) Differentially accessible chromatin at the PDCD1 locus in naïve CD8+ HA-28z and CD19-28z CAR T cells, shown as normalized ATAC-seq tracks. The element highlighted shows increased, sustained accessibility in HA-28z CAR T cells. (D) PCA of all samples by global chromatin accessibility profile. PC1 (32% variance) separate HA-28z CAR T cells from CD19-28z CAR T cells. (E) K-means clustering of the top 5,000 most variable chromatin accessibility peaks, colored by mean Z score within each cluster. (F) TF motif enrichment within each cluster from E.

Throughout the 14-d in vitro maturation, the HA-28z CAR T cells more strongly expressed markers of exhaustion, including the inhibitory surface receptors PD-1, TIM-3, LAG3, PD-L1, and CD39, and the secretion of activating cytokine IL-2, as compared to the CD19-28z CAR T cells (). These phenotypic differences can be observed for some markers by day 7, were most evident by day 10, and persist to day 14. We collected HA-28z and CD19-28z CAR T cells starting on the day of purification and naïve/central memory sorting from a donor blood draw (day 0) and subsequently on days 7, 10, and 14 after CAR transduction on day 2. On each collection day, cells were sorted into CD4+/naïve, CD4+/central memory, CD8+/naïve, and CD8+/central memory populations. RNA and chromatin were purified from each population, and RNA sequencing (RNA-seq) and an improved assay for transposase-accessible chromatin by sequencing (Omni-ATAC-seq) were performed on each sample, respectively (Fig. 1 and ).

Aggregating the ATAC-seq data reveals loci with differential chromatin accessibility both across time and across CAR type, including loci proximal to known T cell exhaustion marker gene PDCD1 (Fig. 1). We found global differences in chromatin accessibility associated with the CAR type that exhibit dynamic genome-wide organization throughout the 14-d maturation time course, with a notable contraction of differentially accessible loci in CD19-28z CAR T cells by the later time points, most prominently on days 10 and 14 CD19-28z CAR T cells (). The chromatin profiles of the two CAR types exhibit the highest level of differential accessibility at day 7 of the time course, before the greatest differential increase in surface expression of exhaustion markers, such as PD-1, TIM-3, and CD39, and dysfunctional IL-2 secretion, manifests. Although thousands of loci are differentially accessible throughout the time course, many are restricted to specific time points (). Thus, the greatest magnitude of global chromatin and gene regulatory changes precede the strongest expression of the exhaustion phenotype, suggesting a causal role, while other loci that become accessible at later time points may play a role in maintaining the exhausted gene expression profile.

Clustering all ATAC-seq profiles reveals that CD4+, CD8+, naïve-derived, and central memory-derived HA-28z CAR T cells share high correlation among their global chromatin profiles across all time points (). The day 0 (prior to CAR transduction and activation) T cells form a cluster most closely associated with the late (days 10 and 14) naïve CD19-28z CAR T cells. Principal component analysis (PCA) of global Omni-ATAC-seq profiles implicate the CAR type as the primary driver of differences in chromatin accessibility, as the PC1 axis separates the profiles of CAR T cells sharing a time point but differing in CAR type (Fig. 1). Genomic regions enrichment of annotations tool (GREAT) analysis (34) of the combined 500 most positive and 500 most negative contributing accessible loci to PC1 reveal a genetic network in which 9 of the top 20 biological process gene ontology (GO) terms relate to leukocyte differentiation (e.g., “positive regulation of leukocyte differentiation,” “T cell differentiation,” etc.). GREAT analysis of the PC2-associated loci shows enrichment for leukocyte activation and cytokine signaling terms (11 of 20, e.g., “regulation of T cell activation,” “regulation of cytokine production,” etc.). In line with the expectation that both pathways are down-regulated following acute T cell stimulation, PCA suggests a cyclical pattern of global chromatin response to the in vitro maturation in CD19-28z CAR T cells, as their chromatin profiles approach the nontransduced day 0 populations by days 10 and 14. This pattern appears disrupted in exhausted HA-28z CAR T cells, as their PCA trajectory leaves them more distant from the day 0 populations.

Global Patterns of Chromatin Accessibility Highlight the Role of AP-1 in Exhaustion.

Performing hierarchical clustering on the 5,000 ATAC-seq peaks, or accessible loci, with the greatest variance across normalized samples resulted in clusters of peaks differentially accessible in CD19-28z or HA-28z CAR T cells, many of which appear proximal to genes associated with T cell exhaustion and function, such as PDCD1, CTLA4, HAVCR2, TIGIT, LAG3, and TOX (28, 35, 36), while also revealing that differential clusters of peaks associated with CD4+/CD8+ and naïve/central memory phenotypes remain a small fraction of the driving differences between cell populations ( and Dataset S1).

Aggregating the accessibility signal within each cluster, we find 10 clusters of peaks with unique temporal dynamics of enrichment within the T cells across the time course (Fig. 1). Transcription factor motif analysis on the peaks within each cluster showed high enrichment in the CAR type-associated clusters for certain families of transcription factor motifs. Notably, the AP-1 motif is enriched in HA-28z CAR T cell peaks starting at day 7 (clusters 1 to 3), suggesting a role for members of this TF family in the establishment of exhaustion, in accordance with our previous findings that inhibitory AP-1 family transcription factors promote exhaustion in human CAR T cells (3), while RUNX and ETS family motifs are enriched in peaks that appear starting at day 10 or 14 in HA-28z CAR T cells (clusters 4 and 5), suggesting a role in maintenance of exhaustion (Fig. 1).

Integration of Transcriptome and Epigenome Identify Coregulated Regulatory Element-Gene Modules.

We used the PECA2 and EnrichTF methods (37, 38) to integrate our multiomic dataset, incorporating RNA expression, chromatin accessibility, and time to identify coregulated modules of genetic regulatory elements, transcription factors, and target genes. Across the CAR T cells collected at day 7 of ex vivo culture, 6,501 DNA elements were identified as differentially accessible in HA-28z CAR T cells, and 698 genes were identified as differentially expressed in HA-28z CAR T cells over CD19-28z CAR T cells, both meeting the criteria of P < 0.05 and a fold change of ≥1.2. A number of AP-1 family TFs are represented in the list of the 20 TFs with the highest day 7 HA-28z-specific regulatory scores as called by EnrichTF, with NFE2L1, a bZIP TF known to interact with the AP-1 complex at gene promoters (39, 40), and its close homolog NFE2L3 (41), as the highest-scoring TFs ().

Across day 10 CAR T cells, we found 4,054 DNA elements and 458 target genes to be differentially accessible or expressed, respectively, in HA-28z CAR T cells, as compared to CD19-28z CAR T cells. Similarly to day 7 CAR T cells, NFE2L1 and NFE2L3 were the highest-scoring HA-28z-specific regulatory TFs at day 10, with other AP-1 TFs in the highest-scoring list as well ().

We compared the TFs called by EnrichTF as differentially regulatory in HA-28z CAR T cells across time points to identify the TFs that play roles in gene regulation specifically at an early time point (day 7) or later time point (day 10) during the development of exhaustion in HA-28z CAR T cells. TFs such as NFE2L1, FOXP4, and JUND were found to be both significant regulators in the HA-28z CAR T cells and higher-scoring in the day 7 TF analysis compared to the day 10 TF analysis (). Conversely, EGR3, BATF3, and ARID3B were among the TFs found to be significant regulators in the HA-28z CAR T cells and higher scoring in the day 10 TF analysis compared to the day 7 TF analysis ().

The association of transcription factors to their target genes allows us to cluster these associations into regulatory modules by the strength of their associations in each sample using PECA2. The HA-28z CAR T cells, across CD4+/CD8+ and naïve/central memory subtypes, show two distinct and mutually exclusive regulator-gene modules in their significantly coregulated TF/target gene pairs (Fig. 2 and ). The CD19-28z CAR T cells, however, do not show a strong breakdown of their TF/target gene associations into distinct modules (Fig. 2 and ). Known exhaustion marker genes are found across both regulatory modules in HA-28z CAR T cells, including PDCD1 (module 1) and CTLA4, LAG3, JUN, and JUNB (module 2) (, Dataset S2). Transcription factors in the two regulatory modules include AP-1 TFs JUN, JUNB, and ATF3 (module 1), and BATF and BATF3 (module 2) (, Dataset S3). NFE2L1 was found to be the highest-scoring TF at the day 7 time point in module 2, again suggesting its importance in the regulation of the exhaustion genetic program.

Fig. 2.

Integration of RNA-seq and ATAC-seq reveals independent regulatory gene-TF modules. (A) Representative transcription factor by gene regulatory score matrix as calculated by PECA2 for day 7/CD4+/central memory/HA-28z CAR T cells and (B) matching day 7/CD4+/central memory/CD19-28z CAR T cells. (C) The Pearson correlation between the mRNA expression of PD-1 and mRNA expression of the transcription factors BATF and (D) BATF3. (E) Pearson correlation between the mRNA expression of PD-1 and the normalized chromatin accessibility at one locus 76 kb upstream of the PDCD1 transcription start site, identified by PECA2 as regulatory to PD-1. (F) PECA2 inferred gene regulatory network for day 7/CD4+/CM CAR T cells, showing TFs whose log2 fold change between HA-28z and CD19-28z CAR T cells is >1. Node color presents log2 fold change and node size represents enrichment score, or the geometric mean of fold change and −log10(P value) on HA-specific genes and HA-specific regulatory elements.

As an example of PECA2 identifying regulatory element/TF/target gene associations, we found a locus ∼76 kb upstream of the PDCD1 promoter, which contained a BATF/BATF3 consensus motif sequence. BATF mRNA expression correlated with PD-1 mRNA expression (Fig. 2), more closely than BATF3 does (Fig. 2), while the chromatin accessibility of this locus also correlates with the expression of PD-1 (Fig. 2). Together, this suggests a role for BATF regulating expression of PD-1 by binding at the PDCD1-76 kb locus. We can infer a broader gene regulatory network through PECA2 for each CD19-28z/HA-28z matched pair, such as the network that maintains the differentially expressed genes between the exhausted and nonexhausted day 7/CD4+/CM CAR T cells, which highlights a central role for BATF, BATF3, several AP-1 factors, and NFE2L3 (Fig. 2).

Human vs. Mouse T Cell Exhaustion Exhibit Distinct Chromatin Landscapes.

The noncoding genome is under far less stringent constraint over evolutionary time, and DNA regulatory elements can rewire to produce a conserved gene expression outcome even among related mammalian species (42–44). We compared our chromatin accessibility profiles in human T cells to the previously published mouse exhaustion datasets. It has been previously reported that loci proximal to the mouse Pdcd1 gene are differentially accessible under prolonged, chronic stimulation, many of which share sequence homology with loci proximal to PDCD1 in our CD19-28z and HA-28z CAR T cells (Fig. 3) (19). In addition, the genomic features underlying the ATAC-seq peak sequences are quite similar across species and T cell exhaustion state ().

Fig. 3.

Exhausted human and mouse CAR T cells both share regulatory circuits and rely on species-restricted regulatory circuits. (A) Chromatin accessibility at the PDCD1 locus in mouse- and human-derived CAR T cells under nonexhaustive (acute stimulation/CD19-28z CAR) and exhaustive (chronic stimulation/HA-28z CAR) conditions. Arrow highlights mouse and human elements that are accessible in the exhausted state. (B) Proportions of accessible chromatin loci determined to be accessible at conserved sequences between species (“conserved”), accessible only in one species with sequence homology in the other species (“species specific, mapped”), or accessible only in one species without sequence homology in the other species (“species specific, unmapped”). (C) Enrichment of TF binding motifs across all the accessible chromatin loci in mouse- and human-derived CAR T cells, highlighting motifs that are common between species, have higher enrichment and variability in human, and have higher enrichment and variability in mouse. (D) Hierarchical clustering of normalized chromatin accessibility in individual loci across mouse- and human-derived CAR T cells, broken into clusters by the species in which the underlying sequence appears, and in which sample accessibility is enriched and differential between exhausted and nonexhausted cells. (E) The enrichment of TF binding motif sequences in each cluster of accessible loci.

However, the interspecies sequence conservation of the underlying genome varies with genomic features. We found that the majority of peaks at gene promoters were found in both species, while the majority of peaks found at intergenic, intronic, and regions overlapping with the Fantom5 CD4+ T cell enhancer list (45) were largely unique to either human or mouse, both with and without homologous sequences across species (Fig. 3).

To address whether the differences in chromatin profiles associated with T cell exhaustion between species is associated with specific transcription factors, we used transcription factor motif enrichment to identify the motifs enriched in peaks that are either shared across species, such as members of the AP-1 family, or are biased toward accessibility in only one of the species, such as the mouse-biased POU family and the human-biased STAT family (Fig. 3). We then performed differential analysis to identify exhaustion-associated accessible chromatin loci across mouse and human samples by their normalized accessibility and identified clusters representing peaks differentially accessible in exhausted or nonexhausted T cells, with shared homology or not (Fig. 3). The majority of exhaustion-associated peaks were found to be species specific—that is, most loci accessible in exhausted T cells were accessible in either (but not both) mouse or human T cells; most loci accessible in nonexhausted cells were similarly accessible in either (but not both) species; and few loci were accessible in both species in the corresponding cell types (). A small fraction of peaks had reversed accessibility dynamics, with high accessibility in the exhausted T cells of one species and the nonexhausted T cells of the other species (). As we observed in exhausted human CAR T cells (), the thousands of loci associated with exhaustion in mouse T cells demonstrate an early, strong wave of increased accessibility followed by a gradual decrease in accessibility in a return to a more basal state (Fig. 3). Together, this analysis illustrates a high level of divergence in the regulatory programs underpinning T cell exhaustion in mouse and human through the utilization of largely unique enhancer sequences and, in some instances, evolving contrasting function for the same locus, while maintaining a characteristic chromatin signature and dynamic profile over time.

We performed transcription factor enrichment on each cluster of peaks, and notably observed AP-1 TF motifs most significantly enriched in peaks associated with the exhausted HA-28z human CAR T cells (cluster 8) (Fig. 3). While several peak clusters represent loci accessible primarily in one species, the TF motifs enriched in those clusters are often shared between species. For example, clusters 4 and 8 contain peaks accessible in exhausted mouse and human T cells, respectively, yet both have high enrichment for such TF motifs as RARa, STAT4, and SMAD3. Several TF motifs within these clusters do demonstrate species bias, however, such as high enrichment of HNF6 and MEF2B in the human exhaustion-biased cluster 8 and high enrichment of Hoxa9 and ZNF416 motifs in the mouse exhaustion-biased cluster 4. We also note that clusters composed of peaks without homologs in the other species have distinct TF motif enrichment patterns, suggesting significant higher-level divergence in the T cell exhaustion regulation program between mouse and human.

Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis on GREAT annotations of exhaustion-associated clusters of peaks (clusters 2, 4, 6, 8, and 10) reveal largely shared pathways for the accessible chromatin loci in exhausted T cells across species, including apoptosis, T cell receptor signaling, and JAK-STAT signaling, along with several mouse-biased pathways, such as the PD-L1/PD1 pathway, and the human-biased NF-κB signaling pathway (). Gene ontology analysis of the annotations of exhaustion-associated clusters indicates that the human and mouse shared exhaustion-related genomic loci are enriched for an immune cell activation and differentiation signal, while the human-specific exhaustion ATAC signature is notably enriched for regulation of cytokines IL-2 and IL-4 ().

HiChIP Identifies Distal Enhancers Associated with CAR T Cell Exhaustion.

Many DNA regulatory elements are located tens to hundreds of kilobases away from its target genes, and three dimensional (3D) chromosome folding brings regulatory elements and target genes in spatial proximity while excluding intervening genes. To understand the 3D chromosome landscape, we performed chromatin conformation capture coupled with immunoprecipitation and sequencing (HiChIP) directed against histone H3 lysine 27 acetylation (H3K27ac), a histone mark associated with active enhancers, to construct the enhancer connectome (46). We performed H3K27ac HiChIP and identified differential 3D chromatin looping in the day 0 starting T cell population, day 10 exhausted CD8+/naïve HA-28z CAR T cells, and day 10 nonexhausted CD8+/naïve CD19-28z CAR T cells (, Dataset S4). A number of gene promoters associated with T cell exhaustion lie within the top 5,000 most variable 10-kb HiChIP loop anchor windows, including PDCD1, HAVCR2, LAG3, TIGIT, and TOX (Fig. 4).

Fig. 4.

Gene-distal loci enriched in AP-1 TF motif sequences are differentially proximal to exhaustion-regulating genes. (A) The 5,000 most variable 10-kb HiChIP loop anchors across day 0 nontransfected primary T cells, day 10 CD19-28z, and day 10 HA-28z CAR T cells, clustered by interaction Z scores, and labeled with select genes whose transcription start sites lie within the loop anchor. (B) HiChIP interaction matrix of differentially proximal chromatin loci between day 10 CD19-28z and HA-28z CAR T cells at the PDCD1, (C) CTLA4, and (D) HAVCR2 loci. (E) Virtual-4C HiChIP interaction map anchored at the transcription start site and one-dimensional (1D) sequencing coverage for both ATAC-seq and HiChIP at PDCD1, (F) CTLA4, and (G) HAVCR2. (H) The enrichment of TF motif binding sequences within each loop anchor, ranked by enrichment significance, in day 10 CD19-28z and (I) HA-28z CAR T cells.

Interaction maps at genes associated with T cell exhaustion, such as PDCD1, CTLA4, and HAVCR2, demonstrate extensive distal chromosomal contact, often reaching hundreds of kilobases from the gene promoter (Fig. 4 ). The H3K27ac HiChIP signal is substantially stronger for these three loci associated with exhaustion in HA-28z than nonexhausted CD19 CAR T cells (Fig. 4 ). Virtual chromosome conformation capture-on-chip (4C) analysis anchored at gene promoters highlights loci of differential chromosomal proximity that show both differential and similar levels of chromatin accessibility as measured by ATAC-seq (Fig. 4 ). Globally, we find a total of 13,788 promoters accessible across all ATAC-seq samples, with 192 (1.4%) showing both differential ATAC-seq accessibility and HiChIP looping between matching CD19-28z and HA-28z CAR T cell samples, 615 (4.5%) exhibiting differential accessibility only, 1,053 (7.6%) showing differential HiChIP looping only, and 11,928 (86.5%) promoters showing no difference in either data type (). Together, this indicates that 3D chromosome looping provides an additional level of control beyond enhancer accessibility.

Transcription factor motif enrichment of accessible chromatin identified in ATAC-seq peaks is independently enriched in HiChIP loop anchors and identified the AP-1 family motif as specifically enriched in the exhausted HA-28z CAR T cell loop anchors (Fig. 5 ). In contrast, the motifs of the general DNA looping factors CTCF and its homolog BORIS, involved in enhancer-promoter contacts (47, 48), are equally enriched in HiChIP anchor points in HA-28z and CD19-28z CAR T cells (Fig. 5 ). These results suggest that AP-1 family factors likely drive both the chromatin accessibility and higher order 3D genome architecture for T cell exhaustion.

Fig. 5.

Genome editing at putative enhancer loci causes reduced expression of PD-1 in Jurkat cells. (A) Virtual-4C interaction map and 1D ATAC-seq and HiChIP sequencing coverage at the PDCD1 locus in day 10 CD19-28z and HA-28z CAR T cells highlighting proximal (5 kb upstream) and distal (98 kb downstream) targets for CRISPR-Cas9 excision. (B) PD-1 expression 7 d following stimulation with CD3/CD28 beads in unstimulated Jurkat cells, mock-transfected cells, Jurkat cells transfected with nontargeting Cas9-RNPs, and Jurkat cells transfected with Cas9-RNPs targeting the 5-kb upstream or 98-kb downstream regulatory loci. (C) PD-1 expression following 14 d of stimulation. (D) PD-1 expression in CD19-28z and HA-28z CAR T cells following 10 d of ex vivo culture and transfection with Cas9-RNPs targeting the 5-kb upstream regulatory locus for excision. (E) The proportion of PD-1high and PD-1low HA-28z CAR T cells following Cas9-RNP transfection and nontargeting Cas9 transfection, broken down by CD4+/CD8+ T cell type.

Genome Editing of Regulatory Loci Reduces Expression of PD-1.

To support our identification of regulatory loci using multiomic assays, we used CRISPR-Cas9 gene editing in immortalized Jurkat T cells to remove putative regulatory sequences from the genome. We selected two loci near PDCD1, one 5 kb upstream and one 98 kb downstream of the transcriptional start site, each exhibiting differential chromatin accessibility between CD19-28z and HA-28z CAR T cells during the exhaustion time course, as well as 3D proximity with the promoter as reported by HiChIP (Figs. 4 and 5). Both loci have been previously demonstrated to be differentially accessible in the Jurkat T cell line following prolonged in vitro stimulation, a similarly dysfunctional phenotypic state to canonical in vivo T cell exhaustion () (49–51). CRISPR-Cas9 guide RNA pairs were designed to flank the ATAC-seq peak at each locus, removing 700 to 1,000 bp (Fig. 5).

We transfected Jurkat cells with purified guide RNA-Cas9 ribonucleoprotein (RNP) pairs targeting one of the exhaustion-associated peaks, purified RNP incorporating a nontargeting control guide RNA, or a mock transfection without RNP. Following 7 d of stimulation, all samples were consistently PD-1low, at levels similar to unstimulated Jurkat cells (Fig. 5 and ). Following 14 d of stimulation, the mock and nontargeting RNP transfections resulted in cells predominantly PD-1high . In contrast, PD-1 expression was bimodal in cells that received RNPs targeting the exhaustion-associated enhancers, and the majority of cells remained PD-1low (Fig. 5 and ). Bimodality of PD-1 expression is likely due to inefficiency of CRISPR-Cas9 transfection and editing combined with incomplete penetrance of ablating single enhancers.

We then targeted the exhaustion-associated locus 5 kb upstream of PDCD1 for deletion in HA-28z CAR T cells. The CAR T cells were transfected with pairs of guide RNA-Cas9 RNPs concurrent with CAR transfection. Following 10 d of ex vivo culture, a greater proportion of cells were PD-1low, compared to CAR-only HA-28z CAR T cells (Fig. 5). Both CD4+ and CD8+ CAR T cells demonstrated an increased fraction of PD-1low when the 5 kb upstream enhancer was targeted for editing (Fig. 5). Together, these results demonstrate the feasibility of targeting exhaustion-specific enhancers in human CAR T cells to potentially manipulate the exhaustion program.

Discussion

Here we describe the integrated landscape of accessible chromatin, histone acetylation, and long-range chromatin looping in a dynamic model of human CAR T cell exhaustion. These data extend our previously reported transcriptional and epigenetic features that differentiate nonexhausted CAR T cells and CAR T cells that have acquired the functional hallmarks of exhaustion (3). Assaying the developing CAR T cells through a time course of their maturation, we uncovered the epigenetic and transcriptional regulatory features that mark different stages in the acquisition of the exhausted phenotype in HA-28z CAR T cells, as compared to nonexhausting CD19-28z CAR T cells. Exhausted CAR T cells show global alterations in their chromatin accessibility profiles at day 7, which precedes the manifestation of phenotypic and functional hallmarks of exhaustion observed later in the time course. The global patterns of exhaustion-associated ATAC-seq peaks and the underlying TF motifs we identified are very similar across CD4+/CD8+ and naïve/central memory T cell subtypes (Fig. 1 and ), suggesting the mechanisms we describe are general to multiple T cell types.

We found putative enhancer elements and transcription factors such as NFE2L1 associated with the exhausted T cells at the time point preceding significant phenotypic divergence between CD19-28z and HA-28z CAR T cells, indicating their functional role in the development of exhaustion, perhaps setting up a regulatory landscape more conducive to a subsequent transcriptional program associated with the exhausted phenotype. NFE2L1, a bZIP family TF, has been previously shown to exhibit increased expression in activated primary T lymphocytes (52) and association with the expression, protein abundance, and nuclear binding activity of several AP-1 factors (53), including c-Jun, a TF we previously identified as a mediator of dysfunction in exhausted CAR T cells (3). Our previously reported single-cell RNA-seq analysis of exhausted GD2-28z CAR T cells (3) demonstrated that NFE2L1 and its homolog NFE2L3 are more highly expressed in a GD2-28z–specific subset of cells that also express exhaustion markers PD-1, CTLA4, and LAG3 ().

When comparing the exhausted human CAR T cell chromatin state to published mouse models of T cell exhaustion, we identified distinct chromatin signatures marking exhaustion in each species. The majority of accessible chromatin loci we found in exhausted human CAR T cells were not accessible in dysfunctional mouse T cells, either because the locus in mouse was not accessible or there was no homologous sequence in mouse (and vice versa). Even though much of the dysfunctional chromatin signature was unique to each species, the transcription factor motif enrichment and associated-gene pathway enrichment for those signatures were largely shared.

Cancer immunotherapies must find a balance between mounting a strong antitumor immune response and overstimulation of the immune system. Immune checkpoint inhibitors, namely anti–PD-1 and anti–PD-L1 therapeutics, have been associated with cytokine-release syndrome (54), a potentially fatal condition characterized by hyperactivation and massive release of cytokines by the immune system (55). CD19-targeting CAR T cells have themselves been associated with cytokine-release syndrome in patients (1, 56). One challenge lies in tuning the expression levels of key regulators of dysfunction during specific cell states to prevent dysfunction while maintaining robust antitumor activity. A notable example is the Tox gene, found to be a driver of exhaustion in a mouse model of T cell exhaustion, yet deletion of the gene results in poor T cell persistence in tumors (57).

Enhancers are a promising target for genome editing in the search for effective treatment for genetic diseases (58, 59). The benefits of editing enhancer elements, rather than gene deletion or other methods of editing the gene promoter or body, are twofold: enhancers are cell type specific and gene expression remains possible under different cell states under the control of different enhancer elements. We found that deletion of PDCD1-associated enhancers reduced the expression of PD-1 in an in vitro model of T cell exhaustion. The enhancers were identified using a combination of chromatin accessibility profiling, chromosome conformation capture, and transcriptome analysis, highlighting the importance of multiomic profiling in order to capture a fuller picture of the genetic regulatory networks. For example, one exhaustion-biased enhancer of PDCD1 that we identified was 98 kb from the PDCD1 promoter, spanning a region that includes several genes and residing within another gene’s intron. Many current methods of enhancer analysis rely on nearest-gene associations to assign regulatory function, but evidence is mounting that this is insufficient without including orthogonal techniques, such as chromosome conformation.

We were also able to delete the exhaustion-associated enhancer most proximal to the PDCD1 promoter in HA-28z CAR T cells, and we found a marked reduction in PD-1high CAR T cells in the population following ex vivo maturation. This suggests that the epigenetic landscape of CAR T cells may be modulated to improve their efficacy prior to reinfusion into patients. Importantly, the PDCD1 gene itself was not deleted, leaving open the possibility of its expression in nonexhaustion contexts and cell states. More investigation is required to explore the functional consequences and potential of the deleted enhancer, especially in regard to the cells’ antitumor response, their in vivo ability to avoid exhaustion, and the prevention of overstimulation and cytokine-storm syndrome.

Materials and Methods

T Cell Isolation.

Buffy coats from healthy donors were obtained from the Stanford Blood Center. T cells were purified using the RosetteSep Human T Cell Enrichment Kit (Stem Cell Technologies) with the Lymphoprep density gradient medium and SepMate-50 tubes. The 2 × 107 T cells were aliquoted into CryoStor cryopreservation media (Stem Cell Technologies) for long-term storage.

CAR Vector Production.

The 293GP cell line was used to produce murine stem cell virus-based splice-gag vector retroviral vectors encoding CD19-28z and HA-28z CARs, as previously described. At 24 and 48 h following plasmid transfection, the culture media were replaced. Supernatants were collected and replaced at 48 and 72 h following transfection. Cell debris was removed by centrifugation, and the supernatants were frozen at −80 °C.

CAR T Cell Production.

Primary T cells were thawed and cultured in AIM V media + 5% fetal bovine serum (FBS), 10 mM Hepes, 2 mM GlutaMAX, 100 U/mL penicillin, and 10 µg/mL streptomycin (Gibco). T cells were activated the same day using Human T-Expander CD3/CD28 Dynabeads (Gibco) in a 3:1 beads:cells ratio and 100 U/mL recombinant human IL-2 (Peprotech). Nontissue culture-treated 12-well plates were prepared by incubating 1 mL phosphate-buffered saline (PBS) + 25 µg/mL RetroNectin (Takara) overnight at 4 °C. Wells were subsequently washed with PBS and blocked with 2% bovine serum albumin (BSA) for 15 min. A 1-mL thawed retroviral supernatant was added to each well and centrifuged for 2 h at 3,200 rpm at 32 °C. On day 2 following T cell activation, cells were added to the prepared plates and maintained at 0.5 to 1 × 106 cell/mL.

Flow Cytometry Antibodies.

The following antibodies were used (from BioLegend): CD4-APC-Cy7 (clone OKT4), CD8-PerCp-Cy5.5 (clone SK1), TIM-3-BV510 (clone F38-2E2), CD39-FITC or APC-Cy7 (clone A1), CD95-PE (clone DX2), CD3-PacBlue (clone HIT3a), PD-1-BV421 (clone EH12.2H7), and CD3-PE (clone HIT3a); (from eBioscience): PD-1-PE-Cy7 (clone eBio J105), LAG-3-PE (clone 3DS223H), CD45RO-PE-Cy7 (clone UCHL1), and CD45-PerCp-Cy5.5 (clone HI30); (from BD): CD45RA-FITC or BV711 (clone HI100), CCR7-BV421 (clone 150503), CD122-BV510 (clone Mik-β3), CD62L-BV605 (clone DREG-56), CD4-BUV395 (clone SK3), and CD8-BUV805 (clone SK1).

The anti–CD19-28z CAR antibody was provided by B. Jena and L. Cooper, Division of Pediatrics, The University of Texas MD Anderson Cancer Center, Houston, TX. The anti–HA-28z CAR antibody was provided by the National Cancer Institute, Frederick, MD. Each was conjugated with Dylight288 and/or 650 antibody labeling kits (Thermo Fisher).

Cytokine Production.

A total of 1 × 105 CAR+ T cells were cultured for 24 h in 96-well flat-bottom plates in 200 µL media in triplicate. Enzyme-linked immunosorbent assay (BioLegend) for IFNγ and IL-2 was conducted on collected supernatants.

Omni-ATAC-seq.

Omni-ATAC-seq was conducted as previously described. Briefly, 100,000 CAR T cells were sorted into complete media, centrifuged at 4 °C at 500 × g, and resuspended in ATAC-seq resuspension buffer (RSB) (10 mM Tris⋅HCl, 10 mM NaCl, 3 mM MgCl2, 0.1% Nonidet P-40, 0.1% Tween-20, and 0.01% digitonin). Cell suspensions were divided into two replicates. Samples were placed on ice for 3 min and washed with RSB-2 (10 mM Tris⋅HCl, 10 mM NaCl, 3 mM MgCl2, and 0.1% Tween-20). Centrifugation at 500 × g for 10 min at 4 °C made a pellet of nuclei, which was resuspended in 50 µL transposition mix (25 μL 2× TD buffer, 2.5 μL transposase (Illumina), 16.5 μL PBS, 0.5 μL 1% digitonin, 0.5 μL 10% Tween-20, and 5 μL H2O). Transposition was conducted at 37 °C for 30 min in a shaking thermomixer at 1,000 rpm. DNA was purified with the MinElute PCR Purification Kit (Qiagen). Libraries were amplified by PCR using the NEBNext Hi-Fidelity PCR Master Mix and custom primers (IDT) and quantified with the KAPA Library Quantification Kit. Libraries were sequenced on the Illumina NextSeq at the Stanford Functional Genomics Facility with paired-end 75-bp reads.

Adapter sequences were trimmed from the reads and sequences aligned to the hg19 genome with bowtie2. Reads were filtered with Picard tools against mitochondrial sequences, quality scores of <20, and PCR duplicates. Peaks were called using MACS2 on Tn5-corrected insertion sites. A union peak set was compiled by extending peak summits to 500 bp, merging all summits, running bedtools cluster, selecting summits with the highest MACS2 score, and filtering by the ENCODE hg19 blacklist. Peaks were annotated with HOMER (60). TF motifs were found by chromVAR (61) motifmatchr with the chromVARmotifs HOMER set and HOMER. Genome coverage tracks were made by pileups binning every 100 bp with rtracklayer. Differential peak calling was done with DESeq2 (62) using transcription start site-annotated peaks as control loci or consistent “housekeeping peaks.” Clustering was performed using the regularized-log transform values from DESeq2. TF motif enrichment was found using a hypergeometric test on the motifmatchr representation in peak subsets.

RNA-Seq.

Total RNA was collected from 2 × 106 CAR T cells with the RNEasy Plus Mini isolation kit (Qiagen). Library preparation and RNA-seq was performed by BGI America on the BGISEQ-500 as single-end 50-bp reads.

Human and Mouse ATAC-Seq Data Comparison.

Mouse ATAC-seq fastq files were downloaded from GEO (GSM2337202) (19). Mouse ATAC-seq data were aligned and filtered based on the mm9 reference genome and analyzed with the same pipeline and tools as human Omni-ATAC-seq. Peaks within ±3,000 bp were categorized as promoter; peaks that overlap with the Fantom5 (63) CD4+ T cell enhancer list were defined as enhancer. Each ATAC-seq peak was annotated with its nearby genes using GREAT (34) by the basal plus extension default setting. The University of California, Santa Cruz tool liftOver module was used to “lift over” peaks between mm9 and hg19. Bedtools intersect module with default settings was used to define the species-specific peaks and conservative peaks. Peaks that cannot be “lifted over” between mm9 and hg19 were defined as species-specific peaks without homology. Mouse peaks which could be lifted over from mm9 to hg19, but do not overlap with any human peaks, were defined as mouse-specific peaks with homology. Human-specific peaks with homology were defined with the same method. Conserved peaks were peaks that could be lifted over, and the corresponding genomic region is accessible in both human and mouse samples. The MultiBamCov module from bedtools was used to generate peak read count matrices from bam files. Accessibility Z scores of motifs were generated by the chromVar R package using Jaspar motifs, and the top motifs with the highest variability were reported. DESeq2 was used to call differentially accessible peaks between exhausted T cells and nonexhausted T cells at day 7 in human samples, and day 8 in mouse samples, with false discovery rate (FDR) <0.05 and fold change of >2. Motif enrichment analysis based on the differential ATAC-seq peaks was conducted by Homer with FDR <0.05.

HiChIP.

A total of 500,000 cells per replicate were pelleted at 300 × g for 5 min and resuspended in 500 µL 1% methanol-free formaldehyde (Thermo Fisher). After a 10-min room temperature incubation with rotation, cross-linking was quenched with glycine to a final concentration of 125 mM. Cross-linked cells were washed in PBS and pellets stored at −80 °C until subsequent next steps. Pellets were resuspended in 500 µL cold Hi-C lysis buffer (10 mM Tris⋅HCl pH 7.5, 10 mM NaCl, 0.2% Nonidet P-40, 1× protease inhibitor [Roche]) and rotated for 30 min at 4 °C. Samples were centrifuged at 2,500 × g for 5 min at 4 °C to pellet nuclei. Pelleted nuclei were washed in cold Hi-C lysis buffer and resuspended in 100 μL of 0.5% sodium dodecyl sulfate (SDS) at 62 °C for 10 min without rotation. SDS was quenched by adding 285 μL of water and 50 μL of 10% Triton X-100, with rotation at 37 °C for 30 min. Fragmentation was carried out by adding 50 µL NEB buffer 2 and 8 μL of 25 U/μL MboI restriction enzyme (New England Biolabs). Samples were incubated at 37 °C for 2 h with rotation, then 62 °C for 20 min. A total of 52 µL of fill-in master mix was added, consisting of 37.5 µL 0.4 mM biotin-deoxyadenosine triphosphate (Thermo Fisher), 1.5 µL each 10 mM deoxycytidine/-guanosine/-thymidine triphosphate (NEB), and 10 µL 5 U/µL DNA Pol I large Klenow fragment (NEB), then incubated at 37 °C for 1 h. A total of 948 µL ligation master mix was added, consisting of 150 µL 10× NEB T4 DNA ligase buffer with 10 mM ATP (NEB), 125 µL 10% Triton X-100, 3 µL 50 mg/mL BSA, 10 µL 400 U/µL T4 DNA ligase (NEB), and 660 µL water, then incubated at room temperature for 4 h. Nuclei were then pelleted at 2,500 × g for 5 min and supernatant was discarded and resuspended in 880 µL nuclear lysis buffer (50 mM Tris⋅HCl, 10 mM ethylenediaminetetraacetic acid (EDTA), 1% SDS, 1× protease inhibitor). DNA was sonicated in a Covaris E220 and Covaris millitube with the following parameters: fill level 10, duty cycle 5, PIP 140, cycles/burst 200, and time 4 min. Samples were centrifuged for 15 min at 4 °C at 16,000 × g. A total of 2× volume of ChIP dilution buffer (0.01% SDS, 1.1% Triton X-100, 1.2 mM EDTA, 16.7 mM Tris⋅HCl pH 7.5, 167 mM NaCl) was added, and each total sample was split across two 1.5-mL tubes. A total of 21 µL Protein A beads (Thermo Fisher), for our total of 3.5 million cells, were washed in ChIP dilution buffer and resuspended in 50 µL ChIP dilution buffer per tube (100 µL total). Beads were added to sample tubes and incubated at 4 °C for 1 h with rotation. Samples were positioned on a magnet, and the supernatant transferred to new tubes. A total of 2.625 µg anti-H3K27ac antibody (Abcam) was diluted 10× in ChIP dilution buffer (for 3.5 million cells total), distributed equally across tubes, and incubated overnight at 4 °C with rotation. Samples were put on a magnet and supernatants transferred to new tubes. A total of 21 µL Protein A beads were washed in ChIP dilution buffer, resuspended in 50 µL ChIP dilution buffer per tube (100 µL total), added to each sample, and incubated for 2 h at 4 °C with rotation. Samples were placed on a magnet and the beads washed three times each with low-salt wash buffer (0.1% SDS, 1% Triton X-100, 2 mM EDTA, 20 mM Tris⋅HCl pH 7.5, 150 mM NaCl), high-salt wash buffer (0.1% SDS, 1% Triton X-100, 2 mM EDTA, 20 mM Tris⋅HCl pH 7.5, 500 mM NaCl), and LiCl wash buffer (10 mM Tris⋅HCl pH 7.5, 250 mM LiCl, 1% Nonidet P-40, 1% sodium deoxycholate, 1 mM EDTA). Beads were resuspended in 100 μL DNA elution buffer (50 mM sodium bicarbonate pH 8.0, 1% SDS) and incubated at room temperature for 10 min with rotation, then at 37 °C for 3 min with shaking. Samples were placed on a magnet and supernatants transferred to new tubes. Elution was repeated with another 100 µL elution buffer, transferring to the same tube. A total of 10 µL Proteinase K (Thermo Fisher) was added to each sample and incubated at 55 °C for 45 min with shaking, followed by 65 °C for 1.5 h with shaking. Samples were purified with Zymo Research DNA Clean and Concentrator columns and eluted in 10 μL of water. A total of 5 μL streptavidin C-1 beads (Thermo Fisher) were washed with Tween wash buffer (5 mM Tris⋅HCl pH 7.5, 0.5 mM EDTA, 1 M NaCl, 0.05% Tween-20) and resuspended in 10 μL 2× biotin binding buffer (10 mM Tris⋅HCl pH 7.5, 1 mM EDTA, 2M NaCl). Beads were added to samples and incubated at room temperature for 15 min with rotation. Samples were placed on a magnet and supernatants discarded. Beads were washed twice in 500 µL Tween wash buffer with 2 min shaking at 55 °C. Beads were washed in 100 μL 1× TD buffer (2× TD buffer: 20 mM Tris⋅HCl pH 7.5, 10 mM magnesium chloride, 20% dimethylformamide). Beads were resuspended in 25 μL 2× TD buffer. A total of 2.5 μL Tn5 transposase (Illumina) was added for each 50 ng of DNA in each sample, and water was added up to 50 µL. Samples were incubated at 55 °C for 10 min with shaking. Samples were placed on a magnet and supernatants removed. Beads were resuspended in 50 mM EDTA and incubated for 30 min at 50 °C. Samples were placed on a magnet and supernatants removed. Samples were washed twice more with 50 mM EDTA at 50 °C for 3 min. Samples were washed twice with Tween wash buffer at 55 °C for 2 min. Samples were washed in 10 mM Tris. Beads were resuspended in 50 µL PCR master mix (25 µL 2× Phusion HF buffer (NEB), 1 µL 12.5 µM Nextera Ad1 primer (Illumina), 1 µL 12.5 µM Nextera Ad2 barcoded primer (Illumina), 23 µL water. PCR was run as follows: 72 °C for 5 min, 98 °C for 1 min, and eight cycles of 98 °C for 15 s, 63 °C for 30 s, and 72 °C for 1 min. Amplified libraries were placed on a magnet and transferred to new tubes. Libraries were Zymo column purified and eluted together into the same 10 µL water. Libraries were polyacrylamide gel electrophoresis purified and size selected for 300 to 700 bp, then column cleaned into 10 µL water. Libraries were sequenced on an Illumina NextSeq at the Stanford Functional Genomics Facility with paired-end 75-bp reads.

HiChIP data were analyzed as previously described (47). Paired-end reads were aligned to hg19 with the Hi-C Pro pipeline, interactions called with Juicer HICCUPS (64) and Fit-HiC (65), and interaction maps were made with Juicebox (66).

Cell Lines.

The Jurkat T cell line was obtained from ATCC (clone E6-1) and cultured in RPMI + 10% FBS, 10 mM Hepes, 100 U/mL penicillin, and 100 µg/mL streptomycin (Gibco). Cells were stimulated with Human T-Expander CD3/CD28 Dynabeads (Gibco). At sorting, propidium iodide was used to stain and exclude dead and dying cells.

CRISPR-Cas9 Editing.

Cas9-RNP complexes were delivered into Jurkat cells using the Alt-R CRISPR-Cas9 system (IDT) and the Neon Transfection System (Thermo Fisher). Guide-RNA sequences were designed using the Broad Institute’s online sgRNA design tool and CRISPETa paired guide designer (67). RNA oligonucleotides were synthesized by IDT:

PDCD1-5 kb-upstream-1–5′: /AltR1/rGrUrArUrUrGrCrCrArArGrGrGrCrArCrCrCrGrGrGrUrUrUrUrArGrArGrCrUrArUrGrCrU/AltR2/

PDCD1-5 kb-upstream-1–3′: /AltR1/rUrCrGrUrUrGrGrGrCrUrGrUrUrUrGrUrUrArArGrGrUrUrUrUrArGrArGrCrUrArUrGrCrU/AltR2/

PDCD1-5 kb-upstream-2–5′: /AltR1/rArUrGrCrGrUrArCrUrGrCrArCrCrUrUrCrCrCrArGrUrUrUrUrArGrArGrCrUrArUrGrCrU/AltR2/

PDCD1-5 kb-upstream-2–3′: /AltR1/rArArUrCrGrUrUrGrGrGrCrUrGrUrUrUrGrUrUrArGrUrUrUrUrArGrArGrCrUrArUrGrCrU/AltR2/

PDCD1-98 kb-downstream-1–5′: /AltR1/rArUrCrUrGrUrCrUrGrCrUrArArGrGrUrCrCrArCrGrUrUrUrUrArGrArGrCrUrArUrGrCrU/AltR2/

PDCD1-98 kb-downstream-1–3′: /AltR1/rArUrGrArArCrArUrArUrUrCrUrGrUrGrGrArUrGrGrUrUrUrUrArGrArGrCrUrArUrGrCrU/AltR2/

PDCD1-98 kb-downstream-2–5′: /AltR1/rArUrCrUrGrArCrGrGrCrCrGrArUrCrArCrArCrArGrUrUrUrUrArGrArGrCrUrArUrGrCrU/AltR2/

PDCD1-98 kb-downstream-2–3′: /AltR1/rArCrGrArUrCrArCrUrCrUrArCrUrGrCrCrCrUrCrGrUrUrUrUrArGrArGrCrUrArUrGrCrU/AltR2/.

Pairs of RNPs were assembled and transfected into Jurkat cells according to manufacturer protocol. Briefly, 220 pmol of each guide-RNA in a pair and 440 pmol tracrRNA were combined in 10 µL Tris-EDTA buffer, heated to 95° for 5 min, then placed at room temperature to cool. A total of 10 µL of 36 µM Cas9 enzyme was added to each crRNA-tracrRNA pair and incubated at room temperature for 10 min. For each transfection, 5 × 106 Jurkat cells were collected, centrifuged, washed in PBS, and resuspended in 80 µL buffer R. The 80-µL cell suspension was mixed with the 20 µL RNP complex and electroporated in the Neon Transfection System at 1,600 V, 10 ms pulse width, and three pulses. Cells were transferred to prewarmed media in six-well plates. At 24 h, cells were stimulated with Human T-Expander CD3/CD28 Dynabeads (Gibco) in a 3:1 beads:cells ratio and 100 U/mL recombinant human IL-2 (Peprotech).