Detection of mitochondrial DNA mutations in T cells following 5-FU or cisplatin exposure.

Introduction
T cells have been in the limelight of cancer immunotherapy for their roles in checkpoint blockade immunotherapies, adoptive T cell therapies, and neoantigen vaccines. Despite early successes, immunotherapy remains ineffective for the majority of patients with solid tumors. T cell dysfunction has been considered one of the major barriers to successful immunotherapy, and studies have shown an association between the decreased efficacy of cancer immunotherapies and T cell dysfunction.1,2,3,4 Unmasking sources of T cell dysfunction is crucial to address the current limitations of cancer immunotherapy. In addition to chronic antigen stimulation and tumor microenvironment that lead to T cell dysfunction, recent mouse studies suggested that mitochondrial function might be linked to T cell dysfunction.5 For example, one study showed that the persistent antigen stimulation could impair mitochondrial oxidative phosphorylation (OXPHOS) in mouse T cells, driving the expression of genes associated with T cell dysfunction.6 Another study found that the chronic stimulation combined with hypoxia further disrupted mitochondrial integrity, pushing T cells toward dysfunction in mice.7 The dysfunctional mitochondria generated intolerable levels of reactive oxygen species (ROS) and promoted T cell exhaustion. Lastly, mitochondrial oxidative damage was linked to the dysfunction of regulatory T cells, in turn leading to heightened autoimmune responses in an experimental autoimmune encephalitis model.8 However, the specific sources of mitochondrial dysfunction that result in T cell dysfunction remain unknown.
Mammalian mitochondrial DNA (mtDNA) contains a very high density of genes, and it encodes 22 transfer RNAs (tRNAs), two ribosomal RNAs (rRNAs), and 13 proteins that are necessary for OXPHOS.9 However, mtDNA has an estimated 20 times higher mutation rate than genomic DNA.10 Unlike the nuclear genome, mtDNA lacks protective histones. Additionally, unlike the nuclear DNA polymerase, the mtDNA polymerase has been found to lack efficiency and proofreading capabilities. Although individual mtDNA mutations may be insignificant, they can accumulate over time in cells and their daughter cells. As a result, a mixture of wild-type and mutated mtDNA may be present in individual cells. This observation is known as mitochondrial heteroplasmy.11 A recently developed single-cell technique enables the detection of mtDNA mutation and heteroplasmy in individual cells.12,13,14,15
Two common chemotherapeutic agents were selected for the investigation in this study. 5-Fluorouracil (5-FU) is a thymidylate synthase inhibitor, and it has been used to treat breast, head and neck, and colorectal cancer, as well as several other types of cancer. 5-FU works to both inhibit nucleotide synthesis and incorporate its metabolites into the growing DNA or RNA strand to inhibit synthesis.16 5-FU acts on actively replicating cells, therefore it also damages healthy cellular DNA, which is presumed to cause the side effects associated with 5-FU treatment.17 Another common platinum-based drug, cisplatin, is used to treat lung, ovarian, bladder, and several other cancers. The cytotoxicity of cisplatin results from its ability to crosslink purines to inhibit DNA synthesis and induce tumor cell apoptosis.18 Cisplatin, like 5-FU, is designed to target tumor cells but also damages healthy cells to cause side effects.19 5-FU and cisplatin thereby indiscriminately target DNA to induce cellular distress.
Mitochondrial integrity is central to T cell metabolism and function, yet the impact of chemotherapy-induced mtDNA mutations on immunotherapy efficacy remains poorly understood. Here, we aimed to study the impact of chemotherapeutic agents on T cells, in addition to the anti-proliferative properties of these agents. Because mtDNA is more susceptible to mutations than genomic DNA, we hypothesize that 5-FU and cisplatin may induce mtDNA mutations in T cells and impact T cell function. In this study, we exposed healthy T cells to 5-FU and cisplatin and used a sensitive single-cell sequencing technique to detect mtDNA mutations. We discovered that common chemotherapeutics, 5-FU and cisplatin, induced mtDNA mutations in T cells, with potential consequences for mitochondrial activity and T cell phenotypes.

Results

5-FU inhibits T cell proliferation and survival at high doses but has a minimal impact on T cell transcriptome profiles
To investigate the potential impact of 5-FU on T cells, peripheral blood T cells from four healthy donors were stimulated with anti-CD3/CD28 beads on day 0. Proliferating T cells were then treated with 5-FU at various concentrations from day 3 to day 7. Treated T cells received a second round of CD3/CD28 stimulation on day 10, followed by the 5-FU treatment. A total of three rounds of 5-FU treatments were performed (Figure 1A). The total numbers of T cells at the end of each cycle were counted and normalized before the next CD3/CD28 stimulation (Figure 1B). As expected, the results indicated that 5-FU strongly inhibited T cell proliferation and survival at high doses. After three cycles, T cells could not survive beyond the 1.56 μM dose.
In addition to T cell proliferation, we sought to understand the impact of 5-FU on the gene expression profiles of T cells at the single-cell level. T cells treated with three cycles of 5-FU (1.56 μM), together with the negative vehicle control, were subjected to single-cell RNA sequencing (scRNA-seq) analysis, followed by clustering and trajectory analyses. The 5-FU-treated and untreated T cells from three donors were integrated for clustering analysis based on gene expression profiles (Figure 1C). Clustering analysis identified 10 clusters across the donors and treatments (Figure 1D). For donor D2008, we found that cluster 7 had more untreated T cells compared to 5-FU-treated T cells (7.48% vs. 2.77%) (Figure 1D). On the other hand, for donor D6053, clusters 3, 9, and 10 had more 5-FU-treated T cells compared to untreated T cells (15.10% vs. 6.47%, 3.20% vs. 0.50%, 3.39% vs. 1.27%, respectively) (Figure 1D). Donor D1007 had no significant difference in cluster size between the two treatment groups. Therefore, we found no consistent effect of 5-FU across all three donors in clustering analysis. Similarly, we identified multiple T cell-related genes, MAL, OX40, and RUNX3, that were differentially expressed in the 5-FU-treated T cells, compared to the negative control (Figure S1). However, those differences were not observed across all three donors. Lastly, trajectory analysis was performed to study T cell differentiation in these samples (Figure 1E). No difference was found when we compared the treated vs. untreated samples from individual donors, but differences were found between donors (Figure 1E). When examining MKI67, TIGIT, and CCR6 gene expression across inferred pseudotime, again, no difference in expression between 5-FU-treated cells and untreated cells was found in any donor (Figure S2). Taken together, 5-FU inhibits T cell proliferation and survival in a dose-dependent manner. However, we did not observe any significant difference between 5-FU-treated T cells and the negative control across all donors based on the scRNA-seq analysis and trajectory analysis.

A single-cell mtDNA sequencing technique reveals that 5-FU induces mtDNA mutations in T cells
Because mtDNA is more susceptible to mutations compared to genomic DNA, we hypothesize that 5-FU can induce mtDNA mutations in T cells. To test this, we utilized a single-cell technique to sequence mtDNA and perform ATAC-seq (assay for transposase-accessible chromatin using sequencing) concurrently (Figure 2A).12,13,14,15 In this technique, cells were fixed before being lysed, so mitochondria and mtDNA were retained prior to the single-cell ATAC-seq sample preparation. This allowed us to obtain an mtDNA library, as well as a chromatin accessibility library for clustering analysis. Importantly, this single-cell approach was significantly more sensitive than conventional bulk sequencing, because it could help us to detect subtle differences in mtDNA at the single-cell level. Here, the same 5-FU-treated T cell samples described in Figure 1 were subjected to single-cell mtDNA/ATAC-seq (Figure 2A).
For the donor D1007, 4,566 high-quality cells were obtained from the 5-FU-treated sample. Similarly, 5,586 high-quality cells were obtained from the untreated sample. The 5-FU-treated and untreated cells were merged and integrated for clustering analysis based on chromatin accessibility (Figure 2B). Twenty-one distinct clusters of cells were obtained (Figure 2C). Next, we attempted to identify mtDNA mutations induced by 5-FU. Considering the complexity of mitochondrial heteroplasmy, the analysis was streamlined and only focused on mtDNA mutations that were found in the 5-FU-treated T cells but were absent from the untreated T cells, similar to the previous approach for identifying tumor-specific mtDNA mutations.20 Unlike mtDNA heteroplasmy, these de novo mutations were induced by the 5-FU treatment. In this pair of samples, a total of 1,000 mtDNA mutations were detected in the 5-FU sample but not in the negative control (Figure 2D). Next, we followed the original publication and applied the restricted quality-control parameters to filter out low-quality mtDNA mutations.14 A total of 78 mtDNA mutations were identified after the filtering. Lastly, we focused our investigation on mtDNA mutations that induced non-synonymous substitutions and that were detected in 10 or more cells with 10% or more mtDNA heteroplasmy. A total of five mtDNA mutations were identified in the 5-FU-treated T cells (Figure 2D; Table S1). The top non-synonymous mtDNA mutations with the highest numbers of cells carrying mtDNA mutations are shown (Figures 2E–2G). T cells with mtDNA mutations are highlighted in red in the UMAP plots.
5-FU-treated T cells with 10775G>A (ND4:V6I) mutation were predominately clustered in cluster #14, a cluster that was mostly contributed by 5-FU-treated T cells (Figures 2B and 2E). Notably, an elevation of T cell exhaustion marker EOMES was found in cluster #14, followed by another exhaustion marker TIGIT (Figures S3A and S3B).21 Additionally, a minimum of interleukin (IL)-7R signal was found in cluster #14, suggesting a lack of self-renewal, memory-like capacity21 (Figure S3C). Furthermore, cluster #14 also shows a significant decrease in CXCR6 accessibility (Figure S3D). Lastly, Gene Ontology (GO) enrichment analysis suggested that the top 100 down-regulated genes were associated with the regulation of T cell responses. (Figure S3E). Therefore, the gene profile of cluster #14 suggested that this cluster was associated with a T cell exhaustion phenotype. We also detected other mutations unique to the 5-FU-exposed T cells, including 9047T>C (ATPase6:I174T) and 11592G>A (ND4:R278Q), but they were not concentrated in clusters (Figures 2F and 2G). To investigate whether these mtDNA mutations were carried by CD4+ or CD8+ T cells, 5-FU-treated T cells from D1007 were subsequently divided into CD4+ and CD8+ subpopulations, based on the CD4 and CD8A chromatin accessibilities. After division, cluster analysis was performed separately (Figures S4A and S4B). Cells with these three mtDNA mutations were found in the CD8+ T cell populations but not in the CD4+ T cell populations (Figures S4C and S4D).
Using the same approach, samples from donor D6053 were analyzed by single-cell mtDNA sequencing. 2,942 high-quality cells were obtained from the 5-FU-treated sample, and 3,475 high-quality cells were obtained from the untreated sample. 5-FU-treated and untreated T cells were merged and integrated for cluster analysis (Figure S5A). Nineteen distinct clusters of T cells were identified (Figure S5B). A total of eight unique mtDNA mutations in the 5-FU-treated T cells were identified by mtDNA genotyping analysis (Figure S5C). The non-synonymous mtDNA mutations 4632G>A (ND2:A55T), 7785T>C (COII:I67T), and 8154G>A (COII:G190E) are highlighted in red in UMAP plots (Figure S5D). Although non-synonymous mtDNA mutations were detected, these mtDNA mutations were not concentrated in any clusters. Analysis on D2008 and D0010 similarly found that 5-FU-specific non-synonymous mtDNA mutations were present in some cells but were scattered in UMAP plots (Figures S6 and S7). The same as mutations identified in D1007, these mutations were found in the CD8+ T cell subset. Taken together, these results indicate that 5-FU treatment can induce mtDNA mutations in T cells from four out of four donors, which can be detected by the single-cell mtDNA sequencing approach (Table S1). However, one caveat is that mtDNA mutations might induce severe damage, leading to T cell death, which cannot be detected by this experimental approach.

Cisplatin induces high numbers of mtDNA mutations in T cells
Since we observed that 5-FU could induce mtDNA mutations in T cells, we further tested whether cisplatin, another common chemotherapeutic agent, could also induce mtDNA mutations in T cells. Following the same approach, peripheral blood T cells isolated from three healthy donors were stimulated with anti-CD3/CD28 beads for three cycles and treated with cisplatin at various concentrations (Figure 3A). Similar to 5-FU, the results showed that cisplatin strongly inhibited T cell proliferation and survival at higher doses (Figure 3B). T cells could not survive the dose higher than 0.313 μM at the end of cycle 3. T cells treated with cisplatin (0.313 μM), together with the negative control, were subjected to the same single-cell mtDNA/ATAC sequencing. The same restricted criteria for identifying non-synonymous mtDNA mutations were applied.
From donor D1007, 10,353 high-quality cells were obtained from the cisplatin-treated sample. Similarly, 8,900 high-quality cells were obtained from the untreated sample. Cisplatin-treated and untreated T cells were merged and integrated for downstream analysis. Clustering analysis from the chromatin accessibility data identified 24 clusters in the analyzed T cells (Figures 3C and 3D). A total of 36 mtDNA mutations unique to cisplatin-treated T cells were detected after quality control and filtering criteria (Figure 3E). The top identified mtDNA mutations are highlighted in red in the UMAP plots (Figure 3F). Additionally, 6451T>C (COI:L183P), 7378T>C (COI:L492P), and 8584G>A (ATPase6:A20T) were associated with cluster #20, which was characterized by elevated EOMES and decreased CCR5 and CTLA4 accessibilities (Figures S8A–S8C). T cells carrying 3850G>A (ND1:A182T) were in part associated with cluster #14, which was characterized by the lack of SELL (CD62L) and elevated TNFSF11 (RANK ligand) and TIGIT accessibilities (Figures S8A and S8B). Additionally, cluster #14 is the only cluster with a lack of SELL accessibility (Figure S8D). The mutations 12359C>T (ND5:T8I) and 5382C>T (ND2:L305F) concentrated across a few clusters (Figure 3F). Lastly, cisplatin-treated T cells from D1007 were subsequently divided into CD4+ and CD8+ subsets (Figures S9A and S9B), and the six mtDNA mutations described above were found in the CD8+ subset but not in the CD4+ subset (Figures S9C and S9D).
Using the same approach, mtDNA mutations specific to cisplatin-treated T cells were also identified in donor D2008. In this sample, 8,040 high-quality cells were obtained from the cisplatin-treated sample and were merged and integrated with untreated T cells (Figure S10A). Clustering analysis of the integrated samples identified 23 distinct clusters of cisplatin-exposed T cells (Figure S10B). After filtering workflows, a total 31 total cisplatin-induced non-synonymous mtDNA mutations were detected (Figure S10C). One of the top mutations from the mtDNA genotyping analysis was 3413G>C (ND1:G36A), which was not detected in the untreated T cells (Figure S10D). Two other mutations unique to the cisplatin-treated T cells, 10753T>C (ND4L:L95P) and 15698C>A (Cytb:R318S), are also shown in Figure S10D. Similar to the previous method, cisplatin-treated T cells from D2008 were also divided into CD4+ and CD8+ subsets and clustered separately (Figures S11A and S11B). In contrast to D1007, there were two mutations detected for D2008, 13552G>A and 13808T>C, that were mostly found in the CD4+ T cell population (Figure S11C). The 3413G>C mutation was found in the CD8+ T cell population (Figure S11D). Analysis of D6053 also led to the detection of 16 non-synonymous mtDNA mutations that were not present in untreated T cells, but the mutations were not concentrated in any clusters (Figure S12). Taken together, the results indicate that cisplatin can induce mtDNA mutations in T cells from three out of three donors.

Distinct mtDNA mutations are detected in the fraction of T cells with reduced mitochondrial activity
Mitochondrial membrane potential is a key indicator of mitochondrial activity,22 and it can be measured directly by a fluorescent dye TMRM (tetramethylrhodamine, methyl ester). TMRM has been used previously to enrich T cell populations with different mitochondrial membrane potentials.23 To test whether mtDNA mutations could damage mitochondrial activity, cisplatin-treated T cells were stained with TMRM and incubated at 37°C for 30 min to measure the mitochondrial activity. The TMRM low and high fractions, about 10%–15% of the T cell populations, were sorted by fluorescence-activated cell sorting (FACS), followed by the downstream single-cell analysis (Figure 4A). We started the analysis with the D1007 cisplatin-treated T cells because of the high numbers of cells with mutations. The TMRM low fraction, 12.5% of the cisplatin-treated T cells from donor D1007, was sorted for single-cell mtDNA/ATAC-seq analysis (Figure 4B). The TMRM-high (11.0%) T cell population was also obtained. To identify mtDNA mutations linked to mitochondrial dysfunction, we searched for mutations present in both the TMRM low fraction and the original-cisplatin-treated T cells. A total of 5,808 high-quality cells were obtained from this TMRM low fraction. Likewise, 4,762 high-quality cells were obtained from the TMRM-high fraction. The TMRM low fraction was merged and integrated with the original-cisplatin-treated T cells and untreated T cells for downstream analysis. Cluster analysis of the combined dataset identified 29 distinct clusters (Figures 4C and 4D).
Next, detected mtDNA mutations were highlighted in red for the original-cisplatin-treated T cells and in purple for the TMRM low fraction in UMAP plots, respectively (Figures 4E and 4F). Six aforementioned mtDNA mutations detected in the original-cisplatin-treated T cells were also detected in the TMRM low fraction of treated T cells (Figures 3F, 4E, and 4F). Five of the six originally detected mutations were detected only in the TMRM low fraction and not the TMRM high fraction. Cisplatin-treated T cells harboring the mutation 3850G>A (ND1:A182T) were detected and clustered in the TMRM low fraction. Eleven cells (0.01%) in the TMRM high fraction harbored the 3850G>A (ND1:A182T) mutation (Table S1). The mutations 5382C>T (ND2:L305F) and 12359C>T (ND5:T8I) were detected exclusively in the original-cisplatin-treated cells and the TMRM low fraction of cisplatin-treated T cells. The other three originally detected mutations, 6451T>C (COI:L183P), 7378T>C (COI:L492P), and 8584G>A (ATPase6:A20T) were also detected in the TMRM low fraction but not in the TMRM high fraction (Table S1). Next, we used the same approach to study the TMRM high fraction. Single-cell datasets from TMRM high and TMRM low fractions were merged with the original-cisplatin-treated and nontreated T cell datasets (Figures S13A and S13B). Three unique mtDNA mutations, 12560A>G (ND5:Q75R), 10158T>C (ND3:S34P), and 15570T>C (Cytb:L275P), were identified in the TMRM high fraction that were not detected in the TMRM low fraction (Figures S13C–S13E; Table S1).
Similarly, TMRM staining and sorting were performed for another cisplatin-treated donor, D2008. A total of 3,004 high-quality cells were obtained from the TMRM low fraction (12.8%). Similarly, 2,409 high-quality cells were obtained from the TMRM high fraction (12.5%). The original-cisplatin-treated and untreated T cells were merged and integrated with the TMRM low T cell fraction for downstream analysis (Figure S14A). Cluster analysis from the merged scATAC-seq dataset yielded 25 clusters (Figure S14B). Three of the originally detected cisplatin-induced mtDNA mutations were also detected in the TMRM low T cell fractions (Figures S14C and S14D). Next, the datasets from TMRM high and low fractions were merged with the original D2008 cisplatin-treated and nontreated T cell datasets (Figure S15A). Clustering analysis of the combined dataset identified 28 clusters of cells (Figure S15B). Eight unique mutations were detected in both the original-cisplatin-treated cells and the TMRM high fraction but not in the TMRM low fraction (Table S1). Three of them are shown in Figures S15C–S15E. The 14087T>C (ND5:I584T) was detected in both TMRM high and low fractions, suggesting that it had no impact on mitochondrial activity. Collectively, by using TMRM staining to sort TMRM low T cells, we identified several mutations associated with decreased mitochondrial activity. These results suggested that some mtDNA mutations induced by cisplatin might reduce mitochondrial activity in T cells. Conversely, mutations that were only found in the TMRM high fraction might potentially improve mitochondrial health.

mtDNA mutations in T cells with reduced mitochondrial activity are also detected in T cells with reduced mitochondrial mass
To further study the relationship between cisplatin-induced mtDNA mutations and mitochondrial function, we used a fluorescent dye MitoTracker Green FM as an indicator of mitochondrial mass. Similar to the TMRM staining, cisplatin-treated T cells from D1007 were stained with MitoTracker Green FM. The mitochondrial mass low and high fractions, about 10%–15% of T cell populations, were sorted by FACS followed by single-cell mtDNA/ATAC-seq analysis (Figure 5A). A total of 7,427 high-quality cells were obtained from the mitochondria mass low fraction, and a total of 3,474 cells were obtained from the mitochondrial mass high fraction. For the purpose of this analysis, the original untreated and cisplatin-treated cells, together with the mitochondria mass low fraction and the TMRM low fraction, were merged and integrated, as shown in Figure 5B. Cluster analysis based on chromatin accessibility identified 34 total clusters (Figure 5C). Detected mtDNA mutations are highlighted in red for the original-cisplatin-treated T cells, purple for the TMRM low fraction, and green for the mitochondria mass low fraction (Figures 5D–5F). Notably, four of the six mtDNA mutations in both the original-cisplatin-treated cells and the TMRM low fraction were also detected in the mitochondria mass low fraction, which included 3850G>A, 5382C>T, 12359C>T, and 6451T>C (Figure 5F). The mutations 3850G>A, 12359C>T, and 6451T>C were detected in the mitochondria mass low fraction but not the mitochondria mass high fraction. 5382C>T was detected in 15 cells (0.03%) in the mitochondria mass high fraction. Of the six original-cisplatin-treated identified mutations, 7378T>C and 8584G>A were not detected in the mitochondrial mass fraction at all. These results suggested that these four mtDNA mutations were associated with both low mitochondrial mass and low mitochondrial activity.
Lastly, the mitochondria mass high and low fractions were merged and integrated with the TMRM high and low fractions, together with the original-cisplatin-treated T cells, and untreated T cells, followed by cluster analysis (Figures S16A and S16B). The three previously identified TMRM high mutations were also detected in the mitochondria mass high fraction, but they were not found in the TMRM low or mitochondria mass low fractions (Figures S16C–S16E). This once again suggested that some cisplatin-induced mtDNA mutations might not damage mitochondrial function.

mtDNA mutations are detected in tumor-infiltrating T cells from colorectal cancer patients after exposure to chemotherapeutic agents
Our results showed that 5-FU and cisplatin could induce mtDNA mutations in T cells. We then hypothesized that T cells residing in tumors might carry mtDNA mutations after patients were treated with chemotherapy. To test this, we utilized the same single-cell mtDNA/ATAC-seq technique to study three tumors from patients with metastatic colorectal cancer who had been exposed to FOLFOX (5-FU and oxaliplatin) chemotherapy before surgery to remove their recurrent tumors. Briefly, after resection, tumor specimens were dissociated and cryopreserved. On the day of the experiment, tumor specimens were thawed, and CD3+ TILs were then sorted and subjected to single-cell mtDNA/ATAC sequencing (Figure 6A).
A patient with colorectal cancer received 12 cycles of FOLFOX therapy and radiation therapy 3 months prior to the resection of liver metastasis. After sorting and single-cell sequencing, 6,413 high-quality TILs were obtained from the tumor specimen. The cluster analysis of TILs identified 10 distinct clusters of T cells (Figure 6B). Additionally, 4,849 high-quality cells were obtained from the peripheral blood mononuclear cell (PBMC) sample that also underwent cluster analysis (Figure 6C). In total, 701 mutations were detected in the TIL sample that were not found in the PBMC sample (Figure 6D). After applying restricted quality-control parameters and other criteria, a total of 24 non-synonymous mtDNA mutations were detected in the TIL sample. The predominant mtDNA mutation identified from the TIL sample was 6516G>A (COI:G250S), which was concentrated in cluster #3, followed by clusters #1 and #5 (Figure 6E); 12014 > T (ND4:L419F) and 4675T>C (ND2:I69T) mutations were mostly found in cluster #1. An mtDNA mutation 9554G>A that was detected in both PBMCs and TILs is shown as a control (Figures 6F and 6G). For identifying gene markers, since clusters #1, #3, and #5 were clustered near each other on the UMAP, the gene profiles for these clusters were similar. We then compared gene markers between cluster #3 and cluster #2 (also #3 vs. #4 and #3 vs. #6). Several T cell exhaustion markers, ENTPD1 (CD39), CXCL13, PDCD1 (PD-1), TOX, LAG3, and CTLA4, were elevated in cluster #3, compared to #2 (Figure S17). Data for statistical analysis are shown in Table S2. We also found elevated cytotoxic/effector T cell markers, including PFR1 (perforin 1), GZMA (granzyme A), IFNG (interferon gamma [IFN-γ]), and ZNF683 (Hobit). Additionally, cluster 3 was associated with diminished T cell markers correlated with younger (naive/memory) T cells, such as IL-7R and transcription factor 7 (TCF7) (Figure S17). Similar observations were also found in clusters #1 and #5 (Figure S17; Table S2). Taken together, mtDNA mutations were identified in TILs associated with increased exhaustion, cytotoxicity, and decreased younger T cell phenotypes.
To further characterize the phenotypes of TILs carrying mtDNA mutations, we performed a projection analysis using scATAC-seq reference clusters, which were previously established by Satpathy et al. using TILs isolated from basal cell carcinoma tumor biopsies (Figure 6H).24 The projection of TILs onto the reference clusters led to a large population of TILs clustering within the intermediate and terminal CD8+ exhausted T cell (Tex) clusters (Figure 6H). T cells with the three identified mtDNA mutations also clustered within the CD8+ Tex clusters (Figure 6J). Mutation 6516G>A predominately was clustered within the terminal Tex cluster. This analysis further validated the exhausted phenotypes of TILs with selected mtDNA mutations.
Similarly, the second patient with colorectal cancer had undergone four cycles of FOLFOX chemotherapy and radiation. A tumor specimen from liver metastasis and blood samples were obtained approximately 6 months after chemotherapy and radiation treatments. A total of 2,873 high-quality cells from TILs were obtained, and 19 main clusters were identified (Figure S18A). Additionally, 5,968 high-quality cells from PBMCs were obtained (Figure S18B). After filtering, a total of nine TIL-specific mtDNA mutations were identified (Figure S18C). Three TIL-specific mtDNA mutations are shown (Figure S18D). The third patient with colorectal cancer was treated with 12 cycles of FOLFOX therapy. A tumor specimen from the recurrent liver metastasis, together with a blood sample, was obtained 15 months after completion of FOLFOX therapy. A total of 2,870 high-quality cells from TILs were obtained, and chromatin accessibility clustering identified 14 distinct clusters of TILs (Figure S18E). Likewise, 4,681 high-quality cells from PBMCs were obtained (Figure S18F). A total of five TIL-specific mutations were identified after filtering (Figure S18G). The detected mtDNA mutations included 3424G>A (ND1:V40M), 9229A>G (COII:Y8C), and 14093T>C (ND5:L586), but we found fewer cells carrying non-synonymous mtDNA mutations compared to the other two donors. (Figure S18H). In summary, we detected TILs with mtDNA mutations in all three colorectal cancer patients who underwent FOLFOX chemotherapy prior to the analysis. A higher number of TILs with mtDNA mutations was detected in patient 1 but less in patient 2 and even less in patient 3. This observation loosely correlated with the time between chemotherapy and the time of tumor sampling, since tumors were obtained 3, 6, and 15 months after the last dose of chemotherapy, respectively. We hypothesized that TILs were eventually replaced with fresh T cells from peripheral blood, and therefore, very few mtDNA mutation-carrying TILs were identified from the third patient. Importantly, a large cohort of patients with various post-chemotherapy time points is required to study the potential correlation between them.

Discussion
The inherent dependence of T cell activation and effector function on the metabolic bioenergetics of the cell has been well established. However, it remains unclear how mitochondrial health may lead to successful cancer immunotherapy. Studies have been looking into mitochondria to potentially remove hurdles and improve T cell function.25,26,27 In an aforementioned study, Rho0 T cells were generated by treating OT-I T cells with low-dose ethidium bromide, leading to depletion of mtDNA.7 Rho0 T cells had high levels of ROS and high expression levels of exhaustion markers PD-1, Tim-3, and TIGIT. In another study, Baldwin JG et al. demonstrated that T cells could receive healthy mitochondria from bone marrow stromal cells through nanotubular connections. T cells with donated mitochondria showed stronger anti-tumor responses and less exhaustion phenotypes, providing a novel approach to revive aging T cells for the use of adoptive cell therapy.28 In a parallel study, a mtDNA mutation at tRNA-Leu (UUA/G) (3290T>C) was identified from a TIL cell line generated from a melanoma specimen. The 3290T>C mutation in mitochondria could be transferred from the melanoma tumor cell line to TIL. Importantly, TIL with transferred mitochondria showed significant reductions of T cell functions in vitro and in vivo.29 Herein, we found that exposing T cells to 5-FU and cisplatin could induce mtDNA mutations. Additionally, several mtDNA mutations in T cells were associated with low mitochondrial activity or changes in gene markers linked to T cell functions. It is critical to validate these identified mtDNA mutations through functional assays. However, gene editing for mtDNA remains technically challenging, suffering from low efficiency and low specificity.30 When the technique of mtDNA editing becomes mature and widely available, it will be interesting to study these identified mtDNA mutations in T cells.
Chemotherapies are employed for their disruption of replication and induction of tumor cell apoptosis, and chemotherapeutic agents can induce characteristic mutational signatures in genomic DNA (gDNA) in tumor cells.31 For instance, both base excision repair and mismatch repair pathways have been implicated in the resistance of tumor cells to 5-FU,32 indicating that 5-FU directly acts on gDNA. Like 5-FU, cisplatin’s effect on gDNA is demonstrated through the resistance pathways that have been identified in cisplatin-exposed tumor cells. Nucleotide excision repair is the most commonly employed repair pathway in cisplatin-treated tumor cells, in addition to mismatch repair and base excision repair.33 Notably, as mtDNA is more susceptible to mutation than gDNA, it is plausible that 5-FU and cisplatin also target mtDNA to induce DNA defects.
Although 5-FU is designed to target tumor cells, 5-FU may also damage healthy cells. Patients treated with 5-FU had higher levels of single-and double-strand breaks in the DNA of their healthy peripheral blood lymphocytes after 5-FU treatment.34 A more recent study found that 5-FU induced gDNA mutations in small intestinal organoid cultures, which were detectable using a whole-genome sequencing technique.35 Furthermore, in a genome-wide mutation analysis of cancer patients treated with 5-FU, a colorectal cancer patient was found to have healthy colon stem cells that had a 5-FU-specific mutational footprint.36 Because 5-FU can induce damage in gDNA, it is possible that 5-FU may also impact the more vulnerable mtDNA. Our findings of 5-FU-treated T cells harboring mtDNA mutations not found in the non-exposed cells align with previous findings of 5-FU inducing specific gDNA mutations. A similar study performing scRNA-seq of intestinal organoids found that one of the most prominent pathways affected by 5-FU treatment was mitochondrial ATP synthesis.37 Our findings suggested the involvement of mtDNA mutations in the OXPHOS pathway, for example, NADH dehydrogenase, cytochrome c oxidase, and ATP synthase, and their potential impacts on chemotherapeutically treated T cells. However, with the current available data and publications, we could not exclude the possibility that 5-FU and cisplatin might reduce OXPHOS through other pathways, independent of mtDNA mutations.
In this study, we identified numerous mtDNA mutations that were only found in T cells treated with 5-FU or cisplatin. This indicated that some chemotherapeutics might induce mtDNA mutations and affect T cell activities. It is possible that 5-FU or cisplatin induced an mtDNA mutation in a single T cell, and this T cell might expand and become a clonal population carrying this specific mtDNA mutation. As a result, TCR clonotypes might link to mtDNA mutations identified in this study. However, none of the current single-cell techniques enable us to investigate this possibility. To address this, we plan to develop a single-cell technique that can enable us to obtain TCR clonotype, ATAC data, and mtDNA sequence from the same cell.38,39 Alternatively, it is also possible that exhausted T cells are more sensitive to mtDNA damage. This hypothesis also needs to be tested in the future. Furthermore, to understand the impact of mtDNA mutations on cells, a study examining heteroplasmic shifts in an epithelial cell population found that selection pressure for a non-synonymous mtDNA mutation depended on the cellular environment, whether that be nutrient deprivation, hypoxia, or ATP synthase inhibition.40 Therefore, the effect of an mtDNA mutation, whether it was harmful, neutral, or beneficial to the cell, depended on their environment. Similarly, another study proposed that T cells might undergo purifying selection and reduce the heteroplasmy for certain non-synonymous mtDNA mutations.41 Additionally, T cell subpopulation and metabolic state might influence the purifying selection of non-synonymous mtDNA mutations.42
To improve the efficacy of immunotherapy, combining chemotherapy or radiation with immunotherapy has been developed in recent years.43 The Food and Drug Administration (FDA) has approved nivolumab with platinum-doublet chemotherapy as neoadjuvant treatment for NSCLC based on the results of the phase III study CHECKMATE-77T (NCT04025879).44 Additionally, pembrolizumab has been approved in combination with chemotherapy, paclitaxel plus carboplatin, as a therapy for advanced endometrial cancer, based on a clinical trial showing a significantly longer progression-free survival compared to chemotherapy alone (13.1 vs. 8.7 months).45 In this trial, chemotherapy and pembrolizumab were administered concurrently on day 1 of a 21-day cycle, for up to six cycles of treatment, followed by a maintenance treatment of pembrolizumab only. Similarly, the FDA has approved two PD-L1 inhibitors, atezolizumab and durvalumab, plus chemotherapy for first-line treatment of extensive-stage small-cell lung cancer.46,47 In both studies, chemotherapy was administered on day 1 concurrently with PD-L1 inhibition, followed by a maintenance treatment of PD-L1 inhibition only. Importantly, although atezolizumab plus carboplatin and etoposide chemotherapy showed a statistically significant overall survival compared to chemotherapy alone, the improvement was only 2 months.46 Lastly, a clinical trial assessed the efficacy of pembrolizumab plus etoposide and platinum (carboplatin or cisplatin) compared to chemotherapy alone as a first-line therapy for extensive-stage small cell lung cancer.48 Similar to the studies described above, chemotherapy was administered concurrently with pembrolizumab for four cycles of treatment, followed by maintenance cycles of pembrolizumab only. A significant improvement in 12-month progression-free survival was observed (13.6% vs. 3.1%); however, there was not a significant difference in median overall survival (10.8 vs. 9.7 months). These moderate increases in overall survival suggest the potential for improvement in the combination therapy.
In all of the above-mentioned studies, chemotherapy and ICB therapy were administered concurrently, followed by a maintenance treatment regimen of ICB therapy alone. Some evidence of clinical efficacies was shown, which led to the FDA approvals. However, it remains largely unknown about the impact of chemotherapy on ICB therapy, other than the additive effects on anti-tumor responses. Multiple studies have raised concerns that chemotherapy may damage T-cell-mediated responses. For example, Mariniello et al. demonstrated that chemotherapy (cisplatin + pemetrexed) could attenuate the CD8+ T cell responses to concomitant PD-1 blockade in a chronic LCMV infection model.49 The experimental results suggested that sequential administration of chemotherapy and PD-1 blockade might be beneficial for viral control in vivo. Additionally, a previous study showed that T cells exposed to the chemotherapeutic agents, cyclophosphamide, doxorubicin, or cytarabine, could have lingering effects on T cells, such as impaired proliferation and mitochondrial functions. These results suggested the potential impact on T cell manufacture and clinical efficacy for the adoptive cell therapy.50 In this study, we demonstrated that common chemotherapeutics 5-FU and cisplatin could have a direct effect on mtDNA, which might impact T cell function and potentially alter the effectiveness of immunotherapy. Our findings suggest using caution when administering certain chemotherapeutics together with checkpoint immunotherapy. Avoiding chemotherapy prior to immunotherapy may optimize immunotherapy efficacy as an avenue to improve clinical outcomes. Taken together, our results underline an underappreciated consequence of chemotherapy on T cell mtDNA, with implications for the design of combination therapies. Future studies should evaluate whether alternative dosing schedules or mitochondria-targeted protectants could mitigate these effects, ultimately enhancing the efficacy of immunotherapy.

Limitations of the study
This study has multiple limitations due to the biological and technical barriers. First, although mitochondrial heteroplasmy is a fascinating phenomenon in mitochondrial biology, it creates significant complexities and challenges for standardized data analysis and interpretation. Second, direct functional studies of mtDNA mutations cannot be performed because direct mtDNA editing in T cells is currently unavailable. Future studies will be needed to validate the impacts of these mtDNA mutations in T cells. Third, it is possible that certain T cell phenotypes, such as exhaustion, are more sensitive to mtDNA damage. Future investigations will be required to further study the link between mtDNA mutations and T cell phenotypes. Fourth, the relationship between mtDNA mutations and TCR clonotypes cannot be investigated due to the lack of a suitable single-cell technique. Once the technique becomes available, the relationship between them needs to be studied and validated. Lastly, a large-scale study is required to study the correlation between the timing of chemotherapy treatments and the quantity of mtDNA mutations in TILs.

Resource availability

Lead contact
Requests for further information or resources should be directed to the lead contact, Yong-Chen William Lu (ylu@uams.edu).

Materials availability
This study did not generate new unique reagents.

Data and code availability

•Single-cell RNA-seq data and mtDNA/ATAC-seq data are publicly available at the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) with the GSE number (GSE310683).

•This paper does not report original code. The bioinformatics for the single-cell data was performed using open-source tools, as listed in the key resources table.

•Any additional information required to reanalyze data reported in this work is available from the lead contact upon request.

Acknowledgments
The authors would like to thank Bing Guan, Andrea Harris, Sue L Johnson, Donald J Johann Jr., and Emilio Murillo for their suggestions and technical support. Next-generation sequencing was conducted at the UAMS Genomics Core and was analyzed by the computational resources of UAMS High Performance Computing cluster. This work was supported by the Winthrop P. Rockefeller Cancer Institute. C.K. was also supported by the 10.13039/100008519UAMS MD/PhD program. R.P.M.D was supported by 10.13039/100000054National Cancer Institute, 10.13039/100000002NIH, United States (R01 CA245083). Y.L. was also supported by the Arkansas Breast Cancer Research Program (ABCRP), United States; the 10.13039/100015317Translational Research Institute, 10.13039/100008519UAMS, United States (UL1 TR003107 and KL2 TR003108); and 10.13039/100000054National Cancer Institute, NIH, United States (R01 CA143130).

Author contributions
Conceptualization, C.K. and Y.-C.W.L.; data curation, C.K.; formal analysis, C.K. and Y.-C.W.L.; funding acquisition, Y.-C.W.L.; investigation, C.K. and Y.-C.W.L.; methodology, C.K., Y.-K.L., R.P.M.D., and C.M.P.; project administration, Y.-C.W.L.; resource, C.M.Q., S.R.P., Y.-K.L., R.P.M.D., C.M.P., and L.B.; supervision, Y.-C.W.L.; validation, C.K. and Y.-C.W.L.; visualization, C.K.; writing (original draft, review and editing), C.K., C.M.P., and Y.-C.W.L.

Declaration of interests
The authors declare no competing interests.

Declaration of generative AI and AI-assisted technologies in the writing process
No AI-assisted technologies have been used in this study.

STAR★Methods

Key resources table

Experimental model and study participant details

Healthy human donor cells
Peripheral mononuclear cells (PBMCs) from four healthy human donors were purchased from Stemcell Technologies. The age range of donors was 23–34 years old, with a median age of 28.5 years. There were 2 female and 2 male donors. No known health, immune status, or previous treatments are available for the healthy donors.

Colorectal cancer patient samples
The tumor and blood specimens from patients with metastatic colorectal cancer were obtained through the University of Arkansas for Medical Sciences Tissue Biorepository and Procurement Service (TBAPS), with the Institutional Review Board not human subject research determination (IRB#262500 and 273430). Written informed consent was obtained from patients. The NIH guidelines and the Declaration of Helsinki were also followed. Three patient samples were obtained, from 1 female and 2 males. The age of the patients were 49, 61, and 65, respectively. All patients had undergone treatment cycles of FOLFOX chemotherapy regimens prior to surgical resection of tumor samples. Because this is not a large-scale study, the generalizability for gender is not discussed, which is part of the limitations for the scope of this study.

Method details

Inducing mtDNA mutations in T cells with chemotherapeutic agents
2 x 106 PBMCs were thawed and stimulated with 25 μL of Dynabeads human T-activator CD3/CD28 (Gibco) in 1 mL RPMI medium containing 5% human serum and IL-2 (600 IU/mL) in a 24-well plate, according to manufacturer’s protocol. On day 3, 1 mL RPMI-5 medium containing 600 IU/mL IL-2 and various concentrations of 5-FU or cisplatin were added to the culture until day 7. If necessary, proliferating T cells were fed and split with RPMI medium containing IL-2 and 5-FU or cisplatin. On day 7, culture medium was carefully removed and replaced with fresh RPMI medium containing IL-2, without 5-FU or cisplatin. On day 10, T cells were harvested and counted. 1 × 106 T cells were re-stimulated with CD3/CD28 Dynabeads, followed by the same 5-FU or cisplatin treatment from day 13 to day 17. On day 20, treated T cells were harvested and re-stimulated. The T cell stimulation and 5-FU/cisplatin treatments were performed for a total of 3 cycles. T cells were cryopreserved for the downstream analysis.

Library preparation for single-cell mtDNA and ATAC (assay for transposase-accessible chromatin) sequencing
Treated T cells were thawed and stained with 7-AAD (BD Bioscience), followed by FACS sorting for live T cells (FACSAria II, BD Bioscience). For TMRM and MitoTracker staining, treated T cells were thawed and recovered two days in RPMI medium containing 5% human serum and IL-2 (600 IU/mL). T cells were stained with TMRM (20 nM final concentration) or MitoTracker Green (50 nM) for 30 min at 37°C, according to the manufacturer’s protocol (Molecular Probes/Invitrogen). After washing, stained cells were stained with 7-AAD, followed by FACS sorting for the TMRM or MitoTracker Green FM high and low fractions (approximately 10–15% of the T cell populations). After sorting, T cells were fixed and lysed, for the preparation for single-cell ATAC and mtDNA sequencing. This method has been established by Lareau CA et al. and published previously,14 with minor modifications in our laboratory.20 Briefly, to modify the scATAC-seq protocol to include mtDNA, T cells were fixed with 1% formaldehyde/PBS for 10 min. Fixed cells were lysed, washed, and counted, followed by performing the standard scATAC-seq sample preparation (10X Genomics). In the last step, the mtDNA/ATAC libraries were enriched based on DNA size from 180 b.p. to 600 b.p. (Pippin Prep instrument, Sage Science). The final libraries were subjected to deep sequencing (NovaSeq 6000, illumina).

Library preparation for single-cell RNA sequencing (scRNA-seq)
Treated T cells were thawed and stained with 7-AAD (BD Bioscience), followed by FACS sorting for live T cells (FACSAria II, BD Bioscience). After sorting, T cells were prepared for single-cell analysis and sequencing. A single-cell 5′ reagent kit was used by following the manufacturer’s protocol (10X Genomics). Briefly, 10,000 T cells per channel were loaded on the Chromium Controller, with a targeted cell recovery of 6000 single cells. The pooled single-cell cDNA samples were first universally amplified by a 15-cycle PCR reaction, followed by fragmentation and barcoding. The gene expression libraries were then sequenced by an Illumina Novaseq 6000 sequencer.

Isolation of tumor-infiltrating T cells
Tumor specimens were dissociated by a gentleMACS dissociator (Miltenyi Biotec) and cryopreserved.52 Briefly, the tumor specimens, approximately 1–3 g, were cut into small pieces (2–3 mm) and transferred to gentleMACS C Tubes (Miltenyi Biotec, Germany). The digestion medium contained 10 mL RPMI medium, 5% human serum, 0.13 U/mL Liberase TM (Roche/Sigma-Aldrich), and 600U/mL DNase I (Roche/Sigma-Aldrich). The tumor specimens were dissociated by selecting the following programs: h_tumor_01, h_tumor_02, and h_tumor_03 programs. The C Tubes were incubated at 37°C for 30 min between each program. At the end of the programs, the single-cell suspensions were passed through a 40 μm cell strainer and washed once with PBS containing 5 mM EDTA and then cryopreserved. On the day of the experiment, tumor specimens were thawed and stained with Fc Block and 7-AAD, followed by anti-CD3 and anti-CD8 antibody staining (BD Biosciences). CD8 staining was included to ensure the correct gating of CD3+ T cells. Total CD3+ T cells, including CD8+ and CD8− cells, were sorted and subjected to the same single-cell mtDNA/ATAC sequencing analysis.

Quantification and statistical analysis

Single-cell mtDNA/ATAC-seq
The combination of chromatin accessibility and mitochondrial genome sequencing data was first processed by a Cell Ranger pipeline (CellRanger-ATAC, v2.1.0; 10X Genomics). The chromatin accessibility sequencing data was mapped to a hard-mask modified reference genome database (GRCh38-r107) in order to mask regions of nuclear genome that share consistent homology with the mtDNA, by following the instruction provided by mgatk.15 In some datasets, the optional “force-cells” command was used to remove low-quality cells. On average, 669,865,207 sequenced read pairs were obtained (range 182,592,167 to 1,503,707,473), with 16,897 median high-quality fragments per cell (range 4,213 to 32,095). The average library complexity was at 66% (percent duplicates range 44.9%–84.6%).
The single-cell mtDNA genotyping and heteroplasmy data were analyzed by mgatk (v0.6.6). Considering the complexity of mtDNA heteroplasmy, we streamlined the analysis and defined treatment-induced mtDNA mutations as detected in the treatment samples but not in the untreated negative controls (Figure 2D).20 Similarly, TIL-specific mtDNA mutations were defined as they were detected in the TIL, but not in the autologous PBMC samples. Next, restricted quality-control parameters, the same as the Lareau CA et study, were applied to remove low-quality mtDNA mutations and potential sequencing artifacts.14 For example, mtDNA mutations with strand concordance values <0.65 and VMR (variance mean ratio) < 0.01 were removed. The filtering was performed using the default filter settings in the mitochondrial genotyping workflow in Signac. Lastly, we selected non-synonymous mtDNA mutations with ≥10% heteroplasmy in ≥10 cells for the downstream analysis and visualization.
The scATAC-seq data were analyzed by Signac (v1.14) and Seurat (v5.3) in R (v4.5.0). The standard mitochondrial genotyping workflow was used to integrate scATAC-seq and mtDNA datasets. Low-quality cells and potential doublets were filtered out using the following criteria: nCount_peaks >10000 & nCount_peaks <50000, pct_reads_in_peaks >15, TSS enrichment >2, mtDNA depth ≥ 10.The filter for nCount_peaks was set at < 40000 for the following sample pairs: D1007-5-FU, D1007-cisplatin, D2008-cisplatin, TILs (P1, P2, P3). Additionally, the TSS enrichment was set at >4 for TILs (P1, P2, P3). The average mtDNA coverage after filtering was approximately 60 (range 30–100). Next, the standard merging objects workflow in Signac was used for merging and integrating multiple samples. Unified peak sets were created using the Reduce function from the GenomicRanges package. Combined datasets were then normalized and clustered to obtain UMAP plots and violin plots. Clustering analysis was performed at the fixed resolution of 2.0 for all combined datasets. The resolution was set at 0.5 for patient TIL and PBL datasets. Clusters with fewer than 100 cells were excluded from analyses. Due to the higher numbers of cells, clusters with <220 cells were excluded for Figures S10–S12. To analyze CD4+ and CD8+ T cell populations separately, normalization and clustering were performed following the standard pipeline described above. Clusters were classified as CD4+ or CD8+ based on their CD4/CD8A signal ratios. Because CD8A showed higher signal than CD4 in general, clusters with a CD8A/CD4 signal ratio ≥4 were designated as CD8+ clusters, whereas clusters with CD8A/CD4 ratio <1 were designated as CD4+ clusters. Clusters with ratios between 1 and 4 were considered double-positive or ambiguous and excluded from downstream analyses. Additionally, clusters lacking CD3E signals were excluded. The resulting CD4+ and CD8+ clusters were separated into distinct datasets and then underwent re-normalization and re-clustering for visualization in UMAP plots.
Comparing differences in chromatin accessibility across cell clusters was performed through the standard mitochondrial genotyping workflow. p-values for cluster gene comparisons can be found in the supplementary figure legends and Table S2, with significant genes defined as a p-value <0.05 and a log2 fold change >0.58 or −0.58 (equal to 50% increase or decrease).

scRNA-seq
The single-cell RNA sequencing data were first processed by a Cell Ranger pipeline (v7.1.0; 10X Genomics) and mapped to the standard reference genome database GRCh38-2020-A, 10X Genomics). The following parameters were used to filter out potentially low-quality cells and doublets (nFeature_RNA >200 & nFeature_RNA <5000 & nCount_RNA <20000 & percent.mt < 6). After filtering, the standard scRNA-seq integration workflow in Seurat was used to integrate scRNA-seq datasets from 3 donors. SelectIntegrationFeatures, FindIntegrationAnchors, and IntegrateData were used to integrate the datasets, followed by normalization, PCA analysis, and clustering to obtain UMAP and violin plots. Significant genes were defined as a p-value <0.05 and a log2 fold change >0.58 or −0.58.

Trajectory analysis
Trajectory analysis of the scRNA-seq samples was performed following the standard Monocle3 (v.1.4.26) workflow. Samples were first filtered and processed into a Seurat object using the following filtering parameters: nFeature_RNA >200 & nFeature_RNA <6000 & nCount_RNA <20000 & percent.mt < 6. The datasets were subsequently processed, integrated, and clustered using the standard Monocle3 pipeline. For pseudotime analysis, root nodes were defined based on peak MKI67 expression.

Projection analysis
Projection analysis of tumor samples was performed using the reference clusters established by Satpathy et al. based on the processed tumor-infiltrating T cell scATAC-seq dataset obtained from basal cell carcinoma (BCC) samples.24 Because this dataset was aligned to hg19, genomic coordinates were converted to hg38 using the UCSC liftOver tool. This dataset was then processed as described by Satpathy et al., using the same normalization and manual singular value decomposition (SVD) parameters to recreate the reference clusters for our projection analysis. The Signac Reference Mapping workflow was followed with minor modifications due to software updates. Specifically, the FindTransferAnchors function in Signac requires a Seurat-calculated latent somatic indexing (LSI) reduction. To accommodate this, a Seurat-style LSI was built using the same reference features defined by Satpathy et al., therefore, the query dataset could be projected onto the reference’s latent space. This allowed the FindTransferAnchors function to identify shared structure and establish mapping anchors between the reference and query datasets. The original reference clusters, generated using the manual SVD workflow, were preserved for use with MapQuery to project our query cells into the same UMAP embedding as the reference. The resulting projection clusters closely approximated the original coordinates due to minor differences between the Seurat-generated LSI and the manually calculated reference SVD.

Pathway analysis
Pathway analysis was performed using the Gene Ontology (GO) Enrichment Analysis tool (release version 2025-10-10).53,54,55 For each mtscATAC-seq cluster of interest, the top 100 upregulated and the top 100 downregulated genes, ranked by log2 fold change, were input for GO enrichment analysis to identify enriched pathways associated with T cell function.

Mutation comparison
Comparison of mtDNA mutations between samples was performed using Venny 2.1 and visualized using DeepVenn.