Melanoma and SCC have shared etiologies (e.g., both are UV-induced cutaneous conditions) and similar immune-oncology-based therapeutic approaches are adopted in treating them, but these conditions can progress very differently, lending themselves to the interrogation of CD8 T cell exhaustion and analysis for common or distinct mechanisms of dysfunction with applications for the development of improved therapies. CD8⁺ TILs in melanoma and SCC share important common features, including the coordinated expression of immune checkpoint receptors and the development of resistance to immunotherapy (29, 31, 32), yet the mechanisms underlying TIL dysfunction in these tumors remain elusive. To address this, we performed transcriptomic profiling on flow-sorted CD8⁺PD-1⁺TIM-3⁺ T cells from melanoma and SCC (total CD8⁺/P⁺/T⁺ TILs) and compared them to peripheral CD8⁺ T cells (pCD8s) from healthy donors. A total of 3,957 genes were differentially expressed (p ≤ 0.05) in total CD8⁺/P⁺/T⁺ TILs vs. pCD8⁺ T cells ( Figure 1A ). Gene set variation analysis (GSVA), beginning with gene ontology (GO) molecular function annotation, identified the corresponding top upregulated pathways, including chemokine activity and three pathways related to metabolism: carbohydrate kinase activity, oxidoreductase activity II, and aspartic-type peptidase activity ( Figure 1A ). The top downregulated pathways included C2H2 zinc finger domain binding, p53 binding, rRNA binding, ARF guanyl nucleotide exchange factory activity, and androgen receptor binding ( Figure 1A ). Clear differential expression of genes in total CD8⁺/P⁺/T⁺ (SCC and MEL) TILs vs. pCD8s indicated upregulation of immune checkpoint receptors including TIM-3 (HAVCR2), PD-1 (PDCD1), LAG3, and TIGIT ( Figure 1B ).

Moreover, additional pathway enrichment analysis via gene set enrichment analysis (GSEA) using the Hallmark database (MSigDB) identified significant induction of multiple metabolic pathways in total CD8⁺/P⁺/T⁺ TILs, including glycolysis, xenobiotic metabolism, fatty acid metabolism, mTORC1 signaling/PI3K AKT MTOR signaling, adipogenesis, and bile acid metabolism, corroborating our GSVA data ( Figure 1C ). We performed an analysis of pathway interactivity (via shared differentially expressed genes (DEGs) according to Jaccard similarity coefficient) to provide a broader overview of how the enriched pathways interact; as expected, this indicated a close correlation between mTORC1 signaling, PI3K AKT MTOR signaling, and glycolysis, with a less-direct connection to bile acid metabolism (via mTORC1 signaling) ( Figure 1D ). We were especially intrigued by the significant differential enrichment of bile acid metabolism, which has not been proposed as a key regulatory pathway in TIL function. Linear regression analysis of PD-1 (surface expression determined by flow cytometry) with gene expression revealed a strong positive correlation between PD-1 protein expression and upregulation of the bile acid metabolism, peroxisome, and glycolysis pathways in CD8⁺ TILs and peripheral blood mononuclear cells (PBMCs) ( Figure 1E ). These findings suggest that CD8⁺ T cells adopt a unique metabolic program associated with the latest stages of exhaustion when PD-1 expression is at its highest.

Moreover, core biological pathways that have been previously shown to be associated with TGF-β attenuation of the tumor response to immune checkpoint blockade were also significantly differentially regulated within the dysfunctional CD8⁺PD-1⁺TIM-3⁺ melanoma and SCC TILs ( Supplementary Figure 1A ).

To objectively investigate additional exhausted T cell phenotypes and associated metabolic signatures, we integrated the expression of immune checkpoint receptors (ICRs) and activation markers with the transcriptional profiles of CD8⁺ TILs (melanoma and SCC) and peripheral CD8⁺ T cells ( Figure 2 ). We used flow cytometric analysis of a smaller, randomly selected subset of samples to measure the expression of immune checkpoint receptors PD-1, TIM-3, TIGIT, BTLA, LAG3, and CD38, which have all been shown to regulate T cell function (33). Coordinate expression of multiple ICRs is associated with progressive CD8⁺ TIL dysfunction; therefore, we used our flow cytometry data to infer levels of CD8⁺ TIL exhaustion.

We employed uniform manifold approximation and projection (UMAP), a non-linear dimensionality-reduction machine learning approach, to analyze our multiparametric flow data and objectively determine the heterogeneity of CD8⁺ T cells based on surface expression of BTLA, CD3, CD8, CD27, CD38, CD45RA, CCR7, LAG3, PD-1, TIGIT, and TIM-3. UMAP analysis revealed 17 unique clusters, nine of which showed significantly different frequencies in CD8⁺ TILs compared to PBMCs ( Figure 2A ). A two-way hierarchically clustered heatmap of MFI (geometric mean fluorescence intensity) based on surface expression shows the distinct surface signature profiles of the 17 clusters identified by UMAP analysis ( Figure 2B ). Notably, the frequencies of lymphocytes co-expressing multiple immune checkpoints (CD38, TIGIT, LAG3, PD-1, and TIM-3), the negative regulators of T cell function, were enriched in CD8⁺ TILs compared to peripheral cells, especially in clusters 6, 15, and 16.

Of note, cluster 6 represents what is likely to be the most highly dysfunctional CD8 population, characterized by relative upregulation of CD38, TIGIT, PD-1, LAG3, and TIM-3, and is significantly enriched in the CD8⁺ TIL population compared to the peripheral CD8⁺ T cells ( Figures 2B, C ). Similar to cluster 6, clusters 10, 15, and 16 have similar ICR expression patterns with increased coordinate immune checkpoint expression levels ( Figure 2B ). Interestingly, cluster 16 is significantly enriched in CD8⁺ TILs of both melanoma and SCC with high expression of PD-1 and TIM-3 ( Figures 2B, C ). Likewise, cluster 15 is also enriched in CD8⁺ TILs, but only in melanoma, with moderate immune checkpoint expression ( Figures 2B, C ). Our findings indicate that clusters 6, 10, 15, and 16, with upregulation of ICRs, are likely to represent progressively more exhausted CD8⁺ T cell clusters ( Figure 2B ).

We then performed linear regression modeling of gene expression with frequencies of cells from each flow UMAP cluster in each donor to identify transcriptomic pathway enrichment that correlates significantly with frequencies of CD8⁺ T cell clusters identified by UMAP. Interestingly, higher frequencies of cells in cluster 6 were strongly correlated with positive enrichment of the bile acid metabolism and peroxisome pathways ( Figure 2C ). Likewise, cluster 16 frequency was also positively correlated with enrichment of the peroxisome, oxidative phosphorylation, and glycolysis pathways ( Figure 2C ). In summary, our data demonstrate a significant correlation between the upregulation of immune checkpoint-expressing dysfunctional CD8⁺ TILs and activation of metabolic functions, particularly the bile acid metabolism, peroxidase, and glycolysis pathways, at the transcriptomic level. Targeting this metabolic-associated dysfunctional CD8⁺ Tex/act cell population could reinvigorate exhausted CD8⁺ T cells, thus enhancing immunotherapy responses in cancer patients.

To determine whether these metabolic pathways were specifically dysregulated in melanoma, we performed differential gene expression analysis of the melanoma CD8⁺ TILs versus healthy peripheral CD8⁺ T cells, resulting in a distinct gene expression signature as shown by the heatmap in Figure 3A . A total of 4,280 genes were differentially expressed in MEL CD8⁺ TILs versus peripheral CD8⁺ T cells, while 2,057 of those genes were commonly differentially regulated in MEL and SCC TILs (p ≤ 0.05, Supplementary Figure 1B ). Statistically, the most upregulated pathways based on GO molecular function annotation (via GSVA) were the oxidoreductase, structural constituent of muscle, and carbohydrate-binding pathways ( Figure 3A ). The top downregulated pathways included C2H2 zinc finger domain binding, N-acetyltransferase activity, rRNA binding, peptide N-acetyltransferase activity, and transcription coactivator activity. The top 50 differentially regulated genes between MEL TILs and pCD8s by statistical significance are shown in Figure 3B , and the ICRs PD-1 (PDCD1), TIM-3 (HAVCR2), CTLA4, and LAG3 were among these. Also significantly upregulated were FABP5 (epidermal fatty acid binding protein), which can act as an intracellular receptor that binds to free fatty acids (FFA) in the cytosol, and Ki-67 (MKI67), indicating enhanced proliferation of the CD8⁺ TILs.

To further evaluate the biology underlying differentially expressed genes found in CD8⁺ TILs and PBMCs in melanoma, we performed pathway enrichment analysis of key canonical pathways. Notably, the peroxisome, bile acid, and fatty acid metabolism pathways were coordinately regulated in melanoma TILs compared to pCD8s ( Figure 3C ). Not only are peroxisome proliferator-activated receptor (PPAR) pathways associated with FABP5, which was upregulated at the RNA level in dysfunctional CD8⁺ TILs, but they also regulate bile acid synthesis (34–36). Our data showed strong positive enrichment of bile acid metabolism pathways in dysfunctional CD8⁺ T cells from melanoma tumors. CYP27A1, a key regulatory enzyme of the bile acid metabolism pathway, was significantly upregulated in the CD8⁺ T cells of melanoma TILs as compared to PBMCs ( Figure 3A ). Upregulation of additional genes from the bile acid and peroxisome signaling pathways was also evident in TIL gene and pathway expression analysis, including upregulation of SOD1, ACSL5, FADS2, CAT, DHCR24, SLC27A2, and IDH2 in melanoma TILs compared to pCD8s ( Supplementary Figure 2 ).

As shown in Figure 3C , pathway enrichment analysis showed that the mTORC1 signaling pathway was significantly and positively enriched in CD8⁺ melanoma TILs as compared to pCD8s. The glycolysis pathway was highly integrated with many other pathways, particularly the mTORC1 and peroxisome pathways, which interact with bile acid metabolism ( Figure 3D ). Compared to pCD8 T cells from control individuals, the glycolysis pathway was the most significantly upregulated pathway under pathway interaction network analysis in the Hallmark canonical pathway database in CD8⁺ TILs from melanoma patients ( Figure 3D ). The differential regulation of multiple metabolic pathways, including bile acid metabolism, mTORC1, and glycolysis ( Figure 3D ), may reflect the fact that the CD8⁺ TIL cells are hosted in a hypoxic, acidic, and nutrient-depleted TME where the cellular metabolism has been reprogrammed as an adaptive mechanism to facilitate their proliferation and survival in melanoma. Our western blot analysis confirmed increased expression of CYP27A1, which was upregulated in the bile acid metabolic pathway, and three additional glycolytic genes (PKM2, PHGDH, and PCK2) in melanoma CD8⁺ TILs compared to PBMCs ( Figure 3E ).

A subset of our cohort of melanoma patients were treated with anti-PD-1 immunotherapy. There were 1,072 differentially expressed genes (p ≤ 0.05) in patients receiving therapy (n = 5) as compared to those who did not receive treatment (n = 8), as summarized in the volcano plot in Figure 3F . CXCR6 was upregulated in the melanoma anti-PD-1 immunotherapy signature, indicating a potential shift toward a proinflammatory tumor microenvironment, which could lead to subsequent metastasis. CTTNBP2NL (CTTNBP2 N-terminal-like protein) was also upregulated; this is known to be associated with STRIPAK complexes, which have been broadly linked to metabolism, immune regulation, and cancer tumorigenesis (37). On the other hand, CD160 and HRH2 (histamine receptor H2) were downregulated in the anti-PD-1 therapy signature. To identify common features of immunotherapy treatment across multiple cancer types, we performed a similar analysis including patients with SCC and melanoma cancer types and multiple immunotherapy treatments. The comparison of patients receiving immunotherapy (n = 16) and those who did not (n = 14) revealed 1,118 differentially expressed genes (p ≤ 0.05). CYP27A1 and ABCA1 were among those genes significantly downregulated in patients receiving immunotherapy. Consistent with the anti-PD-1 therapy signature, we also saw downregulation of CD160 and HRH2 ( Supplementary Figure 3A ). Notably, pathways involved in fatty acid metabolism, mTORC1 signaling, glycolysis, and xenobiotic metabolism were all downregulated in the immunotherapy signature ( Supplementary Figure 3B ).

The mTORC1 signaling pathway was highly significantly enriched in tumors as compared to peripheral blood ( Figures 1C , 3C ). Using intracellular flow cytometry, we conducted further analysis to discover that two of the key targets of mTORC1 complex activation, S6 ribosomal protein (target of S6 kinase, which is a direct target of mTOR phosphorylation) and 4EBP1 (a direct target of mTOR phosphorylation), were significantly more phosphorylated in the MEL TILs compared to peripheral CD8 T cells, validating our findings at the phosphoprotein level ( Figure 4A ). Furthermore, treatment with rapamycin (sirolimus), which specifically blocks MTORC1 signaling during TCR (anti-CD3/CD28) stimulation of melanoma TILs, resulted in consistent downregulation of PD-1 expression, particularly in the CD8⁺CD45RA⁻ memory T cells and CD8⁺CD45RA-CCR7⁻CD27⁻ effector memory T cells in the tumors ( Figure 4B ).

Moreover, to validate the differential enrichment of the mTOR pathway, we analyzed the single-cell transcriptomes of CD8⁺ TILs from a new cohort of metastatic melanoma patients collected prior to anti-PD-1 immunotherapy. Single-cell analysis of CD8 TILs also allowed us to compare distinct T cell populations within the tumor environment. UMAP dimensional reduction of the single-cell RNA sequencing (scRNAseq) dataset revealed a total of nine clusters ( Figure 4C ). We observed distinct clusters of cells with varying levels of T cell exhaustion. We observed high levels of T cell exhaustion markers PDCD1 (PD-1), HAVCR2 (TIM3), LAG3, CXCL13, and TOX in clusters 2, 4, and 8. Cells in these clusters expressed low levels of genes associated with non-exhausted T cells (TCF7 (TCF1), CCR7, SELL, and IL7R) ( Supplementary Figure 4A ). Cluster 8 demonstrated positive enrichment of the mTORC1 and PI3K/AKT/MTOR signaling pathways at the single-cell resolution, validating our bulkseq findings ( Figure 4D ). Violin plots of several genes associated with cluster 8 (PHGDH, E2F1, PSMA4, BUB1, CDK2, MKI67, IDH2, and FABP5) are shown in Figure 4E (mTOR related) and Supplementary Figure 4B , all of which were also differentially regulated in our bulk RNAseq signature of dysfunctional CD8⁺ TILs. We also confirmed PHGDH expression at the protein level ( Figure 3E ). Cluster 4 also displayed an exhausted T cell phenotype as well as enrichment of key metabolic pathways, including peroxisome, cholesterol homeostasis, and fatty acid metabolism ( Supplementary Figure 4C ). These results support our transcriptional analysis of flow-sorted MEL and SCC CD8⁺PD-1⁺TIM-3⁺ TILs and pCD8s ( Figure 1 ) as well as our findings comparing MEL CD8+ TILs and pCD8s ( Figure 3 ). Overall, our data suggest that CD8+ TIL exhaustion is associated with dysregulation of multiple pathways that warrant further investigation, including peroxisome, bile acid, and fatty acid metabolism, cholesterol homeostasis, and mTORC1 signaling.

Samples were obtained at the Cleveland Clinic under a protocol approved by the Cleveland Clinic’s institutional review board, and written informed consent was obtained from each patient. Peripheral blood lymphocytes (PBLs) and tumor specimens were obtained from patients with cutaneous melanoma (MEL, n = 23), SCC (n = 8), and control individuals without skin disease (n = 8), as previously described (29). Of the MEL patients, 60.9% were male, with a mean age of 58.3 ( ± 15.8); three patients presented with a primary tumor, while 20 patients had metastatic disease. Fourteen of the MEL patients received some type of immunotherapy (IO); out of these patients, eight received anti-PD-1/PD-L1 directed IO (either pembrolizumab or nivolumab). Of the SCC patients, 85.7% were male, with a mean age of 73 ( ± 7.9); three patients presented with a primary tumor, while five patients had metastatic disease. Age was unknown for one SCC patient. Tissue for scRNAseq libraries was obtained from patients with cutaneous melanoma (n = 8), of whom 37.5% were male, with a mean age of 60 ( ± 12.7). Two patients presented with a primary tumor, while six patients had metastatic disease. All tumor specimens came from different individuals. Clinical data are summarized in Supplementary Table 1 .

After surgical resection, tumor specimens were washed with antibiotic-containing media and minced with crossed scalpels under sterile conditions. Tissue was dissociated via enzymatic digestion using 1,500 U/ml collagenase IV (Gibco, Grand Island, NY, USA/Life Technologies, Carlsbad, CA, USA), 1,000 U/ml hyaluronidase (Sigma, St. Louis, MO, USA), and 0.05 mU/ml DNase IV (Gibco) in Roswell Park Memorial Institute (RPMI) medium for 1 h at 37°C followed by mechanical agitation. Centrifugation over a Ficoll–Hypaque gradient was used to separate debris from the single-cell suspension. Finally, cells were cryopreserved in 10% dimethyl sulfoxide (DMSO) + bovine serum. Similarly, peripheral blood mononuclear cells were purified from buffy coats by centrifugation over a Ficoll–Hypaque gradient, followed by cryopreservation.

Flow cytometric analysis of immune checkpoint receptor and activation marker expression and T cell memory subset distribution was carried out as follows. For memory subset phenotyping and assessment of negative regulator and ligand expression, a cocktail of the following monoclonal antibodies was used: anti-CD3 Alexa700 (BD, San Jose, CA, USA), anti-CD4 Q.605 (Invitrogen, Carlsbad, CA, USA), anti-CD8 PerCP (BioLegend, San Diego, CA, USA), anti-CD45RA BV650 (BioLegend), anti-CD27 APC-eFluor 780 (eBioscience, San Diego, CA, USA), anti-CCR7 PE-CF594 (BD), anti-CD14 V500 (BD), anti-CD19 BV510 (BioLegend), Live/Dead Amcyan (Invitrogen), anti-BTLA (BD), anti-TIM3 BV421 (BioLegend), anti-PD-1 PE-Cy7 (BioLegend), anti-CTLA-4 APC (BD), anti-TIGIT PE (eBioscience), and anti-LAG3 FITC (Novus, Centennial, CO, USA). After washing, cells were resuspended in staining buffer and sorted on an ARIA-SORP. The following antibodies were used for intracellular flow cytometry to assess the phosphorylation of mTOR targets: anti-S6 (S235/S236) V450 (BD) and anti-p4E-BP1 (T36/46) PE (BD). For the rapamycin assay, cells were stimulated for 3 days with CD3/CD28 beads and then treated with 25 nM of rapamycin overnight. Cells were stained and analyzed for PD-1 surface expression. The following antibodies were used to identify memory T cells and effector memory T cells: anti-CD8 BV711 (BD), anti-CD45RA BV650 (BioLegend), anti-CCR7 FITC (R&D Systems, Minneapolis, MN, USA), and anti-CD27 APC eFluor 780 (eBioscience). Data were analyzed using FlowJo software (TreeStar) and our custom pipeline for dimensionality reduction.

Isolated CD8 cells (50,000) from TILs were lysed in 50 μl of radioimmunoprecipitation assay (RIPA) buffer (150 mM of NaCl, 1 mM of EDTA, 0.5% Triton X-100, 0.5% deoxycholic acid, 0.5% sodium dodecyl sulfate (SDS), and 100 mM of Tris (pH 7.5) containing protease and phosphatase inhibitors). An equal sample volume (100 μl) was used for reducing gel electrophoresis (4%–12%), followed by immunodetection with antibodies toward CYP27A1 (Abcam, Cambridge, UK), Actin (Abcam), PHGDH (Cell Signaling, Danvers, MA, USA), PKM2 (Cell Signaling), and PCK2 (Cell Signaling). Three biological replicates were performed.

RNA was purified from CD8 T cells using RNeasy Micro Kits (Qiagen, Hilden, Germany), followed by low-input RNASeq library generation using Takara SMART-Seq v4 Ultra Low/Nextera XT with Nextera Index v2 Set A. Paired-end sequencing reactions were run on an Illumina NextSeq 550 High Output platform (25M total reads per sample). Raw demultiplexed fastq paired-end read files were trimmed of adapters and filtered using the program skewer to remove reads with an average Phred quality score of less than 30, or trimmed to a length of less than 36 (80). Trimmed reads were then aligned using the HISAT2 aligner to the Homo sapiens NCBI reference genome assembly version GRCh38 and sorted using SAMtools (81, 82). Aligned reads were counted and assigned to gene meta-features using the program featureCounts as part of the Subread package (83). These count files were imported into the R programming environment and were assessed for quality control, normalized, and analyzed via the limma-trend method (84) for differential gene expression testing as well as regression modeling and GSVA (85). Statistical power was guided by the detection of a twofold gene expression change in ≥8 samples per group at a minimum power of 80% at p = 0.05, assuming a sequencing depth between 5× and 8× and a coefficient of variation of 0.4 (commonly used for human studies) (86). Linear regression modeling was performed using the limma framework. RNAseq data supporting this study have been deposited in the Gene Expression Omnibus (GEO) public database with the accession number pending.

All cells were resuspended in Dulbecco’s phosphate-buffered saline (DPBS) with 0.04% bovine serum albumin (BSA) and immediately processed for scRNAseq, as follows. Cell count and viability were determined using trypan blue on a Countess FL II, and approximately 12,000 cells were loaded for capture onto the Chromium System using the v2 single-cell reagent kit according to the manufacturer’s protocol (10X Genomics, Pleasanton, CA, USA). Following capture and lysis, cDNA was synthesized and amplified (12 cycles) as per the manufacturer’s protocol (10X Genomics). The amplified cDNA from each channel of the Chromium System was used to construct an Illumina sequencing library and was sequenced on an Illumina HiSeq 2500 with 150-cycle sequencing (asymmetric reads per 10X Genomics). Illumina basecall files (*.bcl) were converted to FASTQs using CellRanger v3.0, which uses bcl2fastq v2.17.1.14. FASTQ files were then aligned to the GRCh38 human reference genome and transcriptome using the CellRanger v3.0 software pipeline with default parameters, as reported previously (87); this demultiplexes the samples, generates a gene-versus-cell expression matrix based on the barcodes, and assigns UMIs that enable determination of the individual cell from which the RNA molecule originated. Overall, an estimated total of 71,238 cells were analyzed.

Bioinformatic analysis of cells based on whole transcriptomes was performed using the R package Seurat (version 3.0) (88). The dimensionality of gene–barcode matrices was first reduced to 18 principal components using principal component analysis (PCA). PCA-reduced data were further reduced to two-dimensional space using the UMAP method and visualized. Graph-based clustering of cells was conducted in the PCA space; a sparse nearest-neighbor graph of the cells was built first, and Louvain modularity optimization was then applied. The number of nearest neighbors was logarithmic in accordance with the number of cells. In the last step, repeated cycles of hierarchical clustering and merging of cluster pairs that had no significant differential expression was performed, until no more cluster pairs could merge. Differential gene expression analyses of each cluster were conducted using the Wilcoxon rank sum test. The log2 fold-change in expression of a certain gene (UMIs) in one cluster vs. all other clusters, and the corresponding adjusted p-values, were calculated for each cluster, with pathway enrichment scores generated using the GSEA method (89).

Unless otherwise indicated, the Student’s t-test was used, with p ≤ 0.05 chosen as the level of significance, and analyses were performed using GraphPad Prism 8.