DOI: 10.5061/dryad.8kprr4xx9
title: Continuous expression of TOX safeguards exhausted CD8 T cell epigenetic fate

description: CD8 T cell exhaustion is a major barrier limiting anti-tumor therapy.
 Though checkpoint blockade temporarily improves exhausted CD8 T cell (Tex)
 function, the underlying epigenetic landscape of Tex remains largely
 unchanged, preventing their durable “reinvigoration.” Whereas the
 transcription factor (TF) TOX has been identified as a critical initiator
 of Tex epigenetic programming, it remains unclear whether TOX plays an
 ongoing role in preserving Tex biology after cells commit to exhaustion.
 Here, we decoupled the role of TOX in the initiation versus maintenance of
 CD8 T cell exhaustion by temporally deleting TOX in established Tex.
 Induced TOX ablation in committed Tex resulted in apoptotic-driven loss of
 Tex, reduced expression of inhibitory receptors including PD-1, and a
 pronounced decrease in terminally differentiated subsets of Tex cells.
 Simultaneous gene expression and epigenetic profiling revealed a critical
 role for TOX in ensuring ongoing chromatin accessibility and
 transcriptional patterns for key Tex gene modules in committed Tex cells.
 Moreover, when exposed to effector-driving conditions, inducibly
 TOX-deleted established Tex acquired an altered chromatin landscape with
 increased accessibility at cytotoxic genes typically accessible in Teff
 cells, thus undergoing partial reprogramming into a more functional state.
 Together, these findings suggest that continuous TOX expression in
 established Tex acts as a durable epigenetic barrier to reinforce the Tex
 developmental fate by simultaneously maintaining Tex epigenetic commitment
 while restraining differentiation into Teff. Manipulation of TOX even
 after Tex establishment could therefore provide a therapeutic opportunity
 to rewire Tex biology in settings of chronic infection or cancer. The
 secondary goal of this dissertation was to develop a novel Tex
 fate-mapping mouse model driven to track the fate of developing Tex and
 manipulate Tex in a lineage-restricted fashion. Given the selectively high
 expression of TOX in Tex, compared to other peripheral non-Tex CD8
 lineages, we used the Tox locus to drive this model (termed ToxTREx),
 which was engineered with a T2A-hmKO2-P2A-CreERT2 cassette knocked into
 the Tox locus after the last exon. We confirmed intact TOX function, hmKO2
 reporter detection and TOX-driven Cre recombinase activity in this model
 and identified further optimization that will be necessary to improve Cre
 efficiency and specificity of this model for the Tex lineage.
 Nevertheless, the ToxTREx model could enable insightful studies that
 address existing and emerging questions in Tex ontogeny, differentiation
 and function.
description: Cells from inducible-Cre (Rosa26CreERT2/+Toxfl/fl P14) mice where
 TOX was temporally deleted from mature populations of LCMV-specific T
 exhausted cells after establishment of chronic LCMV infection 5 days post
 infection were subjected to scRNA and scATACseq coassay,naive cells and WT
 cells were used as controls. Analysis pipeline
 developed by Josephine Giles and vignettes published by Satija and Stuart
 labs.Transcript count and peak accessibility matrices deposited in
 GSE255042,GSE255043. Seurat/Signac was used to process
 the scRNA and scATACseq coassay data The processed
 Seurat/Signac object above was subsequently used for downstream RNA and
 ATAC analyses as described below: DEGs between TOX WT
 and iKO cells within each subset were identified using FindMarkers
 (Seurat, Signac), with a log<sub>2</sub>-fold-change threshold
 of 0, using the SCT assay. DACRs were identified using FindMarkers using
 the "LR" test, with a log<sub>2</sub>-fold-change
 threshold of 0.1, a min.pct of 0.05, and included the number of counts as
 a latent variable. DEGs and DACRs were filtered with a false discovery
 rate (FDR) of less than 0.05, using the Benjamini–Hochberg method to
 adjust <em>P</em>-values. Gene ontology (GO) analysis of DEGs
 used Metascape (https://metascape.org/) with all expressed genes as the background gene list. AddModuleScore (Seurat, Signac) was used to calculate per-cell gene set enrichment scores and peak set enrichment scores. Gene set enrichment analysis (GSEA) was performed and visualized using clusterProfiler (4.8.1) and enrichplot (1.20.0). T<sub>ex-prog</sub> and T<sub>ex-term</sub> transcriptional signatures were generated from the Giles et al dataset (NI 2022), using FindMarkers to identify DEGs between the "Exh-Prog" cluster and the "Exh-Term" and "Exh-TermGzma" clusters, filtered with a FDR of less than 0.05 and log<sub>2</sub> fold-change threshold of 0.25. Peak sets of subset transition ACRs were identified using FindMarkers to perform pairwise comparisons between TOX WT T<sub>ex-prog</sub> and T<sub>ex-int</sub>, T<sub>ex-klr</sub>, and T<sub>ex-term</sub> respectively (using the "LR" test, with a min.pct of 0.05, included the number of counts as a latent variable and cut-off at top 5000 ACRs by log<sub>2</sub> fold-change threshold). Peak sets of naïve, T<sub>eff</sub>, T<sub>mem</sub> and T<sub>ex</sub>-specific ACRs were generated as previously described (Khan et al, Nature, 2019). Bedtools intersect was used to identify overlapping peaks between different datasets. Genome coverage tracks were generated in Signac using CoveragePlot, PeakPlot, and AnnotationPlot. Genomic locations of peaks were determined using ClosestFeature (Signac) and nearestTSS (edgeR), with promoter-TSS defined as -1kB to +100bp. The average local chromatin accessibility at each gene was determined using GeneActivity (Signac). TF motif enrichment was calculated using the JASPAR2020 function getMatrixSet (species 9606) and Signac functions CreateMotifMatrix, CreateMotifObject, and FindMotifs. Open peaks in the clusters of interest were used as background, using AccessiblePeaks and MatchRegionStats (Signac).
description: # Continuous expression of TOX safeguards exhausted CD8 T cell epigenetic
 fate
 [https://doi.org/10.5061/dryad.8kprr4xx9](https://doi.org/10.5061/dryad.8kprr4xx9) Seurat/Signac pipeline for multiomic scRNA-seq and scATAC-seq dataset, generated following inducible TOX deletion in LCMV-Cl13 Author Yinghui Jane Huang ## Script information **Purpose:** Generate and process Seurat/Signac object for downstream analyses **Written:** Nov 2021 through Oct 2022 **Adapted from:** Analysis pipeline developed by Josephine Giles and vignettes published by Satija and Stuart labs **Input dataset:** Transcript count and peak accessibility matrices deposited in GSE255042,GSE255043 # Signac Object Generation 1\) Create individual signac objects for each sample from the raw 10x cellranger output. 2\) Merge individual objects to create one seurat object. 3\) Add metadata to merged seurat object. Following are the steps in the attached html file for analysis of the paired data (ATAC+RNA) * Load fragments data(***atac_fragments.tsv.gz***) and 10x cellranger output counts data file(***filtered_feature_bc_matrix.h5***) * Get gene annotations * Create Seurat object with RNA and ATAC data *  filter out low quality cells * call peaks using MACS2 with the CallPeaks() function. Here we call peaks on all cells together * combine peaks and merge individual seurat/signac objects to one Object * Perform standard Signac normalization steps Further downstream analyses The processed Seurat/Signac object above was subsequently used for downstream RNA and ATAC analyses as described below: DEGs between TOX WT and iKO cells within each subset were identified using FindMarkers (Seurat, Signac), with a log2-fold-change threshold of 0, using the SCT assay. DACRs were identified using FindMarkers using the "LR" test, with a log2-fold-change threshold of 0.1, a min.pct of 0.05, and included the number of counts as a latent variable. DEGs and DACRs were filtered with a false discovery rate (FDR) of less than 0.05, using the Benjamini–Hochberg method to adjust *P*-values. Gene ontology (GO) analysis of DEGs used Metascape ([https://metascape.org/](https://metascape.org/)) with all expressed genes as the background gene list. AddModuleScore (Seurat, Signac) was used to calculate per-cell gene set enrichment scores and peak set enrichment scores. Gene set enrichment analysis (GSEA) was performed and visualized using clusterProfiler (4.8.1) and enrichplot (1.20.0). Tex-prog and Tex-term transcriptional signatures were generated from the Giles et al dataset (NI 2022), using FindMarkers to identify DEGs between the "Exh-Prog" cluster and the "Exh-Term" and "Exh-TermGzma" clusters, filtered with a FDR of less than 0.05 and log2 fold-change threshold of 0.25. Peak sets of subset transition ACRs were identified using FindMarkers to perform pairwise comparisons between TOX WT Tex-prog and Tex-int, Tex-klr, and Tex-term respectively (using the "LR" test, with a min.pct of 0.05, included the number of counts as a latent variable and cut-off at top 5000 ACRs by log2 fold-change threshold). Peak sets of naïve, Teff, Tmem and Tex-specific ACRs were generated as previously described (Khan et al, Nature, 2019). Bedtools intersect was used to identify overlapping peaks between different datasets. Genome coverage tracks were generated in Signac using CoveragePlot, PeakPlot, and AnnotationPlot. Genomic locations of peaks were determined using ClosestFeature (Signac) and nearestTSS (edgeR), with promoter-TSS defined as -1kB to +100bp. The average local chromatin accessibility at each gene was determined using GeneActivity (Signac). TF motif enrichment was calculated using the JASPAR2020 function getMatrixSet (species 9606) and Signac functions CreateMotifMatrix, CreateMotifObject, and FindMotifs. Open peaks in the clusters of interest were used as background, using AccessiblePeaks and MatchRegionStats (Signac). Processed Seurat/Signac R objects are attached ***20241111_ProcessedSignacObject_InducibleTOXDeletion_WithoutLCMVArmRechallenge.Rds*** -Processed Seurat/Signac R data object for the samples described under GEO [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE255042](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE255042) ***20241111_ProcessedSignacObject_InducibleTOXDeletion_WithLCMVArmRechallenge.Rds*** Processed Seurat/Signac R data object for the samples described under GEO [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE255043](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE255043) Following commands can be used in R to load the r objects for exploration of the completely processed data ``` ProcessedSignacObject_InducibleTOXDeletion_WithoutLCMVArmRechallenge=readRDS('20241111_ProcessedSignacObject_InducibleTOXDeletion_WithoutLCMVArmRechallenge.Rds') ProcessedSignacObject_InducibleTOXDeletion_WithLCMVArmRechallenge.Rds=readRDS('20241111_ProcessedSignacObject_InducibleTOXDeletion_WithLCMVArmRechallenge.Rds') ``` ## Code/software Code with comments is in the attached html documents ,***SeuratSignacPipelineForGSE255042.html *** & ***SeuratSignacPipelineForGSE255043.html*** have the complete end to end R script to process raw data from 10x cellranger output shared on GEO to create a normalized processed Signac/Seurat Object shared as Rds files here. # Access information Publicly accessible locations of the raw data: * [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE255042](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE255042) * [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE255043](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE255043)