Integrative transcriptomic and single-cell analyses identify KLRD1 enrichment in exhausted CD8⁺ T cells in cutaneous squamous cell carcinoma

Integrative transcriptomic and single-cell analyses identify KLRD1 enrichment in exhausted CD8⁺ T cells in cutaneous squamous cell carcinoma | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Integrative transcriptomic and single-cell analyses identify KLRD1 enrichment in exhausted CD8⁺ T cells in cutaneous squamous cell carcinoma Dan-dan Zou, Meng-en Li, Yi-chao Jin, Li-tao Suo, Yuan Chen, Zhen Guan, and 3 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-8726237/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Cutaneous squamous cell carcinoma (cSCC) is characterized by complex interactions between tumor cells and the immune microenvironment. Here, we performed integrated transcriptomic and single-cell analyses to dissect immune-related molecular features and identify key regulators of CD8⁺ T cell exhaustion in cSCC. Differential expression analysis of GSE42677 identified 2,149 genes significantly dysregulated in cSCC compared with normal skin, with enrichment in immune- and inflammation-related pathways. Immune cell infiltration profiling revealed distinct patterns between tumor and normal tissues, and weighted gene co-expression network analysis (WGCNA) highlighted modules negatively correlated with activated CD8⁺ T cells and other immune populations. Protein–protein interaction network analysis and hierarchical clustering defined two molecular subtypes of cSCC with divergent immune landscapes. Single-cell RNA sequencing analysis further revealed expansion of exhausted CD8⁺ T (CD8_Exh) cells in tumors. Subpopulation-specific pathway enrichment and hub gene analysis demonstrated distinct functional programs across CD8⁺ T cell subsets, including antigen processing, apoptosis, MAPK and NF-κB signaling, and cytokine interactions. Pseudotime and trajectory analyses using CytoTRACE, Monocle3, and Slingshot identified differentiation trajectories from naive to exhausted and effector CD8⁺ T cells, highlighting dynamic expression of key genes such as CCL5, CD3D, CTSD, KLRD1, and PTPN6 in the CD8_Exh subset. Spearman correlation with pseudotime indicated a strong positive association for KLRD1, suggesting its role in driving CD8⁺ T cell exhaustion. Experimental validation using quantitative PCR and immunofluorescence confirmed elevated expression of KLRD1 in cSCC tissues. Collectively, these results provide a comprehensive landscape of immune regulation in cSCC and identify KLRD1 as a candidate marker of exhausted CD8⁺ T cells, offering a potential target for immunotherapeutic strategies. Cutaneous squamous cell carcinoma (cSCC) CD8⁺ T cell exhaustion KLRD1 Immune microenvironment Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Introduction Cutaneous squamous cell carcinoma (cSCC) is a common skin cancer, and while most lesions are manageable with surgery, a subset progresses to more aggressive disease 1 . Despite its prevalence, the molecular and immune mechanisms underlying cSCC progression are not fully understood. Tumor infiltrating lymphocytes (TILs), particularly CD8⁺ T cells, are generally considered important mediators of antitumor immunity 2 . However, high immune infiltration does not always correlate with effective tumor control 3 . Tumors can induce functional restraint or exhaustion in T cells, leading to heterogeneous immune responses within the tumor microenvironment 4 , 5 . Such functional heterogeneity has been reported in other solid tumors, where immune infiltration coexists with inhibitory signaling and immune regulation 6 . In cSCC, chronic inflammation and dysregulated immune signaling are common features of the tumor microenvironment 7 . Studies have reported infiltration of CD8⁺ T cells, regulatory T cells, and natural killer cells, alongside evidence of immune suppression and variable cytotoxic activity 8 . Despite this, the diversity of CD8⁺ T cell states including effector, memory, and exhausted populations and the transcriptional programs associated with functional restraint remain incompletely characterized 9 . Moreover, whether molecular subtypes with distinct immune microenvironments exist in cSCC, and how they relate to potential responses to immunotherapy, has not been systematically investigated. Although CD8⁺ T cell exhaustion has received increasing attention as an important factor shaping antitumor immunity and immunotherapy response, its cellular heterogeneity and molecular characteristics in cSCC remain poorly defined 9 . Recent studies have further suggested that exhausted CD8⁺ T cells represent a heterogeneous and dynamically regulated population, rather than a uniform dysfunctional state 10 , 11 . However, how such exhaustion-related programs are organized and maintained within the cSCC tumor microenvironment remains unclear. To address these questions, we performed an integrative analysis combining bulk transcriptomics, immune infiltration profiling, and single-cell RNA sequencing. Our study reveals intertumoral immune heterogeneity, identifies molecular subtypes with distinct immune profiles, and characterizes transcriptional programs associated with CD8⁺ T cell exhaustion. Among these, KLRD1 emerges as a candidate marker linked to exhaustion states. These results provide insight into the immune landscape of cSCC and may inform strategies to improve immune-based interventions in this cancer type. Methods Data Acquisition and Preprocessing Public transcriptomic data (GSE42677 and GSE45164) were downloaded from the GEO database using the GEOquery R package (v2.70.0). Raw expression matrices were extracted via the exprs function, and probe IDs were annotated with gene symbols using the platform annotation file. For probes mapping to multiple gene symbols, only the first symbol was retained. Duplicate gene symbols and missing values were removed using dplyr (v1.1.4) and na.omit. The final gene expression matrix, with genes as rows and samples as columns, was log2-transformed (log2(x + 1)) to stabilize variance and exported for downstream analyses. Differential Expression Analysis Differentially expressed genes (DEGs) between cutaneous squamous cell carcinoma (cSCC) and normal samples were identified using the limma (v3.58.1). A design matrix representing the sample groups was constructed, and linear models were fitted. Contrasts were defined for cSCC vs normal, followed by empirical Bayes moderation. DEGs were defined as those with |log2 fold change| > 0.3 and adjusted p-value < 0.05. Volcano plots and heatmaps were generated using ggplot2 (v3.5.2), EnhancedVolcano (v1.20.0), and pheatmap (v1.0.13). The top significantly up- and down-regulated genes were selected for visualization. Hierarchical Clustering of Key Genes Expression values of predefined key genes (including FN1, MMP9, CXCL8, ESR1, PXDN, TIMP1, CXCL10, TLR2, CXCL1, MMP3, SPP1, CD36, GZMB, SOX9, MMP1, KIT, GATA3, NT5E, BMP2, CXCL9, PLAUR, S100A12, COL3A1, MMP13, CXCL13) were converted to numeric format. Sample-wise Euclidean distances were calculated, and hierarchical clustering was performed using the Ward.D2 method via the hclust function. Clusters were visualized using factoextra (v1.0.7) and FactoMineR (v2.11), and dendrograms were annotated with cluster assignments. Immune Cell Infiltration Analysis (ssGSEA) Gene sets representing immune cell types were constructed from CSV files of marker genes, with gene identifiers as row names. Each column was converted into a list element for single-sample Gene Set Enrichment Analysis (ssGSEA). Normalized gene expression matrices (genes × samples) were log2-transformed (log2(TPM + 1) or log2-CPM + 1). ssGSEA was performed using the GSVA (v1.50.5) with the ssgsea method and kcdf = "Gaussian". Absolute gene ranking was applied (abs.ranking = TRUE), generating enrichment scores for each immune cell type per sample. Scores were normalized to a 0–1 range using Min-Max scaling for comparability across samples. Raw and normalized ssGSEA scores were saved for downstream analyses. Immune cell infiltration patterns were visualized using heatmaps, with rows scaled to highlight relative differences and columns grouped according to clinical classification. For detailed group comparisons, half-violin, half-point, and boxplot overlays were generated using ggplot2 and gghalves (v0.1.4). Data were reshaped into long format, and statistical significance between groups (e.g., C1 vs C2) was assessed using two-sided t-tests, with significance indicated on the plots. Custom colors were applied to violin fills and points. Functional Enrichment Analysis KEGG pathway enrichment analysis of differentially expressed genes (DEGs) or selected gene modules was first performed using KOBAS ( http://bioinfo.org/kobas/genelist/ ) to identify significantly enriched pathways. Genes were mapped to KEGG pathway terms, and enrichment significance was assessed using the hypergeometric test with multiple testing correction (FDR < 0.05). The resulting pathway lists, including the number of input genes per pathway and associated p-values, were subsequently visualized using WeightedTreemaps (v0.1.4) in R. Voronoi-based treemaps were constructed, with cell sizes proportional to the number of genes mapped to each pathway, and colors representing distinct pathway categories. Gradient color schemes were applied to indicate the level of enrichment, facilitating intuitive interpretation of pathway significance and gene involvement. Protein–Protein Interaction (PPI) Network Analysis To explore potential interactions among differentially expressed genes (DEGs), protein–protein interaction (PPI) networks were constructed using Cytoscape (v3.10.4). DEGs were uploaded to the STRING database ( https://cn.string-db.org/ ) to obtain predicted and experimentally validated protein interactions, with a confidence score cutoff of 0.4 (medium confidence). The resulting network data were imported into Cytoscape for visualization and further analysis. Key modules or clusters within the network were identified with the MCODE, and functional enrichment of the hub genes or modules was annotated using KEGG and Gene Ontology databases. Graphical layouts and node colors were adjusted in Cytoscape to reflect interaction strength and gene significance. Single Cell RNA-seq Data Processing and Quality Control Raw count data from GSE144236 were imported into R (v5.1.0) and converted to a sparse matrix for memory efficiency. A Seurat (v5.2.1) object was created, and patient metadata were integrated. Cells with low gene counts or high mitochondrial gene percentage were filtered out to ensure high quality single cell data. Normalization was performed using the LogNormalize method with a scale factor of 10,000. Highly variable genes (n = 2,000) were identified using the vst method, and all genes were scaled. Principal component analysis (PCA) was performed on the top variable genes, and the first 10 PCs were used to construct the K-nearest neighbor (KNN) graph (FindNeighbors) and cluster cells (FindClusters, resolution = 0.5). Dimensionality reduction was visualized using UMAP and t-SNE. T Cell Extraction and Subpopulation Identification T cells were extracted from the full dataset based on level1_celltype == "Tcell" and filtered for cells with ≥ 500 detected genes. PCA and clustering were repeated on this subset, and UMAP was used for visualization. Subpopulations were annotated, and the proportion of each T cell subtype per group (Tumor vs Normal) was calculated. Differences in proportions were assessed using t-tests implemented via ggsignif (v0.6.4). CD8⁺ T cells (CD8_EM, CD8_EMRA, CD8_Exh, CD8_Naive) were further subsetted and reanalyzed. PCA (npcs = 15), KNN graph construction, and clustering (resolution = 0.6) were repeated. UMAP and t-SNE visualizations confirmed the expression of canonical markers CD8A and CD8B using FeaturePlot and VlnPlot. Proportions of CD8⁺ T cell subtypes in Tumor and Normal samples were visualized with stacked bar plots. Differential Gene Expression and KEGG Pathway Enrichment For each CD8⁺ T cell subpopulation, differential expression analysis was conducted using Seurat’s FindMarkers function with the Wilcoxon rank-sum test. Genes with adjusted p-value 0.3 were considered significant. Significantly up- and down-regulated genes were subjected to KEGG pathway enrichment using clusterProfiler (v4.10.1). Enriched pathways were visualized as bubble plots in ggplot2, where bubble size represents the number of input genes and color indicates enrichment significance (-log10 p-value). Pathways were manually ordered according to biological relevance. Single-Cell Trajectory and Pseudotime Analysis CD8⁺ T cell differentiation was analyzed using Monocle3 (v1.3.7), CytoTRACE (v0.3.3), and Slingshot (v2.10.0). In Monocle3, Seurat RNA counts were converted to a cell_data_set object, cells were clustered, and pseudotime trajectories were inferred; dynamically expressed genes along pseudotime were identified with Moran’s I and visualized. CytoTRACE estimated differentiation potential, visualized on UMAP with gene-specific patterns. Slingshot reconstructed pseudotime trajectories on a SingleCellExperiment object using Seurat cluster labels, and trajectories were visualized on UMAP coordinates. Gene Expression Dynamics along Pseudotime To examine dynamic expression patterns, gene expression values were extracted from the Seurat RNA assay and matched with Monocle3-calculated pseudotime. Spearman correlation coefficients were computed to quantify the association between gene expression and pseudotime. Expression trends were visualized using scatter plots with pseudotime on the x-axis and gene expression on the y-axis, with LOESS curves fitted to highlight overall trajectories. Statistical significance of correlations was reported, and results were saved for downstream integration. Quantitative PCR (qPCR) Analysis Total RNA was extracted from cells or tissues using TRIzol reagent (Invitrogen) according to the manufacturer’s protocol. RNA quality and concentration were assessed with a NanoDrop spectrophotometer (Thermo Fisher). Complementary DNA (cDNA) was synthesized using the PrimeScript RT reagent kit (Takara). qPCR was performed on a QuantStudio 6 Flex system (Applied Biosystems) using SYBR Green Master Mix (Takara). Relative gene expression was calculated using the 2^−ΔΔCt method with GAPDH as the internal control. Primers used for each target gene are listed in Supplementary Table S1 . Experiments were performed in triplicate. Immunofluorescence (IF) Staining Cells or tissue sections were fixed with 4% paraformaldehyde for 15 min and permeabilized with 0.1% Triton X-100 for 10 min. Samples were blocked with 5% bovine serum albumin (BSA) for 1 h at room temperature, followed by incubation with primary antibodies overnight at 4°C. After washing, samples were incubated with appropriate fluorophore-conjugated secondary antibodies for 1 h at room temperature in the dark. Nuclei were counterstained with DAPI. Images were captured using a confocal laser scanning microscope and analyzed with ZEN (v3.3). Results Differential gene expression and immune cell infiltration patterns in cSCC Differential expression analysis was performed on the GSE42677 dataset to compare cutaneous squamous cell carcinoma (cSCC) tissues with normal skin samples. Using the criteria of |log2FC| ≥ 0.3 and p < 0.05, a total of 2,149 differentially expressed genes (DEGs) were identified, including 1,161 upregulated and 988 downregulated genes in cSCC (Fig. 1 A). The expression patterns of the top DEGs were visualized using a heatmap. Genes such as CDHR2, PKD2L1, SLC30A3, MORN1, SPOCK2, VIPR2, LIMK1, LHX5, and KIAA1614 were markedly upregulated in cSCC tissues, whereas CP, PTH1R, FLT4, DSPP, PCDHA9, PIP5K1B, RELN, S1PR1, and CMAHP showed relatively lower expression compared with normal skin samples (Fig. 1 B). To characterize the immune microenvironment of cSCC, immune cell infiltration scores were evaluated and visualized as a heatmap. Increased infiltration scores were observed in cSCC samples for CD56dim natural killer cells, plasmacytoid dendritic cells, activated CD8⁺ T cells, macrophages, type 1 T helper cells, regulatory T cells, myeloid-derived suppressor cells, natural killer cells, natural killer T cells, γδ T cells, CD56bright natural killer cells, type 17 T helper cells, neutrophils, T follicular helper cells, and activated dendritic cells. In contrast, higher infiltration levels of activated B cells, mast cells, central memory CD4⁺ T cells, memory B cells, activated CD4⁺ T cells, type 2 T helper cells, effector memory CD8⁺ T cells, eosinophils, central memory CD8⁺ T cells, monocytes, immature B cells, effector memory CD4⁺ T cells, and immature dendritic cells were predominantly observed in normal skin tissues (Fig. 1 C). Weighted gene co-expression network analysis (WGCNA) was further conducted to explore gene modules associated with immune infiltration patterns. The MEblue module exhibited significant negative correlations with the majority of immune cell infiltration scores, including activated CD8⁺ T cells (r = − 0.76, p = 1×10⁻⁴), activated CD4⁺ T cells (r = − 0.70, p = 6×10⁻⁴), and natural killer T cells (r = − 0.86, p = 2×10⁻⁶) (Fig. 1 D). Intersection analysis between DEGs and genes within the MEblue module identified 331 overlapping genes (Fig. 1 E). Functional enrichment analysis of these intersecting genes revealed significant enrichment in immune- and inflammation-related pathways, including cytokine–cytokine receptor interaction, complement and coagulation cascades, JAK–STAT signaling pathway, leukocyte transendothelial migration, NF-κB signaling pathway, TGF-β signaling pathway, Th17 cell differentiation, NOD-like receptor signaling pathway, Th1 and Th2 cell differentiation, TNF signaling pathway, Toll-like receptor signaling pathway, chemokine signaling pathway, viral protein interaction with cytokine and cytokine receptor, and IL-17 signaling pathway (Fig. 1 F). Identification of immune-related hub genes and T cell–enriched molecular subtypes in cSCC Protein–protein interaction (PPI) network analysis was performed based on the intersecting genes, and hub genes were identified using a degree threshold greater than 15. A total of 25 hub genes were screened, including PXDN, PLAUR, CXCL13, CXCL9, MMP1, MMP13, NT5E, SOX9, GZMB, MMP3, CXCL1, GATA3, KIT, TIMP1, COL3A1, S100A12, SPP1, CXCL8, CXCL10, BMP2, and FN1. Among these, genes such as FN1, MMP9, CXCL8, MMP3,PXDN, CXCL10, TLR2, TIMP1, ESR1, and CXCL1 occupied central positions within the PPI network, suggesting their potential key roles in regulating the tumor immune microenvironment of cSCC (Fig. 2 A). Based on the expression patterns of these hub genes, hierarchical clustering analysis was applied to cSCC samples from the GSE45164 dataset, resulting in two distinct molecular subtypes, termed C1 and C2 (Fig. 2 B). Differential expression analysis between the two subtypes identified 113 differentially expressed genes, including 58 upregulated and 55 downregulated genes (Fig. 2 C). Genes such as KRT34, FMO2, ZBED2, DUOX1, CIDEA, DLX2, IL18, MAP2, FETUB, and IL37 were upregulated in the C1 subtype, whereas IL1B, CXCL1, G0S2, SRGN, ADAMDEC1, CXCL6, SELL, EVI2B, FCGR3B, and TDO2 were downregulated (Fig. 2 D). To further characterize the immune features of the two molecular subtypes, immune infiltration levels were compared between C1 and C2. The immune infiltration heatmap revealed that type 2 T helper cells, CD56bright natural killer cells, plasmacytoid dendritic cells, and mast cells exhibited higher infiltration scores in the C1 subtype. In contrast, enhanced infiltration of activated B cells, effector memory CD8⁺ T cells, activated CD8⁺ T cells, and activated CD4⁺ T cells was observed in the C2 subtype (Fig. 2 E). Functional enrichment analysis of the subtype-specific differentially expressed genes highlighted significant enrichment in immune- and inflammation-related pathways, including cytokine–cytokine receptor interaction, TNF signaling pathway, IL-17 signaling pathway, NF-κB signaling pathway, chemokine signaling pathway, NOD-like receptor signaling pathway, Th17 cell differentiation, natural killer cell–mediated cytotoxicity, JAK–STAT signaling pathway, and Toll-like receptor signaling pathway (Fig. 2 F). Furthermore, a PPI network was constructed based on the subtype-specific differentially expressed genes, in which immune-associated genes such as IL18, MMP3, CXCL2, IL1B, FCGR3A, SELL, FCGR3B, CXCL1, CD69, IL15, CCL18, CXCL13, LTF, IL6, PTGS2, SPP1, and CCL20 were identified as central nodes (Fig. 2 G), suggesting the presence of key immune regulatory modules distinguishing the two cSCC subtypes. Differential Immune Infiltration and Expansion of Exhausted CD8⁺ T Cells in cSCC To establish a comprehensive cellular landscape, we first performed UMAP visualization of the entire dataset, revealing 14 distinct cell types and clear separation between normal and tumor-derived cells (Fig. 3 A–B). Focusing on T cells, nine T cell subtypes were identified and their distribution across normal and tumor tissues was characterized (Fig. 3 C–D). The proportion of CD8⁺ T cells was markedly higher in tumor tissues compared to normal tissues. Quantitative comparison showed that in tumors, the proportions of CD8_Exh, CD8_Trge, CD8_EMRA, and CD8_EM were notably higher than in normal tissue, whereas CD8_Naïve and CD4_RGCC/CD4_Naïve were more abundant in normal tissue (Fig. 3 E). For a refined analysis, CD8⁺ T cells were extracted and clustered into five clusters, comprising CD8_Exh, CD8_Naïve, CD8_EM, and CD8_EMRA subpopulations (Fig. 3 F–G). Their spatial distribution revealed that CD8_Naïve cells predominantly resided in normal tissue, while the other three clusters were mainly found in tumor tissue (Fig. 3 H). Quantitative proportion analysis confirmed the UMAP-based observations, with CD8_Naïve cells enriched in normal samples and CD8_Exh, CD8_EMRA, and CD8_EM enriched in tumors (Fig. 3 I–J). Functional characteristics and key gene expression of CD8⁺ T cell subsets in cSCC Differential expression analysis was first performed among the four CD8⁺ T cell subsets (Naive, EM, EMRA, and Exh) with a threshold of |log2FC| ≥ 0.3 and adjusted p-value < 0.05, yielding subset specific differentially expressed genes (DEGs) (Fig. 4 A). Each subset’s DEGs were subjected to KEGG pathway enrichment analysis, and only significantly enriched pathways (p < 0.05) are shown to highlight the main functional features of each subset. STRING database was used to construct protein–protein interaction (PPI) networks for the DEGs, with hub genes identified based on degree values. UMAP based expression visualization was then performed to map key pathway related genes onto single-cell data, illustrating their spatial distribution across subsets and thereby reflecting functional distinctions. In the Exh subset, enriched pathways included antigen processing and presentation, apoptosis, carbon metabolism, cellular senescence, chemokine signaling, C-type lectin receptor signaling, cytokine–cytokine receptor interaction, glycolysis/gluconeogenesis, Jak-STAT signaling, NF-kappa B signaling, oxytocin signaling, p53 signaling, T cell receptor signaling, TNF signaling, and Toll-like receptor signaling (Fig. 4 B). Hub genes in Exh included GAPDH, UBA52, RPS27A, EEF2, and RPS2 (Fig. 4 C). Notably, PTPN6, CD3D, TNFRSF18, CTSD, CCL5, KLRD1, CSF1, and TNFRSF1B were highly expressed in the CD8_Exh subset (Fig. 4 D). For the CD8_EMRA subset, enriched pathways included natural killer cell mediated cytotoxicity, cytokine–cytokine receptor interaction, chemokine signaling, antigen processing and presentation, T cell receptor signaling, MAPK signaling, NF-kappa B signaling, TNF signaling, cell cycle, and apoptosis (Fig. 4 E). Hub genes comprised RPL3, RPL11, RPL5, RPS16, and RPS14 (Fig. 4 F), with KLRD1, DUSP2, FCGR3A, and FGR showing higher expression within this subset (Fig. 4 G). In the CD8_EM subset, significant pathways included MAPK signaling, chemokine signaling, apoptosis, cytokine–cytokine receptor interaction, lysosome, natural killer cell-mediated cytotoxicity, and antigen processing and presentation (Fig. 4 H). Hub genes comprised GAPDH, CXCR3, GZMB, SELL, GZMK, and RPL13A (Fig. 4 I), with CXCR4, CXCR3, DUSP2, and HSPA1B particularly highly expressed (Fig. 4 J). For the CD8_Naive subset, enriched pathways included cell adhesion molecules (CAMs), chemokine signaling, hematopoietic cell lineage, intestinal immune network for IgA production, oxytocin signaling, PD-L1 expression and PD-1 checkpoint pathway in cancer, primary immunodeficiency, riboflavin metabolism, and T cell receptor signaling (Fig. 4 K). Hub genes included GAPDH, IFNG, CD4, CD44, CD8A, and GZMB (Fig. 4 L), with CD44, CD4, CCR6, and IFNGR2 highly expressed in this subset (Fig. 4 M). Trajectory Analysis Reveals Potential Key Genes CCL5, CD3D, CTSD, KLRD1, and PTPN6 in CD8⁺_Exh Cells in cSCC CytoTRACE analysis was used to predict the differentiation potential of CD8⁺ T cells, showing that naive cells exhibited the highest scores, followed by Exh, EM, and EMRA cells in descending order (Fig. 5 A–C). UMAP visualization of differentiation potential revealed a clear gradient across subsets, providing a basis for subsequent pseudotime analysis. Monocle3 pseudotime analysis was performed to construct single-cell trajectories, validating the differentiation order predicted by CytoTRACE. Two main trajectories were identified: one progressing from naive → Exh, and the other from naive → EM → EMRA (Fig. 5 D). This analysis allowed quantification of dynamic cellular states along the differentiation path. Expression of key Exh pathway genes was examined across subsets and along pseudotime. Genes including CCL5, CD3D, CTSD, KLRD1, PTPN6, TNFRSF18, and TNFRSF1B were most highly expressed in CD8_Exh cells (Fig. 5 E). Visualization of these genes along the Monocle3 pseudotime trajectories reflected their dynamic changes, highlighting their regulatory roles during the exhaustion and effector differentiation process (Fig. 5 F). Slingshot trajectory analysis was conducted as an independent validation of the Monocle3 results. UMAP visualization showed the spatial distribution and expression gradients of key genes along the trajectories, emphasizing the involvement of Exh-related genes at the terminal stage of differentiation. Two main trajectories consistent with Monocle3 were observed (Fig. 5 G). UMAP mapping of Exh key genes along the Slingshot trajectories further confirmed that these genes increased in expression at the end of the differentiation trajectory corresponding to CD8_Exh cells. In particular, CCL5, CD3D, CTSD, KLRD1, and PTPN6 exhibited elevated expression in the CD8⁺_Exh subset (Fig. 5 H). Integrated validation identifies KLRD1 as a leading candidate gene in exhausted CD8⁺ T cells To validate whether the key genes identified by single-cell analysis exhibited significant expression differences in CD8⁺ T cells and to assess the consistency of their expression patterns with bulk transcriptomic data, we performed integrative analyses across multiple levels. At the single-cell level, the expression of seven candidate genes (CCL5, CD3D, CTSD, KLRD1, PTPN6, TNFRSF18, and TNFRSF1B) was quantified in CD8⁺ T cells from tumor and normal tissues. All seven genes showed significantly higher expression in tumor-derived CD8⁺ T cells compared with those from normal skin (Fig. 6 A). At the bulk transcriptomic level, differential expression analysis between cSCC and normal skin tissues demonstrated that CCL5, CD3D, KLRD1, and TNFRSF1B exhibited expression patterns consistent with the single-cell results, with significantly elevated expression in cSCC samples (Fig. 6 B), further supporting the robustness of these candidate genes. To further evaluate the association of Exh-related genes with CD8⁺ T cell differentiation and exhaustion, Monocle3 pseudotime analysis was performed, focusing on CTSD, KLRD1, and PTPN6. Spearman correlation analysis between gene expression and pseudotime revealed that KLRD1 showed a strong positive correlation with pseudotime (rho = 0.504, p = 1.75 × 10⁻⁴¹), indicating progressively increased expression during CD8⁺ T cell exhaustion and suggesting its suitability as a primary candidate for subsequent functional studies. CTSD displayed a weak but significant positive correlation (rho = 0.091, p = 0.0228), supporting its potential role as a secondary candidate, whereas PTPN6 showed no significant correlation with pseudotime (rho = 0.061, p = 0.126) (Fig. 6 C). Collectively, these analyses provide a dynamic view of key Exh-related genes during CD8⁺ T cell differentiation and exhaustion, offering a rational basis for prioritizing candidate genes for downstream functional validation and for exploring their roles in the immune microenvironment of cutaneous squamous cell carcinoma. KLRD1 is enriched in exhausted CD8⁺ T cells in cSCC First, quantitative PCR was performed on paired cutaneous squamous cell carcinoma (cSCC) and adjacent normal skin tissues to validate the expression changes of KLRD1, PTPN6, and CTSD, thereby assessing the consistency between single cell transcriptomic findings and bulk tissue samples. All three genes were significantly upregulated in cSCC compared with normal tissues (Fig. 7 A). Next, triple immunofluorescence (IF) staining was performed to simultaneously examine the expression and spatial distribution of CD8 (cyan), PD-1 (red), and KLRD1 (yellow) at the tissue level. This analysis aimed to determine whether KLRD1 was localized within CD8⁺ T cells in the tumor microenvironment and to evaluate its spatial association with the exhaustion marker PD-1. Representative images are shown in the order of merged view, CD8, PD-1, and KLRD1 channels, with normal skin tissues displayed in the upper panels and cSCC tissues in the lower panels. In normal skin tissues, CD8⁺ T cells were sparsely distributed, with weak or barely detectable PD-1 and KLRD1 signals. In contrast, cSCC tissues exhibited a marked accumulation of CD8⁺ T cells (cyan), accompanied by pronounced PD-1 (red) and KLRD1 (yellow) staining. Importantly, merged images revealed frequent co-localization of KLRD1 signals with CD8⁺ cells, and substantial overlap between KLRD1 and PD-1 within the same cellular regions. This spatial pattern suggests that KLRD1 is preferentially expressed in PD-1⁺ CD8⁺ T cells within the tumor microenvironment. Overall, these IF results demonstrate that CD8, PD-1, and KLRD1 are all enriched in cSCC tissues compared with normal skin and display a coordinated spatial distribution consistent with an exhausted CD8⁺ T cell phenotype (Fig. 7 B). Quantitative analysis of IF images was further performed to assess the densities of KLRD1⁺, CD8⁺KLRD1⁺, CD8⁺PD-1⁺, and CD8⁺KLRD1⁺PD-1⁺ cells in tumor and normal tissues. Although statistical significance was not reached due to limited sample size and the rarity of certain cell subsets, all four metrics consistently exhibited higher levels in cSCC tissues. This trend was concordant with the single-cell transcriptomic findings indicating elevated KLRD1 expression in exhausted CD8⁺ T cells, thereby providing supportive evidence for the involvement of KLRD1 in the cSCC immune microenvironment (Fig. 7 C). Discussion Bulk transcriptomic analyses reveal immune activation and molecular heterogeneity in cSCC Our integrative analyses of bulk transcriptomic data and immune infiltration profiles indicate that cSCC is far from an “immune cold” tumor; instead, it exhibits pronounced immune and inflammatory activation alongside substantial molecular heterogeneity 12 – 14 . The extensive dysregulation of inflammatory and immune related pathways including TNF, NF-κB, IL-17, and JAK–STAT signaling underscores the interplay between tumorigenesis and local immune remodeling 14 , 15 . These findings are consistent with reports, where chronic inflammation not only drives tumor progression but also shapes the immune contexture and therapeutic responsiveness 16 – 18 . Immune infiltration analyses revealed notable intertumoral variability: tumor tissues were enriched in activated CD8⁺ T cells, regulatory T cells, and NK cells, whereas memory B and effector memory CD8⁺ T cells in normal skin retained cytotoxic and immune surveillance related transcriptional programs 12 , 19 – 21 . This observation highlights a critical distinction: immune infiltration in tumors does not automatically translate into effective antitumor activity, as tumor-infiltrating lymphocytes display transcriptional features of functional restraint or exhaustion, in contrast to the cytotoxic and immune surveillance programs retained by effector memory CD8⁺ T cells in normal skin. In this context, infiltrating T cells may be functionally restrained by intrinsic transcriptional programs, such as the MEblue module identified in our weighted gene coexpression network analysis. Notably, the MEblue module exhibited significant negative correlations with the majority of immune cell infiltration scores, including activated CD8⁺ T cells, activated CD4⁺ T cells, and natural killer T cells, suggesting that this module may represent a broad immunoregulatory or compensatory transcriptional program. Such inverse associations imply the presence of regulatory circuits that constrain immune activation despite substantial immune infiltration. Similar patterns have been reported in other tumor types, where high immune cell abundance coexists with inhibitory signaling and immune exclusion, including non–small cell lung cancer 22 . Stratification of cSCC samples based on immune related hub genes revealed two molecular subtypes with distinct immune landscapes. Subtype C1 was enriched for immunoregulatory populations, including Th2 cells, plasmacytoid dendritic cells, and mast cells, suggesting a more tolerogenic or modulatory microenvironment. In contrast, C2 exhibited heightened cytotoxic and inflammatory signatures. This mirrors the “immune contexture” concept in other solid tumors, where composition, activation state, and spatial distribution of tumor-infiltrating lymphocytes predict immunotherapy outcomes and disease progression 23 . Notably, the existence of these subtypes in cSCC may partially explain variable clinical responses to immune checkpoint blockade observed in clinical cohorts. Dynamic expansion and exhaustion programming of CD8⁺ T cells in cSCC Single cell profiling further revealed that CD8⁺ T cell expansion in cSCC is accompanied by pronounced functional heterogeneity. Tumor infiltrating CD8⁺ T cells largely comprise exhausted (Exh), effector memory (EM), and terminal effector (EMRA) populations, whereas naïve cells remain largely confined to normal skin. The preferential accumulation of exhausted and terminal effector cells suggests sustained antigenic stimulation and immunosuppressive cues from the tumor microenvironment 24 , 25 . Functional analyses showed that CD8_Exh cells activate TNF, NF-κB, TCR, and p53 pathways while undergoing metabolic and inflammatory reprogramming, whereas CD8_EM and EMRA populations are enriched for cytotoxic, NK-like, MAPK, and cell cycle–associated pathways. CD8_Naïve cells predominantly express adhesion molecules and baseline TCR/PD-1 signaling. Collectively, these data suggest that exhaustion is an actively regulated, rather than passive, state a perspective increasingly recognized in tumor immunology, where partial effector function persists to balance immune surveillance with tissue protection 26 – 28 . Trajectory analyses using CytoTRACE and Monocle3 revealed that exhaustion develops progressively along two primary differentiation axes: naïve → Exh and naïve → EM → EMRA, with exhausted cells occupying terminal branches. Key genes, including CCL5, CD3D, CTSD, KLRD1, and PTPN6, exhibit gradual upregulation along pseudotime, highlighting a coordinated transcriptional program driving exhaustion. These findings support the emerging view that exhaustion markers should be interpreted in a temporal and regulatory context rather than as static indicators of dysfunction, reflecting dynamic metabolic and transcriptional remodeling within exhausted CD8⁺ T cells 29 . KLRD1 as a marker and potential modulator of exhaustion in cSCC Among exhaustion-associated genes, KLRD1 stands out as a robust candidate. Its expression increases progressively along the exhaustion trajectory, spatially colocalizes with CD8⁺ and PD-1⁺ cells, and is detectable at both bulk and tissue levels. Prior studies in HNSC and other squamous cell carcinomas have linked KLRD1to terminally exhausted CD8⁺ and NK cells, where it participates in inhibitory signaling to prevent overactivation 30 . In this context, KLRD1 may function not only as a biomarker of exhaustion but also as a regulatory element within a transcriptional module that integrates chronic antigenic stimulation, effector differentiation, and metabolic adaptationa mechanism increasingly recognized across tumor immunology studies. The potential translational significance of KLRD1 is multifold. Beyond serving as a biomarker of exhaustion, its dynamic expression profile and tissue-level detectability make it a candidate for stratifying patients or guiding immunotherapeutic interventions. For example, combinatorial therapies targeting PD-1 alongside modulation of KLRD1 mediated pathways could enhance cytotoxic T cell function while limiting deleterious overactivation, although this remains speculative and requires functional validation. Limitations and future perspectives Despite the multi-layered evidence supporting KLRD1 as an exhaustion-associated marker, several limitations exist. Immunofluorescence relied on a single marker, limiting resolution of functional states at the single-cell level, and sample size constraints prevented statistically robust tissue quantification. Furthermore, the mechanistic role of KLRD1 in regulating exhaustion remains untested in vitro or in vivo. Future work should assess its functional contribution to T cell exhaustion, interactions with canonical checkpoints (e.g., PD-1, LAG-3), and predictive value for immunotherapy response. Expanding this integrative framework to larger cohorts and other tumor types will further refine strategies to exploit immune heterogeneity for clinical benefit. Conclusion Together, our integrative analyses reveal a highly dynamic and heterogeneous immune landscape in cSCC, identify KLRD1 as a key exhaustion-associated gene, and provide mechanistic and translational insights into CD8⁺ T cell dysfunction. These findings offer a foundation for future studies aimed at modulating exhaustion pathways to improve immunotherapeutic efficacy in cSCC and potentially other solid tumors. Declarations Author Contributions: Dan-dan Zou conceived and designed the study. Meng-en Li contributed to methodology development. Yi-chao Jin performed software implementation and provided resources. Li-tao Suo carried out experimental validation. Li Zhang conducted formal data analysis. Yuan Chen was responsible for investigation and data curation. Dan-dan Zou and Meng-en Li wrote the original draft. Xiao-chuan Wang and Zhi-qiong Wang contributed to writing—review and editing. Zhen Guan prepared the visualizations. Zhi-qiong Wang supervised the study. Xiao-chuan Wang managed the project. Dan-dan Zou acquired funding for the study. Acknowledgments: This work was supported by Kunming University of Science and Technology & the First People's Hospital of Yunnan Province Joint Special Project on Medical Research (KUST-KH2023023Y), Yunnan Province Science and Technology Department & Kunming Medical University Joint Special Project for Applied Basic Research (202501AY070001-261), Yunnan Fundamental Research Projects (grant NO. 202401AU070030), High-level Talents Scientific Research Project of Yunnan Provincial Health Commission (2023-KHRCBZ-B13) and The First People's Hospital of Yunnan Province (KHYJ-2023-6-08). Data Availability Statement: All the data analyzed in this study were obtained from publicly available repositories (GEO database). Conflicts of Interest: The authors have no relevant financial or non-financial interests to dis-close. Ethics and Consent to Participate declarations: Not applicable. References Winge MCG, et al. Advances in cutaneous squamous cell carcinoma. Nat Rev Cancer. 2023;23:430–49. https://doi.org/10.1038/s41568-023-00583-5 . Li F, et al. The association between CD8 + tumor-infiltrating lymphocytes and the clinical outcome of cancer immunotherapy: A systematic review and meta-analysis. EClinicalMedicine. 2021;41:101134. https://doi.org/10.1016/j.eclinm.2021.101134 . Pagès F, et al. Immune infiltration in human tumors: a prognostic factor that should not be ignored. Oncogene. 2010;29:1093–102. https://doi.org/10.1038/onc.2009.416 . Heintzman DR, Fisher EL, Rathmell JC. Microenvironmental influences on T cell immunity in cancer and inflammation. Cell Mol Immunol. 2022;19:316–26. https://doi.org/10.1038/s41423-021-00833-2 . Jiang W, et al. Exhausted CD8 + T Cells in the Tumor Immune Microenvironment: New Pathways to Therapy. Front Immunol. 2020;11:622509. https://doi.org/10.3389/fimmu.2020.622509 . Ouyang P, et al. Infiltration characteristics and regulatory mechanisms of CD8(+) T lymphocytes in solid tumors: spatial distribution, biological functions, and interactions with the immune microenvironment. Front Immunol. 2025;16:1661545. https://doi.org/10.3389/fimmu.2025.1661545 . Buruiană A, et al. The Tumor Stroma of Squamous Cell Carcinoma: A Complex Environment That Fuels Cancer Progression. Cancers (Basel). 2024;16. https://doi.org/10.3390/cancers16091727 . Farhood B, Najafi M, Mortezaee K. CD8(+) cytotoxic T lymphocytes in cancer immunotherapy: A review. J Cell Physiol. 2019;234:8509–21. https://doi.org/10.1002/jcp.27782 . Baessler A, Vignali D. A. T cell exhaustion. Annu Rev Immunol. 2024;42:179–206. Philip M, Schietinger A. CD8 + T cell differentiation and dysfunction in cancer. Nat Rev Immunol. 2022;22:209–23. Li C, Yuan Y, Wang Q, Jiang X. Epigenetic Regulation of CD8⁺ T Cell Exhaustion: Recent Advances and Update. Front Immunol. 2025;16:1700039. Rajesh A, et al. Characterisation of the immune microenvironment of cutaneous squamous cell carcinoma in immunosuppression. Exp Dermatol. 2022;31:1720–8. https://doi.org/10.1111/exd.14650 . Luci C, et al. Cutaneous Squamous Cell Carcinoma Development Is Associated with a Temporal Infiltration of ILC1 and NK Cells with Immune Dysfunctions. J Invest Dermatol. 2021;141:2369–79. https://doi.org/10.1016/j.jid.2021.03.018 . Dong Q, Zhang Z, Li S, Liang L. Mechanisms of immunotherapy in cutaneous squamous cell carcinoma in the tumor microenvironment. Front Immunol. 2025;16:1660272. https://doi.org/10.3389/fimmu.2025.1660272 . Saeidi V, Doudican N, Carucci JA. Understanding the squamous cell carcinoma immune microenvironment. Front Immunol. 2023;14:1084873. https://doi.org/10.3389/fimmu.2023.1084873 . Greten FR, Grivennikov SI. Inflammation and Cancer: Triggers, Mechanisms, and Consequences. Immunity. 2019;51:27–41. https://doi.org/10.1016/j.immuni.2019.06.025 . Nishida A, Andoh A. The Role of Inflammation in Cancer: Mechanisms of Tumor Initiation, Progression, and Metastasis. Cells. 2025;14. https://doi.org/10.3390/cells14070488 . Wen Y, et al. Chronic inflammation, cancer development and immunotherapy. Front Pharmacol. 2022;13:1040163. https://doi.org/10.3389/fphar.2022.1040163 . Strobl J, Haniffa M. Functional heterogeneity of human skin-resident memory T cells in health and disease. Immunol Rev. 2023;316:104–19. https://doi.org/10.1111/imr.13213 . Tokura Y, Phadungsaksawasdi P, Kurihara K, Fujiyama T, Honda T. Pathophysiology of Skin Resident Memory T Cells. Front Immunol. 2020;11:618897. https://doi.org/10.3389/fimmu.2020.618897 . Chen X, Zheng Y, Man X, Li W. Tissue-resident memory T cells and their function in skin diseases. Chin Med J (Engl). 2025;138:1175–83. https://doi.org/10.1097/cm9.0000000000003499 . Anichini A, Perotti VE, Sgambelluri F, Mortarini R. Immune Escape Mechanisms in Non Small Cell Lung Cancer. Cancers (Basel). 2020;12. https://doi.org/10.3390/cancers12123605 . Wu T, Wu X, Wang HY, Chen L. Immune contexture defined by single cell technology for prognosis prediction and immunotherapy guidance in cancer. Cancer Commun (Lond). 2019;39:21. https://doi.org/10.1186/s40880-019-0365-9 . Xie H, Xi X, Lei T, Liu H, Xia Z. CD8(+) T cell exhaustion in the tumor microenvironment of breast cancer. Front Immunol. 2024;15:1507283. https://doi.org/10.3389/fimmu.2024.1507283 . Zhang B, et al. CD8(+) T cell exhaustion and its regulatory mechanisms in the tumor microenvironment: key to the success of immunotherapy. Front Immunol. 2024;15:1476904. https://doi.org/10.3389/fimmu.2024.1476904 . Chow A, Perica K, Klebanoff CA, Wolchok JD. Clinical implications of T cell exhaustion for cancer immunotherapy. Nat Rev Clin Oncol. 2022;19:775–90. https://doi.org/10.1038/s41571-022-00689-z . Baessler A, Vignali DA. A. T Cell Exhaustion. Annu Rev Immunol. 2024;42:179–206. https://doi.org/10.1146/annurev-immunol-090222-110914 . Jenkins E, Whitehead T, Fellermeyer M, Davis SJ, Sharma S. The current state and future of T-cell exhaustion research. Oxf Open Immunol. 2023;4:iqad006. https://doi.org/10.1093/oxfimm/iqad006 . Zhong T, Sun S, Zhao M, Zhang B, Xiong H. The mechanisms and clinical significance of CD8(+) T cell exhaustion in anti-tumor immunity. Cancer Biol Med. 2025;22:460–80. https://doi.org/10.20892/j.issn.2095-3941.2024.0628 . Dong C, Lin Z, Hu Y, Lu Q. KLRD1 (CD94): A Prognostic Biomarker and Therapeutic Candidate in Head and Neck Squamous Cell Carcinoma. J Cancer. 2025;16:982–95. https://doi.org/10.7150/jca.104762 . Additional Declarations No competing interests reported. Supplementary Files TableS1.xlsx Table S1. Primer sequences used for quantitative PCR (qPCR). The table lists forward and reverse primer sequences for genes analyzed in this study. Primers were designed to specifically amplify target transcripts, and qPCR was performed to assess gene expression levels in cSCC and paired normal skin samples. Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-8726237","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":584055790,"identity":"6c3a1784-b943-426c-895c-52d7daace043","order_by":0,"name":"Dan-dan Zou","email":"","orcid":"","institution":"First People's Hospital of Yunnan Province","correspondingAuthor":false,"prefix":"","firstName":"Dan-dan","middleName":"","lastName":"Zou","suffix":""},{"id":584055791,"identity":"50646142-7be1-487c-820b-6e8ffdd613d7","order_by":1,"name":"Meng-en Li","email":"","orcid":"","institution":"Yunnan Cancer Hospital","correspondingAuthor":false,"prefix":"","firstName":"Meng-en","middleName":"","lastName":"Li","suffix":""},{"id":584055792,"identity":"b8283a08-e9d0-40fd-a825-9aadd136f686","order_by":2,"name":"Yi-chao Jin","email":"","orcid":"","institution":"First People's Hospital of Yunnan Province","correspondingAuthor":false,"prefix":"","firstName":"Yi-chao","middleName":"","lastName":"Jin","suffix":""},{"id":584055793,"identity":"afb82a4c-51b7-405d-8395-5c9027a515ad","order_by":3,"name":"Li-tao Suo","email":"","orcid":"","institution":"First People's Hospital of Yunnan Province","correspondingAuthor":false,"prefix":"","firstName":"Li-tao","middleName":"","lastName":"Suo","suffix":""},{"id":584055794,"identity":"b4f95aa0-45ca-42ce-91ef-2003caea19b5","order_by":4,"name":"Yuan Chen","email":"","orcid":"","institution":"First People's Hospital of Yunnan Province","correspondingAuthor":false,"prefix":"","firstName":"Yuan","middleName":"","lastName":"Chen","suffix":""},{"id":584055795,"identity":"81907f8e-1002-47d6-9606-e158d3e3a6a3","order_by":5,"name":"Zhen Guan","email":"","orcid":"","institution":"First People's Hospital of Yunnan Province","correspondingAuthor":false,"prefix":"","firstName":"Zhen","middleName":"","lastName":"Guan","suffix":""},{"id":584055796,"identity":"f88b4e20-93f9-4c99-b831-106427e8db90","order_by":6,"name":"Li Zhang","email":"","orcid":"","institution":"First People's Hospital of Yunnan Province","correspondingAuthor":false,"prefix":"","firstName":"Li","middleName":"","lastName":"Zhang","suffix":""},{"id":584055797,"identity":"d2654de4-b51a-44dd-9c34-a088cb081acd","order_by":7,"name":"Xiao-chuan Wang","email":"","orcid":"","institution":"First People's Hospital of Yunnan Province","correspondingAuthor":false,"prefix":"","firstName":"Xiao-chuan","middleName":"","lastName":"Wang","suffix":""},{"id":584055798,"identity":"7995b04e-7f35-4a42-9d92-3f2f76946454","order_by":8,"name":"Zhi-qiong Wang","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAnElEQVRIiWNgGAWjYBACAyCWYKiQkJMnTcuBMxbGhg0kaTnYVpHIcIBYLeYSyQ9vf5wnkcDYwPzw0Q1itFjOSDO2OLhNIo+dgc3YOIcoh91IMJMAailmbOBhkyZSS/o3iYNzJBIbDhCvJQdoSwNJWs68KbY4c0zC2LCZaL8cT994o6KmTk6evfnhY6K0IAAzacpHwSgYBaNgFOADAI7IMOISU1sDAAAAAElFTkSuQmCC","orcid":"","institution":"First People's Hospital of Yunnan Province","correspondingAuthor":true,"prefix":"","firstName":"Zhi-qiong","middleName":"","lastName":"Wang","suffix":""}],"badges":[],"createdAt":"2026-01-29 02:38:42","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-8726237/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-8726237/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":101675526,"identity":"0fd6e452-3b51-43c0-9396-a63e4252bbf3","added_by":"auto","created_at":"2026-02-02 13:29:12","extension":"jpg","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":1220214,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eScreening immune related hub genes in cutaneous squamous cell carcinoma (cSCC).\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eA. Volcano plot of differentially expressed genes (DEGs) between cSCC and normal samples in GSE42677. Thresholds: |log2FC| ≥ 0.3, p \u0026lt; 0.05.\u003c/p\u003e\n\u003cp\u003eB. Heatmap showing expression of top DEGs. Orange represents high expression; cyan represents low expression.\u003c/p\u003e\n\u003cp\u003eC. Heatmap of immune infiltration scores across samples. Orange indicates high infiltration; cyan indicates low infiltration.\u003c/p\u003e\n\u003cp\u003eD. Correlation heatmap between WGCNA modules and immune infiltration phenotypes. Orange represents positive correlation; cyan represents negative correlation.\u003c/p\u003e\n\u003cp\u003eE. Venn diagram showing overlap between WGCNA blue module genes and DEGs.\u003c/p\u003e\n\u003cp\u003eF. KEGG pathway enrichment analysis of intersected genes; only significantly enriched pathways are shown (p \u0026lt; 0.05).\u003c/p\u003e","description":"","filename":"fig1.jpg","url":"https://assets-eu.researchsquare.com/files/rs-8726237/v1/717921056cb6dced881dc1a5.jpg"},{"id":101675521,"identity":"59959d81-d8a8-4d3a-9586-cbef8f57713a","added_by":"auto","created_at":"2026-02-02 13:29:12","extension":"jpg","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":837861,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eHub gene selection and immune related molecular subtype analysis.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eA. Protein protein interaction (PPI) network of intersected genes. Hub genes with degree \u0026gt; 15 are highlighted as potential key regulators in the immune microenvironment.\u003c/p\u003e\n\u003cp\u003eB. Hierarchical clustering of 10 cSCC samples from GSE45164 based on intersected genes in the PPI network, separating patients into two subtypes (C1 and C2).\u003c/p\u003e\n\u003cp\u003eC. Volcano plot of DEGs between C1 and C2 Thresholds: |log2FC| ≥ 0.3, p \u0026lt; 0.05.\u003c/p\u003e\n\u003cp\u003eD. Heatmap of DEGs between C1 and C2 Orange indicates high expression; cyan indicates low expression.\u003c/p\u003e\n\u003cp\u003eE. Heatmap of immune infiltration scores in C1 vs. C2 Orange indicates high scores; cyan indicates low scores.\u003c/p\u003e\n\u003cp\u003eF. KEGG enrichment analysis of DEGs between subtypes; only significant pathways are shown (p \u0026lt; 0.05).\u003c/p\u003e\n\u003cp\u003eG. PPI network of DEGs between C1 and C2.\u003c/p\u003e","description":"","filename":"fig2.jpg","url":"https://assets-eu.researchsquare.com/files/rs-8726237/v1/437bd54ec6a1d7f683de262f.jpg"},{"id":101753364,"identity":"e088e1ed-05fc-4f65-ad59-a090119508aa","added_by":"auto","created_at":"2026-02-03 10:39:52","extension":"jpg","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":641876,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eT cell landscape in tumor and normal tissues from GSE144236.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eA. UMAP plot showing overall cell type distribution in GSE144236.\u003c/p\u003e\n\u003cp\u003eB. UMAP plots comparing cell distributions between normal and tumor tissues.\u003c/p\u003e\n\u003cp\u003eC. UMAP of T cell subpopulations.\u003c/p\u003e\n\u003cp\u003eD. UMAP showing T cell distributions in tumor vs. normal tissues.\u003c/p\u003e\n\u003cp\u003eE. Left: Stacked bar plot of CD4⁺ and CD8⁺ T cell proportions in normal and tumor tissues; Right: stacked bar plot of all T cell subpopulations within CD4⁺ and CD8⁺ compartments.\u003c/p\u003e\n\u003cp\u003eF. UMAP of CD8⁺ T cell clusters based on Seurat analysis.\u003c/p\u003e\n\u003cp\u003eG. UMAP showing distributions of CD8⁺ T cell subpopulations.\u003c/p\u003e\n\u003cp\u003eH. UMAP comparing CD8⁺ T cell distributions in tumor vs. normal tissues.\u003c/p\u003e\n\u003cp\u003eI. Stacked bar plot of CD8⁺ T cell subpopulation proportions in tumor and normal tissues.\u003c/p\u003e\n\u003cp\u003eJ. Proportions of the four CD8⁺ T cell subpopulations in normal vs. tumor tissues.\u003c/p\u003e","description":"","filename":"fig3.jpg","url":"https://assets-eu.researchsquare.com/files/rs-8726237/v1/606b71b2bf8e766f723fb57f.jpg"},{"id":101754353,"identity":"0c73244e-837f-432c-a5e7-3be4e0a307f3","added_by":"auto","created_at":"2026-02-03 10:42:14","extension":"jpg","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":999236,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eFunctional analyses of CD8⁺ T cell subpopulations based on DEGs.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eA. DEG analysis of four CD8⁺ T cell subpopulations (Naïve, EM, EMRA, Exh); thresholds: |log2FC| ≥ 0.3, padj \u0026lt; 0.05.\u003c/p\u003e\n\u003cp\u003eB. KEGG pathway enrichment of DEGs in Exh subpopulation; only significant pathways are shown.\u003c/p\u003e\n\u003cp\u003eC. PPI network of Exh DEGs; hub genes selected based on degree.\u003c/p\u003e\n\u003cp\u003eD. Expression distribution of key Exh pathway genes in UMAP (green: high, purple: low).\u003c/p\u003e\n\u003cp\u003eE. KEGG enrichment of EMRA DEGs.\u003c/p\u003e\n\u003cp\u003eF. PPI network of EMRA DEGs; hub genes highlighted.\u003c/p\u003e\n\u003cp\u003eG. Expression distribution of key EMRA pathway genes in UMAP.\u003c/p\u003e\n\u003cp\u003eH. KEGG enrichment of EM DEGs.\u003c/p\u003e\n\u003cp\u003eI. PPI network of EM DEGs; hub genes highlighted.\u003c/p\u003e\n\u003cp\u003eJ. Expression distribution of key EM pathway genes in UMAP.\u003c/p\u003e\n\u003cp\u003eK. KEGG enrichment of Naïve DEGs.\u003c/p\u003e\n\u003cp\u003eL. PPI network of Naïve DEGs; hub genes highlighted.\u003c/p\u003e\n\u003cp\u003eM. Expression distribution of key Naïve pathway genes in UMAP.\u003c/p\u003e","description":"","filename":"fig4.jpg","url":"https://assets-eu.researchsquare.com/files/rs-8726237/v1/fe3cb60fc6347bc9499f903a.jpg"},{"id":101675536,"identity":"b40bc919-5015-4a19-b6f8-1266c4f7c665","added_by":"auto","created_at":"2026-02-02 13:29:12","extension":"jpg","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":945864,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eDifferentiation potential and pseudotime trajectory analysis of CD8⁺ T cells.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eA. CytoTRACE predicted cell ordering overlaid on UMAP; color indicates differentiation potential.\u003c/p\u003e\n\u003cp\u003eB. CytoTRACE UMAP showing spatial distributions of different CD8⁺ T cell subpopulations.\u003c/p\u003e\n\u003cp\u003eC. Average CytoTRACE scores of each CD8⁺ T cell subpopulation, ordered from high to low.\u003c/p\u003e\n\u003cp\u003eD. Monocle3 pseudotime UMAP; cells colored by pseudotime.\u003c/p\u003e\n\u003cp\u003eE. Scatter plots of key Exh pathway gene expression across four CD8⁺ T cell subpopulations.\u003c/p\u003e\n\u003cp\u003eF. Scatter plots of key Exh pathway gene expression along Monocle3 pseudotime.\u003c/p\u003e\n\u003cp\u003eG. Slingshot pseudotime trajectory UMAP.\u003c/p\u003e\n\u003cp\u003eH. UMAP showing expression of key Exh genes along trajectory (green: high, purple: low).\u003c/p\u003e","description":"","filename":"fig5.jpg","url":"https://assets-eu.researchsquare.com/files/rs-8726237/v1/ba58c7247ce9a6c3f5386f12.jpg"},{"id":101675532,"identity":"edef474f-8073-4ecb-b2c5-e3ebc8048bd6","added_by":"auto","created_at":"2026-02-02 13:29:12","extension":"jpg","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":510550,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eComparison of key gene expression at singlecell and bulk transcriptome levels.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eA. Expression of six key genes (CCL5, CD3D, CTSD, KLRD1, PTPN6, TNFRSF18, TNFRSF1B) in CD8⁺ T cells at the single cell level. Gene expression was compared between tumor (Tumor) and normal (Normal) tissues. Statistical significance was assessed using a t-test and indicated as follows: *p \u0026lt; 0.05, **p \u0026lt; 0.01, ***p \u0026lt; 0.001, ****p \u0026lt; 0.0001.\u003c/p\u003e\n\u003cp\u003eB. Expression of six key genes in cSCC vs normal skin at bulk transcriptome level.\u003c/p\u003e\n\u003cp\u003eC. Correlation between pseudotime trajectory and key Exh gene expression in CD8⁺ T cells.\u003c/p\u003e","description":"","filename":"fig6.jpg","url":"https://assets-eu.researchsquare.com/files/rs-8726237/v1/afcb044fce3913a0fd7823c0.jpg"},{"id":101675535,"identity":"4393f74a-5f96-4bdd-b9dc-394452bfb3ac","added_by":"auto","created_at":"2026-02-02 13:29:12","extension":"jpg","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":649242,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eValidation of KLRD1, PTPN6, and CTSD expression in cSCC and localization of KLRD1 in exhausted CD8⁺ T cells.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eA. qPCR validation of KLRD1, PTPN6, and CTSD expression in cSCC vs paired normal skin; significance indicated (*p \u0026lt; 0.05, **p \u0026lt; 0.01).\u003c/p\u003e\n\u003cp\u003eB. Representative immunofluorescence (IF) images showing CD8, PD-1, and KLRD1 expression and colocalization in cSCC tissues. From left to right: merged image (Merge), CD8, PD-1, and KLRD1. Scale bars maintained across images.\u003c/p\u003e\n\u003cp\u003eC. Quantitative analysis of IF staining comparing tumor and paired normal tissues for KLRD1⁺ cells, CD8⁺KLRD1⁺ cells, CD8⁺PD-1⁺ cells, and CD8⁺KLRD1⁺PD-1⁺ cells. Despite non-significant differences, all cell types show a trend of increased abundance in tumor tissues.\u003c/p\u003e","description":"","filename":"fig7.jpg","url":"https://assets-eu.researchsquare.com/files/rs-8726237/v1/615d1ce352a423867efae9ed.jpg"},{"id":102320103,"identity":"73d25b37-b5a2-41a3-83f4-5663d2e5485e","added_by":"auto","created_at":"2026-02-10 13:27:20","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":7188702,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-8726237/v1/2eaa8d1e-f557-44ea-81a6-d5590c8f8321.pdf"},{"id":101675537,"identity":"280b538d-e0fc-4cce-8052-4c41b1061c4e","added_by":"auto","created_at":"2026-02-02 13:29:13","extension":"xlsx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":9781,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eTable S1. Primer sequences used for quantitative PCR (qPCR).\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe table lists forward and reverse primer sequences for genes analyzed in this study. Primers were designed to specifically amplify target transcripts, and qPCR was performed to assess gene expression levels in cSCC and paired normal skin samples.\u003c/p\u003e","description":"","filename":"TableS1.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-8726237/v1/44c572a14fe313725b827e1c.xlsx"}],"financialInterests":"No competing interests reported.","formattedTitle":"Integrative transcriptomic and single-cell analyses identify KLRD1 enrichment in exhausted CD8⁺ T cells in cutaneous squamous cell carcinoma","fulltext":[{"header":"Introduction","content":"\u003cp\u003eCutaneous squamous cell carcinoma (cSCC) is a common skin cancer, and while most lesions are manageable with surgery, a subset progresses to more aggressive disease\u003csup\u003e\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e\u003c/sup\u003e. Despite its prevalence, the molecular and immune mechanisms underlying cSCC progression are not fully understood.\u003c/p\u003e \u003cp\u003eTumor infiltrating lymphocytes (TILs), particularly CD8⁺ T cells, are generally considered important mediators of antitumor immunity\u003csup\u003e\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e\u003c/sup\u003e. However, high immune infiltration does not always correlate with effective tumor control\u003csup\u003e\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u003c/sup\u003e. Tumors can induce functional restraint or exhaustion in T cells, leading to heterogeneous immune responses within the tumor microenvironment\u003csup\u003e\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e,\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u003c/sup\u003e. Such functional heterogeneity has been reported in other solid tumors, where immune infiltration coexists with inhibitory signaling and immune regulation\u003csup\u003e\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eIn cSCC, chronic inflammation and dysregulated immune signaling are common features of the tumor microenvironment\u003csup\u003e\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e\u003c/sup\u003e. Studies have reported infiltration of CD8⁺ T cells, regulatory T cells, and natural killer cells, alongside evidence of immune suppression and variable cytotoxic activity\u003csup\u003e\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e\u003c/sup\u003e. Despite this, the diversity of CD8⁺ T cell states including effector, memory, and exhausted populations and the transcriptional programs associated with functional restraint remain incompletely characterized\u003csup\u003e\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e\u003c/sup\u003e. Moreover, whether molecular subtypes with distinct immune microenvironments exist in cSCC, and how they relate to potential responses to immunotherapy, has not been systematically investigated. Although CD8⁺ T cell exhaustion has received increasing attention as an important factor shaping antitumor immunity and immunotherapy response, its cellular heterogeneity and molecular characteristics in cSCC remain poorly defined\u003csup\u003e\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eRecent studies have further suggested that exhausted CD8⁺ T cells represent a heterogeneous and dynamically regulated population, rather than a uniform dysfunctional state\u003csup\u003e\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e,\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e\u003c/sup\u003e. However, how such exhaustion-related programs are organized and maintained within the cSCC tumor microenvironment remains unclear.\u003c/p\u003e \u003cp\u003eTo address these questions, we performed an integrative analysis combining bulk transcriptomics, immune infiltration profiling, and single-cell RNA sequencing. Our study reveals intertumoral immune heterogeneity, identifies molecular subtypes with distinct immune profiles, and characterizes transcriptional programs associated with CD8⁺ T cell exhaustion. Among these, KLRD1 emerges as a candidate marker linked to exhaustion states. These results provide insight into the immune landscape of cSCC and may inform strategies to improve immune-based interventions in this cancer type.\u003c/p\u003e"},{"header":"Methods","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003eData Acquisition and Preprocessing\u003c/h2\u003e \u003cp\u003ePublic transcriptomic data (GSE42677 and GSE45164) were downloaded from the GEO database using the GEOquery R package (v2.70.0). Raw expression matrices were extracted via the exprs function, and probe IDs were annotated with gene symbols using the platform annotation file. For probes mapping to multiple gene symbols, only the first symbol was retained. Duplicate gene symbols and missing values were removed using dplyr (v1.1.4) and na.omit. The final gene expression matrix, with genes as rows and samples as columns, was log2-transformed (log2(x\u0026thinsp;+\u0026thinsp;1)) to stabilize variance and exported for downstream analyses.\u003c/p\u003e \u003c/div\u003e\n\u003ch3\u003eDifferential Expression Analysis\u003c/h3\u003e\n\u003cp\u003eDifferentially expressed genes (DEGs) between cutaneous squamous cell carcinoma (cSCC) and normal samples were identified using the limma (v3.58.1). A design matrix representing the sample groups was constructed, and linear models were fitted. Contrasts were defined for cSCC vs normal, followed by empirical Bayes moderation. DEGs were defined as those with |log2 fold change| \u0026gt; 0.3 and adjusted p-value\u0026thinsp;\u0026lt;\u0026thinsp;0.05. Volcano plots and heatmaps were generated using ggplot2 (v3.5.2), EnhancedVolcano (v1.20.0), and pheatmap (v1.0.13). The top significantly up- and down-regulated genes were selected for visualization.\u003c/p\u003e\n\u003ch3\u003eHierarchical Clustering of Key Genes\u003c/h3\u003e\n\u003cp\u003eExpression values of predefined key genes (including FN1, MMP9, CXCL8, ESR1, PXDN, TIMP1, CXCL10, TLR2, CXCL1, MMP3, SPP1, CD36, GZMB, SOX9, MMP1, KIT, GATA3, NT5E, BMP2, CXCL9, PLAUR, S100A12, COL3A1, MMP13, CXCL13) were converted to numeric format. Sample-wise Euclidean distances were calculated, and hierarchical clustering was performed using the Ward.D2 method via the hclust function. Clusters were visualized using factoextra (v1.0.7) and FactoMineR (v2.11), and dendrograms were annotated with cluster assignments.\u003c/p\u003e\n\u003ch3\u003eImmune Cell Infiltration Analysis (ssGSEA)\u003c/h3\u003e\n\u003cp\u003eGene sets representing immune cell types were constructed from CSV files of marker genes, with gene identifiers as row names. Each column was converted into a list element for single-sample Gene Set Enrichment Analysis (ssGSEA). Normalized gene expression matrices (genes \u0026times; samples) were log2-transformed (log2(TPM\u0026thinsp;+\u0026thinsp;1) or log2-CPM\u0026thinsp;+\u0026thinsp;1). ssGSEA was performed using the GSVA (v1.50.5) with the ssgsea method and kcdf = \"Gaussian\". Absolute gene ranking was applied (abs.ranking\u0026thinsp;=\u0026thinsp;TRUE), generating enrichment scores for each immune cell type per sample. Scores were normalized to a 0\u0026ndash;1 range using Min-Max scaling for comparability across samples. Raw and normalized ssGSEA scores were saved for downstream analyses.\u003c/p\u003e \u003cp\u003eImmune cell infiltration patterns were visualized using heatmaps, with rows scaled to highlight relative differences and columns grouped according to clinical classification. For detailed group comparisons, half-violin, half-point, and boxplot overlays were generated using ggplot2 and gghalves (v0.1.4). Data were reshaped into long format, and statistical significance between groups (e.g., C1 vs C2) was assessed using two-sided t-tests, with significance indicated on the plots. Custom colors were applied to violin fills and points.\u003c/p\u003e\n\u003ch3\u003eFunctional Enrichment Analysis\u003c/h3\u003e\n\u003cp\u003eKEGG pathway enrichment analysis of differentially expressed genes (DEGs) or selected gene modules was first performed using KOBAS (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://bioinfo.org/kobas/genelist/\u003c/span\u003e\u003cspan address=\"http://bioinfo.org/kobas/genelist/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) to identify significantly enriched pathways. Genes were mapped to KEGG pathway terms, and enrichment significance was assessed using the hypergeometric test with multiple testing correction (FDR\u0026thinsp;\u0026lt;\u0026thinsp;0.05). The resulting pathway lists, including the number of input genes per pathway and associated p-values, were subsequently visualized using WeightedTreemaps (v0.1.4) in R. Voronoi-based treemaps were constructed, with cell sizes proportional to the number of genes mapped to each pathway, and colors representing distinct pathway categories. Gradient color schemes were applied to indicate the level of enrichment, facilitating intuitive interpretation of pathway significance and gene involvement.\u003c/p\u003e \u003cdiv id=\"Sec8\" class=\"Section2\"\u003e \u003ch2\u003eProtein\u0026ndash;Protein Interaction (PPI) Network Analysis\u003c/h2\u003e \u003cp\u003eTo explore potential interactions among differentially expressed genes (DEGs), protein\u0026ndash;protein interaction (PPI) networks were constructed using Cytoscape (v3.10.4). DEGs were uploaded to the STRING database (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://cn.string-db.org/\u003c/span\u003e\u003cspan address=\"https://cn.string-db.org/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) to obtain predicted and experimentally validated protein interactions, with a confidence score cutoff of 0.4 (medium confidence). The resulting network data were imported into Cytoscape for visualization and further analysis. Key modules or clusters within the network were identified with the MCODE, and functional enrichment of the hub genes or modules was annotated using KEGG and Gene Ontology databases. Graphical layouts and node colors were adjusted in Cytoscape to reflect interaction strength and gene significance.\u003c/p\u003e \u003c/div\u003e\n\u003ch3\u003eSingle Cell RNA-seq Data Processing and Quality Control\u003c/h3\u003e\n\u003cp\u003eRaw count data from GSE144236 were imported into R (v5.1.0) and converted to a sparse matrix for memory efficiency. A Seurat (v5.2.1) object was created, and patient metadata were integrated. Cells with low gene counts or high mitochondrial gene percentage were filtered out to ensure high quality single cell data. Normalization was performed using the LogNormalize method with a scale factor of 10,000. Highly variable genes (n\u0026thinsp;=\u0026thinsp;2,000) were identified using the vst method, and all genes were scaled. Principal component analysis (PCA) was performed on the top variable genes, and the first 10 PCs were used to construct the K-nearest neighbor (KNN) graph (FindNeighbors) and cluster cells (FindClusters, resolution\u0026thinsp;=\u0026thinsp;0.5). Dimensionality reduction was visualized using UMAP and t-SNE.\u003c/p\u003e\n\u003ch3\u003eT Cell Extraction and Subpopulation Identification\u003c/h3\u003e\n\u003cp\u003eT cells were extracted from the full dataset based on level1_celltype == \"Tcell\" and filtered for cells with \u0026ge;\u0026thinsp;500 detected genes. PCA and clustering were repeated on this subset, and UMAP was used for visualization. Subpopulations were annotated, and the proportion of each T cell subtype per group (Tumor vs Normal) was calculated. Differences in proportions were assessed using t-tests implemented via ggsignif (v0.6.4).\u003c/p\u003e \u003cp\u003eCD8⁺ T cells (CD8_EM, CD8_EMRA, CD8_Exh, CD8_Naive) were further subsetted and reanalyzed. PCA (npcs\u0026thinsp;=\u0026thinsp;15), KNN graph construction, and clustering (resolution\u0026thinsp;=\u0026thinsp;0.6) were repeated. UMAP and t-SNE visualizations confirmed the expression of canonical markers CD8A and CD8B using FeaturePlot and VlnPlot. Proportions of CD8⁺ T cell subtypes in Tumor and Normal samples were visualized with stacked bar plots.\u003c/p\u003e \u003cdiv id=\"Sec11\" class=\"Section2\"\u003e \u003ch2\u003eDifferential Gene Expression and KEGG Pathway Enrichment\u003c/h2\u003e \u003cp\u003eFor each CD8⁺ T cell subpopulation, differential expression analysis was conducted using Seurat\u0026rsquo;s FindMarkers function with the Wilcoxon rank-sum test. Genes with adjusted p-value\u0026thinsp;\u0026lt;\u0026thinsp;0.05 and |log2 fold change| \u0026gt; 0.3 were considered significant.\u003c/p\u003e \u003cp\u003eSignificantly up- and down-regulated genes were subjected to KEGG pathway enrichment using clusterProfiler (v4.10.1). Enriched pathways were visualized as bubble plots in ggplot2, where bubble size represents the number of input genes and color indicates enrichment significance (-log10 p-value). Pathways were manually ordered according to biological relevance.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec12\" class=\"Section2\"\u003e \u003ch2\u003eSingle-Cell Trajectory and Pseudotime Analysis\u003c/h2\u003e \u003cp\u003eCD8⁺ T cell differentiation was analyzed using Monocle3 (v1.3.7), CytoTRACE (v0.3.3), and Slingshot (v2.10.0). In Monocle3, Seurat RNA counts were converted to a cell_data_set object, cells were clustered, and pseudotime trajectories were inferred; dynamically expressed genes along pseudotime were identified with Moran\u0026rsquo;s I and visualized. CytoTRACE estimated differentiation potential, visualized on UMAP with gene-specific patterns. Slingshot reconstructed pseudotime trajectories on a SingleCellExperiment object using Seurat cluster labels, and trajectories were visualized on UMAP coordinates.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec13\" class=\"Section2\"\u003e \u003ch2\u003eGene Expression Dynamics along Pseudotime\u003c/h2\u003e \u003cp\u003eTo examine dynamic expression patterns, gene expression values were extracted from the Seurat RNA assay and matched with Monocle3-calculated pseudotime. Spearman correlation coefficients were computed to quantify the association between gene expression and pseudotime. Expression trends were visualized using scatter plots with pseudotime on the x-axis and gene expression on the y-axis, with LOESS curves fitted to highlight overall trajectories. Statistical significance of correlations was reported, and results were saved for downstream integration.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec14\" class=\"Section2\"\u003e \u003ch2\u003eQuantitative PCR (qPCR) Analysis\u003c/h2\u003e \u003cp\u003eTotal RNA was extracted from cells or tissues using TRIzol reagent (Invitrogen) according to the manufacturer\u0026rsquo;s protocol. RNA quality and concentration were assessed with a NanoDrop spectrophotometer (Thermo Fisher). Complementary DNA (cDNA) was synthesized using the PrimeScript RT reagent kit (Takara). qPCR was performed on a QuantStudio 6 Flex system (Applied Biosystems) using SYBR Green Master Mix (Takara). Relative gene expression was calculated using the 2^\u0026minus;ΔΔCt method with GAPDH as the internal control. Primers used for each target gene are listed in Supplementary Table \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e. Experiments were performed in triplicate.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec15\" class=\"Section2\"\u003e \u003ch2\u003eImmunofluorescence (IF) Staining\u003c/h2\u003e \u003cp\u003eCells or tissue sections were fixed with 4% paraformaldehyde for 15 min and permeabilized with 0.1% Triton X-100 for 10 min. Samples were blocked with 5% bovine serum albumin (BSA) for 1 h at room temperature, followed by incubation with primary antibodies overnight at 4\u0026deg;C. After washing, samples were incubated with appropriate fluorophore-conjugated secondary antibodies for 1 h at room temperature in the dark. Nuclei were counterstained with DAPI. Images were captured using a confocal laser scanning microscope and analyzed with ZEN (v3.3).\u003c/p\u003e \u003c/div\u003e"},{"header":"Results","content":"\u003cdiv id=\"Sec17\" class=\"Section2\"\u003e \u003ch2\u003eDifferential gene expression and immune cell infiltration patterns in cSCC\u003c/h2\u003e \u003cp\u003eDifferential expression analysis was performed on the GSE42677 dataset to compare cutaneous squamous cell carcinoma (cSCC) tissues with normal skin samples. Using the criteria of |log2FC| \u0026ge; 0.3 and \u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05, a total of 2,149 differentially expressed genes (DEGs) were identified, including 1,161 upregulated and 988 downregulated genes in cSCC (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eA). The expression patterns of the top DEGs were visualized using a heatmap. Genes such as CDHR2, PKD2L1, SLC30A3, MORN1, SPOCK2, VIPR2, LIMK1, LHX5, and KIAA1614 were markedly upregulated in cSCC tissues, whereas CP, PTH1R, FLT4, DSPP, PCDHA9, PIP5K1B, RELN, S1PR1, and CMAHP showed relatively lower expression compared with normal skin samples (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eB). To characterize the immune microenvironment of cSCC, immune cell infiltration scores were evaluated and visualized as a heatmap. Increased infiltration scores were observed in cSCC samples for CD56dim natural killer cells, plasmacytoid dendritic cells, activated CD8⁺ T cells, macrophages, type 1 T helper cells, regulatory T cells, myeloid-derived suppressor cells, natural killer cells, natural killer T cells, γδ T cells, CD56bright natural killer cells, type 17 T helper cells, neutrophils, T follicular helper cells, and activated dendritic cells. In contrast, higher infiltration levels of activated B cells, mast cells, central memory CD4⁺ T cells, memory B cells, activated CD4⁺ T cells, type 2 T helper cells, effector memory CD8⁺ T cells, eosinophils, central memory CD8⁺ T cells, monocytes, immature B cells, effector memory CD4⁺ T cells, and immature dendritic cells were predominantly observed in normal skin tissues (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eC). Weighted gene co-expression network analysis (WGCNA) was further conducted to explore gene modules associated with immune infiltration patterns. The MEblue module exhibited significant negative correlations with the majority of immune cell infiltration scores, including activated CD8⁺ T cells (r\u0026thinsp;=\u0026thinsp;\u0026minus;\u0026thinsp;0.76, \u003cem\u003ep\u003c/em\u003e\u0026thinsp;=\u0026thinsp;1\u0026times;10⁻⁴), activated CD4⁺ T cells (r\u0026thinsp;=\u0026thinsp;\u0026minus;\u0026thinsp;0.70, \u003cem\u003ep\u003c/em\u003e\u0026thinsp;=\u0026thinsp;6\u0026times;10⁻⁴), and natural killer T cells (r\u0026thinsp;=\u0026thinsp;\u0026minus;\u0026thinsp;0.86, \u003cem\u003ep\u003c/em\u003e\u0026thinsp;=\u0026thinsp;2\u0026times;10⁻⁶) (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eD). Intersection analysis between DEGs and genes within the MEblue module identified 331 overlapping genes (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eE). Functional enrichment analysis of these intersecting genes revealed significant enrichment in immune- and inflammation-related pathways, including cytokine\u0026ndash;cytokine receptor interaction, complement and coagulation cascades, JAK\u0026ndash;STAT signaling pathway, leukocyte transendothelial migration, NF-κB signaling pathway, TGF-β signaling pathway, Th17 cell differentiation, NOD-like receptor signaling pathway, Th1 and Th2 cell differentiation, TNF signaling pathway, Toll-like receptor signaling pathway, chemokine signaling pathway, viral protein interaction with cytokine and cytokine receptor, and IL-17 signaling pathway (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eF).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec18\" class=\"Section2\"\u003e \u003ch2\u003eIdentification of immune-related hub genes and T cell\u0026ndash;enriched molecular subtypes in cSCC\u003c/h2\u003e \u003cp\u003eProtein\u0026ndash;protein interaction (PPI) network analysis was performed based on the intersecting genes, and hub genes were identified using a degree threshold greater than 15. A total of 25 hub genes were screened, including PXDN, PLAUR, CXCL13, CXCL9, MMP1, MMP13, NT5E, SOX9, GZMB, MMP3, CXCL1, GATA3, KIT, TIMP1, COL3A1, S100A12, SPP1, CXCL8, CXCL10, BMP2, and FN1. Among these, genes such as FN1, MMP9, CXCL8, MMP3,PXDN, CXCL10, TLR2, TIMP1, ESR1, and CXCL1 occupied central positions within the PPI network, suggesting their potential key roles in regulating the tumor immune microenvironment of cSCC (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eA). Based on the expression patterns of these hub genes, hierarchical clustering analysis was applied to cSCC samples from the GSE45164 dataset, resulting in two distinct molecular subtypes, termed C1 and C2 (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eB). Differential expression analysis between the two subtypes identified 113 differentially expressed genes, including 58 upregulated and 55 downregulated genes (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eC). Genes such as KRT34, FMO2, ZBED2, DUOX1, CIDEA, DLX2, IL18, MAP2, FETUB, and IL37 were upregulated in the C1 subtype, whereas IL1B, CXCL1, G0S2, SRGN, ADAMDEC1, CXCL6, SELL, EVI2B, FCGR3B, and TDO2 were downregulated (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eD). To further characterize the immune features of the two molecular subtypes, immune infiltration levels were compared between C1 and C2. The immune infiltration heatmap revealed that type 2 T helper cells, CD56bright natural killer cells, plasmacytoid dendritic cells, and mast cells exhibited higher infiltration scores in the C1 subtype. In contrast, enhanced infiltration of activated B cells, effector memory CD8⁺ T cells, activated CD8⁺ T cells, and activated CD4⁺ T cells was observed in the C2 subtype (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eE). Functional enrichment analysis of the subtype-specific differentially expressed genes highlighted significant enrichment in immune- and inflammation-related pathways, including cytokine\u0026ndash;cytokine receptor interaction, TNF signaling pathway, IL-17 signaling pathway, NF-κB signaling pathway, chemokine signaling pathway, NOD-like receptor signaling pathway, Th17 cell differentiation, natural killer cell\u0026ndash;mediated cytotoxicity, JAK\u0026ndash;STAT signaling pathway, and Toll-like receptor signaling pathway (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eF). Furthermore, a PPI network was constructed based on the subtype-specific differentially expressed genes, in which immune-associated genes such as IL18, MMP3, CXCL2, IL1B, FCGR3A, SELL, FCGR3B, CXCL1, CD69, IL15, CCL18, CXCL13, LTF, IL6, PTGS2, SPP1, and CCL20 were identified as central nodes (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eG), suggesting the presence of key immune regulatory modules distinguishing the two cSCC subtypes.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec19\" class=\"Section2\"\u003e \u003ch2\u003eDifferential Immune Infiltration and Expansion of Exhausted CD8⁺ T Cells in cSCC\u003c/h2\u003e \u003cp\u003eTo establish a comprehensive cellular landscape, we first performed UMAP visualization of the entire dataset, revealing 14 distinct cell types and clear separation between normal and tumor-derived cells (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eA\u0026ndash;B). Focusing on T cells, nine T cell subtypes were identified and their distribution across normal and tumor tissues was characterized (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eC\u0026ndash;D). The proportion of CD8⁺ T cells was markedly higher in tumor tissues compared to normal tissues. Quantitative comparison showed that in tumors, the proportions of CD8_Exh, CD8_Trge, CD8_EMRA, and CD8_EM were notably higher than in normal tissue, whereas CD8_Na\u0026iuml;ve and CD4_RGCC/CD4_Na\u0026iuml;ve were more abundant in normal tissue (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eE). For a refined analysis, CD8⁺ T cells were extracted and clustered into five clusters, comprising CD8_Exh, CD8_Na\u0026iuml;ve, CD8_EM, and CD8_EMRA subpopulations (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eF\u0026ndash;G). Their spatial distribution revealed that CD8_Na\u0026iuml;ve cells predominantly resided in normal tissue, while the other three clusters were mainly found in tumor tissue (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eH). Quantitative proportion analysis confirmed the UMAP-based observations, with CD8_Na\u0026iuml;ve cells enriched in normal samples and CD8_Exh, CD8_EMRA, and CD8_EM enriched in tumors (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eI\u0026ndash;J).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec20\" class=\"Section2\"\u003e \u003ch2\u003eFunctional characteristics and key gene expression of CD8⁺ T cell subsets in cSCC\u003c/h2\u003e \u003cp\u003eDifferential expression analysis was first performed among the four CD8⁺ T cell subsets (Naive, EM, EMRA, and Exh) with a threshold of |log2FC| \u0026ge; 0.3 and adjusted p-value\u0026thinsp;\u0026lt;\u0026thinsp;0.05, yielding subset specific differentially expressed genes (DEGs) (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eA). Each subset\u0026rsquo;s DEGs were subjected to KEGG pathway enrichment analysis, and only significantly enriched pathways (p\u0026thinsp;\u0026lt;\u0026thinsp;0.05) are shown to highlight the main functional features of each subset. STRING database was used to construct protein\u0026ndash;protein interaction (PPI) networks for the DEGs, with hub genes identified based on degree values. UMAP based expression visualization was then performed to map key pathway related genes onto single-cell data, illustrating their spatial distribution across subsets and thereby reflecting functional distinctions. In the Exh subset, enriched pathways included antigen processing and presentation, apoptosis, carbon metabolism, cellular senescence, chemokine signaling, C-type lectin receptor signaling, cytokine\u0026ndash;cytokine receptor interaction, glycolysis/gluconeogenesis, Jak-STAT signaling, NF-kappa B signaling, oxytocin signaling, p53 signaling, T cell receptor signaling, TNF signaling, and Toll-like receptor signaling (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eB). Hub genes in Exh included GAPDH, UBA52, RPS27A, EEF2, and RPS2 (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eC). Notably, PTPN6, CD3D, TNFRSF18, CTSD, CCL5, KLRD1, CSF1, and TNFRSF1B were highly expressed in the CD8_Exh subset (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eD). For the CD8_EMRA subset, enriched pathways included natural killer cell mediated cytotoxicity, cytokine\u0026ndash;cytokine receptor interaction, chemokine signaling, antigen processing and presentation, T cell receptor signaling, MAPK signaling, NF-kappa B signaling, TNF signaling, cell cycle, and apoptosis (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eE). Hub genes comprised RPL3, RPL11, RPL5, RPS16, and RPS14 (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eF), with KLRD1, DUSP2, FCGR3A, and FGR showing higher expression within this subset (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eG). In the CD8_EM subset, significant pathways included MAPK signaling, chemokine signaling, apoptosis, cytokine\u0026ndash;cytokine receptor interaction, lysosome, natural killer cell-mediated cytotoxicity, and antigen processing and presentation (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eH). Hub genes comprised GAPDH, CXCR3, GZMB, SELL, GZMK, and RPL13A (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eI), with CXCR4, CXCR3, DUSP2, and HSPA1B particularly highly expressed (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eJ). For the CD8_Naive subset, enriched pathways included cell adhesion molecules (CAMs), chemokine signaling, hematopoietic cell lineage, intestinal immune network for IgA production, oxytocin signaling, PD-L1 expression and PD-1 checkpoint pathway in cancer, primary immunodeficiency, riboflavin metabolism, and T cell receptor signaling (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eK). Hub genes included GAPDH, IFNG, CD4, CD44, CD8A, and GZMB (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eL), with CD44, CD4, CCR6, and IFNGR2 highly expressed in this subset (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eM).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003cb\u003eTrajectory Analysis Reveals Potential Key Genes CCL5, CD3D, CTSD, KLRD1, and PTPN6 in CD8⁺_Exh Cells in cSCC\u003c/b\u003e \u003c/p\u003e \u003cp\u003eCytoTRACE analysis was used to predict the differentiation potential of CD8⁺ T cells, showing that naive cells exhibited the highest scores, followed by Exh, EM, and EMRA cells in descending order (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eA\u0026ndash;C). UMAP visualization of differentiation potential revealed a clear gradient across subsets, providing a basis for subsequent pseudotime analysis. Monocle3 pseudotime analysis was performed to construct single-cell trajectories, validating the differentiation order predicted by CytoTRACE. Two main trajectories were identified: one progressing from naive \u0026rarr; Exh, and the other from naive \u0026rarr; EM \u0026rarr; EMRA (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eD). This analysis allowed quantification of dynamic cellular states along the differentiation path. Expression of key Exh pathway genes was examined across subsets and along pseudotime. Genes including CCL5, CD3D, CTSD, KLRD1, PTPN6, TNFRSF18, and TNFRSF1B were most highly expressed in CD8_Exh cells (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eE). Visualization of these genes along the Monocle3 pseudotime trajectories reflected their dynamic changes, highlighting their regulatory roles during the exhaustion and effector differentiation process (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eF). Slingshot trajectory analysis was conducted as an independent validation of the Monocle3 results. UMAP visualization showed the spatial distribution and expression gradients of key genes along the trajectories, emphasizing the involvement of Exh-related genes at the terminal stage of differentiation. Two main trajectories consistent with Monocle3 were observed (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eG). UMAP mapping of Exh key genes along the Slingshot trajectories further confirmed that these genes increased in expression at the end of the differentiation trajectory corresponding to CD8_Exh cells. In particular, CCL5, CD3D, CTSD, KLRD1, and PTPN6 exhibited elevated expression in the CD8⁺_Exh subset (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eH).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec21\" class=\"Section2\"\u003e \u003ch2\u003eIntegrated validation identifies KLRD1 as a leading candidate gene in exhausted CD8⁺ T cells\u003c/h2\u003e \u003cp\u003eTo validate whether the key genes identified by single-cell analysis exhibited significant expression differences in CD8⁺ T cells and to assess the consistency of their expression patterns with bulk transcriptomic data, we performed integrative analyses across multiple levels. At the single-cell level, the expression of seven candidate genes (CCL5, CD3D, CTSD, KLRD1, PTPN6, TNFRSF18, and TNFRSF1B) was quantified in CD8⁺ T cells from tumor and normal tissues. All seven genes showed significantly higher expression in tumor-derived CD8⁺ T cells compared with those from normal skin (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eA). At the bulk transcriptomic level, differential expression analysis between cSCC and normal skin tissues demonstrated that CCL5, CD3D, KLRD1, and TNFRSF1B exhibited expression patterns consistent with the single-cell results, with significantly elevated expression in cSCC samples (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eB), further supporting the robustness of these candidate genes. To further evaluate the association of Exh-related genes with CD8⁺ T cell differentiation and exhaustion, Monocle3 pseudotime analysis was performed, focusing on CTSD, KLRD1, and PTPN6. Spearman correlation analysis between gene expression and pseudotime revealed that KLRD1 showed a strong positive correlation with pseudotime (rho\u0026thinsp;=\u0026thinsp;0.504, p\u0026thinsp;=\u0026thinsp;1.75 \u0026times; 10⁻⁴\u0026sup1;), indicating progressively increased expression during CD8⁺ T cell exhaustion and suggesting its suitability as a primary candidate for subsequent functional studies. CTSD displayed a weak but significant positive correlation (rho\u0026thinsp;=\u0026thinsp;0.091, p\u0026thinsp;=\u0026thinsp;0.0228), supporting its potential role as a secondary candidate, whereas PTPN6 showed no significant correlation with pseudotime (rho\u0026thinsp;=\u0026thinsp;0.061, p\u0026thinsp;=\u0026thinsp;0.126) (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eC). Collectively, these analyses provide a dynamic view of key Exh-related genes during CD8⁺ T cell differentiation and exhaustion, offering a rational basis for prioritizing candidate genes for downstream functional validation and for exploring their roles in the immune microenvironment of cutaneous squamous cell carcinoma.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec22\" class=\"Section2\"\u003e \u003ch2\u003eKLRD1 is enriched in exhausted CD8⁺ T cells in cSCC\u003c/h2\u003e \u003cp\u003eFirst, quantitative PCR was performed on paired cutaneous squamous cell carcinoma (cSCC) and adjacent normal skin tissues to validate the expression changes of KLRD1, PTPN6, and CTSD, thereby assessing the consistency between single cell transcriptomic findings and bulk tissue samples. All three genes were significantly upregulated in cSCC compared with normal tissues (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eA). Next, triple immunofluorescence (IF) staining was performed to simultaneously examine the expression and spatial distribution of CD8 (cyan), PD-1 (red), and KLRD1 (yellow) at the tissue level. This analysis aimed to determine whether KLRD1 was localized within CD8⁺ T cells in the tumor microenvironment and to evaluate its spatial association with the exhaustion marker PD-1. Representative images are shown in the order of merged view, CD8, PD-1, and KLRD1 channels, with normal skin tissues displayed in the upper panels and cSCC tissues in the lower panels. In normal skin tissues, CD8⁺ T cells were sparsely distributed, with weak or barely detectable PD-1 and KLRD1 signals. In contrast, cSCC tissues exhibited a marked accumulation of CD8⁺ T cells (cyan), accompanied by pronounced PD-1 (red) and KLRD1 (yellow) staining. Importantly, merged images revealed frequent co-localization of KLRD1 signals with CD8⁺ cells, and substantial overlap between KLRD1 and PD-1 within the same cellular regions. This spatial pattern suggests that KLRD1 is preferentially expressed in PD-1⁺ CD8⁺ T cells within the tumor microenvironment. Overall, these IF results demonstrate that CD8, PD-1, and KLRD1 are all enriched in cSCC tissues compared with normal skin and display a coordinated spatial distribution consistent with an exhausted CD8⁺ T cell phenotype (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eB). Quantitative analysis of IF images was further performed to assess the densities of KLRD1⁺, CD8⁺KLRD1⁺, CD8⁺PD-1⁺, and CD8⁺KLRD1⁺PD-1⁺ cells in tumor and normal tissues. Although statistical significance was not reached due to limited sample size and the rarity of certain cell subsets, all four metrics consistently exhibited higher levels in cSCC tissues. This trend was concordant with the single-cell transcriptomic findings indicating elevated KLRD1 expression in exhausted CD8⁺ T cells, thereby providing supportive evidence for the involvement of KLRD1 in the cSCC immune microenvironment (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eC).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e"},{"header":"Discussion","content":"\u003cdiv id=\"Sec24\" class=\"Section2\"\u003e \u003ch2\u003eBulk transcriptomic analyses reveal immune activation and molecular heterogeneity in cSCC\u003c/h2\u003e \u003cp\u003eOur integrative analyses of bulk transcriptomic data and immune infiltration profiles indicate that cSCC is far from an \u0026ldquo;immune cold\u0026rdquo; tumor; instead, it exhibits pronounced immune and inflammatory activation alongside substantial molecular heterogeneity\u003csup\u003e\u003cspan additionalcitationids=\"CR13\" citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e\u003c/sup\u003e. The extensive dysregulation of inflammatory and immune related pathways including TNF, NF-κB, IL-17, and JAK\u0026ndash;STAT signaling underscores the interplay between tumorigenesis and local immune remodeling\u003csup\u003e\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e,\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e\u003c/sup\u003e. These findings are consistent with reports, where chronic inflammation not only drives tumor progression but also shapes the immune contexture and therapeutic responsiveness\u003csup\u003e\u003cspan additionalcitationids=\"CR17\" citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eImmune infiltration analyses revealed notable intertumoral variability: tumor tissues were enriched in activated CD8⁺ T cells, regulatory T cells, and NK cells, whereas memory B and effector memory CD8⁺ T cells in normal skin retained cytotoxic and immune surveillance related transcriptional programs\u003csup\u003e\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e,\u003cspan additionalcitationids=\"CR20\" citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e\u003c/sup\u003e. This observation highlights a critical distinction: immune infiltration in tumors does not automatically translate into effective antitumor activity, as tumor-infiltrating lymphocytes display transcriptional features of functional restraint or exhaustion, in contrast to the cytotoxic and immune surveillance programs retained by effector memory CD8⁺ T cells in normal skin. In this context, infiltrating T cells may be functionally restrained by intrinsic transcriptional programs, such as the MEblue module identified in our weighted gene coexpression network analysis. Notably, the MEblue module exhibited significant negative correlations with the majority of immune cell infiltration scores, including activated CD8⁺ T cells, activated CD4⁺ T cells, and natural killer T cells, suggesting that this module may represent a broad immunoregulatory or compensatory transcriptional program. Such inverse associations imply the presence of regulatory circuits that constrain immune activation despite substantial immune infiltration. Similar patterns have been reported in other tumor types, where high immune cell abundance coexists with inhibitory signaling and immune exclusion, including non\u0026ndash;small cell lung cancer\u003csup\u003e\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e\u003c/sup\u003e. Stratification of cSCC samples based on immune related hub genes revealed two molecular subtypes with distinct immune landscapes. Subtype C1 was enriched for immunoregulatory populations, including Th2 cells, plasmacytoid dendritic cells, and mast cells, suggesting a more tolerogenic or modulatory microenvironment. In contrast, C2 exhibited heightened cytotoxic and inflammatory signatures. This mirrors the \u0026ldquo;immune contexture\u0026rdquo; concept in other solid tumors, where composition, activation state, and spatial distribution of tumor-infiltrating lymphocytes predict immunotherapy outcomes and disease progression\u003csup\u003e\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e\u003c/sup\u003e. Notably, the existence of these subtypes in cSCC may partially explain variable clinical responses to immune checkpoint blockade observed in clinical cohorts.\u003c/p\u003e \u003cdiv id=\"Sec25\" class=\"Section3\"\u003e \u003ch2\u003eDynamic expansion and exhaustion programming of CD8⁺ T cells in cSCC\u003c/h2\u003e \u003cp\u003eSingle cell profiling further revealed that CD8⁺ T cell expansion in cSCC is accompanied by pronounced functional heterogeneity. Tumor infiltrating CD8⁺ T cells largely comprise exhausted (Exh), effector memory (EM), and terminal effector (EMRA) populations, whereas na\u0026iuml;ve cells remain largely confined to normal skin. The preferential accumulation of exhausted and terminal effector cells suggests sustained antigenic stimulation and immunosuppressive cues from the tumor microenvironment\u003csup\u003e\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e,\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eFunctional analyses showed that CD8_Exh cells activate TNF, NF-κB, TCR, and p53 pathways while undergoing metabolic and inflammatory reprogramming, whereas CD8_EM and EMRA populations are enriched for cytotoxic, NK-like, MAPK, and cell cycle\u0026ndash;associated pathways. CD8_Na\u0026iuml;ve cells predominantly express adhesion molecules and baseline TCR/PD-1 signaling. Collectively, these data suggest that exhaustion is an actively regulated, rather than passive, state a perspective increasingly recognized in tumor immunology, where partial effector function persists to balance immune surveillance with tissue protection\u003csup\u003e\u003cspan additionalcitationids=\"CR27\" citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eTrajectory analyses using CytoTRACE and Monocle3 revealed that exhaustion develops progressively along two primary differentiation axes: na\u0026iuml;ve \u0026rarr; Exh and na\u0026iuml;ve \u0026rarr; EM \u0026rarr; EMRA, with exhausted cells occupying terminal branches. Key genes, including CCL5, CD3D, CTSD, KLRD1, and PTPN6, exhibit gradual upregulation along pseudotime, highlighting a coordinated transcriptional program driving exhaustion. These findings support the emerging view that exhaustion markers should be interpreted in a temporal and regulatory context rather than as static indicators of dysfunction, reflecting dynamic metabolic and transcriptional remodeling within exhausted CD8⁺ T cells\u003csup\u003e\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec26\" class=\"Section3\"\u003e \u003ch2\u003eKLRD1 as a marker and potential modulator of exhaustion in cSCC\u003c/h2\u003e \u003cp\u003eAmong exhaustion-associated genes, KLRD1 stands out as a robust candidate. Its expression increases progressively along the exhaustion trajectory, spatially colocalizes with CD8⁺ and PD-1⁺ cells, and is detectable at both bulk and tissue levels. Prior studies in HNSC and other squamous cell carcinomas have linked KLRD1to terminally exhausted CD8⁺ and NK cells, where it participates in inhibitory signaling to prevent overactivation\u003csup\u003e\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e\u003c/sup\u003e. In this context, KLRD1 may function not only as a biomarker of exhaustion but also as a regulatory element within a transcriptional module that integrates chronic antigenic stimulation, effector differentiation, and metabolic adaptationa mechanism increasingly recognized across tumor immunology studies.\u003c/p\u003e \u003cp\u003eThe potential translational significance of KLRD1 is multifold. Beyond serving as a biomarker of exhaustion, its dynamic expression profile and tissue-level detectability make it a candidate for stratifying patients or guiding immunotherapeutic interventions. For example, combinatorial therapies targeting PD-1 alongside modulation of KLRD1 mediated pathways could enhance cytotoxic T cell function while limiting deleterious overactivation, although this remains speculative and requires functional validation.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec27\" class=\"Section3\"\u003e \u003ch2\u003eLimitations and future perspectives\u003c/h2\u003e \u003cp\u003eDespite the multi-layered evidence supporting KLRD1 as an exhaustion-associated marker, several limitations exist. Immunofluorescence relied on a single marker, limiting resolution of functional states at the single-cell level, and sample size constraints prevented statistically robust tissue quantification. Furthermore, the mechanistic role of KLRD1 in regulating exhaustion remains untested in vitro or in vivo. Future work should assess its functional contribution to T cell exhaustion, interactions with canonical checkpoints (e.g., PD-1, LAG-3), and predictive value for immunotherapy response. Expanding this integrative framework to larger cohorts and other tumor types will further refine strategies to exploit immune heterogeneity for clinical benefit.\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e"},{"header":"Conclusion","content":"\u003cp\u003eTogether, our integrative analyses reveal a highly dynamic and heterogeneous immune landscape in cSCC, identify KLRD1 as a key exhaustion-associated gene, and provide mechanistic and translational insights into CD8⁺ T cell dysfunction. These findings offer a foundation for future studies aimed at modulating exhaustion pathways to improve immunotherapeutic efficacy in cSCC and potentially other solid tumors.\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eAuthor Contributions:\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eDan-dan Zou conceived and designed the study. Meng-en Li contributed to methodology development. Yi-chao Jin performed software implementation and provided resources. Li-tao Suo carried out experimental validation. Li Zhang conducted formal data analysis. Yuan Chen was responsible for investigation and data curation. Dan-dan Zou and Meng-en Li wrote the original draft. Xiao-chuan Wang and Zhi-qiong Wang contributed to writing\u0026mdash;review and editing. Zhen Guan prepared the visualizations. Zhi-qiong Wang supervised the study. Xiao-chuan Wang managed the project. Dan-dan Zou acquired funding for the study.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAcknowledgments:\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThis work was supported by Kunming University of Science and Technology \u0026amp; the First People\u0026apos;s Hospital of Yunnan Province Joint Special Project on Medical Research (KUST-KH2023023Y), Yunnan Province Science and Technology Department \u0026amp; Kunming Medical University Joint Special Project for Applied Basic Research (202501AY070001-261), Yunnan Fundamental Research Projects (grant NO. 202401AU070030), High-level Talents Scientific Research Project of Yunnan Provincial Health Commission (2023-KHRCBZ-B13) and The First People\u0026apos;s Hospital of Yunnan Province (KHYJ-2023-6-08).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eData Availability Statement:\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eAll the data analyzed in this study were obtained from publicly available repositories (GEO database).\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eConflicts of Interest:\u003c/strong\u003e\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eThe authors have no relevant financial or non-financial interests to dis-close.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eEthics and Consent to Participate declarations:\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eNot applicable.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\u003cli\u003e\u003cspan\u003eWinge MCG, et al. Advances in cutaneous squamous cell carcinoma. Nat Rev Cancer. 2023;23:430\u0026ndash;49. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.1038/s41568-023-00583-5\u003c/span\u003e\u003cspan address=\"10.1038/s41568-023-00583-5\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLi F, et al. The association between CD8\u0026thinsp;+\u0026thinsp;tumor-infiltrating lymphocytes and the clinical outcome of cancer immunotherapy: A systematic review and meta-analysis. EClinicalMedicine. 2021;41:101134. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.1016/j.eclinm.2021.101134\u003c/span\u003e\u003cspan address=\"10.1016/j.eclinm.2021.101134\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePag\u0026egrave;s F, et al. Immune infiltration in human tumors: a prognostic factor that should not be ignored. Oncogene. 2010;29:1093\u0026ndash;102. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.1038/onc.2009.416\u003c/span\u003e\u003cspan address=\"10.1038/onc.2009.416\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHeintzman DR, Fisher EL, Rathmell JC. Microenvironmental influences on T cell immunity in cancer and inflammation. Cell Mol Immunol. 2022;19:316\u0026ndash;26. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.1038/s41423-021-00833-2\u003c/span\u003e\u003cspan address=\"10.1038/s41423-021-00833-2\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eJiang W, et al. Exhausted CD8\u0026thinsp;+\u0026thinsp;T Cells in the Tumor Immune Microenvironment: New Pathways to Therapy. Front Immunol. 2020;11:622509. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.3389/fimmu.2020.622509\u003c/span\u003e\u003cspan address=\"10.3389/fimmu.2020.622509\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eOuyang P, et al. Infiltration characteristics and regulatory mechanisms of CD8(+) T lymphocytes in solid tumors: spatial distribution, biological functions, and interactions with the immune microenvironment. Front Immunol. 2025;16:1661545. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.3389/fimmu.2025.1661545\u003c/span\u003e\u003cspan address=\"10.3389/fimmu.2025.1661545\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBuruiană A, et al. The Tumor Stroma of Squamous Cell Carcinoma: A Complex Environment That Fuels Cancer Progression. Cancers (Basel). 2024;16. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.3390/cancers16091727\u003c/span\u003e\u003cspan address=\"10.3390/cancers16091727\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eFarhood B, Najafi M, Mortezaee K. CD8(+) cytotoxic T lymphocytes in cancer immunotherapy: A review. J Cell Physiol. 2019;234:8509\u0026ndash;21. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.1002/jcp.27782\u003c/span\u003e\u003cspan address=\"10.1002/jcp.27782\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBaessler A, Vignali D. A. T cell exhaustion. Annu Rev Immunol. 2024;42:179\u0026ndash;206.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePhilip M, Schietinger A. CD8\u0026thinsp;+\u0026thinsp;T cell differentiation and dysfunction in cancer. Nat Rev Immunol. 2022;22:209\u0026ndash;23.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLi C, Yuan Y, Wang Q, Jiang X. Epigenetic Regulation of CD8⁺ T Cell Exhaustion: Recent Advances and Update. Front Immunol. 2025;16:1700039.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eRajesh A, et al. Characterisation of the immune microenvironment of cutaneous squamous cell carcinoma in immunosuppression. Exp Dermatol. 2022;31:1720\u0026ndash;8. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.1111/exd.14650\u003c/span\u003e\u003cspan address=\"10.1111/exd.14650\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLuci C, et al. Cutaneous Squamous Cell Carcinoma Development Is Associated with a Temporal Infiltration of ILC1 and NK Cells with Immune Dysfunctions. J Invest Dermatol. 2021;141:2369\u0026ndash;79. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.1016/j.jid.2021.03.018\u003c/span\u003e\u003cspan address=\"10.1016/j.jid.2021.03.018\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDong Q, Zhang Z, Li S, Liang L. Mechanisms of immunotherapy in cutaneous squamous cell carcinoma in the tumor microenvironment. Front Immunol. 2025;16:1660272. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.3389/fimmu.2025.1660272\u003c/span\u003e\u003cspan address=\"10.3389/fimmu.2025.1660272\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSaeidi V, Doudican N, Carucci JA. Understanding the squamous cell carcinoma immune microenvironment. Front Immunol. 2023;14:1084873. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.3389/fimmu.2023.1084873\u003c/span\u003e\u003cspan address=\"10.3389/fimmu.2023.1084873\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGreten FR, Grivennikov SI. Inflammation and Cancer: Triggers, Mechanisms, and Consequences. Immunity. 2019;51:27\u0026ndash;41. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.1016/j.immuni.2019.06.025\u003c/span\u003e\u003cspan address=\"10.1016/j.immuni.2019.06.025\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eNishida A, Andoh A. The Role of Inflammation in Cancer: Mechanisms of Tumor Initiation, Progression, and Metastasis. Cells. 2025;14. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.3390/cells14070488\u003c/span\u003e\u003cspan address=\"10.3390/cells14070488\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWen Y, et al. Chronic inflammation, cancer development and immunotherapy. Front Pharmacol. 2022;13:1040163. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.3389/fphar.2022.1040163\u003c/span\u003e\u003cspan address=\"10.3389/fphar.2022.1040163\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eStrobl J, Haniffa M. Functional heterogeneity of human skin-resident memory T cells in health and disease. Immunol Rev. 2023;316:104\u0026ndash;19. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.1111/imr.13213\u003c/span\u003e\u003cspan address=\"10.1111/imr.13213\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eTokura Y, Phadungsaksawasdi P, Kurihara K, Fujiyama T, Honda T. Pathophysiology of Skin Resident Memory T Cells. Front Immunol. 2020;11:618897. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.3389/fimmu.2020.618897\u003c/span\u003e\u003cspan address=\"10.3389/fimmu.2020.618897\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eChen X, Zheng Y, Man X, Li W. Tissue-resident memory T cells and their function in skin diseases. Chin Med J (Engl). 2025;138:1175\u0026ndash;83. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.1097/cm9.0000000000003499\u003c/span\u003e\u003cspan address=\"10.1097/cm9.0000000000003499\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eAnichini A, Perotti VE, Sgambelluri F, Mortarini R. Immune Escape Mechanisms in Non Small Cell Lung Cancer. Cancers (Basel). 2020;12. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.3390/cancers12123605\u003c/span\u003e\u003cspan address=\"10.3390/cancers12123605\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWu T, Wu X, Wang HY, Chen L. Immune contexture defined by single cell technology for prognosis prediction and immunotherapy guidance in cancer. Cancer Commun (Lond). 2019;39:21. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.1186/s40880-019-0365-9\u003c/span\u003e\u003cspan address=\"10.1186/s40880-019-0365-9\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eXie H, Xi X, Lei T, Liu H, Xia Z. CD8(+) T cell exhaustion in the tumor microenvironment of breast cancer. Front Immunol. 2024;15:1507283. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.3389/fimmu.2024.1507283\u003c/span\u003e\u003cspan address=\"10.3389/fimmu.2024.1507283\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhang B, et al. CD8(+) T cell exhaustion and its regulatory mechanisms in the tumor microenvironment: key to the success of immunotherapy. Front Immunol. 2024;15:1476904. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.3389/fimmu.2024.1476904\u003c/span\u003e\u003cspan address=\"10.3389/fimmu.2024.1476904\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eChow A, Perica K, Klebanoff CA, Wolchok JD. Clinical implications of T cell exhaustion for cancer immunotherapy. Nat Rev Clin Oncol. 2022;19:775\u0026ndash;90. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.1038/s41571-022-00689-z\u003c/span\u003e\u003cspan address=\"10.1038/s41571-022-00689-z\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBaessler A, Vignali DA. A. T Cell Exhaustion. Annu Rev Immunol. 2024;42:179\u0026ndash;206. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.1146/annurev-immunol-090222-110914\u003c/span\u003e\u003cspan address=\"10.1146/annurev-immunol-090222-110914\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eJenkins E, Whitehead T, Fellermeyer M, Davis SJ, Sharma S. The current state and future of T-cell exhaustion research. Oxf Open Immunol. 2023;4:iqad006. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.1093/oxfimm/iqad006\u003c/span\u003e\u003cspan address=\"10.1093/oxfimm/iqad006\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhong T, Sun S, Zhao M, Zhang B, Xiong H. The mechanisms and clinical significance of CD8(+) T cell exhaustion in anti-tumor immunity. Cancer Biol Med. 2025;22:460\u0026ndash;80. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.20892/j.issn.2095-3941.2024.0628\u003c/span\u003e\u003cspan address=\"10.20892/j.issn.2095-3941.2024.0628\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDong C, Lin Z, Hu Y, Lu Q. KLRD1 (CD94): A Prognostic Biomarker and Therapeutic Candidate in Head and Neck Squamous Cell Carcinoma. J Cancer. 2025;16:982\u0026ndash;95. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.7150/jca.104762\u003c/span\u003e\u003cspan address=\"10.7150/jca.104762\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"Cutaneous squamous cell carcinoma (cSCC), CD8⁺ T cell exhaustion, KLRD1, Immune microenvironment","lastPublishedDoi":"10.21203/rs.3.rs-8726237/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-8726237/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eCutaneous squamous cell carcinoma (cSCC) is characterized by complex interactions between tumor cells and the immune microenvironment. Here, we performed integrated transcriptomic and single-cell analyses to dissect immune-related molecular features and identify key regulators of CD8⁺ T cell exhaustion in cSCC. Differential expression analysis of GSE42677 identified 2,149 genes significantly dysregulated in cSCC compared with normal skin, with enrichment in immune- and inflammation-related pathways. Immune cell infiltration profiling revealed distinct patterns between tumor and normal tissues, and weighted gene co-expression network analysis (WGCNA) highlighted modules negatively correlated with activated CD8⁺ T cells and other immune populations. Protein\u0026ndash;protein interaction network analysis and hierarchical clustering defined two molecular subtypes of cSCC with divergent immune landscapes.\u003c/p\u003e \u003cp\u003eSingle-cell RNA sequencing analysis further revealed expansion of exhausted CD8⁺ T (CD8_Exh) cells in tumors. Subpopulation-specific pathway enrichment and hub gene analysis demonstrated distinct functional programs across CD8⁺ T cell subsets, including antigen processing, apoptosis, MAPK and NF-κB signaling, and cytokine interactions. Pseudotime and trajectory analyses using CytoTRACE, Monocle3, and Slingshot identified differentiation trajectories from naive to exhausted and effector CD8⁺ T cells, highlighting dynamic expression of key genes such as CCL5, CD3D, CTSD, KLRD1, and PTPN6 in the CD8_Exh subset. Spearman correlation with pseudotime indicated a strong positive association for KLRD1, suggesting its role in driving CD8⁺ T cell exhaustion.\u003c/p\u003e \u003cp\u003eExperimental validation using quantitative PCR and immunofluorescence confirmed elevated expression of KLRD1 in cSCC tissues. Collectively, these results provide a comprehensive landscape of immune regulation in cSCC and identify KLRD1 as a candidate marker of exhausted CD8⁺ T cells, offering a potential target for immunotherapeutic strategies.\u003c/p\u003e","manuscriptTitle":"Integrative transcriptomic and single-cell analyses identify KLRD1 enrichment in exhausted CD8⁺ T cells in cutaneous squamous cell carcinoma","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2026-02-02 13:29:07","doi":"10.21203/rs.3.rs-8726237/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"c3f42ea6-8a34-4356-8704-77e0e55d5f36","owner":[],"postedDate":"February 2nd, 2026","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[],"tags":[],"updatedAt":"2026-03-24T17:53:53+00:00","versionOfRecord":[],"versionCreatedAt":"2026-02-02 13:29:07","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-8726237","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-8726237","identity":"rs-8726237","version":["v1"]},"buildId":"XKTyCvWXoU3ODBz1xrDgd","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}