The human CRC tissue microarray (TMA, catalog: HColA180Su16) was sourced from Shanghai Outdo Biotech Co., Ltd. (Shanghai, China). Specifically, the CRC TMA consisted of samples from 104 patients, namely 58 males and 46 females, ranging in age from 24 to 90 years. One patient was excluded due to the absence of clinical data. Survival data for all 104 patients were compiled and employed in the survival analysis. This study, which involved human participants, was reviewed and received approval from the Clinical Research Ethics Committee at Outdo Biotech (Shanghai, China, SHYJS-CP-1707004).

The mIHC staining was conducted using the Opal 7-color Manual IHC kit (Catalog No. NEL801001KT, PerkinElmer, USA), in conjunction with automated quantitative analysis tools from PerkinElmer, USA, according to the manufacturer’s instructions, consistent with our previously reported methodology (12). The primary objective of our current mIHC staining study was to identify cytokeratin (CK), CD8⁺T cells (CD8), activated T cells (PD-1), naïve/memory cells (TCF1), and proliferating cells (KI67) within CRC tissues and their adjacent normal counterparts. Nuclei were stained using 4’,6-diamidino-2-phenylindole (DAPI). The CRC TMA was first dried at 63°C for 1 hour. Subsequently, it was dewaxed in xylene and then rehydrated through a series of graded ethanol solutions. For antigen retrieval, we used a 1 mM EDTA solution at pH 9.0. After this, the tissue sections were allowed to cool at room temperature for 1 hour. Each section was treated with 100-200 µL of blocking solution and left to sit at room temperature for 15 minutes. We employed the following primary antibodies: anti-CD8 (catalog No. PA067, BioDot, USA), anti-TCF1 (1:50 dilution, catalog No. 2203, CST, USA), anti-PD-1 (catalog No. PA153, BioDot, USA), anti-KI67 (catalog No. GM724007, Gene Tech, China), and anti-CK (1:10 dilution, catalog No. PA125, BioDot, USA). The CRC TMA slide was then incubated with HRP-conjugated secondary antibodies (PerkinElmer, USA) in an Opal working solution (PerkinElmer, USA). The slide was mounted with ProLong Diamond Antifade Reagent with DAPI (Thermofisher, USA).

Firstly, panoramic multispectral scanning of the slides was conducted using the TissueFAXS system (Tissue Gnostics Asia Pacific Limited, Austria). The images obtained were then processed with the StrataQuest analysis software (Version No.7.0.1.165, Tissue Gnostics Asia Pacific Limited, Austria). In this software, each fluorophore was spectrally separated into individual channels and archived as distinct files. The DAPI staining served as a guide, ensuring precise differentiation between each cell’s nucleus, cytoplasm, and cell membrane. The expressions of CK, CD8, TCF1, PD-1, and KI67 were integrated with DAPI staining, resulting in binary masks highlighting cells that exhibited these specified biomarkers. Ultimately, the CK binary mask was analyzed to enumerate local tumor cells. For tissue recognition, Pan-CK was employed to segregate stromal nuclei from epithelial tumor tissue regions.

For the TCGA-COAD cohort, the RNA-seq and clinical data were downloaded using the TCGAbiolinks R package. The consensus molecular subtypes (CMS) were predicted from the count matrix using the CMScaller package. Only patients with assigned CMS groups underwent further analysis. The Gene Set Variation Analysis (GSVA) for the CD8⁺T cell subsets was executed using the GSVA package, with the enrichment analysis results visualized through ggplot2 package. The GSVA also performed enrichment analysis for selected chemokines (CXCL9/10/11) and their receptor (CXCR3). To determine the cutoff point for feature genes from the enrichment analysis, the surv_cutpoint function from the survminer package was applied. Based on this cutoff, patients were categorized into two groups: “High score” and “Low score”. Kaplan-Meier survival analysis was subsequently conducted using the survival package, and results were visualized with the ggsurvplot function.

Regarding data from the Synapse platform (Synapse: syn2623706) (16), gene IDs were mapped to symbols via the idmap2 R package. To mitigate potential batch effects among datasets, the ComBat function from the sva package was deployed. CMS predictions for each sample were executed using the CMSclassifier package. Patients without CMS assignments were excluded. The ensuing analyses mirrored those performed for the TCGA-COAD cohort. For rectal cancer patients with incomplete or complete responses to neo-adjuvant therapy (17), the GSVA was utilized to score CD8⁺T cell subsets, and the results were visualized with ggplot2.

Whole tumor data files for human CRC scRNA-seq were sourced from the Figshare website (https://figshare.com/articles/dataset/CRC_CD45_rda/14820318). We employed an unsupervised graph-based clustering algorithm through Seurat v4 (18) to categorize each cell by its gene expression profile. Diverging from the official FindVariableFeatures function, the R package DUBStepR (with default parameters) was employed to pinpoint highly variable genes. Cell cycling scores were determined via the CellCycleScore function, leveraging the gene set of the S and G2M phases from Seurat. The ScaleData function was then implemented to offset cell cycle effects in ensuing analyses. Principal component analysis (PCA) on the distinguished highly variable genes was undertaken, with the Harmony package being employed to mitigate batch effects. Dimensionality was reduced using the RunUMAP function on 30 principal components at a resolution of 0.1. Cells were then annotated in alignment with the clustering outcomes outlined in the article. Similarly, CD45⁺ data files for mouse MC38 scRNA-seq were retrieved from the Figshare website (https://figshare.com/articles/dataset/CRC_CD45_rda/14820318). The clustering methodology was analogous to that applied to the human scRNA-seq dataset. Data files pertaining to CD8⁺T cells for human CRC scRNA-seq were accessed from GEO DataSets (GSE108989). Following official protocols, we conducted unsupervised clustering. Cell annotation was based on the clustering results described in the article. The entire tumor data set for human CRC scRNA-seq, classified by CMS, was procured from the Synapse platform (Synapse: syn26844071) upon request. Clustering adhered to standard protocols, and cell annotation was reliant on the R package scibetR, using pre-processed scRNA-seq data as a reference.

Scores for CD8⁺T cells in treatment-related scRNA-seq datasets, specifically from colorectal and melanoma tumors (GSE205506 and GSE120575), were deduced via the AddModuleScore function. Subsequent visualizations were achieved using the R package ggplot2.

We conducted a transcriptional regulatory network analysis on CD8⁺T cells extracted from the mouse MC38 scRNA-seq dataset. Initially, the co-expression network was computed using GRNBoost2, available in the Python package pySCENIC. Subsequently, regulons were discerned utilizing the R package RcisTarget. The activity of each regulon was quantified for every cell through the R package AUCell. For dimensionality reduction of regulon activity, we adopted t-distributed stochastic neighbor embedding (t-SNE), executed via the runSCENIC_3_scoreCells function in the SCENIC R package (19). Visualization of the area under the curve (AUC) values for the regulons was facilitated using the DoHeatmap function within Seurat.

The scRNA-seq data from human CRC was classified by CMS and annotated using the R package, scibetR. These data were then employed to infer cell-cell interactions. Notably, the TAM-C1QC and CD8-GZMK cells were singled out for further analysis. To identify potential ligand/receptor interactions, we utilized the statistical_analysis function from the cellphonedb method, keeping the default parameters. The interactions identified were subsequently visualized using the dot_plot function from cellphonedb.

Spatial transcriptome sequencing data were downloaded according to the article (20). The data files were processed following the tutorial created by Åsa Björklund & Paulo Czarnewski (https://nbisweden.github.io/workshop-scRNAseq/labs/compiled/seurat/seurat_07_spatial.html

#Deconvolution). Then, we deconvoluted the spatially-indexed datasets by R package SCDC (21) and predicted single-cell distribution according to previously processed single-cell data.

According to the median value, the proportions of each subset in CRC TMA were divided into two groups: “Expression low” and “Expression high”. Subsequently, Kaplan-Meier survival curves and univariate Cox’s regression was calculated. Additionally, the correlation between CK subsets and CD8⁺T cell subsets was performed using GraphPad Prism 9 software. To further investigate the relationship among these variables and clinical data, multivariate Cox’s regression analysis was conducted using the coxph function from the R package survival. The result was visualized by R package forestplot.

Statistical analyses for Kaplan-Meier survival curves and univariate Cox’s regression were conducted using GraphPad Prism 9 software. Multivariate Cox’s regression, bulk RNA-seq, scRNA-seq, and ST data analyses were performed with R version 4.2.2 and RStudio software. The pySCENIC tool, a lightning-fast Python implementation of GRNBoost2, was executed using Python 3.7.6. Specific details about the significance tests can be found in the accompanying figure legends.

T_pex cells represent a unique subset of CD8⁺T cells, displaying traits reminiscent of both naïve and effector/exhausted cells. At the transcriptional level, T_pex cells preserve certain “stemness” attributes seen in naïve cells while concurrently showcasing partial effector functionalities typical of effector and exhausted cells (22). To delve deeper into the role of T_pex cells within the TME of CRC, we utilized single-cell transcriptome data from publicly accessible databases, covering both human and mouse samples (23, 24). This enabled us to scrutinize the gene expression and transcriptional regulation profile of T_pex cells. By leveraging on these scRNA-seq datasets, our aim was to identify the molecular signatures and controlling mechanisms governing T_pex cells in the CRC TME.

Initiating our investigation with human CRC scRNA-seq data, we observed from clustering results highlighted in the referenced article that CD8⁺T cells could be segmented into eight distinct subgroups ( Figure 1A ). Significantly, hallmark genes of T_pex cell (22), such as TCF7 and PDCD1, were chiefly co-expressed within the CD8_C04-GZMK subgroup ( Figures 1B, C ). This observation reaffirmed that the CD8_C04-GZMK cluster epitomized the T_pex cell subset. This finding aligned well with earlier studies indicating that T_pex cells predominantly express the granule-associated gene GZMK over GZMB (25).

Furthermore, we expanded our inquiry using scRNA-seq data from the mouse MC38 model. Relying on the clustering results outlined in the article, CD8⁺T cells were categorized into six distinct subgroups ( Figure 1D ). Importantly, the genes Tcf7, Pdcd1, and Gzmk were predominantly co-expressed within the mT09_CD8_Tem-Ccl5 subgroup, suggesting that this group exemplified the T_pex cell subset. In addition, the transcription factor Eomes was highly expressed in the T_pex cell subset, whereas Tbx21 had increased expression in the T_ex subset (25). Analyzing chemokine receptor expression revealed that the receptor Cxcr3, linked to T-cell recruitment, was up-regulated in the T_pex cell subset. Simultaneously, the receptor Cxcr4, associated with circulation, experienced a similar up-regulation. In contrast, the receptor Cxcr6, tied to T-cell residency (26), manifested a reduced expression within the T_pex cell subset ( Figure 1E ).

SCENIC analysis showed marked discrepancies in transcriptional regulation between lymph node (LN) and TILs. Moreover, specific T-cell populations presented diverse transcriptional landscapes, with the T_pex cell subset (mT09_CD8_Tem-Ccl5) exhibiting a transitional state between LN and TILs ( Figure 1F ). The heatmap detailing AUC values of regulons indicated that T_pex cells bore resemblances to naïve CD8⁺T cells concerning naive attributes, while it mirrored the T_ex cell subset regarding exhaustion characteristics ( Figure 1G ). This finding underscored that T_pex cells were discernible in both human and mouse CRC TMEs.

To bolster our findings, we performed mIHC staining on both CRC tissues and normal colorectal tissues. This enabled us to discern the spatial distribution and proportions of various CD8⁺T cell subsets: overall CD8⁺T cells, TCF1⁺CD8⁺T cells (indicative of naïve/memory cells), PD-1⁺CD8⁺T cells (representing activated cells), and KI67⁺CD8⁺T cells (denoting proliferating cells) ( Figures 2A, B ; Supplementary Figures 1 , 2A , 3A ).

Contrasting the cancerous tissues with their normal counterparts, we observed a reduction in the proportion of CD8⁺T cells within the cancer tissues. Among the CD8⁺T cell subsets, the frequency of TCF1⁺CD8⁺T cells exhibited a decline, while that of KI67⁺CD8⁺T cells experienced an uptick. The proportion of PD-1⁺CD8⁺T cells, however, remained relatively unchanged ( Supplementary Figure 3A ).

Our survival analysis revealed that the presence ratios of CD8⁺T cells and TCF1⁺CD8⁺T cells in tumor tissues were not significantly associated with patient outcomes ( Figures 2C, D ). Intriguingly, a heightened presence of PD-1⁺CD8⁺T cells and KI67⁺CD8⁺T cells was positively correlated with improved patient prognosis ( Figure 2E ; Supplementary Figure 2B ). Thus, it could be inferred that increased infiltration of activated and proliferating CD8⁺T cells augmented the prognosis for CRC patients.

Distinct from other CD8⁺T cell subsets, the spatial distribution of T_pex cells may be intrinsically tied to their biological role. To delve into this, we employed mIHC staining on CRC tissues and their adjacent normal tissues ( Supplementary Figure 1 ). This allowed us to scrutinize the spatial distribution and ratios of TCF1⁺PD-1⁺CD8⁺T cells (referred to as T_pex cells) and KI67⁺TCF1⁺PD-1⁺CD8⁺T cells (proliferating T_pex cells) ( Figure 3A ; Supplementary Figure 2C ).

Upon comparison, both T_pex cell and proliferating T_pex cell proportions in CRC tissues were broadly consistent with those in normal colorectal tissues. A notable exception was a marginal decrease in the proportion of proliferating T_pex cell within the overall T_pex cell population ( Supplementary Figures 3B, C ). Furthermore, survival analysis highlighted a positive correlation between the proportions of T_pex and proliferating T_pex cell, suggesting an association with improved patient outcomes ( Figures 3B–D ; Supplementary Figures 2D–F ).

To delve deeper into the distribution of CD8⁺T cell populations within tumor tissues, we compared the proportions of these populations in both epithelial and stromal regions. Our findings showed that the stromal region had a markedly higher overall proportion of CD8⁺T cells than the epithelial region. Yet, within the epithelial region itself, there was a significantly elevated proportion of activated, memory, and proliferating CD8⁺T cells when compared to their counterparts in the stromal region ( Figure 4A ).

As highlighted in earlier research (27), quiescent tumor cells often cluster together, and an increased proportion of these cells may inhibit immune cell infiltration. In contrast, a heightened presence of proliferating tumor cells can enhance immune cell infiltration. Notably, the balance between proliferating and quiescent tumor cells can be indicative of tumor immunogenicity.

To corroborate this, we delved into the spatial distribution and proportions of CK⁺, TCF1⁺CK⁺, and KI67⁺CK⁺ epithelial cells using mIHC staining ( Supplementary Figure 4A ). In alignment with prior findings, our analysis indicated that the ratios of CK⁺ and TCF1⁺CK⁺ epithelial cells didn’t considerably affect patient outcomes ( Supplementary Figures 4B, C ). However, patients with a greater proportion of KI67⁺CK⁺ epithelial cells presented with more favorable prognoses ( Supplementary Figure 4D ). Additionally, we explored the correlation between the proportions of TCF1⁺CK and KI67⁺CK epithelial cells and various CD8⁺T cell subsets. Interestingly, KI67⁺CK epithelial cells exhibited stronger correlations with most CD8⁺T cell subsets compared to TCF1⁺CK epithelial cells. The sole exception was the naïve/memory TCF1⁺CD8⁺T cells ( Figure 4B ).

Multivariate Cox’s regression analysis revealed several pivotal factors significantly impacting patient survival rates ( Figure 4C ). First, gender stood out as a critical determinant; male patients demonstrated a decreased survival rate. Additionally, the size of the tumor, specifically the maximum diameter, was inversely associated with survival: smaller tumors corresponded to better survival outcomes. Clinical stage was another influential factor, patients diagnosed at earlier tumor stages exhibited superior survival rates. Delving into cellular dynamics, varying levels of infiltrating immune cells emerge as significant. Specifically, a heightened proliferation of CD8⁺T cells and a greater proportion of activated CD8⁺T cells expressing TCF1 (percentage of T_pex cells within the PD-1⁺CD8⁺T cell cohort) positively correlated with improved survival, especially when juxtaposed against the proportion of TCF1⁺CK cells in the overall tumor cell population. In essence, our multivariate Cox’s regression analysis illuminated the key variables affecting patient survival rates, encompassing factors like gender, tumor dimensions, stage of tumor development, and cellular infiltration patterns. These insights bear significant value, offering clinicians a compass for tailoring treatment strategies and making well-informed patient management choices.

To explore the interactions between T_pex cells and other cell types within the TME, we turned to ST sequencing data from CRC patients (20). Analyzing proximity and co-localization patterns granted us a deeper understanding of potential cellular interactions and the specific micro-environmental niches T_pex cells occupies.

To systematically assess correlations across various cell types, we deconvoluted spatially-indexed datasets (21) and predicted the distribution of individual cells (24). Notably, our findings unveiled a positive correlation between T_pex cells and TAM1 (specifically, the TAM-C1QC subset), hinting at a potential interplay between these cell groups ( Figure 5A ). To further elucidate this relationship, we projected the derived cell clusters onto a transcriptome map ( Figure 5B ; Supplementary Figure 5 ). This graphical representation enriched our grasp on the spatial interplay between T_pex cells and TAM1 within the TME.

CRC has been categorized into four distinct CMS based on gene expression profiles and their molecular characteristics (16). To delve deeper into the relationship between T_pex cells and TAM1, we analyzed their correlation using bulk RNA-seq data (16, 28) and computed immune infiltration scores. This analysis unveiled a marked correlation between T_pex cells and TAM1 ( Figures 5C–F ). Intriguingly, the association between T_pex cells and TAM1 appeared more robust in the CMS1 and CMS4 subtypes than in the CMS2 and CMS3 subtypes. This pointed towards a particularly potent interaction between T_pex cells and TAM1 in immunogenic tumors. In essence, our research highlighted the co-localization and potential interplay between T_pex cells and TAM1 within the TME. Furthermore, when viewed through the lens of CMS subtypes, it’s evident that the synergy between T_pex cells and TAM1 was especially pronounced in immunogenic tumors.

Through an integrated analysis of both spatial and bulk transcriptome sequencing data, we uncovered strong evidence of the interaction between T_pex cells and TAM1 within the TME. This emphasized their potential collaborative role in this context. Yet, the precise mechanisms driving their interaction remain to be elucidated. To delve deeper into this dynamic, we turned to scRNA-seq data from CRC, classified according to the CMS (29). Our aim in analyzing this dataset was to further understand the molecular intricacies and functional interplay between T_pex cells and TAM1. By undertaking this comprehensive analysis, we hope to illuminate the driving forces behind their interaction, enriching our understanding of the complex cellular dynamics at play within the TME.

Analysis of cell-cell communication revealed a pivotal role for the co-stimulatory molecules CD80/86-CD28 and ICOS-ICOSL in the interaction between T_pex cells and TAM1 ( Figure 6A ). Notably, the immune checkpoint molecules PD-1-PD-L2 and TIGIT-Nectin2 also exhibited significant enrichment. Within the different CMS groups, CD86, PDCD1LG2 (PD-L2), and NECTIN2 were over-expressed in the CMS1 group ( Figure 6B ). Moreover, IFN-γ, originating from T_pex cells, modulated TAM1 through its receptor ( Figure 6A ). The transcription factor STAT1, which is linked to the IFN-γ pathway, also showed heightened expression in the CMS1 group ( Figure 6C ). TAM1’s immune-regulatory function was not limited to merely producing inflammatory cytokines. It secreted chemokines like CXCL9, CXCL10, and CXCL11, all downstream effectors of the IFN-γ pathway, to attract T cells. Additionally, IFN-γ can induce the up-regulation of CXCR3 on CD8⁺T cells. Considering prior findings that highlight an increased expression of CXCR3 on T_pex cells ( Figure 1D ), we postulated a potential interaction between TAM1 and T_pex cells via the CXCL9/10/11-CXCR3 axis. Our analyses also underscored a marked enrichment of the CXCL9/10-CXCR3 interaction between these cell populations ( Figure 6D ). Notably, the expression of CXCL9 within the CMS1 group was significantly higher compared to other subtypes ( Figure 6E ). This aligned with the previously observed amplification of IFN-γ signaling. Importantly, heightened expression of these chemokines and receptors between TAM1 and T_pex cells was correlated with better prognoses in CRC patients ( Figure 6F ). In essence, the dynamic relationship between TAM1 and T_pex cells encompassed a blend of co-stimulatory molecules, immune checkpoint proteins, and downstream chemokines and their receptors, all integral to the IFN-γ signaling cascade.

CD8⁺T cells are cornerstone players in tumor immunotherapy. To delve deeper into the functionality of T_pex cells and its differentiation from other CD8⁺T cell subsets in the context of immunotherapy, we dissected characteristic gene sets derived from scRNA-seq data of CD8⁺T cells obtained from CRC patients. By integrating insights from both scRNA-seq (30, 31) and bulk RNA-seq (17) datasets related to immunotherapy, we computed scores for each CD8⁺T cell subset. Our initial focus was on analyzing CD8⁺T cell subsets from scRNA-seq datasets of CRC and melanoma patients. This involved evaluating specific scores for each subset ( Figure 7A ; Supplementary Figure 6 ). Notably, our statistical scrutiny unveiled significantly elevated T_pex cell scores in CRC and melanoma patients who showed positive responses to immunotherapy compared to those who did not. In stark contrast, the T_ex scores in responders from both patient groups were markedly reduced relative to non-responders. Furthermore, we navigated through bulk RNA-seq data sourced from rectal cancer patients undergoing neoadjuvant therapy ( Figure 7B ). Corroborating our earlier findings, the data highlighted that patients demonstrating a complete response to immunotherapy had noticeably heightened T_pex cell scores. In contrast, the T_ex cell score was markedly depressed in those showcasing an incomplete response. Collectively, these revelations underscored the potential of T_pex cell proportions as pivotal indicators for prognosticating the efficacy of immunotherapy targeting CD8⁺T cells in CRC.