scRNA-seq data for ESCC were retrieved from the GEO database (accession number: GSE145370), comprising tumor samples from three stage II and four stage III ESCC patients. We followed the data processing workflow as described by Li H et al. (39), and the identification of canonical markers was also guided by original publication (40). T-cell subclustering was conducted using a resolution parameter of 1.0, followed by the extraction of CD8⁺T cells for downstream computational interrogation.

To explore the lineage trajectories of CD8⁺T-cell subsets, we conducted pseudotime analyses using both Monocle (version 2) and Slingshot, applying default parameters. The Seurat object was converted into a CellDataSet object for Monocle analysis. Dimensionality reduction was performed using the DDRTree algorithm, and cells were ordered along inferred developmental trajectories based on the HVGs identified earlier. For statistical modeling of gene expression variance, estimateSizeFactors, estimateDispersions, and dispersionTable functions were applied with default settings to characterize dynamic gene expression changes during differentiation.

A human ESCC tissue microarray (TMA; Catalog no. HEsoS180Su05) was provided by Shanghai Outdo Biotech Co., Ltd. (Shanghai, China). The array included 105 primary tumor tissues and 75 adjacent normal tissues from patients aged 48 to 82 years. Survival status was determined through clinical follow-up. Following the exclusion of cases with incomplete survival information or tissue core loss, the final analysis incorporated 80 tumor specimens and 75 matched adjacent normal tissues. Survival data from these 80 patients were used for subsequent survival analysis. All procedures involving human specimens were reviewed and approved by the Clinical Research Ethics Committee of Outdo Biotech (Approval No. SHYJS-CP-1807012), ensuring compliance with institutional and national ethical requirements.

mIHC was performed following the manufacturer’s protocol and as detailed in prior publications (41–43). Quantitative image acquisition and analysis were conducted using automated tools provided by PerkinElmer.

Whole-slide multispectral imaging and subsequent image analysis, including software and version details, were performed as previously described in our earlier study (43). Within this platform, fluorescent signals were spectrally unmixed, and each fluorophore was isolated into a separate channel and archived individually. DAPI staining was employed as a nuclear reference to guide accurate segmentation of cellular compartments, including nuclei, cytoplasm, and membranes.

Signal intensities for CK, KI67, CD8, TCF1, and PD-1 were co-registered with DAPI to generate binary expression masks for each biomarker. These masks were used to identify and quantify marker-positive cells. Tumor regions were delineated by pan-CK staining, which enabled precise discrimination between epithelial tumor cells and surrounding stromal compartments. Quantification of tumor cells was conducted based on the CK binary mask.

scRNA-seq data from ESCC patients receiving anti-PD-1 therapy were obtained from the Genome Sequence Archive (GSA) at the National Genomics Data Center, China, under accession number HRA003312 from project PRJCA012636 (44). The dataset comprised tumor specimens (baseline and post-intervention) from seven ESCC patients undergoing neoadjuvant chemoimmunotherapy with PD-1 blockade. Among these patients, three exhibited complete responses, while four were classified as non-responders.

We followed the data processing workflow as described by the original authors, including all preprocessing steps and code implementation as provided in the corresponding publication and associated repositories. Cell clustering and subsequent analyses were conducted in Scanpy (version 1.8.1), with annotation informed by canonical markers and validated against reference profiles reported in the original study. Specifically, T cell–enriched populations were reclustered using the Leiden algorithm at a resolution of 1.5 to achieve finer granularity of cellular subtypes.

To investigate the biological functions and signaling mechanisms associated with differential therapeutic responses, we applied R package clusterProfiler (version 4.8.3) to perform GO enrichment, KEGG pathway analyses, and GSEA. Annotations and hallmark gene sets were sourced from the MSigDB database (version 7.5.1). Pathways showing significant differences (P < 0.05) between complete responders (CR) and non-complete responders (NCR) were identified and visualized using bar plots and dot plots to illustrate enrichment trends.

Intercellular communication networks between T_pex cells and other immune cell populations were inferred using the CellChat R package (version 1.1.3) (41). The analysis was conducted on scRNA-seq datasets stratified by treatment response. The number and strength of inferred ligand–receptor interactions were quantified, and the global information flow of each signaling pathway was compared between CR and NCR groups. A minimum threshold of 10 cells per cell type was required for inclusion in the analysis. Significant signaling interactions, including those uniquely enriched in the CR group, were extracted for further investigation. Average expression levels of ligand–receptor pairs in CR-specific pathways were visualized using heatmaps generated by the ComplexHeatmap package (version 2.6.2) (42).

To categorize patients into high- and low-infiltration cohorts according to T_pex cell abundance, the optimal cut-off threshold was identified using the survminer R package via the log-rank test. Kaplan–Meier survival analyses were performed in GraphPad Prism 9, and univariate and multivariate Cox regression analyses were conducted in R (version 4.2.2) within the RStudio environment. For validation purposes, survival analyses based on T_pex-associated gene signatures in ESCC were executed using GEPIA2 (http://gepia2.cancer-pku.cn/#index), employing datasets from The Cancer Genome Atlas (TCGA).

Statistical methods were chosen according to the data distribution: normally distributed datasets were assessed with two-tailed, unpaired Student’s t-tests, whereas the Mann–Whitney U test was utilized for nonparametric comparisons between two groups. Statistical significance was denoted in figures as follows: *P < 0.05; **P < 0.01; ***P < 0.001; ns, not significant.

To elucidate the dynamic landscape and developmental hierarchy of CD8⁺T-cell subsets during ESCC progression, we analyzed a large-scale scRNA-seq dataset encompassing 102,611 cells. After reclustering, 29,073 CD8⁺T cells were identified and classified into six distinct subpopulations: T_pex, terminally exhausted, cytotoxic, memory, and proliferative CD8⁺T cells, each defined by canonical gene expression signatures ( Figures 1A–C ). These CD8⁺T cells were derived from matched tumor and adjacent normal tissues of seven ESCC patients. Notably, T_pex cells were present in both tumor and normal adjacent tissues, with a marginally higher frequency detected in the normal tissue compartment ( Figures 1D, E ).

Pseudotime trajectory inference revealed that T_pex and memory T-cell clusters were positioned at early stages of the differentiation continuum, whereas proliferative and terminally exhausted subsets occupied later pseudotime states ( Figures 1G, H ), supporting the notion that T_pex represents an early exhausted progenitor population. Furthermore, intercellular communication analysis via CellChat indicated that T_pex cells exhibited robust interactions with myeloid populations within the ESCC TME ( Figures 1F, I ), underscoring their potential immunoregulatory role.

mIHC was employed to simultaneously detect CK, CD8, TCF1, and PD-1 expression in ESCC tumor specimens and matched adjacent normal tissues ( Figures 2A, B ). Distinct CD8⁺T cell subpopulations were delineated based on the expression of TCF1 and PD-1, two key markers of differentiation states. TCF1⁺CD8⁺ cells are considered to possess stem-like properties, whereas PD-1⁺CD8⁺ cells are typically regarded as effector or exhausted populations. The TCF1⁺PD-1⁺CD8⁺ subset, co-expressing both markers, is identified as T_pex cells.

Quantitative analysis revealed a significant reduction in T_pex-cell infiltration within tumor tissues relative to adjacent normal tissues (P < 0.05; Figure 2C ). Survival analysis stratified by mIHC-defined T_pex abundance demonstrated that patients exhibiting high T_pex infiltration experienced significantly improved OS compared to those with low infiltration (cut-off = 0.05%, hazard ratio [HR] = 0.581, 95% confidence interval [CI]: 0.353–0.957, P < 0.05; Figure 2D ). Validation in the TCGA-ESCC cohort demonstrated a statistically significant association between increased T_pex abundance and improved disease-free survival (DFS) (P < 0.05; Figure 2E ), establishing its clinical relevance as a favorable prognostic indicator.

Clinicopathological correlations were assessed in a cohort of 80 ESCC patients. Chi-square (χ²) analysis revealed a significant association between T_pex-cell infiltration levels in tumor tissues and tumor size (P < 0.05), while no significant correlations were detected between T_pex-cell infiltration and other clinicopathological parameters. A significant relationship was observed between PD-1⁺CD8⁺T cell infiltration and patient gender. However, infiltration levels of total CD8⁺T cells and TCF1⁺CD8⁺T cells did not demonstrate any substantial correlations with clinicopathological features ( Table 1 ). Both univariate and multivariate Cox regression analyses were performed to evaluate prognostic factors for OS. As is shown, TNM stage remained an independent predictor in the multivariate model, with an HR of 3.410 (95% CI: 1.907–6.096, P < 0.001), confirming its strong association with poor prognosis. ( Table 2 ).

Patients were stratified by TNM stage: those at stage I or II were assigned to the early-stage group, whereas stages III and IV were categorized as advanced-stage. Qualitative assessment of mIHC images showed differential T_pex-cell infiltration between these groups, with representative examples presented ( Figures 3A, B ).

Early-stage patients (I+II) with abundant T_pex-cell infiltration showed a promising trend for prolonged survival relative to those with minimal infiltration in our Kaplan–Meier survival analyses ( Figure 3C ). A similar pattern was observed when patients were stratified by the T_pex to total CD8⁺T-cell ratio, where a higher ratio correlated with marginally better OS ( Figure 3D ).

Within the advanced-stage cohort (stages III and IV), patients exhibiting higher ratios of T_pex to TCF1⁺CD8⁺T cells demonstrated a trend toward extended OS compared to those with lower ratios ( Figure 3E ). Furthermore, an increased T_pex to PD-1⁺CD8⁺T-cell ratio was similarly linked to modestly improved survival outcomes in this group ( Figure 3F ).

Consistent with previous reports, we observed that T_pex infiltration correlated with improved clinical outcomes following ICB treatment in ESCC patients (44). We interrogated scRNA-seq data from project PRJCA012636. After rigorous quality control, 132,482 cells were annotated into 16 major populations using well-established marker gene profiles, encompassing neutrophils, monocytes, macrophages, B cells, endothelial and epithelial cells, plasmacytoid dendritic cells (pDCs), fibroblasts, pericytes, T cells, plasma cells, mast cells, and conventional dendritic cells (cDCs) ( Supplementary Figures 1A, B ).

To characterize immune heterogeneity relative to therapeutic response, we compared the proportions of major immune populations across pre- and post-treatment samples and between response groups ( Supplementary Figure 1C ). The cohort included three CR and four NCR. T cells constituted a predominant immune compartment, comprising over 50% of total cells in CR patients, a markedly higher proportion than in the NCR group.

Unsupervised clustering further subdivided the T-cell compartment into nine distinct subsets defined by their gene signatures, including CD8_T_rm (tissue resident memory), CD8_T_ex (terminally exhausted), CD8_MKi67 (proliferating), CD8_T_em (effector memory), CD8_T_pex, CD4_T_ex, CD4_T_cm (central memory), and regulatory T cells (T_reg) ( Supplementary Figures 1D, E ). Among these, CD8_T_pex emerged as a key transitional population bridging functional and exhausted states.

We visualized CD8_T_pex cells within the CD8⁺T-cell compartment using UMAP embeddings, stratified by treatment time points and clinical response categories ( Figures 4A–D ). The relative abundance of T_pex cells significantly increased post-therapy and was consistently higher in the CR group compared to NCR.

Pseudotime trajectory analysis revealed two primary differentiation pathways originating from the CD8_T_pex cluster, diverging into either terminal exhaustion or proliferation branches ( Figures 4E, F ), consistent with our earlier findings ( Figure 1E ). This underscored the central role of T_pex cells in dictating divergent CD8⁺T-cell fates.

Building on the observed differences in the proportion and trajectory of T_pex, we next assessed whether the microenvironmental signaling milieu differed between CR and NCR groups. Analysis of ligand–receptor interactions revealed that macrophages, plasma cells, and cDCs exhibited substantial differences in incoming signaling, while cDCs and macrophages showed differential outgoing signals ( Figure 5A ). CD8_T_pex cells in CR patients displayed enhanced outgoing interactions with macrophages and cDCs, implying their involvement in modulating the immune microenvironment to favor therapeutic response. Conversely, CD8_T_pex cells in CR patients received attenuated signals from T_reg and NK cells, suggesting a reduced immunosuppressive or cytotoxic regulatory influence.

We further evaluated the overall information flow within the TME. Most signaling pathways enriched in CR were intimately linked to immune activation, including pathways mediating T-cell activation (e.g., CD86, CD137), cytokine signaling (e.g., IL-2, IFN-γ), antigen presentation (e.g., MHC-II), as well as immune cell adhesion and trafficking (e.g., ICAM, SEMA4) ( Figure 5B ). These findings suggested a more immunologically active and coordinated microenvironment in CR patients.

Notably, the CD6 signaling pathway, which modulates T-cell activation, proliferation, and trafficking via interaction with its ligand ALCAM, was significantly upregulated in communications between CD8_T_pex cells and macrophages as well as cDCs (45) ( Figure 5C ). This enhanced CD6 signaling likely reflected an elevated activation status within the CR group. The ICOS pathway, known to regulate T helper cell responses and support adaptive immunity, was selectively activated in the interactions between CD8_T_pex and macrophages, as well as between CD8_T_pex and cDCs in CR patients. ( Figure 5D ).

By quantifying ligand–receptor pair ratios, we observed that distinct tissue regions exhibited unique intercellular signaling patterns. Further cell-type-specific analyses revealed that the interaction probability of ICOS-related ligand–receptor pairs, especially TNFSF9–TNFRSF9, in incoming signals to CD8_T_pex cells was newly induced or significantly upregulated in CR patients relative to NCR. Additionally, chemokine CCL family members were markedly enriched in the outgoing signaling from CD8_T_pex cells in CR patients ( Figure 6A ).

To gain deeper insights into these signaling modifications, pathway enrichment analyses were performed. The findings highlighted key biological processes such as T-cell differentiation, positive regulation of cell adhesion, and histone modification, collectively indicating enhanced immune activation, strengthened intercellular interactions, and possible epigenetic reprogramming of T_pex cells in responders ( Figure 6C ). Moreover, numerous RNA metabolic processes were significantly enriched.

Conversely, genes downregulated in CR patients were enriched for processes such as cytoplasmic translation, intrinsic apoptotic signaling, ribosome biogenesis, and responses to bacterial components or lipopolysaccharide (LPS), indicating a possible reduction in protein synthesis burden and stress-related signaling in T_pex cells from non-responders ( Figure 6D ). These observations were corroborated by gene set enrichment analysis (GSEA), which further supported elevated transcriptional and post-transcriptional activity in the CR group compared to a more suppressed or dysfunctional transcriptional landscape in the NCR group ( Figure 6B ).