This research complies with all relevant ethical regulations.
Approval for vertebrate animal studies was provided by the La Jolla Institute for Immunology (LJI) Animal Care Committee (protocol number AP00001025).
Approval for human studies was provided by the LJI Institutional Review Board (protocol number SCRO_002). Informed consent was obtained from all human blood donors. Donors were compensated per LJI policy.
Mice
C57BL6/J, B6.SJL-PtprcaPepcb/BoyJ (CD45.1) and C57BL/6-Tg(TcraTcrb)1100Mjb/J (OT-I) mice were obtained from The Jackson Laboratory. CD45.1+ OT-I mice were obtained by crossbreeding. Male mice were used for the experiments on antitumor effects in vivo, and both male and female mice were used for the in vitro studies. Six-week-old mice were purchased from The Jackson Laboratory to be used as recipient mice, and rested for at least 1 week after delivery before they were used for the experiments. Mice were age matched and between 7 and 12 weeks old when they were used for the experiments. All mice were bred and/or managed in the animal facility at the LJI. All experiments were performed in compliance with study protocol number AP00001025 approved by the LJI Animal Care Committee. The approved protocol specified a maximum tumor size of 1.77 cm3 or greater for 3 d without signs of regression. This maximum tumor size was not exceeded.
Human peripheral CD8+ T cell isolation
Whole blood samples from healthy subjects were collected by a staff phlebotomist in the Clinical Studies Core at the LJI, and peripheral blood mononuclear cells were isolated using Lymphoprep (STEMCELL Technologies) according to the manufacturer’s protocol. Peripheral CD8+ T cells were negatively isolated using a human CD8+ T cell isolation kit (STEMCELL Technologies) following the manufacturer’s protocol.
Construction of retroviral and lentiviral vectors
CAR expression plasmid
The sequence of the retroviral vector (MSCV-myc-CAR-2A-Thy1.1) encoding the Myc epitope-tagged CAR has been reported previously49,50. It contains the hCD19 single-chain variable fragment49 and the murine CD3ζ and CD28 sequences. The CAR complementary DNA was cloned into an MSCV-puro murine retroviral vector in place of PGK-puro.
Human CD19 retroviral expression plasmid
A PCR-amplified DNA fragment encoding hCD19 was cloned into an MSCV-puro (Clontech) murine retroviral vector as described in previous papers5,6.
Retroviral vectors (MSCV-bZIP-IRES-Thy1.1 and MSCV-bZIP-IRES-eGFP)
To generate pMIG-Batf, the Batf coding sequence was amplified from pMSCV-Batf-IRES-Thy1.1 (C.-W.J.L.; unpublished), derived from pcDNA3.1-Batf (Addgene; 34575), and cloned into pMSCV-IRES-eGFP (Addgene; 27490). DNA fragments encoding Jun, Maff and the Batf HKE mutant were PCR amplified or synthesized as gBlocks (Integrated DNA Technologies) and cloned into MSCV-IRES-eGFP (Addgene plasmid 27490), kindly provided by W. S. Pear (Univ. Pennsylvania). pMIG-IRF4 was purchased from Addgene (58987).
Lentiviral vectors (pTRPE-19.28z-P2A-NGFR, pTRPE-IRES-eGFP and pTRPE-BATF-IRES-eGFP)
The plasmid pTRPE-19.28z, which contains the hCD19 single-chain variable fragment and the human CD3ζ and CD28 sequences, was kindly provided by A. D. Posey Jr (University of Pennsylvania). A fragment containing the P2A and nerve growth factor receptor (NGFR) sequences was PCR amplified and cloned into the pTRPE-19.28z vector to yield pTRPE-19.28z-P2A-NGFR. A fragment containing the internal ribosome entry site (IRES) and enhanced green fluorescent protein (eGFP) sequences was PCR amplified and cloned into the pTRPE-19.28z vector in place of 19.28z to yield pTRPE-IRES-eGFP. DNA fragments encoding human BATF were synthesized as gBlocks (Integrated DNA Technologies) and cloned into pTRPE-IRES-eGFP.
Cloning of NFAT–AP-1 reporter plasmids
A retroviral reporter plasmid containing six tandem NFAT–AP-1 sites driving GFP expression on a self-inactivating retroviral backbone was kindly provided by H. Spits51. Mouse Thy1.1 was cloned into this plasmid, in place of the GFP reporter, using Gibson Assembly. The mouse genes for Jun, Maff, Batf, Batf3, Jund, Fosl2 and Nfil3 were synthesized as gBlocks (Integrated DNA Technologies) and cloned downstream of Thy1.1 with a P2A linker in between using Gibson Assembly.
Cell lines
The B16F0 mouse melanoma cell line was purchased from the American Type Culture Collection. The B16F0-hCD19 cell line was generated by transduction with amphotropic virus encoding hCD19, followed by sorting for cells expressing high levels of hCD19. The B16F10-OVA mouse melanoma cell line was kindly provided by S. Schoenberger (LJI). The Platinum-E Retroviral Packaging Cell Line, Ecotropic (Plat-E) cell line was purchased from Cell Biolabs. All tumor cell lines were tested frequently to be sure they were negative for Mycoplasma contamination and were used at passage 4 after thawing from stock.
Transfections
A total of 3 × 106 Plat-E cells were seeded in 10 cm dishes in media (Dulbecco’s modified Eagle’s medium with 10% fetal bovine serum (FBS), 1% l-glutamine and 1% penicillin/streptomycin) the day before transfection, and the medium was changed just before transfection. For retroviral transduction, we used a mixture of 10 μg retroviral plasmid + 3.4 μg pCL-Eco packaging vectors or PCL10A1. For lentiviral transduction, the mixture contained 10 μg lentiviral plasmid + 7.5 μg Gag pol + 5 μg Rev + 2.5 μg VSV-G packaging vectors. The plasmid mixtures were incubated with 40 μl TransIT-LT1 Transfection Reagent (Mirus Bio) at ~22 °C for 20 min in 1.5 ml Opti-MEM media and then added to the Plat-E cells, after which the cells were incubated at 37 °C in a 10% CO2 incubator for 30–40 h. The supernatant was filtered through a 40-μm filter before being used for transduction of CD8+ T cells.
Tumor experiments
Preparation of B16F0-hCD19 and B16F10-OVA melanoma cells for tumor inoculation
Tumor cells (B16F0-hCD19 or B16F10-OVA) were thawed and cultured in Dulbecco’s modified Eagle’s medium with 10% FBS, 1% l-glutamine and 1% penicillin/streptomycin at 37 °C in a 5% CO2 incubator, and were split and passaged at days 1, 3 and 5 after thawing before inoculation. At the time of tumor inoculation (day 0), cells were trypsinized and resuspended in PBS solution, then injected subcutaneously into 7- to 12-week-old C57BL/6J mice.
Generation and transfer of CAR T cells
Splenic CD8+ T cells from C57BL/6, B6.SJL-PtprcaPepcb/BoyJ, C57BL/6-Tg(TcrαTcrβ)1100Mjb/J or CD45.1xOT-I mice were isolated by negative selection using a CD8 isolation kit (Invitrogen or STEMCELL Technologies), activated with 1 µg ml−1 anti-CD3 and anti-CD28 for 1 d, then removed from the plates and retrovirally transduced using 15 µg ml−1 polybrene at 37 °C, followed by centrifugation at 2,000g for 1–2 h. After transduction, cells were cultured in house-made T cell medium containing 100 U ml−1 human IL-2. A second transduction was performed the next day using the same protocol, after which cells were cultured in T cell media containing 100 U ml−1 human IL-2 for 3 d. On the day of adoptive transfer, the cells were analyzed by flow cytometry to check the transduction efficiency (typically 90% for single retroviral transduction and 80% for double retroviral transductions), and cell counts were obtained using the Accuri flow cytometer. Cells were washed with PBS and resuspended in PBS before adoptive transfer into recipient mice.
Assessing antitumor responses
On day 0, 7- to 12-week-old C57BL/6J mice were injected subcutaneously with 1 × 105 B16F0-hCD19 cells or 2.5 × 105 B16F10-OVA cells. When the tumors were palpable, tumor measurements were recorded with a caliper three to four times a week and the tumor size was calculated in mm2 (length × width). On day 7, 3 × 106 CAR T cells or 1 × 106 OT-I T cells were adoptively transferred into tumor-bearing mice. For all survival experiments, tumor growth was monitored until an experimental endpoint of day 100 after tumor inoculation or until the Institutional Animal Care and Use Committee-approved endpoint of a maximum tumor size measurement exceeding a diameter of 1.77 cm3 for more than 3 d without signs of regression. If mice were pale, had scars or ulcerations or adopted a hunched position, or if their body temperature was low, we euthanized the mice under the guidance of the staff of the Department of Laboratory Animal Care (DLAC) at the LJI. In most cases, tumor sizes were measured in a blinded manner by DLAC staff, except during the holiday season or when the institute was under restricted access due to the COVID-19 shutdown.
Harvesting TILs
On day 0, 7- to 12-week-old C57BL/6J mice were injected subcutaneously with 1 × 105 B16F0-hCD19 or 2.5 × 105 B16F10-OVA cells in PBS. When tumors were palpable, tumor measurements were recorded with a caliper three to four times a week and the tumor size was calculated in mm2 (length × width). On day 12, 1.5 × 106 CAR T cells or 1 × 106 OT-I T cells were adoptively transferred into tumor-bearing mice. On day 20, tumors were collected from the mice and placed into C tubes (Miltenyi Biotec) containing RPMI 1640 with 10% FBS and collagenase D (1 mg ml−1; Roche), hyaluronidase (30 U ml−1; Sigma–Aldrich) and DNase I (100 µg ml−1; Sigma–Aldrich). Tumors were dissociated using the gentleMACS dissociator (Miltenyi Biotec), incubated with shaking at 2,000 r.p.m. for 60 min at 37 °C, filtered through a 70-μm filter and spun down. Lymphocytes were separated using lymphocyte separation medium (MP Biomedicals; 0850494).
NFAT–AP-1 reporter assays
Primary mouse CD8+ T cells were isolated from the spleens of C57BL/6J mice (The Jackson Laboratory; 000664) by negative selection (EasySep; 19853). Up to 5 × 106 freshly isolated CD8+ cells were activated with plate-bound anti-CD3 (145-2C11) and anti-CD28 (37.51) at a final 1 µg ml−1 in T cell media in a six-well plate. After 24 h, cells were transduced with retroviral supernatant at 32 °C for 2 h at 2,000g with 8 µg ml−1 polybrene. After transduction, cells were cultured in T cell media containing 100 U ml−1 IL-2. On day 2, the same transduction was performed. On day 3, cells were surface stained for live CD8+ Thy1.1+ cells as a measure of reporter activity.
Flow cytometry analysis
BD Fortessa, BD LSR III or BD Celesta flow cytometers were used for cell analysis. Cells were resuspended in FACS buffer (PBS, 1% FBS and 2.5 mM ethylenediaminetetraacetic acid (EDTA)) and filtered using a 70-μm filter before running the flow cytometer. Fluorochrome-conjugated antibodies were purchased from BD Biosciences, Thermo Fisher Scientific, Miltenyi Biotec and BioLegend. For surface staining, cells were stained with ~1:100–1:200 dilution of antibodies in FACS buffer (PBS + 1% FBS and 2.5 mM EDTA) for 15 min with Fc block (BioLegend). For cytokine staining, cells were activated with 10 nM phorbol 12-myristate 13-acetate, 500 nM ionomycin and 1 μg ml−1 Golgi plug and/or Golgi Stop in T Cell Media at 37 °C in a 10% CO2 incubator for 4 h. After stimulation, cells were stained for surface markers and resuspended with Fix/Perm (BD Biosciences) buffer for 20 min, washed with FACS buffer twice and stained for cytokines at a final concentration of 1:200 in 1× BD Perm/Wash buffer. For the detection of transcription factors, cells were stained for surface markers first, after which the Foxp3/transcriptional staining kit was used according to the manufacturer’s protocol. All transcription factor antibodies were used at 1:200 dilution. All flow data were analyzed with FlowJo (version 10.6.2).
Mass cytometry by time of flight (CyTOF) analysis
On day 0, 7- to 12-week-old C57BL/6J mice were injected subcutaneously with 1 × 105 B16F0-hCD19 cells. When the tumors were palpable, tumor measurements were recorded with a caliper three to four times a week and the tumor size was calculated in mm2 (length × width). On day 12, 1.5 × 106 CAR T cells were adoptively transferred into tumor-bearing mice. On day 20, tumors were collected from the mice and placed into C tubes (Miltenyi Biotec) containing RPMI 1640 with 10% FBS and collagenase D (1 mg ml−1; Roche), hyaluronidase (30 U ml−1; Sigma–Aldrich) and DNase I (100 µg ml−1; Sigma–Aldrich). The tumors were dissociated using the gentleMACS dissociator (Miltenyi Biotec), incubated with shaking at 2,000 r.p.m. for 60 min at 37 °C, filtered through a 70-μm filter and spun down. Lymphocytes were separated using lymphocyte separation medium (MP Biomedicals; 0850494) and sorted by flow cytometry based on forward versus side scatter gating to obtain highly purified lymphocytes. After sorting, the lymphocytes were rested in T cell media for 4 h. Cells were washed with PBS, centrifuged at 400g for 5 min and the supernatant was discarded by aspiration. Cells were resuspended in PBS with Cell-ID Cisplatin (5 μM), incubated at ~22 °C for 5 min and washed with MACS staining buffer (2 mM EDTA and 2% FBS in PBS) using 5× the volume of the cell suspension. Cells were stained with a cocktail of antibodies to surface proteins with Fc blocking for 15 min at ~22 °C, washed with MACS staining buffer, then fixed and permeabilized using a FoxP3 staining buffer kit (eBioscience) and stained for 1 h at ~22 °C with a cocktail of antibodies to intracellular proteins. Cells were washed twice with Perm/Wash buffer, fixed with 1.6% paraformaldehyde for 10 min at ~22 °C and washed twice again with Perm/Wash buffer. Cells were stained with Cell-ID Intercalator-Ir in Fix/Perm buffer overnight at 4 °C before analysis of the sample using a CyTOF mass spectrometer. All CyTOF data were analyzed with FlowJo (version 10.6.2) or the OMIQ.ai analysis platform.
Cell sorting
Cell sorting was performed by the LJI flow cytometry core using FACSAria I, FACSAria II or FACSAria Fusion (BD Biosciences) flow cytometers. For transcriptional profiling using Smart-seq, 10,000 cells were sorted from the Live/Dead dye-negative CD8+ Thy1.1+ GFP+ population of the isolated TILs or cultured CD8+ T cells. The cells were resuspended in FACS buffer and filtered with a 70-μm filter before sorting. For ATAC-seq, 50,000 live cells were sorted using the same procedure as for Smart-seq. Cells were sorted into 1.5 ml microfuge tubes containing 500 μl 50% FBS. The sorted cells were washed with cold PBS twice before further procedures.
Antibodies
The following antibodies were used: BUV395 rat anti-mouse CD8α (clone 53-6.7; BD Biosciences; 563786); and BV711 anti-rat CD90/mouseCD90.1 (Thy1.1) (clone OX-7; BioLegend 202539).
Primary cell culture
Splenic CD8+ T cells from C57BL/6 mice were isolated using a Dynabeads Untouched Mouse CD8 Cells Kit (Invitrogen) or EasySep Mouse CD8+ T Cell Isolation Kit (STEMCELL Technologies) following the manufacturer’s protocols, following which 3 × 106 CD8+ T cells per well were stimulated with 1 μg ml−1 anti-CD3 and anti-CD28 in T cell media in a six-well plate for 1 d, then removed from the plates and retrovirally transduced using 15 μg ml−1 polybrene at 37 °C followed by centrifugation at 2,000g for 1 h. After transduction, cells were cultured in house-made T cell media containing 100 U ml−1 human IL-2. A second transduction was performed the next day using the same protocol, after which the cells were cultured in T cell media with 100 U ml−1 human IL-2 for 3 d.
Human CAR T cell experiments
Human CD8+ T cells were stimulated with a Dynabeads Human T-Activator CD3/CD28 (Gibco) in X-Vivo (Lonza) medium. After 2 d, Dynabeads were removed from the cells and the cells were lentivirally transduced using retronectin-coated plates (20 μg ml−1) at 32 °C followed by centrifugation at 2,000g for 2 h. Cells were expanded for 2 d in X-Vivo medium with 500 U ml–1 human IL-2. Human CAR T cells were enriched by positive selection for NGFR using MACS columns and beads (Miltenyi Biotec).
In vitro cytotoxicity assay
CAR T cells were labeled with CellTrace Violet dye (Invitrogen) and cocultured with NALM6 tumor cells for 5 h. The percentage cytotoxicity was calculated as 1 − (R5/R0)) × 100, where R5 = (target cells (percentage of total) at 5 h)/(effector cells (percentage of total) at 5 h) and R0 = (target cells (percentage of total) at 0 h)/(effector cells (percentage of total) at 0 h).
In vitro proliferation assay
CellTrace Violet-labeled CAR T cells were cultured in X-Vivo media with 500 U ml−1 human IL-2 for 4 d.
ChIP-seq library preparation
pMIG- or BATF-transduced CD8+ T cells (1 × 106 cells per ml in culture media) were fixed with 1% formaldehyde at ~22 °C for 10 min with nutation. To quench the fixation, 0.5 ml 2.5 M glycine was added per 10 ml, then the cells were incubated on ice for 5 min and washed twice with cold PBS. Fixed cells were transferred to low-binding tubes with 1 ml cold PBS and spun down at 2,000 r.p.m. at 4 °C for 10 min. Cells were pelleted, snap frozen with liquid nitrogen and stored at −80 °C until further processing. To isolate nuclei, cell pellets were thawed on ice and the pellets were resuspended in 1 ml Bioruptor lysis buffer (50 mM HEPES (pH 7.5), 150 mM NaCl, 1 mM EDTA, 10% glycerol, 0.5% NP-40 and 0.25% Triton X-100) and incubated for 10 min at 4 °C with nutation. After centrifugation at 1,700g and 4 °C for 5 min, the resulting nuclear pellets were washed twice with washing buffer (10 mM Tris-HCl (pH 8.0), 200 mM NaCl, 1 mM EDTA and 0.5 mM EGTA). The pellets were resuspended in 100 μl shearing buffer (10 mM Tris-HCl (pH 8.0), 1 mM EDTA and 1% sodium dodecyl sulfate (SDS)) and sonicated using a Bioruptor in 1.5 ml bioruptor tubes (ten cycles; 30 s on and 30 s off). After sonication, the supernatants were transferred to 1.5 ml low-binding tubes, and insoluble debris were removed by centrifugation at 20,000g. Pellets were resuspended in 100 µl shearing buffer, then nine volumes of conversion buffer (10 mM Tris-HCl (pH 7.5), 255 mM NaCl, 1 mM EDTA, 0.55 mM EGTA. 0.11% sodium deoxycholate and 0.11 % Triton X-100) were added. Chromatin was precleared with washed protein A and protein G Dynabeads for 1 h, and the chromatin concentration was measured by Qubit. A total of 5% of the chromatin was saved as input and the chromatin was incubated with anti-BATF (Brookwood Biomedical) or anti-IRF4 (clone D9P5H; Cell Signaling Technology) antibodies and protein A and protein G Dynabeads overnight at 4 °C with rotation. The following day, bead-bound chromatin was washed twice with RIPA buffer (50 mM Tris-HCl (pH 8.0), 150 mM NaCl, 1 mM EDTA, 1% NP-40, 0.1% SDS and 0.5% sodium deoxycholate) and then with high salt buffer (50 mM Tris-HCl (pH 8.0), 500 mM NaCl, 1 mM EDTA, 1% NP-40 and 0.1% SDS), LiCl buffer (50 mM Tris-HCl (pH 8.0), 250 mM LiCl, 1 mM EDTA, 1% NP-40 and 1% sodium deoxycholate) and TE buffer (10 mM Tris-HCl (pH 8.0) and 1 mM EDTA). Chromatin was eluted with 100 μl elution buffer (100 mM NaHCO3, 1% SDS and 1 mg ml−1 RNase A) twice for 30 min at 37 °C using a 1,000 r.p.m. shaking heat block. Some 5 μl proteinase K (20 mg ml−1; Ambion) and 8 μl 5 M NaCl were added to the eluted DNA and the samples were incubated at 65 °C with shaking (1,200 r.p.m.) for de-crosslinking. DNA was purified with a Zymo ChIP DNA Clean & Concentrator (Zymo Research). Libraries were prepared using NEBNext Ultra II Library Prep kits (NEB) following the manufacturer’s instructions, then sequenced using an Illumina NovaSeq 6000 sequencer (50-bp paired-end reads).
ATAC-seq and RNA-seq library preparation
ATAC-seq libraries were prepared following the Omni-ATAC protocol with minor modification52. Some 50,000 cells were collected by sorting and washed twice with cold PBS at 600g for 5 min. Cell pellets were resuspended in 50 μl ATAC lysis buffer (10 mM Tris-HCl (pH 7.4), 10 mM NaCl, 3 mM MgCl2, 0.1% NP-40, 0.1% Tween 20 and 0.01% digitonin) and incubated on ice for 3 min, after which 1 ml washing buffer (10 mM Tris-HCl (pH 7.4), 10 mM NaCl, 3 mM MgCl2 and 0.1% Tween 20) was added and the cells were spun down at 1,000g for 10 min at 4 °C. The supernatant was removed carefully, and the nuclei were resuspended in 50 μl transposition mix (25 μl TD buffer (20 mM Tris-HCl (pH 7.6), 10 mM MgCl2 and 20% dimethylformamide), 2.5 μl 2 μM transposase, 16.5 μl PBS, 0.5 μl 1% digitonin, 0.5 μl 10% Tween 20 and 5 μl water) and incubated at 37 °C for 30 min. DNA was purified with a Qiagen MinElute Kit (Qiagen). Libraries were amplified with KAPA HiFi HotStart Real-Time PCR master mix, and sequenced on an Illumina NovaSeq 6000 sequencer (50-bp paired-end reads). RNA-seq libraries were prepared following the Smart-Seq2 protocol53 modification. Total RNA was extracted from 10,000 sorted cells using the RNeasy Plus Micro kit (Qiagen) and following the Smart-seq2 protocol as described. Libraries were prepared using the Nextera XT Library Prep kit (Illumina) and sequenced on an Illumina NovaSeq 6000 sequencer (50-bp paired-end reads).
ATAC-seq analysis
Genome browser tracks
Paired raw reads were aligned to the Mus musculus genome (mm10) using Bowtie (version 1.0.0 and -X 2000 -m 1 --best --strata -tryhard -S --fr)54. Unmapped reads were trimmed to remove adapter sequences and clipped by one base pair (bp) with Trim_galore (version 0.4.3)55,56 before being aligned again (-X 2500 -m 1 --best --strata -tryhard -S --fr -v 3 -e 100). Sorted alignments from the first and second alignments were merged together with SAMtools (version 1.8)57, followed by removal of reads aligned to the mitochondrial genome using a custom perl script (version 5.18.1). Duplicated reads were removed with Picard tools’ MarkDuplicates (version 1.94)58. Reads aligning to the blacklisted regions (generated by A. Boyle and A. Kundaje as part of the ENCODE and modENCODE projects)59 were removed using bedtools intersect (version 2.27.1)60. Subnucleosomal fragments were defined as mapped paired reads with an insertion distance smaller than 100 bp, obtained from merged mapping results. The Tn5 footprint was obtained by adapting J. Li’s preShift.pl script to take the strand orientation of a given read to take 9 bp around the start or end of the forward and reverse reads ([−4,5] and [−5,4] respectively). The preShift.pl script is available at https://github.com/riverlee/ATAC/blob/master/code/preShift.pl and the adaptation can be found at https://github.com/Edahi/NGSDataAnalysis/blob/master/ATAC-Seq/Tn5_bed9bp_full.pl. For quality control purposes, we used X. Chen’s Fragment_length_density_plot.py Python script. The script is available at https://github.com/Edahi/NGSDataAnalysis/blob/master/ATAC-Seq/Fragment_length_density_plot.py. This program plots the histogram of the distances among the mapped usable reads. Final mapping results were processed using HOMER’s makeTagDirectory program followed by the makeMultiWigHub.pl program (version 4.10.4)61 to produce normalized bigWig genome browser tracks for the whole mapping results, Tn5 footprint and subnucleosomal reads separately.
DARs
We used the complete fragments for peak calling using the MACS2 callpeak function (version 2.1.1.20160309 and -q 0.0001 --keep-dup all --nomodel --call-summits)62. The narrowpeak files from all samples and replicates for in vivo (or in vitro) experiments were merged with bedtools merge (version 2.27.1)60 to generate a universe of peaks, which was used to obtain the Tn5 footprint signal from each sample. After limma-voom normalization63, performed on the Tn5 signal, and a linear model fitted to each region, computation of significance statistics for differential enrichment (accessibility) was done by empirical Bayes moderation of the standard errors, with [−1,1] (lfc) as the interval for the null hypothesis. A region was considered differentially accessible if it had a log2[fold change] ≥ 1 and an adjusted P value (Padj) ≤ 0.05. The Tn5 signals from the in vivo and in vitro experiments were analyzed independent of one another. The MA plots used the merged signal from replicates. The R64 packages used included: IRdisplay65, limma66, edgeR67, Glimma68, Mus.musculus69, RColorBrewer70, ggplot2 (ref. 71), GenomicRanges and GenomicAlignments72 and pheatmap73.
Venn diagrams
DARs from TILs were intersected with bedtools intersect (version 2.27.1)60 with default parameters (a 1-bp overlap was considered an overlap) against the exhaustion- or activation-related regions from Mognol et al.36 (GSE88987). The overlaps were used to plot the Venn diagrams for both BATF and pMIG TILs. A one-tailed Fisher test (Fisher’s exact test on 2 × 2 contingency tables in MATLAB)74 was used to calculate the significance of the overlaps.
Heatmaps
The z score from the limma-voom63 normalized signal from TIL and CD8+ T cell samples in the regions of interest (pMIG or BATF DARs from either TILs or CD8+ T cells) was clustered by the region’s signal (cluster_rows = T) and plotted using the R library pheatmap73.
Quartile boxplots from ChIP regions
The raw Tn5 signal72 from the 2,504 ChIP-seq regions meeting the criterion log2[Tn5 signal in BATF-overexpressing cells/Tn5 signal in pMIG control cells] ≥ 3 was reads per million (RPM) normalized for both BATF and pMIG CD8+ T cells, with the RPM per replicate averaged. The regions were subdivided in quartiles with respect to the pMIG Tn5 RPM signal and the signals for both ATAC-seq and ChIP-seq data were then plotted71 together.
Known motifs analysis
A region was defined as differentially accessible when it had a twofold difference and a Padj (false discovery rate (FDR)) of <0.05. The DARs per condition and per experiment (BATF and pMIG, in vivo and in vitro) were used as input for HOMER’s findMotifsGenome.pl (version 4.10.4)57.
RNA-seq analysis
Genome browser tracks
Paired reads were mapped to STAR75 using the parameters --outFilterMultimapNmax 30 --outReadsUnmapped Fastx --outSAMattributes All --outSAMprimaryFlag OneBestScore --outSAMstrandField intronMotif --outSAMtype BAM SortedByCoordinate --quantMode GeneCounts. Mapping results were processed using HOMER’s makeTagDirectory61 twice: once for the individual replicates and a subsequent time merging them (for a less crowded genome browser session), followed by makeMultiWigHub.pl programs (version 4.10.4) to produce normalized bigWig genome browser tracks.
MA plots of differential gene expression (TILs)
Counts per gene were obtained from STAR’s STAR_gene_counts (version subread-2.0.0-source)75 Differential gene expression was analyzed with R (version 3.5.2) and the following packages: IRdisplay65, limma66, edgeR67, Glimma68, Mus.musculus69, RColorBrewer70 and gplots76. In brief, count reads from STAR were read and voom normalized after both counts per million (CPM) conversion and removal of genes whose CPM was lower than 1 across less than one-third of total samples. After limma-voom normalization, performed on the gene’s signal, and a linear model fitted to each gene, computation of significance statistics for differential gene expression was done by empirical Bayes moderation of the standard errors, without intervals for the null hypothesis ([0,0] lfc). A gene was considered differentially expressed if Padj (FDR) was ≤0.1. Colors in the MA plots for these genes indicate these parameters (red indicates genes more expressed in BATF-transduced compared with control pMIG-transduced cells, whereas blue indicates the opposite and gray indicates genes that were not differentially expressed).
MA plots of differential gene expression (in vitro)
The processing was similar to that for TILs, but [−1.2,1.2] (lfc) was used as the interval for the null hypothesis. A gene was considered differentially expressed if the absolute log2[fold change] was ≥2 and the Padj (FDR) was ≤0.05.
Gene signal heatmaps
The heatmaps are composed of the top 100 most significant (Padj) differentially expressed genes in pMIG control cells between 0 and 6 h after restimulation. The limma-voom-normalized signal for all of the pMIG-, BATF- and HKE-transduced samples was z score transformed gene-wise. The z score-normalized data were then used to plot the heatmaps with the heatmap.2 function from the gplots76 R package.
ChIP-seq analysis
Genome browser tracks
Paired raw reads were aligned to the M. musculus genome (version mm10) using bwa77 mem (version 0.7.15-r1144- dirty). Unmapped reads were trimmed to remove adapter sequences and clipped by 1 bp with Trim_galore (version 0.4.3)55,56 before being aligned again. Sorted alignments from the first and second alignments were merged together with SAMtools (version 1.8), followed by removal of reads aligned to the mitochondrial genome using a custom perl script (version 5.18.1). Duplicated reads were removed with Picard tools’ MarkDuplicates (version 1.94)58. Reads aligning to the blacklisted regions (generated by A. Boyle and A. Kundaje as part of the ENCODE and modENCODE projects) were removed using bedtools60 intersect (version 2.27.1). Final mapping results were processed using the HOMER61 makeTagDirectory program followed by the makeMultiWigHub.pl program (version 4.10.4) to produce normalized bigWig genome browser tracks.
Venn diagram
For each sample, peaks were called using the MACS2 (version 2.1.1.20160309)62 callpeak function, using the sample’s respective input dataset, a q value of 0.05 and the --keep-dup all and --nomodel parameters. The narrowpeak files among replicates were merged using bedtools merge60 (version 2.27.1). To identify overlapping genes by the merged narrowpeak files per condition, we used the University of California, Santa Cruz M. musculus mm10 annotation genes. Called peaks were assigned to a gene if they overlapped with a window containing the body of the gene (the longest transcription unit for the gene locus definition) plus the 20-kilobase (kb) region upstream of the transcription start site (TSS) and the 5-kb region downstream of the 3′ end of the gene. Each gene was considered only once and the whole gene set was used to find shared genes among the samples being compared. The overlap was conducted with the bedtools60 intersect function (version 2.27.1). Venn diagrams of shared overlapping genes were produced using R (version 3.5.2), as well as the libraries VennDiagram78,79 and viridis80.
Probability per base pair of the BATF binding site
Peaks from BATF-transduced CD8+ T cells subjected to ChIP-seq with anti-BATF antibodies were functionally annotated to mm10 using the HOMER61 annotatePeak.pl program. The distance to the nearest TSS and gene name were filtered from the annotation results. A sublist of the genes differentially expressed between BATF- and pMIG-transduced CD8+ T cells, identified by RNA-seq analysis, was used to subset separately the peak annotation results for genes up- and downregulated in BATF-transduced cells. The genomic histograms were generated using R (3.5.2)64 and ggplot2 (ref. 71) with all of the peak results, whereas the up- and downregulated histograms used the subset of genes generated above. The percentage of genes closer than 20 kb was obtaining by taking the absolute value to the closest TSS that was ≤20 kb. The distances were numerically sorted and an empirical cumulative distribution function was generated based on the data.
Removal of spurious peaks
All of the peaks from all of the different conditions and replicates were merged into a singularity table, keeping track of which condition belonged to which region. For the superset of peaks belonging to the αBATF immunoprecipitation, we kept the peaks whose average RPM signal across the pMIG-, BATF- and HKE-transduced input samples was lower than 0.75 times the αBATF immunoprecipitation RPM signal from BATF-overexpressing cells. Similarly, for the aIRF4 immunoprecipitation superset, we kept peaks where said input signal was lower than 0.75 times that of the αIRF4 immunoprecipitation RPM signal from pMIG control cells. These filtered supersets were used for all subsequent analyses.
Normalized αIRF4 ChIP-seq reads report accurately on IRF4 binding
It cannot be taken for granted that a difference in the normalized αIRF4 signal (in RPM) between pMIG and BATF-overexpressing cells reports on a change in IRF4 binding at the peak in question. The general issue is that normalization of the IRF4 signal at a particular peak to total mapped reads introduces a second independent variable into the measurement. If, for example, there was free IRF4 in the nucleus of pMIG cells, and if overexpressed BATF recruited this additional IRF4 to sites in DNA, then a greater total amount of IRF4-bound DNA would be precipitated from BATF-overexpressing cells. For any individual site where exactly the same amount of IRF4-bound DNA was precipitated as from pMIG cells, the normalization would result in an artifactually lower RPM value.
To address this issue, we utilized a subset of nonspecific background DNA regions that were equally represented in the input samples and immunoprecipitated samples from the same cells. The reads mapping to these regions in immunoprecipitated samples—which seem to represent a low fraction of input DNA carried along by nonspecific binding to the protein A/protein G beads—can serve as an internal standard. Specifically, we selected the 20 spurious peak regions with the largest ATAC-seq signal in pMIG cells (see the preceding section), since the high total signal ensured that any fractional contribution to the signal from actual IRF4 binding would be negligible. The spurious αIRF4 ChIP-seq signal from these regions was consistently the same in BATF-overexpressing and pMIG cells (Extended Data Fig. 9), implying that normalization did not distort the comparison between BATF-overexpressing and pMIG samples, and that a decrease in the normalized αIRF4 signal for an individual specific αIRF4 peak meant that there was an actual decrease in IRF4 binding at that peak.
Scatter and contour plots
Each scatterplot is based on the log2 of the RPM immunoprecipitation signal of a subset of regions representing those of interest (for example, the αBATF immunoprecipitation signal from BATF-overexpressing cells versus the αBATF immunoprecipitation signal from pMIG control cells). We took the union of peaks for the illustrated samples and fetched the αIRF4 and/or αBATF average RPM immunoprecipitation signal (as indicated in the graphs) followed by a log2 transformation. These normalized signals were then processed in R using ggplot’s function geom_bin2d(bins = 300) for the scatterplots (density; that is, occurrences of points per region) and geom_density_2d(bins = 30) for the contour plots71.
Overlap measurement as reads-in-peaks percentage
For the Venn diagrams of Fig. 6a and Extended Data Fig. 8b, we took the union of the peaks in the two conditions considered (for example, the union of αBATF ChIP-seq peaks from BATF-overexpressing and pMIG control cells), conserving the information on whether an individual peak was unique to one condition or shared between the two conditions. For each condition, we divided the number of reads that mapped to peaks shared by both conditions by the total reads mapped to peaks in that condition, obtaining the reads-in-peaks percentage in shared peaks and—by complementation—the reads-in-peaks percentage in unique peaks.
Histograms of signal distribution among subsets of peaks
For the histograms of Fig. 6a and Extended Data Fig. 8b, the identities of peaks that were unique to a specified condition or shared with a second condition were used to fetch the RPM-normalized ChIP-seq signals for peaks in each subset. The values were then log2 transformed, the median value was calculated from the log2-transformed data and its distribution was plotted as a histogram using R’s ggplot2 function geom_histogram71.
Heatmaps
We used the deepTools81 computeMatrix function (with the parameters --referencePoint center -a 1000 -b 1000 --binSize 50 --averageTypeBins mean --missingDataAsZero -p 4) to compute the signal matrices across all of the conditions. The regions that were used were the input-corrected peaks (one peakset per condition). The bigWig datasets used to fetch the signal were the HOMER-normalized bigWigs (the same ones used in the genome browser track). We then proceeded to give this program’s output as input to the deepTools plotHeatmap function (with the parameters --averageType mean --plotType se --averageTypeSummaryPlot mean --sortRegions descend --sortUsing mean --sortUsingSamples 6 --refPointLabel Center --missingDataColor yellow).
Statistical analysis
No statistical method was used to predetermine sample size. No data were excluded from the analyses. Tumor-bearing mice were randomly assigned to adoptive-transfer treatment groups. In most cases, tumor sizes were measured in a blinded manner by DLAC staff, except during the holiday season or when the institute was under restricted access due to the COVID-19 shutdown. Investigators were not blinded to sample identity when analyzing T cells recovered from the tumors. Details of the sample sizes, replicates and statistical tests used are provided in the individual figure captions.
Reporting Summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.