Experimental approach and study design

We obtained ten surgical resections of CRCs or GCs from eight different patients (Table 1, Additional file 1: Table S1). They comprised four male and four female patients with an average age of 62 years. The tissue samples included seven resections of primary CRC from seven individual patients. From one patient with GC, we obtained three independent resections—one from the primary tumor and two from metastases to the peritoneum, which is an organ that lines the abdominal cavity.

The samples underwent rapid processing following surgical resection (Methods). From each tumor, we split the tissue, using one portion to generate single-cell suspensions and scRNA-seq and scTCR libraries for sequencing (Fig. 1A, “Methods”). The tissues that were processed immediately into single-cell libraries provide a baseline of the cellular composition of the tumor. We refer to this baseline as time-point zero (“T0”). The other portion was used for TSC culturing. These results were used for determining cellular changes that may occur during the TSC culturing.

For experimental testing of each cancer’s TME, we generated ex vivo tumor slice cultures (TSCs) from the resections. There were four different conditions. The TSCs were treated with either (i) isotype control antibody (“ctrl”), (ii) known T cell activator PMA/Ionomycin (“PMAIono”) as a positive control, (iii) GITR agonist antibody (“GITR”), or (iv) TIGIT antagonist antibody (“TIGIT”). After 24 h of treatment, the cells were harvested and then processed for scRNA-seq and scTCR-seq. The number of experimental conditions tested per sample depended on the available size of each resection. All ten samples had adequate tissue for a baseline T0 and control sample for scRNA-seq. We conducted scRNA-seq on PMA/Ionomycin treatment from eight samples, GITR agonist treatment from nine samples, and TIGIT antagonist treatment from five samples (Table 1). Quality control measures including filtering cells for mitochondrial genes indicative of cell death [39] and doublet identification [24]. Following filtering, our final analysis included a total of 236,483 single cells with an average of 5630 cells per sample (Additional file 1: Table S2).

Baseline immune cell characteristics of the TME from primary gastrointestinal tumors

We determined the baseline cellular composition of the T0 samples. Batch effects were reduced using the Harmony algorithm [25]. Specific cell type clusters were composed of different samples, indicating the elimination of batch effects (Fig. 1B). Using the scRNA-seq data, we made cell type assignments based on canonical marker genes (Methods). Overall, we identified tumor epithelium, myeloid cells, stromal cells, and lymphocytes. We sub-clustered each major lineage to characterize cellular features in greater detail. For these results, we denote the cell type and functional state by listing prominent examples among the associated gene expression markers.

Tumor epithelial cells expressed well-characterized markers of CRC or GC [18, 40] (KRT7, KRT17, ELF3, CEACAM6, FABP1, FABP5, SPINK1, REG4, TFF3) (Additional file 2: Fig. S1A, S1B). Stromal cells included fibroblast subsets with expression of extracellular matrix-related marker genes including MGP, DCN, and Collagen family genes (Fig. 1C, Additional file 2: Fig. S1C). We also detected smooth muscle cells (ACTA2, TAGLN) and pericytes (RGS5, PDGFRB, NOTCH3). Peritoneal metastasis samples (GC_1_2, GC_1_3) also contained mesothelial cells (SLPI, UPK3B, KRT8, KRT18, KRT19) [40, 41]. Endothelial cells expressed known arterial, venous, capillary (PLVAP, VWF, CD320, PECAM1, KDR, ENG, ACKR1, SELE, ICAM2, SRP14, SRGN), or lymphatic markers (CCL21, LYVE1, PROX1, PDPN) [18, 40, 42].

Myeloid lineage cells included macrophages, dendritic cells (DCs), and mast cells (Fig. 1D, Additional file 2: Fig. S1D). Among macrophages we detected previously characterized subsets [21, 43, 44] including infiltrating monocytes (S100A8, S100A9, FCN1, VCAN), proinflammatory (CXCL8, IL1B, IL6, IL8), anti-inflammatory LYVE1+ (LYVE1, FOLR2, PLTP), SPP1+ (SPP1, APOE, TREM2, CTSB, MMP9), and C1QC+ (C1QA, C1QB, C1QC, APOE) macrophages. DCs included conventional DCs (cDC) subsets (CLEC9A, FLT3, IDO1, CD1C, FCER1A, HLA-DQA1, HLA-DQB1, LAMP3, CCR7, CCL22, CCL19) as well as plasmacytoid DCs (pDC) (GZMB, SOX4, JCHAIN, IRF7) [18, 43, 45]. Mast cells highly expressed known marker genes (TPSAB1, TPSB2, KIT, GATA2, CPA3, MS4A2) [43].

Among lymphocytes, we detected B cells expressing B cell markers MS4A1, CD79A, CD79B, CD19, CD83, and CD37 (Additional file 2: Fig. S1E, F) [21]. Plasma cells lacked mature B cell marker genes and expressed known marker genes including SDC1, TNFRSF17, and immunoglobulin genes including JCHAIN, IGKC, IGHG1, IGHM, and IGLC2. We also detected proliferating B cells expressing MKI67, STMN1, and TUBA1B.

We characterized the T and NK functional cell states using a method called cell reference mapping. This method used an established reference from a pan-cancer tumor immune cell atlas [21] and the SingleR algorithm [29]. Each cell’s gene expression is matched to a given reference cell type. This approach provides an unbiased identification of cell subtypes without applying cell clustering methods. We evaluated these cell states based on the expression of lineage markers, transcription factors, surface receptors, cytokine effectors, and other genes that have been extensively described in recent scRNA-seq studies [2, 20, 21].

We detected CD4 naïve cells (CCR7, SELL, LEF1, TCF7, IL7R) (Fig. 1E). We identified regulatory T (Treg) cells with high expression of FOXP3, BATF, IL2RA, co-stimulatory molecules TNFRSF4 and TNFRSF9, and immune checkpoint CTLA4, resembling the profile of intratumoral Tregs identified by previous scRNA-seq studies [21]. Treg cells are immunosuppressive and limit anti-tumor activity through specific effects on cytotoxic CD8 T cells, dendritic cells, and macrophages [46]. As corroborated by other studies [35], we also observed a proliferative subset of Treg cells. These proliferative subsets may reflect a TME response to local tumor antigens [2].

We also detected CXCL13 expressing CD4 T cells with low expression of T helper cell cytokines such as IFNG, GZMA, GZMB, CCL3, and CCL5 [20] (Fig. 1E). These cells expressed several genes associated with T follicular helper (TFh) differentiation including transcription factors NR3C1, TOX2, TOX, TSHZ2, RBPJ, and BHLHE40 but did not express CXCR5 or BCL6 [20]. They also expressed genes associated with CD8 exhaustion including NMB, CD200, and PDCD1. We also detected a proliferating subset of these cells with higher expression of T helper cytokines, possibly reflecting a response to tumor antigens [2]. These cells resemble previous scRNA-seq analysis findings which variously labelled them as CXCL13+ T helper-like cells [19], CD4_CXCL13 cells [47, 48], CD4- CXCL13 with TFh-like features [34], PD-1+CXCR5−CD4+Th-CXCL13 [49], dysfunctional TFh [35], TFh-related cells, or TFh/Th1 cells [20]. We refer to these cells as TFh-like cells. TFh-like cells have been linked to anti-tumor immunity by promoting CD8 and B cell activity [2, 50].

CD8 naïve cells (CCR7, SELL, LEF1, TCF7) lacked effector cytokine or checkpoint expression (Fig. 1F). We also detected NK cells with high expression of various cytotoxic effector genes including NGK7, GNLY, PRF1, CCL4, GZMA, GMZB, and GZMH together with Killer cell lectin-like receptors [21].

Among the CD8 T cells, we identified effector cytotoxic CD8 characterized by high expression of effector cytokine GZMK and low expression of immune checkpoints (Fig. 1F). GZMK expression in CD8 effector cells has been associated with early dysfunction [2, 20, 51]. We observed CXCL13 expressing dysfunctional CD8 T cells with increased expression of inhibitory (LAG3, PDCD1, HAVCR2, CTLA4) and co-stimulatory (TNFRSF9) receptors [2]. These cells also expressed genes (ENTPD1, LAYN) and transcription factors (RBPJ, TOX, PRDM1) linked with exhaustion. Additionally, they continued to express cytotoxic effectors (including GZMA, GZMB, GZMH) reflective of their anti-tumor potential [20]. Dysfunctional CD8 T cells had a subset of proliferating cells (noted by expression of the marker genes MKI67, STMN1, TUBA1B) with intermediate CXCL13 and checkpoint gene expression. Proliferating dysfunctional CD8 T cells have been linked to early dysfunction in a clonal tumor-reactive population [35]. Rare subpopulations of dysfunctional CD8 T cells function as precursor or progenitor cells that give rise to terminally exhausted cells [2]. We examined a set of genes linked to precursors of exhausted T (TPEX) cells (TCF7, CCR7, SELL, IL7R, TNFRSF4, IL6R, IGFL2) [43, 52] (Additional file 2: Fig. S1G). Rare cells in both dysfunctional and dysfunctional proliferating cells expressed TPEX genes. Hence, our analysis identified transitional states of CD8 T cell dysfunction in the TME including progenitor, early, and late dysfunction that have been previously reported [2].

These cell types were identified across all patients in varying proportions (Additional file 1: Table S3). In summary, across all tissue samples, the TME in the baseline T0 resections contained diverse functional T cell states with anti-tumor (cytotoxic CD8, dysfunctional CD8, TFh-like) and immunosuppressive (Treg) properties.

Baseline T cell receptor clonality in the primary tumor TME

To assess clonality of the T cells in the cancer’s TME at the baseline state, we performed scTCR-seq on the baseline T0 samples. We identified TCR chains from an average of 57% of the T cells with matching single-cell gene expression (range 31–78%). Next, we determined whether there was evidence of TCR clonotypes being highly represented within a given sample [19]. This overrepresentation is termed as being “an expansion” for a given T cell clonotype. Moreover, one can assign specific clonotypes to different transcriptional cell states (i.e., Tregs, Tfh) using matched cell barcodes.

To conduct this analysis, the frequency of individual clonotypes was calculated using the Shannon entropy score—this metric quantifies T cell clonotype expansion with a value range of 0 to 1, with 1 indicating high clonality. Cytotoxic and dysfunctional CD8 T cells showed high expansion index of TCR clones (Fig. 2A). High clonality and expansion may be an indicator of tumor antigen-driven expansion in the TME [19]. Next, we examined the frequency distribution of CD8 T cell clonotypes across samples (Fig. 2B). Single-cell clonotypes represented the majority of TCRs, indicating a lack of expansion among these clonotypes. Across the samples, between 22 and 89% of the total cells were represented by expanded clonotypes.

Fig. 2

A TCR expansion index for respective cell types. B Frequencies of clonotypes in CD8 T cells from respective patients together with absolute number of cells and clonotypes examined. C Overlap between TCR clonotypes in cytotoxic and dysfunctional CD8 T cells from respective patients. D, E Frequencies of clonotypes in D TFh-like and E Treg cells from respective patients together with absolute number of cells and clonotypes examined

We investigated the overlap between clonotypes found in cytotoxic and dysfunctional CD8 T cells (Fig. 2C). We detected an average overlap of 36% (range 23 to 52%) between the two cell types across all samples. This analysis excluded tumor CRC_3 with only one clonotype detected in dysfunctional CD8 T cells. Hence, a subset of GZMK+ effector cells were linked to the dysfunctional phenotype in agreement with previous studies [2, 20, 35].

CD4 TFh-like and Treg cells had a low degree of expansion in a subset of samples (Fig. 2A). Examining the clonotype frequency distribution confirmed that three out of ten tumors did not contain expanded clones (Fig. 2D, E). In the remaining seven tumors, expanded clonotypes comprised between 8 and 80% of all cells in TFh-like and 6–55% in Tregs. Although the total number of cells analyzed affected these expansion metrics, our findings resemble previous studies, which detected highest expansion in CD8 T compared to TFh-like and Treg cells [20, 35].

Overall, this analysis identified TCR sequences of expanded clones in infiltrating T cells that may be potentially tumor-reactive in each sample at baseline.

GITR and TIGIT gene and protein expression in the baseline TME

We evaluated the gene expression of the immunotherapy targets TNFRSF18 (encoding protein GITR) and TIGIT among gastrointestinal cancers at the baseline state. In both CRCs and GC tumors, the dysfunctional CD8, TFh-like, and Treg cells had the highest levels TNFRSF18 expression (Fig. 3A, B). The complementary ligand, TNFSF18, encoding the protein GITRL, was expressed by fibroblasts, DCs, and macrophages in CRC. TIGIT expression was highest in cytotoxic CD8, dysfunctional CD8, TFh-like, and Treg cells. The genes PVR and NECTIN2, which encode for TIGIT ligands, were expressed by tumor epithelial, endothelial, fibroblasts, macrophages, and DCs in the TME. These expression patterns are along the lines of other reports [53, 54]. Overall, this result indicated that among all samples, the TME cells expressed genes required for GITR and TIGIT receptor-ligand signaling.

Fig. 3

A, B Scaled average expression of respective genes in various cell types from A all CRC T0 resections and B all GC T0 resections. C, D Immunofluorescence staining for respective proteins or their merged image in an example region of interest from sample CRC-2. Scale bar = 50 μm. E, F Scaled expression of respective genes in various cell types from E CRCs in the publicly available tumor immune atlas dataset and F our previously published GC dataset

We also measured the protein expression of GITR and TIGIT in these baseline tumor tissues (T0 resections) using multiplexed immunofluorescence (mIF) staining. We used two independent antibody panels containing CD8, FOXP3, and TIGIT or CD8, FOXP3, and GITR respectively (Fig. 3C, D). We performed image analysis using a multiplex classifier for detecting single-stain or double-stain positive cells as described previously [18]. From all samples, an average of 37.4% of total CD8 positive cells expressed TIGIT. An average of 53.68% of total FOXP3 positive cells were TIGIT positive (Additional file 2: Fig. S2A, B). Similarly, 42.46% of CD8 cells expressed GITR (Additional file 2: Fig. S2C, D). Among the FOXP3 cells, 74.5% expressed GITR. These results confirmed that among our tumor samples, CD8 T cells and Tregs expressed the TIGIT and GITR protein.

Overall, these results indicated the TME expression of the target receptors and their ligands among the tumors that were used for these experiments. Targeting these receptors has the potential to modify the function of anti-tumor T cell subsets such as cytotoxic CD8, dysfunctional CD8, and TFh-like cells, as well as immunosuppressive Tregs.

Expression of TIGIT and GITR in colorectal and other cancer types

To determine the expression of these two targets among an expanded, independent set of colorectal, gastric, and other tumor types, we analyzed gene expression among a data set of 13 different cancer types [21]. Importantly, this dataset included 25 independent CRCs (Fig. 3E). For gastric cancer, we evaluated our previously published dataset of seven GC samples [18] (Fig. 3F). High TNFRSF18 and TIGIT expression was detected in dysfunctional and cytotoxic CD8 T, TFh-like, Treg, and proliferating cells confirming results from CRC and GC T0 samples.

Other cancer types included breast carcinoma (BC), basal cell (BCC) and squamous cell carcinoma (SCC), endometrial adenocarcinoma (EA), renal cell carcinoma (RCC), intrahepatic cholangiocarcinoma (ICC), hepatocellular carcinoma (HCC), pancreatic ductal adenocarcinoma (PDAC), ovarian cancer (OC), non-small-cell lung cancer (NSCLC), and cutaneous (CM) and uveal melanoma (UM) (Additional File 2: Fig. S2E, F). Across all cancer types, target expression followed similar patterns as GC and CRC. In EA and OC, high TNFRSF18 expression was noted in DCs and NK cells respectively. These results confirmed the expression of these targets among CRC and GC tumors as well as a wide variety of solid tumor types.

Primary tissue slice cultures maintain the native TME composition of gastrointestinal cancers

As noted previously, tissue slice cultures have been demonstrated to maintain a high degree of tissue viability, cellular diversity, and cellular transcriptional profiles [6, 10]. We integrated data from all ex vivo ctrl and treatment experiments and performed cell type identification using marker based and SingleR assignments (Methods).

First, we confirmed that the TSCs cultured for 24 h maintained the cellular characteristics similar to the baseline state of the TME (i.e., T0 cell conditions at time of resection). As noted previously, other groups have used short culture periods to maintain the cellular diversity found in the native TME [4,5,6,7]. We evaluated the TSC cellularity using hematoxylin and eosin (H&E) staining of the cultures. This result showed that cell morphology remained intact with little evidence of necrosis or other signs of overt cell death (Additional file 2: Fig. S2G).

Next, we evaluated the single-cell gene expression across the two conditions for all samples which included: (i) the baseline T0; (ii) TSC following a 24 h incubation with isotype control antibody; (Fig. 4A). We corrected the data for experimental batch but not for the experimental condition, using the Harmony algorithm [25]. Cells belonging to the baseline T0 tissue and control clustered together in the Uniform manifold approximation and projection (UMAP). This result indicated that the cells had similar gene expression profiles. We calculated the Adjusted Rand Index (ARI) to determine the variation in gene expression between the baseline and control culture. Cluster labels compared to the experimental condition had a low ARI value of 0.009 which indicated that clustering was due to the cells having similar gene expression characteristics and not driven by the T0 or TSC experimental condition. Despite the low ARI, we observed shifts in the UMAP embeddings between the two conditions. This indicated that transcriptional profiles in TSCs resemble T0 but are not completely identical.

Fig. 4

A, B UMAP representation of dimensionally reduced data from T0 and 24 h ctrl TSCs following batch-corrected graph-based clustering of all datasets colored by A experimental condition and B cell type. C Quantile-quantile plot comparing the proportion distributions of respective cell lineages across all T0 and ctrl TSCs. D Scatter plot indicating average log expression of marker genes for T0 cell lineages in T0 and ctrl TSC in respective cell lineage, annotated with the number of marker genes examined. Pearson’s co-efficient was calculated using non-log transformed values

We compared cell types in T0 and TSC samples. The TSC samples contained all cell types compared to the matched baseline T0 samples. These cell types included tumor epithelial, macrophages and dendritic cells, NK, B or plasma lymphocytes, mast cells, fibroblasts, and endothelial cells, as well as T cell subsets including naïve, cytotoxic CD8, dysfunctional CD8, Treg, and TFh-like cells (Fig. 4B). The relative proportion of all cell lineages was also maintained in the TSCs compared to the baseline tumor tissue (Fig. 4C).

We identified differentially expressed genes for each cell lineage in the baseline T0 samples (Seurat Wilcoxon test, log2 fold change ≥ 0.4, adjusted p ≤ 0.05). We compared the average gene expression of these genes in each respective cell lineage in T0 to TSCs. Expression was highly correlated across all cell types (Fig. 4D) (Pearson correlation ≥0.72, p ≤ 7.1E-22). Hence, the TSCs maintained the cellular heterogeneity and closely resembled the transcriptional cell states that were present in the original tumor.

General stimulation of T cells and other cell types in the TSC TME

To demonstrate that the TSC cells were functionally responsive, we used an activation control stimulus with phorbol ester 12-myrisate 13-acetate and the calcium ionophore ionomycin (PMA/Ionomycin). In combination, these compounds stimulate downstream pathways associated with T cell activation and have been extensively studied [55]. We evaluated specific subsets of cells from samples treated with ctrl and each respective perturbation. To account for interpatient variability in differential expression (DE) analysis, we utilized model-based analysis of single-cell transcriptomics (MAST) [30] incorporating sample as a random effect in the model [