Introduction

Colon cancer is one of the most common malignancies, with more than 1 million new cases and 500 thousand deaths reported globally in 2020 (Sung et al., 2021). Approximately 19–26% of patients with colon cancer present with synchronous metastatic diseases at the first diagnosis (Siegel et al., 2020), among which 14.5– 17.5% develop liver metastasis (Siegel, Miller, 2020, Wang et al., 2023). Although progress in the treatment of metastatic disease, including improved diagnosis and treatment strategies for liver metastases (Chua et al., 2011, Zampino et al., 2016, Ruers et al., 2017), increased cancer-directed surgery (Wancata et al., 2016), and the development of targeted therapies (Piawah et al., 2019), has greatly improved the survival of these patients in recent decades, the 2-year relative survival rate for patients diagnosed with distant-stage disease was 36% (Siegel, Miller, 2020), significantly lower than those without metastasis. Thus, colon cancer-derived liver metastasis (CCLM) is a clinical challenge that requires urgent development of novel treatment methods.

The tumor microenvironment (TME) is a dynamic environment that governs tumor behavior and is required for tumor cell survival, growth, proliferation, and metastasis (Yao et al., 2023). Unlike previous understanding, it has recently been proposed that certain immune cells present in the TME (such as myeloid-derived suppressor cells (Liu et al., 2023), macrophages (Ngambenjawong et al., 2017), neutrophils (Li et al., 2019), and CD8⁺ T cells (Tiberti et al., 2022)) favor tumor progression. Thus, the contradictory findings of the pro- and anti-tumor effects of tumor-infiltrating immune cells require a better understanding of the immune features and profile heterogeneity of CCLM for the future development of immune- modulatory strategies to stratify and target immune cells.

Natural killer (NK) cells, an important component of tumor-infiltrating immune cells, are cytotoxic lymphocytes that play a key role in recognizing and eliminating malignant cells and are therefore considered the early responders against tumors (Maskalenko et al., 2022). Among various immune cells, tumor-infiltrating NK cells are associated with prognosis in various cancers, such as colorectal cancer (CRC) (Zhong et al., 2022), lung cancer (Villegas et al., 2002), and gastric cancer (Ishigami et al., 2000). Recently, based on their anti-tumor potential, NK cells were exploited for treating malignancies and have proven to be highly promising for the treatment of certain hematologic malignancies (Lamb et al., 2021, Liu et al., 2020) but have limited efficacy for treating solid tumors (Zhang et al., 2024, Hu et al., 2019, Oh et al., 2019). Additionally, intratumor NK cells differ phenotypically or functionally from peripheral NK cells in several malignancies (Bruno et al., 2013, Carrega et al., 2008); a specific subset of tumor-infiltrating NK cells (CD11b⁻ and CD27⁻) was associated with tumor progression in lung and hepatocellular cancers (Zhang et al., 2017, Jin et al., 2013). Notably, the tumor-promoting transformation of NK cells may be educated by cancer cells (Chan et al., 2020, Huergo-Zapico et al., 2018). However, their function remains elusive. In particular, the presence of specific NK cell subsets within colon cancer and their potential relationship to CCLM progression have not been fully characterized to date.

In this study, we combined previously published single-cell RNA sequencing, spatial transcriptomics (ST), and bulk RNA-sequencing data from public datasets, aiming to comprehensively chart the cellular landscape of TME in primary colon cancer and matched liver metastasis. Results of bioinformatics analysis were further validated through an in vitro co-culture experiment (Fig. 1A).

Figure 1.

Results

Integrated scRNA-seq and ST Precisely Quantified Immune Cell Diversity in CCLM

We used a previously published single-cell dataset containing 89 samples from paired samples of colon cancer, adjacent colon, liver metastasis, adjacent liver, lymph nodes along colons, and peripheral blood mononuclear cells from 20 patients to define the heterogeneous immune microenvironment landscapes of CCLM. Subsequently, 178,630 CD45+ cells were integrated.

To further define the main cell type, clustering analysis and SingleR were performed, and monocytes, neutrophil cells, native B cells, plasma cells, T cells, and NK cells were identified from the CCLM samples (Fig. 1B-C; Fig. S1A and B). Notably, T and NK cells were the major components of all immune cells and distinct across different tissue types, with NK cells significantly decreasing in tumor tissues (both primary colon cancer and liver metastasis cancer) compared to corresponding adjacent normal tissues, whereas T cells were the opposite (Fig. 1D and E; Fig. S1C).

To comprehensively analyze the spatial distribution profile of colon primary tumors and CCLM tumors, we collected ST data, which came from eight samples including paired primary colon cancer and liver metastasis. Unsupervised clustering analysis was first performed to cluster similar ST spots, and the annotation of the clusters was further determined according to cell marker genes (Fig. S2A-C). Further, after combining the gene expression features of each sample (Fig. S3) and H&E staining, six morphological regions, including tumor, fibroblast, smooth muscle, B cells, hepatocytes, normal epithelium, and tumor and paratumor areas, were identified. (Fig. 2A-C).

Figure 2.

To integrate the scRNA-seq and ST data, we used the AddModuleScore function in Seurat to quantify the main immune subpopulations. Consistent with scRNA-seq, T cells were remarkably enriched in the tumor region of both primary colon cancer and liver metastasis cancer, whereas NK cells were enriched in the non-tumor area (Fig. 2D).

Biological Relationship Between NK Cells and Metastasis of Colon Cancer Revealed by Bulk RNA Transcriptomics

TME cell composition differences between metastasis and non-metastasis colon cancer were evaluated using multiple tools that could robustly quantify the abundance of cell populations based on transcriptomic datasets, including xCell, EPIC, quanTIseq, and MCPcounter and bulk-seq datasets (TCGA COAD cohort). The percentages of NK cells markedly decreased in metastasis colon cancer using MCPcounter, EPIC, and xCell (Fig. 3A-D).

Figure 3.

To further understand the potential triggers that induce colon cancer metastasis, the DEGs were calculated between the metastasis and non-metastasis colon cancer groups in the TCGA COAD cohort. A total of 1,378 DEGs were determined using the limma package with cutoffs |log fold change|>1.5 and p <0.05, including 817 upregulated and 561 downregulated genes between metastasis and non-metastasis colon cancer (Fig. 3E). Subsequently, GO and KEGG enrichment analyses and GSEA were performed to investigate the correlative functions and pathways. As shown in Fig. 3F, cell killing, NK cell-mediated immunity, NK cell-mediated cytotoxicity, NK cell activation, and MHC-related GO terms were significantly enriched. The KEGG pathways of NK cell-mediated cytotoxicity were also enriched (Fig. 3F), with most of the genes in this pathway markedly downregulated in the metastasis group compared to the non-metastasis group. In contrast, GSEA revealed that NK cell-mediated cytotoxicity was enriched in the metastasis group (Fig. 3G and H).

The Landscape of NK Cells in CCLM Progression

The differences in infiltrated NK cells between tumor-adjacent tissues and tumor tissues suggest that the dynamic remodeling of NK cells may play an important role in CCLM progression. To further understand the landscape of NK cells in the CCLM microenvironments, NK cells were selected, and unsupervised clustering analysis was performed. Eight NK cell subtypes were thus identified (Fig. 4A). Combining the highly variable features of each NK cluster (Fig. 4B) and CIBERSORT-reported canonical NK cell markers (Fig. 4C and D), NK cells were classified into three cell types, including activated NK cells identified by the expression of KIR2DL4, GPR183, GRP171, CD69, and IFNG, resting NK cells marked by GZMK, TTC38, CD160, and PLEKNF1, and the other NK cells of which no characteristic gene was identified. Although all three NK cell subtypes were presented in primary colon tumors, adjacent normal colon tissues, liver metastasis tumors, and lymph nodes (Fig. 4E), the infiltration grade for each of these NK cell subsets was significantly different. At the individual sample level, considerable variability was observed in the NK cell subset composition. Where serial samples were available from individual patients, the expression of resting NK cells increased, whereas that of activated NK cells decreased during disease progression (Fig. 4F and G). Additionally, using the paired primary colon tumor, adjacent normal colon tissues, and liver metastasis tumor samples, we observed that the gradient of the resting NK cells increased, whereas the gradient of the activated NK cells decreased (Fig. 4H).

Figure 4:

To evaluate the clinical relationship between NK cell differences and CCLM, CIBERSORT was performed using the TCGA COAD cohort. The percentages of resting NK cells significantly increased (Fig. 5A), and the activated NK cells decreased in the metastasis group (Fig. 5B), which is consistent with the scRNA-seq results. Also, we observed that the resting NK cells were significantly increased, whereas the activated NK cells decreased in higher T, N, and TNM stages (Fig. S4A). In terms of survival, neutrophils and resting NK cells were associated with a worse prognosis, whereas activated NK cells and M1 macrophages were associated with better outcomes (Fig. S4B). Further, survival analysis showed that colon cancer patients with lower and higher degrees of infiltration of activated NK cells and resting NK cells, respectively, had a significantly worse prognosis in the TCGA COAD and GSE29623 cohorts (Fig. 5C and D). Interestingly, resting NK cells were remarkably enriched in the tumor region, whereas activated NK cells were enriched in the non- tumor area of both primary colon cancer and liver metastasis (Fig. 5E).

Figure. 5

Characterization and Developmental Course of Differential Subsets of NK Cells in CCLM

To further study the heterogeneity of NK cells subgroup, the proportions of all eight clusters of NK cells were calculated (Fig. S5). KIR2DL4+ GPR171+ activated NK cells were lessened in proportion with the progression of the disease from normal colon to colon cancer and liver metastasis cancer, whereas GZMK+ resting NK cells increased (Fig. 6A and B), suggesting that these two NK subsets were associated with metastasis. ST analysis also confirmed that there was more infiltration of the KIR2DL4+ GPR171+ activated NK cell subset in the non-tumor area and GZMK+ resting NK cells in the tumor region, both primary colon cancer and liver metastasis (Fig. 6C).

Figure 6.

Pseudotime cell trajectory analysis of the two NK cell clusters was constructed to investigate the evolutionary dynamics of metastasis-associated NK cells. NK cells from the normal colon group were located at the top right corner of the trajectory curve, suggesting the clear starting point of this evolving trajectory curve. After confirmation of this starting point, developmental routes were determined to begin with the KIR2DL4+ GPR171+ activated NK cells and then develop into GZMK+ resting NK cells (Fig. 6D).

To explore the transitional relationships of marker genes in NK cells in distinct clusters, we examined the expression level of these markers in different clusters through pseudotime analysis. Marker genes of activated NK cells (KIR2DL4 and GPR171) exhibited a gradient descent expression pattern, whereas the expression of resting NK cell marker genes (GZMK) exhibited a gradient rise pattern (Fig. 6E). BEAM analysis was performed to identify the cell cluster marker genes that change as NK cell subtypes move from the activated stage to the resting stage. The branched heatmap showed the expression dynamics of the top 200 significant genes in different cell fate branches. Corresponding GO enrichment analyses of these significant genes further demonstrated that a mass of leukocyte differentiation and immune system process-related GO terms were significantly enriched (Fig. 6F).

Interestingly, patients with colon cancer with higher infiltration of GZMK+ resting NK cells (Fig. 6G) and lower infiltration of KIR2DL4+ GPR171+ activated NK cells (Fig. 6H) exhibited shorter survival in the TCGA COAD cohort.

Tumor Cells-Educated NK Cells Shift Toward Tumor-Promoting Status Depends on Cell-to-Cell Interaction

Given the characteristic spatial distribution of NK cell subpopulations, we hypothesized that the tumor region-enriched NK cells might be educated by tumor cells and thus functionally distinct. To verify the possible effect of tumor cells on the phenotypic switch of NK cells, we set up mixed cell co-culture experiments using colon cell line HCT-116 and NK cell line NK-92 in the ratio 1:1. As controls, NK cells were cultured alone, either in the tumor supernatant or fresh medium (Fig. S6A). After 24 h of co-culturing, NK cells were analyzed using FACS to determine the expression of the NK cells markers above identified as well as functional receptors NKG2A and NKG2D. Upon exposure to tumor cells, NK cells significantly decreased the expression of activated receptor KIR2DL4 (Fig. 7A) and functional receptors, both NKG2A and NKG2D (Fig. 7 B). However, no such effect was observed after co- culturing with tumor supernatant and fresh medium.

Figure. 7

To explore the NK cell-mediated inductive effect on colon cells, we performed co-culture experiments in transwells, and the co-cultured supernatant was collected to incubate fresh HCT-116 (Fig. S6B).

Interestingly, HCT-116 cells in the CNS group underwent a significant increase in proliferation (Fig. 7C-E) and migration (Fig. 7F-H) compared to the CS group, where cancer cells were cultured directly in the supernatant after 24 h of cancer cell culturing. Furthermore, CCK-8 (Fig. 7I) and colony formation assays (Fig. 7J-K) showed that colon cells underwent a significant increase in proliferation when cultured in the supernatant from the co-culture system in which NK and colon cells were in contact (CNS group), but not when co-cultures were performed in the cell supernatant (SNS group) and fresh medium (MNS group). Similarly, the migration and invasion of HCT-116 were significantly increased when cultured in the supernatant from the co-culture system in which NK and colon cells were in contact (CNS group) (Fig. 7L-N).

Further, another colon cancer cell line, DLD-1, was chosen to evaluate the effect of the supernatant in the different co-culture groups. DLD-1 cells in the CNS group did not undergo a prominent increase in proliferation, as HCT-116 and CCK-8 assays showed a minimal proliferation of colon cells in the MNS group (Fig. S7A-C). The same results were obtained for their migration and invasion (Fig. S7D-F). However, no such cell functional changes were induced in the CNS group (Fig. S7G-L).

Resting NK-induced Colon Cancer Malignant Phenotype Promotion Depends on SCF Release

To study the potential mechanism by which the supernatant from the NK and colon cells in direct contact with the co-culture system promoted colon cell migration, Luminex liquid suspension chip detection was used to compare the differential expression of 48 common chemotactic and inflammatory cytokines in the CNS, SNS, and MNS groups. We observed that SCF was upregulated in all three repetitions in the CNS group (Fig. 8A-B). Further, ELISA confirmed the upregulated level of SCF in the CNS group (Fig. 8C).

Figure. 8

Figure. 9.

We further assessed the inhibitory effects of imatinib mesylate on the enhanced proliferation and invasion. CCK-8 (Fig. 8D-E) and colony formation assays (Fig. 8F- I) showed that imatinib mesylate (2 µM) significantly inhibited supernatant-enhanced proliferation of HCT-116. Similarly, the supernatant-enhance migration (Fig. 8J-K) and invasion (Fig. 8L-M) were also inhibited by imatinib.

Discussion

CCLM is a multistep process and is mostly fatal (Manfredi et al., 2006). During disease progression, functional interactions between tumor cells and the surrounding TME are critical for influencing tumor growth, promoting angiogenesis, and finally resulting in metastasis (Quail et al., 2013). Since the genomic divergence between primary and metastatic CRC cells is relatively low (Wei et al., 2017, Yaeger et al., 2018) and accumulating evidence has highlighted the key role of tumor-infiltrating immune cells in dictating the fate of cancer cells (Kumar et al., 2018, Hanahan et al., 2012), we presented a comprehensive cellular and spatial immune landscape of the primary and liver metastatic tumors of CRC by using previously published scRNA- seq and ST data. Consistent with previous studies (Zhou et al., 2021, Wu et al., 2022), we observed that the immune microenvironment of primary and liver metastasis lesions of colon cancer had undergone extensive remodeling with a significantly decreased proportion of NK cells and a strong enrichment of T cells. However, the contrast between the decreased proportion of NK cells in tumor tissue and the enrichment of NK cells in tumor regions may indicate the specific effects of NK cells in CCLM.

Traditionally, NK cells within the TME were known as potential suppressors for tumor growth by directly killing cells and secreting proinflammatory cytokines (Huntington et al., 2007, Freud et al., 2006). However, the function of NK cells highly depends on their maturation status and localization (Cooper et al., 2001). The peripheral blood CD56^dimCD16^high NK cell population predominantly mediates the killing of target cells by secreting perforin and granzymes (Jacobs et al., 2001), whereas the CD56^brightCD16^low NK cells that reside in secondary lymphoid tissues are immature and have reduced cytotoxic capability (Moretta, 2010). Recently, tumor- associated NK cells, which are enriched in tumors with impaired anti-tumor functions, were discovered to be associated with an unfavorable prognosis and resistance to immunotherapy in multiple solid tumors (Tang et al., 2023), such as lung cancer, pancreatic cancer, and esophageal cancer. The potential impact of NK cells on CCML progression is still poorly investigated.

In this study, the proportion of NK cells was higher than in past cognition, similar to a previous study (Xia et al., 2022) in that NK cells occupied the largest immune cell compartment of nearly 50% in cholangiocarcinoma. In contrast, NK cells may not be well distinguished during automatic annotation using SingleR. Consistent with the study that revealed that treatment responders are associated with TME remodeling, including NK cell recruitment (Kim et al., 2022), the proportions of NK cells in PR patients were significantly higher among the LiverP and LiverT groups (Fig. S1D and E), which also increased the overall proportion of NK cells to some extent. Additionally, we observed that NK cells are decreased in primary colon cancer and liver metastatic using scRNA-seq, consistent with previous findings supporting that NK cells play a crucial role in anticancer immunity (Laskowski et al., 2022, Myers et al., 2021). However, we further subdivided NK cells considering the conflicting result, in which NK cells were significantly enriched in the tumor area when using ST. Since this is the first study to comprehensively explore the heterogeneity of NK cells, no previous NK cell subpopulation markers were available; resting and activated NK cell markers in CIBERSORT were thus used to annotate NK cells. Consistent with a recent study that reported that mature and immature NK cells serve different functions (Thacker et al., 2023), our results also provide evidence that resting NK cells in the colon cancer TME possess a potential tumor-promoting activity, whereas activated NK cells play an anti-tumor role as conventional cognition. Bulk RNA-sequencing and the corresponding clinical data further indicated that higher infiltration of resting NK cells and lower infiltration of activated NK cells were correlated with a worse prognosis in patients with CRC. However, since the resting and activated NK cells we identified were only the aggregation of multiple NK cell subpopulations, further studies are required to investigate the key NK cell subpopulations that play a crucial role in CCLM.

Further analysis of NK cell subtypes revealed a dramatic change in KIR2DL4+ GPR171+ activated NK cells and GZMK+ resting NK cells, whereas no other NK cell subpopulations differed in colon cancer progression. KIR2DL4, also known as CD158d, is an unusual member of the killer cell Ig-like receptor family expressed in all NK cells (Andre et al., 2001). Studies revealed that KIR2DL4 activates the cytotoxicity of NK cells (Faure et al., 2002), and NK-92 cells stimulated with KIR2DL4 had higher cytotoxicity against breast cancer cells (Kilic et al., 2023). Consistent with previous studies, we also observed that the expression of the KIR2DL4+ GPR171+ activated NK cell subtype was decreased during CCLM progression, possibly reflecting the anti-tumor activity of this population. However, the GZMK+ resting NK cell subsets were first identified as potential tumor promoters. Moreover, it has been proven that GZMK^high CD8⁺ T effector memory cells, a cell subset that is particularly similar to GZMK+ NK cells, were associated with poor clinical outcomes in patients with colorectal tumors. Our results also showed that GZMK+ resting NK cells highly infiltrate the cancer region, and this predicts a worse prognosis in patients with colon cancer. However, despite the interesting and gratifying clinical relevance, the mechanism remains unknown.

Considering the evolutionary trajectory and the characteristic distribution of these metastasis-associated NK cell subtypes, we predicted that activated NK cells might differentiate into resting NK cells under the action of tumor cells. Consistent with a recent study that indicated that exposure to cancer cells causes NK cells to lose their cytotoxic ability (Chan, Knutsdottir, 2020), we identified the tendency of NK cells to differentiate to an resting state with decreasing expression of activated marker KIR2DL4 as well as typical functional NK cell markers: NKG2D (Yu et al., 2023) and NKG2A which are only triggered by cell-to-cell interactions as validated by FACS. Unfortunately, resting markers GZMK can not be detected in all the three groups since it’s a secreted protein.

Additionally, NK-mediated tumor cell editing appears to partly depend on the release of cytokines in in vitro co-culture experiments. SCF, the natural ligand of c- Kit with membrane-bound (mSCF) and soluble forms (sSCF) (Hsu et al., 1997), enhances the growth and migration of cancer cells such as cervical cancer (Aguilar et al., 2014), Ewing’s sarcoma (Landuzzi et al., 2000) and CRC (Yasuda et al., 2007). Although the dimeric mSCF is remarkably more active and induces a more persistent activation of the receptor (Hsu, Wu, 1997, Langley et al., 1993), the tumor-promoting effect of the supernatant from the co-culture system in the current study and the reversion effect of imatinib for the co-culture supernatant induced pro-oncogenic effects suggests an important role of sSCF on colon cancer cells. However, no up- regulation of KITLG (the SCF encoding gene) was found in HCT-116 after co- culturing with NK cells (Fig. S9A). In addition, when using DLD-1, a colon cancer cell with a nearly 15 times lower expression of KITLG than HCT-116 (Fig. S9B), we also did not observe co-cultured supernatant induced tumor-promoting effects. We thus hypothesized that in our experimental system, NK cells may promote the cleavage and release of the extracellular part of SCF on the tumor cell surface after direct contact with tumor cells, and thus increase the concentration of sSCF, which induces the malignant phenotype of colon cancer cells in a paracrine manner (Fig. 8J). However, the detailed mechanism still needs to be further studied.

In summary, we have unveiled the spatial and cellular profiles of TME from tumors and paratumor tissues of CCLM. Our analysis uncovered the different states, functions, and dynamic nature of NK cells in different CCLM settings, which can be used for further identification of regulatory mechanisms and for developing potential therapeutic targets.

Materials and Methods

Single-cell RNA Sequencing Analysis

The R (v4.0.5) package Seurat (v4.1.0) (Hao et al., 2021) was used to process the scRNA data. Since dataset quality control had been performed and Seurat objects had been created in previous studies, we did not further filter the scRNA-seq data or remove the batch effects. The SCTransform method was used to normalize the data. After selecting 2,000 highly variable genes using the Find Variable Features function in Seurat, principal component analysis was performed using these genes to reduce the data dimensions. Based on the ElbowPlot function in Seurat, we chose the top 30 principal components to run the FindNeighbors function. Next, the cells were clustered using the FindClusters function with a resolution of 0.1 for clustering immune cells. A uniform manifold approximation and projection algorithm was used for data visualization, as previously reported (Peng et al., 2023). Differentially expressed genes (DEGs) of each subset were identified using the FindAllMarkers function in Seurat. SingleR (v1.0.0) was used to name each cluster (Zheng et al., 2021).

Pseudotime Analysis

To analyze the differentiations of NK cells, monocle2 (http://cole-trapnell-lab.github.io/monocle-release), which uses an algorithm to learn the changes in gene expression sequences that each cell must undergo as part of a dynamic biological process (Hong et al., 2022), was used for pseudotime trajectory analysis to identify the transitional relationships among different clusters. The cells were reduced dimensionally using the DDRTree method, sequenced in pseudotime, and finally visualized (Qiu et al., 2017).

Bulk RNA-seq Data Analysis

The bulk transcriptome RNA-seq data and corresponding clinical data were obtained from The Cancer Genome Atlas of colon adenocarcinoma through the UCSC Xena browser (GDC hub) (https://gdc.xenahubs.net)(Goldman et al., 2020). In total, 459 colon cancer samples with survival information were downloaded, and 310 samples were finally enrolled (11 formalin-fixed samples were excluded, and an additional 134 and 4 samples were excluded due to deletion of metastasis information and sequencing matrix, respectively). Transcriptomic data from 65 colon cancer samples in GSE29623 were obtained as the validation cohorts. CIBERSORT (Chen et al., 2018), xCELL (Aran et al., 2017), EPIC (Racle et al., 2017), and MCP- counter (Becht et al., 2016), which use gene expression to infer the proportions of tumor-infiltrating immune cells, were used to analyze the TME cell type. R packages survival (v3.2-10) and survminer (v0.4.9) were used for survival analysis. The Youden index was selected as the cutoff value to differentiate patients into distinct groups (high or low). The Kaplan–Meier survival curve was modeled using the survfit function. Subsequently, a Cox proportional hazards regression model was established to determine the independent risk factors. DEGs between the metastasis and non- metastasis groups were determined with the filtering condition of log2 |fold change|>1 and p-value <0.05, using the R package limma. GO and KEGG pathway functional enrichment analyses were performed using the clusterProfiler R package (v3.18.1) to assign various biological processes, molecular functions, cellular components, and pathways of identified marker genes in the cluster of interest (Yu et al., 2012), and p <0.05 was regarded as statistically enriched. To explore the different KEGG pathways and hallmark gene sets between the metastasis and non-metastasis groups, gene set enrichment analysis (GSEA) was performed using data from The Molecular Signatures Database (c2.cp.kegg.v7.3.symbols) and the fgsea R package (v1.12.0) (Liberzon et al., 2015). Pathways with an adjusted p-value below 0.05 were deemed to be significantly enriched.

ST Data Analysis

Seurat was also used for ST data processing and visualization. We used the SCT method to normalize the ST data; the functions SelectIntegrationFeatures, PrepSCTIntegration, FindIntegrationAnchors, and IntegrateData were used to integrate the ST data. An unsupervised clustering method was subsequently used to cluster similar ST spots. Cell population annotations were based on hematoxylin and eosin (H&E) staining sections and the highly variable genes in each cluster. Scores of cell-specific signatures (top 20 DEGs) from scRNA-seq were calculated using the AddModuleScore function. SpatialDimPlot and SpatialFeaturePlot were combined to visualize the cell expression level in the ST data (Peng, Ren, 2023).

Cell Lines and Co-cultures

Colon carcinoma cell lines HCT-116 and DLD-1 were purchased from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China) and grown in a RPMI- 1640 medium (GIBCO/BRL) supplemented with 10% FBS (GIBCO/BRL) in a humidified chamber at 37 °C and 5% CO₂. The human NK cell line NK-92 was also purchased from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China) and maintained in RPMI-1640 supplemented with 12.5% horse serum (GIBCO), 12.5% FBS, 100 U/mL rhIL-2, 0.1 mmol/L b-mercaptoethanol, and 0.02 mmol/L folic acid.

Colon cancer cells were plated at a density of 1×10⁵ cells/well in 6-well plates for 24 h. Supernatants were collected, and cells and debris were removed by centrifugation. NK cells were co-cultured with cancer cells in a ratio of 1:1 in a fresh mixed medium for a further 24 h. Additionally, NK cells were cultured in the soluble supernatant for 24 h with a supernatant/RPMI-1640 10% FBS and NK medium in the ratio 1:1. When indicated, co-cultures were performed in Transwell devices (JET biofil), maintaining the same ratios and culture times.

Flow Cytometry Analysis

NK cells in each co-cultured group were collected and washed twice in 1 × PBA. Then, they were incubated with PE-labeled KIR2DL4 (Affinity Biosciences), FITC- labeled GZMK (Affinity Biosciences), FITC-labeled NKG2A (Biolegend) and PE- labeled NKG2D antibody (Biolegend) for 30 min, rewashed, and resuspended in PBS. BD flow cytometry was used for detection.

Cell Counting Kit-8 (CCK-8) Cell Viability and Cell Colony Formation Assay

Colon cancer cells (5,000 cells/well) were cultured in 96-well plates with with co-cultured supernatants from different co-culture groups alone or with both co- cultured supernatants and imatinib mesylate (2 µM) as previously reported (Yasuda, Sawai, 2007) for 24 h. At pre-determined time points, 10 μL of CCK-8 reagent (Dojindo, Japan) was added and incubated for 2 h at 37 °C, and then the absorbance was measured using a microplate reader (Thermo Scientific) at 450 nm. All experiments were carried out in triplicate.

For the cell colony formation assay, colon cancer cells (500 cells/well) were seeded in 12-well plates with co-cultured supernatants alone or with both co-cultured supernatants and imatinib mesylate at 37 °C for 1 week. Then, cell colonies were fixed with 4% paraformaldehyde for 10 min and stained with 0.5% crystal violet for 5 min. Cell colonies containing >20 cells were counted. All experiments were carried out in triplicate.

Transwell Migration/Invasion Assays

The polycarbonate membrane in the transwell chambers was coated with Matrigel (Corning, NY, USA). Next, we transferred 1 × 10⁵ cells from the serum-free medium with or without imatinib mesylate into the top chamber, added co-cultured supernatants alone or with both co-cultured supernatants and imatinib mesylate in the bottom chamber, and incubated at 37 °C for 24 h. Then, we removed the non-invading cells on the top side of the membrane by scrubbing, fixed the migrating or invading cells at the bottom side of the membrane with 4% paraformaldehyde, and stained with 0.5% crystal violet. The number of cells was counted under a microscope (Leica, London, UK) from four randomly chosen fields per well to determine the number of cells in each group.

Total RNA Extraction and Quantitative Real-Time PCR (qPCR)

Total RNA was extracted from cells using TRIZOL reagent (Invitrogen, USA). RNA was reverse transcribed into cDNA using a PrimeScript RT Reagent Kit (Takara, Japan). qPCR was performed using QuantStudio™ Test Development Software (Thermo Scientific, Waltham, MA, USA) with SYBR Green qPCR Master Mix (EZBioscience, Roseville, MN, USA). The sequence of KITLG and housekeeping gene GAPDH primers is listed in Table S1. KITLG mRNA data were normalized to that of GAPDH.

Luminex Liquid Suspension Chip Detection and ELISA

Luminex liquid suspension chip detection was performed by Wayen Biotechnologies (Shanghai, China). The Bio-Plex Pro Human Chemokine Panel 48- plex kit was used in accordance with the manufacturer’s instructions. Briefly, supernatants from different co-cultured groups were incubated in 96-well plates embedded with microbeads for 1 h and then incubated with a detection antibody for 30 min. Subsequently, streptavidin-PE was added in each well for 10 min, and values were read using the BioPlex MAGPIX System (Bio-Rad). The Human SCF ELISA kit (Solarbio, China) was used according to the manufacturer’s instructions.

Statistical Analysis and Visualization

All statistical analyses were performed using SPSS (v23.0; IBM SPSS, Chicago, IL, USA) and R (v4.0.5), and data visualization was performed on R packages Seurat (v4.1.0), ggplot2 (v3.3.5), ggsignif (v0.6.1), and pheatmap (v1.0.12).

Data Availability

The raw data of single-cell RNA sequencing and ST data (10X Genomics) were derived from previously published research at the website: http://www.cancerdiversity.asia/scCCLM/)(Wu, Yang, 2022). And two bulk transcriptomics data of colon cancer were downloaded from TCGA cohort COAD at the website: https://xena.ucsc.edu/ and NCBI’s Gene Expression Omnibus with accession number GSE29623. All data generated or analyzed during this study are available from the corresponding author on reasonable request.

Funding

Acknowledgements

We thank Prof. Qiang Gao from Fudan University and Prof. Xiaoming Zhang from Shanghai Institute of Immunity and Infection, CAS for sharing the SC RNA sequencing and ST data. We thank Zhejiang Key Laboratory of New Techniques for Diagnosis and Treatment of Critical Diseases of Pancreas and Liver and Zhejiang International science and Technology Cooperation base for tumor transformation research for supporting this study. We thank Editage Group (https://www.editage.cn/) for polishing the draft of this manuscript.

Authors’ Contributions

CCM and YYC: Conceptualization, experimental design and execution, data curation, method development, formal analysis, statistics, writing original draft, and editing the manuscript. DX and TMZ: conceptualization, data curation and method development. DFM and ZH: experimental design. WKX and CL: experimental design. YXL and JYY: Validation. DX, MDL and XYX: Experimental design, supervision, funding acquisition, and manuscript editing. XS: Conceptualization, supervision, writing manuscript, project administration.