- Research
- Open access
- Published:
Single-cell transcriptomic analysis reveals dynamic changes in the liver microenvironment during colorectal cancer metastatic progression
Journal of Translational Medicine volume 23, Article number: 336 (2025)
Abstract
Background
Metastasis is a leading cause of cancer-related deaths, with the liver being the most frequent site of metastasis in colorectal cancer. Previous studies have predominantly focused on the influence of the primary tumor itself on metastasis, with relatively limited research examining the changes within target organs.
Methods
Using an orthotopic mouse model of colorectal cancer, single-cell sequencing was employed to profile the transcriptomic landscape of pre-metastatic and metastatic livers. The analysis focused on identifying cellular and molecular changes within the hepatic microenvironment, with particular emphasis on inflammatory pathways and immune cell populations.
Results
A neutrophil subpopulation with high Prok2 expression was identified, showing elevated levels in the pre-metastatic and metastatic liver. Increased infiltration of Prok2⁺ neutrophils correlated with poor prognosis in liver metastatic colorectal cancer patients. In the liver metastatic niche (MN), these neutrophils showed high App and Cd274 (PD-L1) expression, suppressing macrophage phagocytosis and promoting T-cell exhaustion.
Conclusion
A Prok2⁺ neutrophil subpopulation infiltrated both pre-metastatic and macro-metastatic liver environments, potentially driving immunosuppression through macrophage inhibition and T-cell exhaustion. Targeting Prok2⁺ neutrophils could represent a novel therapeutic strategy for preventing liver metastasis in colorectal cancer patients.
Background
Colorectal cancer (CRC) is a prevalent malignancy, ranking third in incidence and second in mortality globally [1, 2]. While the 5-year survival rate for early-stage colorectal cancer following radical surgery exceeds 90%, this rate plummets to just 10–20% once distant metastasis occurs [3]. The liver is the primary site for CRC metastasis, with about 50% of CRC patients developing liver metastases during the course of their disease, making liver metastasis the leading cause of death in CRC patients [4]. Unfortunately, colorectal cancer liver metastasis (CRCLM) patients typically show poor response to various treatment options, highlighting the critical need for a deeper understanding of the mechanisms driving CRCLM [5].
The pre-metastatic niche (PMN) plays a vital role in tumor progression, forming prior to metastasis and serving as a key step in the metastatic process [6,7,8]. The liver's rich blood supply, coupled with its immunosuppressive environment, makes it a frequent site of metastasis [9]. Researchers used bulk RNA sequencing to analyze changes in the liver before metastasis following orthotopic cecal modeling in mice, revealing an increase in MDSCs in the pre-metastatic liver [10]. However, bulk RNA-seq is insufficient to analyze transcriptional characteristics of specific cell types. The liver consists of parenchymal cells (hepatocytes) and non-parenchymal cells (such as liver sinusoidal endothelial cells, hepatic stellate cells, and Kupffer cells), all of which contribute to the formation of the PMN [11, 12]. For example, hepatocytes secrete serum amyloid A (SAA), which facilitates the establishment of the PMN. Meanwhile, hepatic stellate cells become activated, M2 macrophages are recruited to release immunosuppressive factors, and myeloid-derived suppressor cells (MDSCs) are attracted to the site [13,14,15]. However, the precise mechanisms underlying PMN formation and the complex interactions among various cell types remain poorly understood.
Previous single-cell transcriptomic studies on the liver microenvironment during metastasis have largely focused on established liver metastases. Although differences between metastatic lesions and primary tumors have been observed, there is a lack of comprehensive transcriptomic analysis tracking the liver's progression from normal tissue, through the pre-metastatic stage, to the formation of metastatic lesions. Understanding the changes in the liver from normal to pre-metastatic states can better help us comprehend the early processes of metastasis formation. Due to the difficulty of obtaining clinical samples of normal and pre-metastatic liver tissue, we utilized a mouse model of colorectal cancer with orthotopic cancer cell injection. Based on time points and liver status, we categorized the liver samples into normal liver, pre-metastatic liver, and early tumor-colonized liver, followed by single-cell transcriptome sequencing.
In this study, we constructed a single-cell atlas of the liver PMN and metastatic niche (MN) in mice, analyzing the changes in various cell types during this process. Among these cell populations, we identified a special group of neutrophils that highly express Prok2 (Prokineticin 2). These neutrophils show increased infiltration in the PMN liver and further accumulation in the MN liver. Additionally, we found that they are associated with exhausted T cells in the liver and are linked to poor survival in patients with colorectal cancer liver metastasis. Studies have shown that infiltrating neutrophils can promote metastasis through the secretion of growth factors such as Prok2 [16]. However, the value of Prok2+ neutrophils within the PMN, as well as their interactions with other cells at different stages of metastasis, has not yet been reported. Prok2+ neutrophils may represent a potential therapeutic target for anti-metastatic immunotherapy.
Materials and methods
Cell line
CT26-GFP cells were cultured in 1640 medium supplemented with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin at 37 °C in a humidified 5% CO2 incubator. GFP was overexpressed in CT26 cells to enable flow cytometric detection of tumor cells in the livers of the orthotopic mouse model of colorectal cancer.
Orthotopic mouse model of colorectal cancer
All experimental animals were obtained from the Animal Center of the Sixth Affiliated Hospital of Sun Yat-sen University (Guangzhou, China) and handled in strict accordance with protocols approved by the Institutional Animal Care and Use Committee of the same institution. Female BALB/c mice, aged 4–6 weeks, were randomly assigned to either the control or tumor groups. Under anesthesia, both groups underwent surgery. A sterilized incision was made just left of the center of the abdomen to expose the cecum. In the tumor group, 1 × 10^6 CT26-GFP cells, suspended in 50 μL of PBS, were slowly injected into the plasma membrane layer of the cecum. The control group received an equal volume of PBS. After the injection, a sterile cotton swab was applied to the injection site for 2 min to prevent leakage of the CT26-GFP cell suspension into the peritoneal cavity. Finally, the incision was sutured using sterile needle and thread. The maximum diameter of all tumors is less than 1.5 cm, in compliance with the ethical standards of the animal center.
Isolation of liver cells
Liver single cell suspension preparation was performed by the laboratory staff of NovelBio Biopharmaceutical Technologies Ltd. Liver lysis was performed as described previously [17], whole mouse livers were removed, cut into small pieces on ice, and then enzymatically lysed with 1 mg/mL Collagenase II (Worthington) for 30 min on a shaker at 37 °C. 70-micron cell filters were sieved and then erythrocytes were lysed in erythrocyte lysis buffer (Miltenyi Biotec), rinsed twice with PBS, and then the MACS Dead Cell Removal Kit (Miltenyi Biotec) was used to further enrich the single cell suspension.
Single-cell library preparation and sequencing
Single cell capture was realized using the BD Rhapsody system. Whole transcriptome libraries were prepared according to the BD Rhapsody Single Cell Whole Transcriptome Amplification Workflow and 150 bp paired-end sequencing was performed using the NovaSeq 6000 and NovaSeq X Plus (Illumina, San Diego, California, USA).
Single-cell sequencing quality control
Matrix data were organized and analyzed using the Seurat R package [18]. Low-quality cells were filtered based on the following criteria: mitochondrial gene expression > 25%, erythroid gene expression > 1%, a minimum of 200 gene tests, and Unique Molecular Identifiers (UMIs) > 1000. Additionally, doublet prediction and filtering were performed using the DoubletFinder R package [19].
Cell type annotation
To eliminate batch effects between samples, we applied the default parameters of the Harmony R package [20] to integrate the data from each sample. To facilitate the comparison of gene expression across samples, we normalized the data and selected the top 4000 highly variable genes for further analysis using a scaling factor of 10,000. The data were log-transformed using the “ScaleData” function from the Seurat R package [18], and the first 50 principal components were reduced using PCA. Dimensionality reduction was performed using PCA, and clustering was done using the first 20 principal components with 'FindNeighbors' and 'FindClusters' functions. Cell types were identified based on known markers, and we categorized the cells into 8 major types: B cells, myeloid cells, T/NK cells, cholangiocytes, endothelial cells, hepatocytes, fibroblasts, and tumor cells. To further investigate liver immune cell subpopulation changes during colorectal cancer liver metastasis, T/NK cells and myeloid cells were isolated using the “Subset” function, normalized again, and cell subpopulations were identified based on known markers. Dimensionality reduction clustering was performed using the Harmony R package [20].
Deconvolution
We used single-cell sequencing data from the control and PMN groups as a reference to predict the immune cell abundance of each sample from the bulk RNA-seq data using the BisqueRNA R package [21]. The datasets included GSE109480 (which collected pre-metastatic liver samples using the LSL-KrasG12D/+; LSL-Trp53R127H/+; Pdx1-cre (KPC) mouse model of pancreatic ductal adenocarcinoma (PDAC)) [22], PRJNA590588 [23] (which collected pre-metastatic niche (PMN) and metastatic niche (MN) liver samples using an intrasplenic injection mouse model) and GSE147044 (which collected pre-metastatic liver samples using an orthotopic colorectal cancer mouse model) [10]. These three datasets from the GEO database were used to validate the cell occupancy of each population in the control and PMN groups in single-cell sequencing, totaling 13 healthy mouse livers and 17 PMN livers. These datasets were selected due to their comprehensive coverage of pre-metastatic and metastatic niches in similar mouse models.
Differential gene analysis and enrichment analysis
Differential gene expression in hepatocytes during liver metastasis of colorectal cancer was analyzed using the “FindMarkers” function from the Seurat R package [18]. Top 20 up-regulated genes were labeled in a volcano plot. The “fgsea” function in the clusterProfiler R [24] package was used for the Gene Set Enrichment Analysis (GSEA) analysis. Differential analysis of bulk RNA-seq data was done using the limma R package [25], and finally ggpubr R package [26] was used to show the top 20 genes that were differentially up-regulated in the PMN group compared to the control group in the volcano plot.
Survival analysis
Data for the survival analysis were obtained from GSE159216 [27], which includes 283 patient samples with colorectal cancer liver metastasis. The RNA sequencing data was derived from these metastases, and patients were monitored for prognosis. The survival R package [28] was employed to conduct the survival analysis. Human gene symbols were converted to mouse symbols using the “bitr” function from the clusterProfiler R package [24]. Subsequently, the BisqueRNA R package [21] was utilized to estimate the N0-N3 cluster neutrophil percentages for each colorectal cancer liver metastasis patient, based on N0-N3 neutrophil single-cell sequencing data through deconvolution. The classification of patients into high or low N1-N3 groups was determined using the “surv_cutpoint” function in the survminer R package [26]. Finally, Kaplan–Meier survival curves were generated to visualize the survival outcomes.
Integrated analysis of human neutrophil single-cell sequencing data
Our research indicates that neutrophils exhibit the most significant changes in proportion between the pre-metastatic niche (PMN) and metastatic niche (MN) in the liver. Neutrophils are primarily generated in the bone marrow [29, 30]. Once they mature, these cells enter the peripheral blood and migrate to the liver as needed. Considering that peripheral blood is easily accessible and that literature reports elevated levels of neutrophils in peripheral blood prior to tumor metastasis [31], we collected single-cell sequencing data from both the peripheral blood and liver tissue of healthy individuals and patients with liver metastasis to further investigate the proportion of neutrophils. First, we collected healthy human peripheral blood single-cell sequencing data from the literature [32], which includes data from 2 healthy human peripheral blood samples. To effectively capture neutrophils, the Microwell-seq platform was utilized for this sequencing. Additionally, we obtained healthy human liver single-cell sequencing data from the GSE134355 dataset [33], which includes 5 liver samples. This sequencing also used the Microwell-seq platform. Furthermore, sequencing data for peripheral blood and liver single cells from patients with liver metastases of colorectal cancer were sourced from Wu et al. [34]. This dataset contains 11 liver metastasis samples along with 2 peripheral blood samples, sequenced using the 10 × Genomics platform. To identify neutrophils, we focused on specific markers: CXCR2, S100A8, G0S2, S100A12, and CSF3R.
To facilitate better comparison of gene expression profiles across different sequencing platforms, we utilized the Harmony R package [20] to integrate data from various platforms based on distinct samples. The calculation of the N1 signature was performed using the “AddModuleScore” function from the Seurat package [18], enabling us to score each cell accordingly.
Gene set score
N1 neutrophils were analyzed differentially in comparison to all other cell types using the “FindMarkers” function. Genes exhibiting high expression levels (p < 0.05, LogFC > 2) were intersected with those associated with poor prognosis in patients with liver metastases from colorectal cancer (GSE159216). This analysis yielded a total of 24 genes identified as N1 signatures. The N1 neutrophil-associated signature, along with M1 and M2 signatures, was compiled from various literature sources and is organized in Table S3.
Trajectory analysis
Pseudotime analysis was conducted using the Monocle3 R package [35]. Based on existing literature and gene expression data, the pseudotime of all cell types was predicted, utilizing cells at early stages of development as the root node. For instance, the developmental onset corresponding to the putative time allowance for neutrophils was identified as N3 neutrophils, characterized by high expression levels of genes such as Camp and Ltf. The root nodes for macrophages were designated as Kuppfer3 and MAC3, respectively.
Cell–cell interactions
To examine the changes in communication between different cells in the liver microenvironment during various stages of colorectal cancer liver metastasis, we employed the CellChat R package [36] for our analysis. This package utilizes a database of receptor pairs known as CellChatDB, which is a manually curated resource that incorporates data from the KEGG database and receptor pairs reported in the literature. Using the cell–cell contact library within the CellChatDB, we predicted the likelihood and p-values of cell–cell interactions. Finally, we presented the receptor ligands that exhibit differences across the various subgroups.
Immunofluorescence
Tissue was firstly fixed in 4% paraformaldehyde, embedded in paraffin and sectioned at 4 μm thickness. Then those paraffin-embedded tissue sections were deparaffinized, rehydrated and microwave antigen retrieved in EDTA buffer. Sections were blocked with 5% BSA for 30 min at room temperature. Then those samples were following incubated with anti-MPO (5ug/ml, AF3667, R&D) and anti-PROK2 (1:50, PS04727S, Abmart) overnight at 4 °C. After rinsing by PBS, fluorochrome conjugated secondary antibodies (1:1000, A11012, Invitrogen; 1:1000, ab150129, Abcam) were added for incubation 1 h at room temperature. After rinsing by PBS, DAPI (D1306, Invitrogen) was then used for staining the nuclei. Observation and photographing were performed with the confocal microscopy Cell Observer (ZEISS, Germany). Briefly, the positive for Prok2 and Mpo were count as Prok2+ neutrophils.
Results
Single-cell atlas of the liver across different stages of colorectal cancer metastasis
To temporally track the alterations of live microenvironment during liver metastasis of colorectal cancer, the GFP-labelled CT26 cells were implanted in the cecum of mice to establish an orthotopic model. Liver samples were collected on day 19 (n = 2) and 38 (n = 2) representing pre-metastatic and metastatic niche respectively determined by flow cytometry to evaluate the metastatic lesions with existence of GFP positive tumor cells. Simultaneously, the liver tissues from the same-age mice without harboring tumor were collected as control (n = 2). As expected, rare GFP positive cells were colonized in the liver with obvious tumor burden in the cecum on day 19, while around 2% GFP positive tumor cells were captured in the liver on day 38, suggesting a successful formation of pre-metastatic and metastatic niche in the murine model (Fig. 1A, B, C and Fig. S1). Next, single-cell suspensions were prepared from the aforementioned livers. To capture the diverse range of cell types consisting the liver microenvironment we avoided cell sorting step, and single-cell capture was achieved by using a BD Rhapsody system (Fig. 1D). Following preliminary quality control and dimensionality reduction clustering, a total of 53,966 high-quality cells were obtained (Table S1). These included immune cells such as T/NK cells, B cells, and myeloid cells; non-immune cells such as endothelial cells and fibroblasts; as well as epithelial cells, including tumor epithelial cells, hepatocytes, and cholangiocytes (Fig. 1E). The cell types belonging to each cluster were annotated using known markers (Fig. 1F) [37,38,39]. Intriguingly, we observed significant changes of the hepatic microenvironment starting from the PMN formation. In the PMN, there was an obvious increase in myeloid cells, while the infiltration of T/NK cells, endothelial cells, and B lymphocytes showed a slight decrease, with no change observed in the abundance of hepatocytes. Additionally, the proportion of myeloid cells further increased, accompanied by significant reduction of T/NK cells, endothelial cells, and B cells in the MN (Fig. 1G). Consistent with FACS analysis, single-cell sequencing data showed that no tumor cells were captured in the liver of PMN, whereas a small number of tumor cells were present in the liver of MN (Fig. 1H). We also validated our findings by analyzing bulk RNA-seq data available in public databases, which included 13 healthy liver samples and 17 PMN liver samples. We performed immune cell deconvolution to predict the proportions of T/NK cells, myeloid cells, and B cells. In line with our sequencing results, the PMN liver tissues displayed increased infiltration of myeloid cells and reduced presence of T/NK cells and B cells compared to the control group (Fig. 1I, J, K).
Single-cell atlas of liver in a mouse model of colorectal cancer liver metastasis. Flow cytometry analysis of GFP signals in liver tissues from control (A), pre-metastatic niche (PMN) (B), and metastatic niche (MN) (C) groups of the colorectal cancer orthotopic mouse model. D Diagram of liver single-cell suspension preparation. E Uniform Manifold Approximation and Projection (UMAP) plot of liver samples from all groups. F Cell type annotation markers. G The proportion of each cell type in different groups. H The proportion of tumor cells in the control, PMN, and MN groups. The proportion of T/NK cells (I), Myeloid cells (J) and B cells (K) among immune cells in the control and PMN groups. P value was calculated using Wilcoxon test
Inflammatory pathways are significantly activated in hepatocytes from both the PMN and MN of colorectal cancer liver metastasis
Given that hepatocytes, account for 80% of all liver parenchymal cells [40, 41] and are responsible for carrying out most liver functions [42, 43], we next sought to explore transcriptomic changes in hepatocytes during metastatic progression. Our findings revealed that inflammatory protein-related genes, such as Saa1, Saa2, Orm1, Orm2, Lcn2 and the granulocyte-associated chemokine Cxcl1, were significantly elevated in hepatocyte prior to the onset of colorectal cancer liver metastases (Fig. 2A, Table S2). In the MN, hepatocytes express granulocyte chemotactic factors such as S100a8 and S100a9, as well as additional chemotactic factors like Cxcl2 and Cxcl10 (Fig. 2B, Table S2). Moreover, compared to hepatocytes in the PMN, the expression of these inflammatory factors is further elevated within the MN (Fig. 2C). Interestingly, the inflammation-related proteins secreted by hepatocytes, such as Saa1, Saa2, Orm1, Orm2, and Lcn2, do not show a further significant increase (Fig. 2C, Table S2). Inflammation-associated pathways, including the IL-6/JAK/STAT3 and TNF-α/NF-κB pathways, were activated in hepatocytes before liver metastases occurred (Fig. 2D). Consistent with our study, activation of IL-6-STAT3-SAA signaling in hepatocytes has been shown to support creating a metastasis-favorable microenvironment [22]. IL-6/JAK/STAT3 and TNF-α/NF-κB were more strongly activated in hepatocytes within the MN compared to the PMN Additionally, the bile acid metabolism pathway is downregulated in the PMN and further downregulated in the MN, suggesting that metabolic changes in the liver also play a role in the formation of this process (Fig. 2E, F).
The differential genes and pathways in hepatocytes between the pre-metastatic niche (PMN) and metastatic niche (MN). A Differential volcano plot of hepatocytes between the PMN group and the control group. B Differential volcano plot of hepatocytes between the MN group and the control group. C Differential volcano plot of hepatocytes between the MN group and the PMN group. Red dots represent the top 20 upregulated genes, orange dots represent upregulated genes with p < 0.05 and Log2fc > 1, and blue dots represent downregulated genes with p < 0.05 and Log2fc < −1. GSEA analysis of differential pathways in hepatocytes between the PMN group vs control (D), MN group vs control (E), and MN group vs PMN group (F). Red indicates upregulated pathways, and blue indicates downregulated pathways. NES represents normalized enrichment score. G Differential volcano plot in the pre-metastatic niche (PMN) liver of the pancreatic cancer orthotopic liver metastasis model compared to the control group. H Differential volcano plot in the PMN liver of the pancreatic cancer spleen metastasis model compared to the control liver. I Differential volcano plot in the PMN liver of the colorectal cancer orthotopic liver metastasis model compared to the control liver, showing the top 19 upregulated genes in PMN compared to control
To further validate our conclusion, we analyzed the bulk RNA-seq data of a pre-metastatic pancreatic cancer liver model [10, 22, 23] and similarly found that inflammation-related protein genes, including Saa1, Saa2, Orm1, Orm2, Lcn2, and Cxcl1, were significantly upregulated in the PMN liver (Fig. 2G, H, Table S2). Moreover, in the bulk RNA-seq data of the PMN liver from a colorectal cancer orthotopic liver metastasis model [10], these inflammatory proteins as well as Cxcl1 were also significantly elevated compared to the sham-operated liver group (Fig. 2I, Table S2). These results suggest that in both PMN and MN, hepatocytes, as the most abundant cells in the liver, undergo inflammatory activation as their primary change. The high expression of inflammatory proteins and myeloid cell-associated cytokines may be the initiating factors for the increased infiltration of myeloid cells in PMN and MN.
N1 neutrophils (Prok2+ neutrophils) are increased in PMN, MN and associated with poor prognosis in patients with colorectal cancer liver metastasis
Increased myeloid cell infiltration is one of the most prominent features in PMN and MN (Fig. 1G). Next, we focus on the variations among different types of cells within the myeloid cell population. we conducted a detailed subpopulation analysis of hepatic myeloid cells using known markers [37] (Fig. 3A, B). Seven distinct subpopulations of myeloid cells were identified after UMAP clustering: macrophages, neutrophils, monocytes, DC1, DC2, pDC cells, and basophils. The proportion of neutrophils progressively increases in both PMN and MN, while the proportion of monocytes also rises in these stages. Conversely, pDCs progressively decrease in PMN and MN, cDCs mainly decrease in MN, and macrophages show no significant change in proportion (Fig. 3C).
N1 neutrophils infiltrate more in the liver pre-metastatic niche (PMN) and metastatic niche (MN) are associated with poor prognosis in colorectal cancer liver metastasis patients. A UMAP plot of myeloid cells. B Myeloid cell marker expression plot. C Proportion of myeloid cell subpopulations in control, PMN, and MN. D UMAP plot of neutrophils. E Proportions of neutrophil subpopulations. F Proportion of N0-N3 neutrophils in the PMN liver of the colorectal cancer liver metastasis model. G Proportion of N0-N3 neutrophils in the PMN liver and MN of the pancreatic cancer liver metastasis model. H Prok2+ neutrophils in the livers of the control, PMN, and MN groups through co-staining Mpo and Prok2 using immunofluorescence (IF). I Proportion of N0-N3 neutrophils across different pseudotime stages. J Proportion of cells at different pseudotimes in control, PMN, and MN. K Heatmap showing neutrophil-related signatures in N0-N3 neutrophils, with values scaled from −1 to 1. L Survival curve of colorectal cancer liver metastasis patients based on the high or low proportion of N1 neutrophils. M Expression levels of the N1 neutrophil signature in peripheral blood and liver tissues of colorectal cancer liver metastasis patients and healthy controls. "Liver normal" and "liver tumor" refer to liver tissues from healthy donors and colorectal cancer liver metastasis patients, respectively. "PBMC normal" and "PBMC tumor" refer to peripheral blood mononuclear cells from healthy donors and colorectal cancer liver metastasis patients, respectively
Among these, neutrophils showed the most significant increase in PMN and MN (Fig. 3C). Thus, we focused on neutrophils and used UMAP to categorize them into four subpopulations (N0, N1, N2, N3) to further explore their potential biological states (Fig. 3D). Neutrophil subtypes N1 and N3 followed the same pattern as the overall neutrophil increase in myeloid cells, with gradual elevation during PMN and MN formation (Fig. 3E). To determine whether this phenomenon is general, predictive deconvolution analysis was performed using bulk RNA-seq data of CRC and PDAC samples from Zeng et al. [10] and Li et al. [23]. Interestingly, in line with our finding, the data also indicated a significant rise in N1 and N3 neutrophil infiltration during the PMN phase (Fig. 3F, G). Notably, N1 neutrophils increased further as metastasis progressed (Fig. 3G). N1 and N3 neutrophils exhibited high expression of early neutrophil markers like Retnlg and Mmp8 compared to other neutrophil subtypes (Fig. S2A, S2B, S2C and Table S3). Immature and early neutrophils are closely linked to tumor metastasis [44, 45]. While both N1 and N3 populations shared certain gene expressions, N3 neutrophils expressed earlier bone marrow-derived neutrophil markers like Ltf and Camp at higher levels than N1 neutrophils [46] (Fig. S2A, Fig. S2D, S2E).
Conversely, N1 neutrophils had higher expression of metastasis-related genes such as Prok2 (Fig. S2F), which is critical in chemotherapy-promoted breast cancer metastasis and serves as a mediator of neutrophil-dependent angiogenesis during early tumor progression [16, 47]. We have confirmed the results by detecting Prok2+ neutrophils in the liver of the control, PMN, and MN groups by co-staining Mpo and Prok2 with immunofluorescence (IF). The result is quite consistent with our single-cell analysis (Fig. 3H). Pseudotime analysis indicated that N3 neutrophils developed earlier than N1 neutrophils (Fig. 3I, Fig. S2G). In contrast, control neutrophils developed later and more maturely than those in PMN and MN liver neutrophils (Fig. 3J).
Gene set scoring revealed that N1 neutrophils expressed high levels of MDSC-related genes, suggesting a potential link between N1 neutrophils and tumor immunosuppression. In comparison, N3 neutrophils expressed genes linked to phagocytosis and granule formation, placing them in a more immature state (Fig. 3K). Meanwhile, N0 and N2 neutrophils exhibited characteristics of senescence and apoptosis (Fig. 3K). Next, we focus on these specific neutrophil subpopulations to explore whether they are associated with the survival of patients with liver metastasis. Interestingly, among four subpopulations of neutrophil, we found that only the proportion of N1 neutrophils was associated with poorer prognosis in patients with colorectal cancer liver metastases. (Fig. 3L, Fig. S2H, S2I, and S2J).
To further investigate the predictive potential of N1 neutrophils, we identified an N1 signature by intersecting genes highly expressed by N1 neutrophils with genes associated with poor prognosis in colorectal cancer liver metastases (Fig. S2K and Table S3). As anticipated, the N1 signature correlated significantly with poor patient outcomes (Fig. S2L). Moreover, sequencing data from healthy individuals and patients with colorectal cancer liver metastases revealed that the N1 signature was significantly elevated in both peripheral blood and liver neutrophils of patients compared to healthy individuals (Fig. 3M). This indicates that the N1 signature is a potential biological marker for predicting the prognosis of liver metastasis in colorectal cancer. The above results suggest that we have identified a neutrophil subpopulation likely playing a key role in PMN and MN formation. This subpopulation expresses higher levels of immunosuppressive and tumor metastasis-related genes, which are associated with poorer prognosis in metastatic tumors.
Increased infiltration of immunosuppressive macrophages was primarily observed in MM
Next, we focused on macrophages, as they played a crucial role in the pre-metastatic lung microenvironment in cancer metastasis [8, 15, 48]. Using UMAP dimensionality reduction clustering, we identified six types of macrophages (Fig. 4A), among which three types expressed Kupffer cell markers such as Cd5l, Timd4, Lipa, and Hmox1 (Fig. 4B) [9]. Kupffer cells (KCs) are resident macrophages in the liver and represent the largest group of mononuclear phagocytes in the body [49]. KCs rapidly differentiate into mature macrophages in response to various injuries and play an important role in controlling liver inflammation and immunity [50]. The three types of KCs showed differences in characteristic gene expression, such as Kupffer2, which highly expressed cytokines like Il6 (Fig. 4C and Table S4). Kupffer2 cells also showed a significant increase in infiltration in MN, while Kupffer3 cells showed decreased infiltration (Fig. 4D).
Increased infiltration of suppressive macrophages in the metastatic niche (MN). A UMAP plot of macrophages. B Bubble plot of Kupffer cell markers. C Display of the top 5 markers for each macrophage cluster. D Proportion of macrophage clusters. E UMAP trajectory plot of macrophages in pseudotime. F Heatmap showing the expression of genes in Kupffer2 and Kupffer3 cells, with the top 40 genes upregulated in Kupffer2 compared to Kupffer3. G GO enrichment bubble plot of differential genes between Kupffer2 and Kupffer3. GO enrichment was performed with p < 0.05 and LogFC > 1, showing the top 10 functions based on p-value. H KEGG pathway enrichment bubble plot. I UMAP plot of CD14 signature expression values. J Diagonal plot of differential genes between MAC2 and MAC1, percentage difference is calculated by subtracting the expression levels between two cell types. K, L UMAP plots of M1 and M2 macrophage signature expression scores
Pseudotime analysis revealed that Kupffer2 cells develop later than Kupffer3 cells (Fig. 4E). We conducted differential analysis between Kupffer2 and Kupffer3 cells to explore how tumors impact macrophages in the MN. Kupffer2 cells in the MN secrete large amounts of cytokines, chemokines, and Il6 (Fig. 4F and Table S4). GO enrichment analysis showed that Kupffer2 cells were significantly enriched for functions related to neutrophil chemotaxis and cell–cell interaction compared to Kupffer3 (Fig. 4G and Table S4). KEGG enrichment analysis further revealed that Kupffer2 cells activated the TNF/NF-kB pathway more than Kupffer3 cells (Fig. 4H).
According to McGinnis CS et al., myeloid cells activated by TLR/NF-kB inflammation infiltrate the lung microenvironment and are closely associated with lung metastasis of breast cancer. Notably, one of the essential subtypes of myeloid cells in this process is CD14+ “activated” cells [51]. Interestingly, our study also found that the Cd14 inflammation score was predominantly expressed in Kupffer2 cells (Fig. 4I).
In addition to Kupffer cells, we identified other groups of macrophages, including MAC1, MAC2, and MAC3. Among them, MAC2 infiltration significantly increased in the MN, while MAC1 infiltration decreased (Fig. 4D). Differential analysis revealed that MAC2 expressed higher levels of M2 macrophage markers such as Arg1 and Ctsd compared to MAC1 (Fig. 4J and Table S4). Additionally, as expected, M1 tumor-suppressive macrophage markers were expressed at low levels in most macrophages (Fig. 4K), while M2 tumor-promoting markers were highly expressed in MAC2 and Kupffer2 cells (Fig. 4L).
In summary, our study demonstrates that two immunosuppressive macrophage populations, Kupffer2 and MAC2, show significantly increased infiltration in the MN.
Exhausted T cells mainly accumulate in the MN
Next, we focus on lymphocytes which primarily function as executor of anti-tumor immune response. The main cytotoxic cells in antitumor immunity are T cells and NK cells [52,53,54]. Using known markers, we classified the T/NK cells into nine subpopulations (Fig. 5A, B). Interestingly, exhausted T cells showed a significant increase in infiltration in the MN, and Tregs were also increased in the MN (Fig. 5C). This suggests that the MN in the liver is an immunosuppressive microenvironment, promoting tumor growth and colonization.
Increased infiltration of exhausted T cells in the metastatic niche (MN). A UMAP plot of T/NK cells. B Marker expression for T/NK cells. C Proportion of T/NK cell subpopulations. D Differential volcano plots of various T/NK cell subpopulations comparing MN with control. E GO enrichment bubble plot showing the results of GO pathway enrichment for differential genes between MN exhausted T cells and control exhausted T cells. F Bar chart showing the correlation between the proportion of N1 neutrophils and exhausted T cells in colorectal cancer liver metastasis patients. G Bar chart showing the correlation between the proportion of Kupffer2 macrophages and exhausted T cells in colorectal cancer liver metastasis patients. H Bar chart showing the correlation between the proportion of N1 neutrophils and Kupffer2 macrophages. I Bar chart showing the correlation between the proportion of N1 neutrophils and MAC2 macrophages. The p-values in these figures were calculated using the Wilcoxon test
Further analysis revealed that exhausted T cells in the MN expressed higher levels of genes related to T cell dysfunction, such as Fosl2 and Crem, compared to exhausted T cells in control samples (Fig. 5D). Additionally, GO enrichment analysis indicated that MN T cells were primarily enriched in pathways related to T cell apoptosis (Fig. 5E and Table S5). Previous studies have shown that myeloid-derived suppressor cells (MDSCs) are closely associated with the suppression of T cell function [44, 55]. MDSCs are divided into G-MDSCs and M-MDSCs, which are derived from neutrophils and monocytes, respectively [56].
We analyzed the correlation between different neutrophil and macrophage subtypes and the infiltration of exhausted T cells in the livers of metastatic patients with colorectal cancer. Interestingly, compared to other neutrophil subtypes, only N1 neutrophils were correlated with the proportion of exhausted T cells (Fig. 5F, Fig. S3A, S3B, and S3C). Patients with a higher proportion of N1 neutrophils had more infiltration of CD8-exhausted T cells (Fig. 5F). This suggests that N1 neutrophils may promote T cell exhaustion. Additionally, patients with a higher proportion of Kupffer2 also had more infiltrated CD8-exhausted T cells (Fig. 5G, Fig. S3D, S3E, S3F, S3G and S3H).
Furthermore, we examined the relationship between neutrophils and macrophages. We found that the proportion of N1 neutrophils was unrelated to the proportion of Kupffer2 (Fig. 5H, Fig. S3I, S3J, S3K, and S3L), but positively correlated with the proportion of MAC2 (Fig. 5I).
Diverse interactions between parenchymal cells and immune cells contribute to the formation of an immunosuppressive microenvironment
Cell–cell interactions play key roles in the formation of PMN and MN. To study the specific interactions between cells, we used CellChat to analyze the membrane protein ligand-receptor pairs in the PMN and MM. Overall, the analysis revealed that compared to the control, the number and intensity of ligand-receptor interactions between cells significantly increased in PMN and MN (Fig. 6A). The changes were more pronounced in the MN, where 787 significantly different ligand-receptor pairs were observed compared to PMN (Fig. 6A).
Cell-to-cell interactions. A The total number and proportion of different receptor-ligand interactions between cells in the control, PMN, and MN groups. B Number of receptor-ligand interactions among hepatocytes, macrophages, neutrophils, NK cells, and exhausted T cells in the control, PMN, and MN groups. C Bubble plot showing the upregulated receptor-ligand interactions in the MN group. D Expression levels of Cd274 in neutrophils, hepatocytes, and macrophages, and Cd6 and Pdcd1 in CD8 exhausted T cells in the control, PMN, and MN groups. E Upregulated receptor-ligand interactions between neutrophils and macrophages in the MN group. F Expression levels of App in neutrophils across different groups. G Expression levels of Cd74 in macrophages across different groups
We performed cell–cell interaction analysis between hepatocytes, macrophages, neutrophils, and CD8-exhausted T cells. The statistical analysis of meaningful ligand-receptor pairs indicated that hepatocytes, macrophages, and neutrophils all sent more signals to CD8-exhausted T cells in the MN (Fig. 6B). The Cd274/Pdcd1 and Alcam/Cd6 signaling axes were the main pathways increased in the MN (Fig. 6C). Specifically, Cd274 expression significantly increased in MN neutrophils, macrophages, and hepatocytes, while Pdcd1 expression was notably higher in MN exhausted T cells (Fig. 6D).
Interestingly, we also found that neutrophils in the MN might further suppress macrophage phagocytic function through the App/Cd74 signaling axis (Fig. 6E, F, G). In summary, in the MN of the liver, hepatocytes, neutrophils, and macrophages may all contribute to T cell exhaustion, leading to the formation of an immunosuppressive microenvironment.
In addition, we separately analyzed the interactions between different subtypes of neutrophils and macrophages with CD8-exhausted T cells. We found that most neutrophil subtypes highly expressed Cd274 in the MN, including N1 neutrophils (Fig. S4A, and S4B), with progressive increases in both the PMN and MN (Fig. S4C). Furthermore, compared to other types of neutrophils, N1 neutrophils exhibited higher expression of the App gene, interacting more significantly with Kupffer2 and MAC2 (Fig. S4D).
Among macrophages, MAC2 expressed higher levels of Cd74 (Fig. S4E), with progressive increases in both the PMN and MN (Fig. S4F), while Kupffer2 expressed higher levels of Cd274 (Fig. S4E).
Discussion
Tumor metastasis is a complex, multi-step process, and the tumor microenvironment plays a critical role in determining metastasis [57, 58]. However, most previous studies have focused on changes within the tumor microenvironment lacked temporal resolution in their analysis. Liver metastasis is one of the primary causes of death in colorectal cancer (CRC) patients [59]. Research has indicated that patients with liver metastases from colorectal cancer have a poor response to immunotherapy [60]. To further elucidate changes in the liver microenvironment at different stages of metastasis, and to analyze the gene changes associated with the poor response to immunotherapy in liver metastases, we used the CT26-GFP cell line to establish a CRC orthotopic model and collected liver samples for single-cell sequencing before and after tumor metastasis. These sequencing data revealed many metastasis-related gene sets and changes in immune and liver parenchymal cells under pre-metastatic niche (PMN) and metastatic niche (MN) conditions, providing some clues as to why immunotherapy is less effective in patients with CRC liver metastasis.
We identified a subpopulation of neutrophils characterized by elevated Prok2 expression. These cells exhibited increased infiltration in the liver prior to and further after metastasis (Fig. 3C). Notably, the enhanced presence of Prok2+ neutrophils in the liver correlated with poor prognosis in patients with colorectal cancer (CRC) liver metastasis (Fig. 3L). Prokineticin 2 (Prok2), an inflammatory cytokine-like molecule predominantly produced by macrophages and neutrophils infiltrating damaged tissues [61, 62], plays a pivotal role in neovascularization associated with various cancers. Research has demonstrated that combining anti-Bv8/PROK2 with anti-VEGF therapies effectively inhibits tumor formation in genetically engineered mice with cancer stem cell characteristics (CIC). Furthermore, a PK2/Bv8/PROK2 antagonist has been shown to suppress tumorigenic processes [63]. This evidence underscores the potential value of targeting Prok2+ neutrophils as a therapeutic approach to improve outcomes for patients suffering from CRC liver metastases.
To investigate this further, we developed an N1 signature based on gene expression profiles of N1 neutrophils and their association with patient prognosis. Given the clinical difficulty of obtaining liver tissue from CRC patients before metastasis, as well as the unclear origins of N1 neutrophils, we analyzed sequencing data from the peripheral blood neutrophils of both CRC liver metastasis patients and healthy individuals. Our findings revealed a significant elevation of the N1 signature in the peripheral blood neutrophils of CRC liver metastasis patients compared to healthy controls. However, the question of whether the N1 signature is elevated in peripheral blood neutrophils before CRC metastasis remains to be experimentally validated. Moving forward, our goal is to monitor early-stage CRC patients to explore the correlation between the N1 signature in peripheral blood neutrophils and the development of CRC liver metastasis. This research may highlight Prok2 as a promising target for intervention in CRC treatment strategies.
Prok2+ neutrophils interact with other cell types as well. These Prok2+ neutrophils also exhibit high expression of Cd274 and App (Fig. 6C), which may lead to T cell exhaustion and inhibit macrophage phagocytosis, respectively. Consistent with this, in samples from human colorectal cancer liver metastases, the proportion of Prok2+ neutrophils is positively correlated with exhausted T cells. However, there is no correlation between the proportion of Prok2+ neutrophils and Kupffer2 cells, which may be related to the characteristic of Kupffer2 cells as resident macrophages in the liver. In most cases, Kupffer2 cells proliferate and increase in number due to stimuli such as tumors and inflammation, rather than being recruited [64, 65]. Although studies have shown that Prok2 can promote the chemotaxis of macrophages [66, 67], the macrophages referred to here are more likely to be monocyte-derived macrophages rather than Kupffer cells. Consistent with this, in our study, the proportion of N1 neutrophils was positively correlated with MAC2 cells (Fig. 5I).
We also observed gene changes in resident liver cells during the progression of CRC liver metastasis, including hepatocytes and Kupffer cells. Before liver metastasis in CRC, hepatocytes expressed large amounts of inflammatory proteins such as Saa1 and the neutrophil chemokine Cxcl1. The increased secretion of Cxcl1 is likely one of the factors driving the infiltration of myeloid cells in the liver during the PMN phase of colorectal cancer liver metastasis. Studies have shown that Saa1 secreted by the liver can act on T cells in pancreatic cancer, promoting T cell exhaustion. The poor response to immunotherapy in CRC liver metastasis patients may be closely related to the cytokines secreted by the liver. Interestingly, we found that Saa1 did not further increase in MN but instead showed a decline. One possible reason is that these inflammation-related proteins belong to the acute-phase proteins and mainly participate in PMN formation, while in the later MN stage, hepatocytes shift to primarily secrete inflammatory cytokines [68,69,70].
Tumor cells also have a significant impact on resident Kupffer cells in the liver. In the MN liver, the TNF/NF-kB pathway in Kupffer2 cells was significantly activated, which may further promote metastasis. This may be closely related to the inflammatory stimuli from the tumor. Research has shown that following liver resection or viral infection, Kupffer2 cells and hepatocytes rapidly produce IL-6, promoting liver regeneration [71, 72]. Additionally, IL-6 promotes T cell exhaustion by facilitating the expression of immune checkpoint receptors (IR) in CD8 + T cells [73, 74]. At the same time, we observed similar results in human specimens. The infiltration of Kupffer2 cells in liver metastatic lesions of colorectal cancer patients was positively correlated with the infiltration of CD8 exhausted T cells. Immunotherapies such as PD-L1 inhibitors primarily target early-stage stem-like exhausted T cells [75]. Once T cells enter the terminal exhaustion phase, they no longer respond to PD-L1 inhibitors [76]. We found that after CRC liver metastasis, various cells, including neutrophils, macrophages, and hepatocytes, highly expressed Cd274, promoting further T cell exhaustion. These are likely to be one of the factors contributing to the lack of response to immunotherapy in patients with colorectal cancer liver metastases.
Conclusions
Our study highlights the dynamic changes in the liver microenvironment during colorectal cancer (CRC) liver metastasis, providing new insights into the mechanisms underlying the poor response to immunotherapy in metastatic cases. We identified key immune cell subpopulations, including N1 neutrophils (Prok2 + neutrophils) and Kupffer cells, that play significant roles in promoting immunosuppression and T cell exhaustion within the metastatic niche (MN). The increased infiltration of neutrophils and macrophages, alongside their high expression of Cd274 and other immunosuppressive markers, contributes to the creation of a tumor-supportive environment, hindering the effectiveness of immunotherapy. Furthermore, we established an N1 signature that correlates with poor prognosis in CRC liver metastasis patients and demonstrated that these findings could be extended to peripheral blood neutrophils, offering a potential biomarker for early detection and monitoring of metastasis. Overall, our study not only deepens the understanding of cellular interactions and gene expression changes in the liver during CRC metastasis but also identifies novel targets for therapeutic intervention and potential biomarkers to improve the treatment of CRC liver metastasis.
Availability of data and materials
The accession number for the processed data from this paper is GEO: GSE284449.
Abbreviations
- PMN:
-
Pre-metastatic niche
- MN:
-
Metastatic niche
- CRC:
-
Colorectal cancer
- CRCLM:
-
Colorectal cancer liver metastases
- MDSC:
-
Myeloid-derived suppressor cells
- GFP:
-
Green fluorescent protein
- FACS:
-
Fluorescence-activated cell sorting
- UMAP:
-
Uniform manifold approximation and projection
- PDAC:
-
Pancreatic ductal adenocarcinoma
- KCs:
-
Kupffer cells
- KEGG:
-
Kyoto encyclopedia of genes and genomes
References
Rawla P, Sunkara T, Barsouk A. Epidemiology of colorectal cancer: incidence, mortality, survival, and risk factors. Prz Gastroenterol. 2019;14:89–103.
Xi Y, Xu P. Global colorectal cancer burden in 2020 and projections to 2040. Transl Oncol. 2021;14:101174.
AlZaabi A, AlHarrasi A, AlMusalami A, AlMahyijari N, Al Hinai K, ALAdawi H, Al-Shamsi HO. Early onset colorectal cancer: Challenges across the cancer care continuum. Ann Med Surg. 2022;82:104453.
Cervantes A, Adam R, Roselló S, Arnold D, Normanno N, Taïeb J, Seligmann J, De Baere T, Osterlund P, Yoshino T, Martinelli E. Metastatic colorectal cancer: ESMO Clinical Practice Guideline for diagnosis, treatment and follow-up. Ann Oncol. 2023;34:10–32.
Jin B, Wu X, Xu G, Xing J, Wang Y, Yang H, Du S, Mao Y. Evolutions of the management of colorectal cancer liver metastasis: a bibliometric analysis. J Cancer. 2021;12:3660–70.
Patras L, Shaashua L, Matei I, Lyden D. Immune determinants of the pre-metastatic niche. Cancer Cell. 2023;41:546–72.
Liu H, Zhang G, Gao R. Cellular and molecular characteristics of the premetastatic niches. Animal Model Exp Med. 2023;6:399–408.
Wang Y, Jia J, Wang F, Fang Y, Yang Y, Zhou Q, Yuan W, Gu X, Hu J, Yang S. Pre-metastatic niche: formation, characteristics and therapeutic implication. Signal Transduct Target Ther. 2024;9:236.
Li W, Yang Y, Yang L, Chang N, Li L. Monocyte-derived Kupffer cells dominate in the Kupffer cell pool during liver injury. Cell Rep. 2023;42:113164.
Zeng D, Wang M, Wu J, Lin S, Ye Z, Zhou R, Wang G, Wu J, Sun H, Bin J, et al. Immunosuppressive microenvironment revealed by immune cell landscape in pre-metastatic liver of colorectal cancer. Front Oncol. 2021;11:620688.
Xu Y, Yao Y, Yu L, Zhang X, Mao X, Tey SK, Wong SWK, Yeung CLS, Ng TH, Wong MY, et al. Clathrin light chain A-enriched small extracellular vesicles remodel microvascular niche to induce hepatocellular carcinoma metastasis. J Extracell Vesicles. 2023;12:e12359.
Grünwald B, Harant V, Schaten S, Frühschütz M, Spallek R, Höchst B, Stutzer K, Berchtold S, Erkan M, Prokopchuk O, et al. Pancreatic premalignant lesions secrete tissue inhibitor of metalloproteinases-1, which activates hepatic stellate cells via CD63 signaling to create a premetastatic niche in the liver. Gastroenterology. 2016;151:1011-1024.e1017.
Stone ML, Lee J, Lee JW, Coho H, Tariveranmoshabad M, Wattenberg MM, Choi H, Herrera VM, Xue Y, Choi-Bose S, et al. Hepatocytes coordinate immune evasion in cancer via release of serum amyloid A proteins. Nat Immunol. 2024;25:755–63.
Lin Q, Ren L, Jian M, Xu P, Li J, Zheng P, Feng Q, Yang L, Ji M, Wei Y, Xu J. The mechanism of the premetastatic niche facilitating colorectal cancer liver metastasis generated from myeloid-derived suppressor cells induced by the S1PR1–STAT3 signaling pathway. Cell Death Dis. 2019;10:693.
Chen XW, Yu TJ, Zhang J, Li Y, Chen HL, Yang GF, Yu W, Liu YZ, Liu XX, Duan CF, et al. CYP4A in tumor-associated macrophages promotes pre-metastatic niche formation and metastasis. Oncogene. 2017;36:5045–57.
Sasaki S, Baba T, Muranaka H, Tanabe Y, Takahashi C, Matsugo S, Mukaida N. Involvement of Prokineticin 2-expressing neutrophil infiltration in 5-fluorouracil-induced aggravation of breast cancer metastasis to lung. Mol Cancer Ther. 2018;17:1515–25.
Meng Y, Ye F, Nie P, Zhao Q, An L, Wang W, Qu S, Shen Z, Cao Z, Zhang X, et al. Immunosuppressive CD10(+)ALPL(+) neutrophils promote resistance to anti-PD-1 therapy in HCC by mediating irreversible exhaustion of T cells. J Hepatol. 2023;79:1435–49.
Hao Y, Stuart T, Kowalski MH, Choudhary S, Hoffman P, Hartman A, Srivastava A, Molla G, Madad S, Fernandez-Granda C, Satija R. Dictionary learning for integrative, multimodal and scalable single-cell analysis. Nat Biotechnol. 2024;42:293–304.
McGinnis CS, Murrow LM, Gartner ZJ. DoubletFinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Syst. 2019;8:329-337.e324.
Korsunsky I, Millard N, Fan J, Slowikowski K, Zhang F, Wei K, Baglaenko Y, Brenner M, Loh P-R, Raychaudhuri S. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods. 2019;16:1289–96.
Jew B, Alvarez M, Rahmani E, Miao Z, Ko A, Garske KM, Sul JH, Pietiläinen KH, Pajukanta P, Halperin E. Accurate estimation of cell composition in bulk expression through robust integration of single-cell information. Nat Commun. 1971;2020:11.
Lee JW, Stone ML, Porrett PM, Thomas SK, Komar CA, Li JH, Delman D, Graham K, Gladney WL, Hua X, et al. Hepatocytes direct the formation of a pro-metastatic niche in the liver. Nature. 2019;567:249–52.
Li Q, Zhang XX, Hu LP, Ni B, Li DX, Wang X, Jiang SH, Li H, Yang MW, Jiang YS, et al. Coadaptation fostered by the SLIT2-ROBO1 axis facilitates liver metastasis of pancreatic ductal adenocarcinoma. Nat Commun. 2023;14:861.
Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, Feng T, Zhou L, Tang W, Zhan L, et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation. 2021;2:100141.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43:e47.
Kassambara A. Practical guide to cluster analysis in R: Unsupervised machine learning. Sthda; 2017.
Moosavi SH, Eide PW, Eilertsen IA, Brunsell TH, Berg KCG, Røsok BI, Brudvik KW, Bjørnbeth BA, Guren MG, Nesbakken A, et al. De novo transcriptomic subtyping of colorectal cancer liver metastases in the context of tumor heterogeneity. Genome Med. 2021;13:143.
Therneau T. A package for survival analysis in S. 2015, 2:2014.
Furze RC, Rankin SM. Neutrophil mobilization and clearance in the bone marrow. Immunology. 2008;125:281–8.
Overbeeke C, Tak T, Koenderman L. The journey of neutropoiesis: how complex landscapes in bone marrow guide continuous neutrophil lineage determination. Blood. 2022;139:2285–93.
Shaul ME, Fridlender ZG. Tumour-associated neutrophils in patients with cancer. Nat Rev Clin Oncol. 2019;16:601–20.
Ma R, Zhou X, Zhai X, Wang C, Hu R, Chen Y, Shi L, Fang X, Liao Y, Ma L, et al. Single-cell RNA sequencing reveals immune cell dysfunction in the peripheral blood of patients with highly aggressive gastric cancer. Cell Prolif. 2024;57:e13591.
Han X, Zhou Z, Fei L, Sun H, Wang R, Chen Y, Chen H, Wang J, Tang H, Ge W, et al. Construction of a human cell landscape at single-cell level. Nature. 2020;581:303–9.
Wu Y, Yang S, Ma J, Chen Z, Song G, Rao D, Cheng Y, Huang S, Liu Y, Jiang S, et al. Spatiotemporal immune landscape of colorectal cancer liver metastasis at single-cell level. Cancer Discov. 2022;12:134–53.
Trapnell C, Cacchiarelli D, Grimsby J, Pokharel P, Li S, Morse M, Lennon NJ, Livak KJ, Mikkelsen TS, Rinn JL. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol. 2014;32:381–6.
Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan C-H, Myung P, Plikus MV, Nie Q. Inference and analysis of cell-cell communication using Cell Chat. Nat Commun. 2021;12:1088.
Guilliams M, Bonnardel J, Haest B, Vanderborght B, Wagner C, Remmerie A, Bujko A, Martens L, Thoné T, Browaeys R, et al. Spatial proteogenomics reveals distinct and evolutionarily conserved hepatic macrophage niches. Cell. 2022;185:379-396.e338.
Zheng L, Qin S, Si W, Wang A, Xing B, Gao R, Ren X, Wang L, Wu X, Zhang J, et al. Pan-cancer single-cell landscape of tumor-infiltrating T cells. Science. 2021;374:6474.
Zhang L, Li Z, Skrzypczynska KM, Fang Q, Zhang W, O’Brien SA, He Y, Wang L, Zhang Q, Kim A, et al. Single-cell analyses inform mechanisms of myeloid-targeted therapies in colon cancer. Cell. 2020;181:442-459.e429.
Namineni S, O’Connor T, Faure-Dupuy S, Johansen P, Riedl T, Liu K, Xu H, Singh I, Shinde P, Li F, et al. A dual role for hepatocyte-intrinsic canonical NF-κB signaling in virus control. J Hepatol. 2020;72:960–75.
Meng Z, Fu X, Chen X, Zeng S, Tian Y, Jove R, Xu R, Huang W. miR-194 is a marker of hepatic epithelial cells and suppresses metastasis of liver cancer cells in mice. Hepatology. 2010;52:2148–57.
Si-Tayeb K, Lemaigre FP, Duncan SA. Organogenesis and development of the liver. Dev Cell. 2010;18:175–89.
Herranz-Itúrbide M, López-Luque J, Gonzalez-Sanchez E, Caballero-Díaz D, Crosas-Molist E, Martín-Mur B, Gut M, Esteve-Codina A, Jaquet V, Jiang JX, et al. NADPH oxidase 4 (Nox4) deletion accelerates liver regeneration in mice. Redox Biol. 2021;40:101841.
Veglia F, Sanseviero E, Gabrilovich DI. Myeloid-derived suppressor cells in the era of increasing myeloid cell diversity. Nat Rev Immunol. 2021;21:485–98.
Talmadge JE, Gabrilovich DI. History of myeloid-derived suppressor cells. Nat Rev Cancer. 2013;13:739–52.
Grieshaber-Bouyer R, Radtke FA, Cunin P, Stifano G, Levescot A, Vijaykumar B, Nelson-Maney N, Blaustein RB, Monach PA, Nigrovic PA. The neutrotime transcriptional signature defines a single continuum of neutrophils across biological compartments. Nat Commun. 2021;12:2856.
Shojaei F, Zhong C, Wu X, Yu L, Ferrara N. Role of myeloid cells in tumor angiogenesis and growth. Trends Cell Biol. 2008;18:372–8.
Morrissey SM, Zhang F, Ding C, Montoya-Durango DE, Hu X, Yang C, Wang Z, Yuan F, Fox M, Zhang H-G, et al. Tumor-derived exosomes drive immunosuppressive macrophages in a pre-metastatic niche through glycolytic dominant metabolic reprogramming. Cell Metab. 2021;33:2040–58.
Woltman AM, Boonstra A, Naito M, Leenen PJM. Kupffer cells in health and disease. In: Biswas SK, Mantovani A, editors. Macrophages: biology and role in the pathology of diseases. New York: Springer; 2014. p. 217–47.
Robinson MW, Harmon C, O’Farrelly C. Liver immunology and its role in inflammation and homeostasis. Cell Mol Immunol. 2016;13:267–76.
McGinnis CS, Miao Z, Superville D, Yao W, Goga A, Reticker-Flynn NE, Winkler J, Satpathy AT. The temporal progression of lung immune remodeling during breast cancer metastasis. Cancer Cell. 2024;42:1018-1031.e1016.
Nutt SL, Huntington ND. 17 - Cytotoxic T lymphocytes and natural killer cells. In: Rich RR, Fleisher TA, Shearer WT, Schroeder HW, Frew AJ, Weyand CM, editors. Clinical immunology (Fifth Edition). London: Elsevier; 2019. p. 247- 259.e241.
Coënon L, Geindreau M, Ghiringhelli F, Villalba M, Bruchard M. Natural Killer cells at the frontline in the fight against cancer. Cell Death Dis. 2024;15:614.
Yang YL, Yang F, Huang ZQ, Li YY, Shi HY, Sun Q, Ma Y, Wang Y, Zhang Y, Yang S, et al. T cells, NK cells, and tumor-associated macrophages in cancer immunotherapy and the current state of the art of drug delivery systems. Front Immunol. 2023;14:1199173.
Zhang H, Li ZL, Ye SB, Ouyang LY, Chen YS, He J, Huang HQ, Zeng YX, Zhang XS, Li J. Myeloid-derived suppressor cells inhibit T cell proliferation in human extranodal NK/T cell lymphoma: a novel prognostic indicator. Cancer Immunol Immunother. 2015;64:1587–99.
Sendo S, Saegusa J, Morinobu A. Myeloid-derived suppressor cells in non-neoplastic inflamed organs. Inflamm Regener. 2018;38:19.
Shi X, Wang X, Yao W, Shi D, Shao X, Lu Z, Chai Y, Song J, Tang W, Wang X. Mechanism insights and therapeutic intervention of tumor metastasis: latest developments and perspectives. Signal Transduct Target Ther. 2024;9:192.
Liu M, Yang J, Xu B, Zhang X. Tumor metastasis: Mechanistic insights and therapeutic interventions. MedComm. 2020;2021(2):587–617.
Ki DH, Jeung HC, Park CH, Kang SH, Lee GY, Lee WS, Kim NK, Chung HC, Rha SY. Whole genome analysis for liver metastasis gene signatures in colorectal cancer. Int J Cancer. 2007;121:2005–12.
Wang C, Sandhu J, Ouyang C, Ye J, Lee PP, Fakih M. Clinical response to immunotherapy targeting programmed cell death receptor 1/programmed cell death ligand 1 in patients with treatment-resistant microsatellite stable colorectal cancer with and without liver metastases. JAMA Netw Open. 2021;4:e2118416.
Watson RP, Lilley E, Panesar M, Bhalay G, Langridge S, Tian SS, McClenaghan C, Ropenga A, Zeng F, Nash MS. Increased prokineticin 2 expression in gut inflammation: role in visceral pain and intestinal ion transport. Neurogastroenterol Motil. 2012;24:65–75.
Lattanzi R, Maftei D, Marconi V, Florenzano F, Franchi S, Borsani E, Rodella LF, Balboni G, Salvadori S, Sacerdote P, Negri L. Prokineticin 2 upregulation in the peripheral nervous system has a major role in triggering and maintaining neuropathic pain in the chronic constriction injury model. Biomed Res Int. 2015;2015:301292.
Curtis VF, Wang H, Yang P, McLendon RE, Li X, Zhou QY, Wang XF. A PK2/Bv8/PROK2 antagonist suppresses tumorigenic processes by inhibiting angiogenesis in glioma and blocking myeloid cell infiltration in pancreatic cancer. PLoS ONE. 2013;8:e54916.
Nguyen-Lefebvre AT, Horuzsko A. Kupffer cell metabolism and function. J Enzymol Metab. 2015;1:101.
Dixon LJ, Barnes M, Tang H, Pritchard MT, Nagy LE. Kupffer cells in the liver. Compr Physiol. 2013;3:785–97.
Tu Q, Yu X, Xie W, Luo Y, Tang H, Chen K, Ruan Y, Li Y, Zhou J, Yin Y, et al. Prokineticin 2 promotes macrophages-mediated antibacterial host defense against bacterial pneumonia. Int J Infect Dis. 2022;125:103–13.
Martin C, Balasubramanian R, Dwyer AA, Au MG, Sidis Y, Kaiser UB, Seminara SB, Pitteloud N, Zhou QY, Crowley WF Jr. The role of the prokineticin 2 pathway in human reproduction: evidence from the study of human and murine gene mutations. Endocr Rev. 2011;32:225–46.
Xu Y, Yamada T, Satoh T, Okuda Y. Measurement of serum amyloid A1 (SAA1), a major isotype of acute phase SAA. Clin Chem Lab Med. 2006;44:59–63.
Sack GH. Serum amyloid A – a review. Mol Med. 2018;24:46.
Thorn CF, Lu ZY, Whitehead AS. Tissue-specific regulation of the human acute-phase serum amyloid A genes, SAA1 and SAA2, by glucocorticoids in hepatic and epithelial cells. Eur J Immunol. 2003;33:2630–9.
Wang M-J, Zhang H-L, Chen F, Guo X-J, Liu Q-G, Hou J. The double-edged effects of IL-6 in liver regeneration, aging, inflammation, and diseases. Exp Hematol Oncol. 2024;13:62.
Naseem S, Hussain T, Manzoor S. Interleukin-6: A promising cytokine to support liver regeneration and adaptive immunity in liver pathologies. Cytokine Growth Factor Rev. 2018;39:36–45.
Huseni MA, Wang L, Klementowicz JE, Yuen K, Breart B, Orr C, Liu LF, Li Y, Gupta V, Li C, et al. CD8(+) T cell-intrinsic IL-6 signaling promotes resistance to anti-PD-L1 immunotherapy. Cell Rep Med. 2023;4:100878.
Soler MF, Abaurrea A, Azcoaga P, Araujo AM, Caffarel MM. New perspectives in cancer immunotherapy: targeting IL-6 cytokine family. J Immunother Cancer. 2023;11:e007530.
Lee JB, Kim HR, Ha S-J. Immune checkpoint inhibitors in 10 years: contribution of basic research and clinical application in cancer immunotherapy. Immune Netw. 2022. https://doiorg.publicaciones.saludcastillayleon.es/10.4110/in.2022.22.e2.
Niborski LL, Gueguen P, Ye M, Thiolat A, Ramos RN, Caudana P, Denizeau J, Colombeau L, Rodriguez R, Goudot C, et al. CD8+T cell responsiveness to anti-PD-1 is epigenetically regulated by Suv39h1 in melanomas. Nat Commun. 2022;13:3739.
Acknowledgements
Not applicable.
Funding
This work was supported by the National Natural Science Foundation of China (82273359, 82473022, 82003163), Guangzhou Science and Technology Projects (2023B01J1004) and the Fundamental Research Funds for the Central Universities, Sun Yat-Sen University (23ptpy145).
Author information
Authors and Affiliations
Contributions
Weidong Pan, Bing Cheng, and Wenyu Wang conceptualized and designed the study. Yue Jiang, Guojie Long, and Xiaoming Huang conducted the bioinformatics analysis and animal experiments and drafted the manuscript.
Corresponding authors
Ethics declarations
Ethics approval and consent to participate
All experimental animals were obtained from the Animal Center of the Sixth Affiliated Hospital of Sun Yat-sen University (Guangzhou, China) and handled in strict accordance with protocols approved by the Institutional Animal Care and Use Committee of the same institution. The maximum diameter of all tumors is less than 1.5 cm, in compliance with the ethical standards of the animal center.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
12967_2025_6351_MOESM1_ESM.pdf
Supplementary material 1: Figure S1. Orthotopic colorectal cancer model, related to Figure 1.Cecum and liver samples from mice at different time points.HE sections of mouse cecal tumor. PMN refers to the pre-metastatic niche, while MN refers to the metastatic niche
12967_2025_6351_MOESM2_ESM.pdf
Supplementary material 2: Figure S2. Characteristics of different neutrophil subtypes, related to Figure 3.Heatmap of top genes in N0-N3 neutrophil subpopulations.Uniform Manifold Approximation and Projectionplot of Retnlg expression in different neutrophil subtypes.UMAP plot of Mmp8 expression in different neutrophil subtypes.UMAP plot of Ltf expression in different neutrophil subtypes.UMAP plot of Camp expression in different neutrophil subtypes.UMAP plot of Prok2 expression in different neutrophil subtypes.UMAP plot of neutrophil pseudotime.Survival curve of prognosis based on high and low N0 neutrophil proportions in colorectal cancer patients.Survival curve of prognosis based on high and low N2 neutrophil proportions in colorectal cancer patients.Survival curve of prognosis based on high and low N3 neutrophil proportions in colorectal cancer patients.Differentially expressed genes between N1 neutrophils and other cells.Survival curve of prognosis based on high and low N1 signature in colorectal cancer patients.
12967_2025_6351_MOESM3_ESM.pdf
Supplementary material 3: Figure S3. Correlation between the proportions of neutrophils, macrophages, and exhausted T cells, related to Figure 4.Proportion of infiltrating exhausted CD8+ T cells in high- and low-infiltration N0, N2, and N3neutrophil colorectal cancer patients with liver metastasis.Proportion of infiltrating exhausted CD8+ T cells in high- and low-infiltration Kupffer1, Kupffer3, MAC1, MAC2and MAC3colorectal cancer patients with liver metastasis.Proportion of infiltrating N1 neutrophil in high- and low-infiltration Kupffer1, Kupffer2, MAC1, and MAC3colorectal cancer patients with liver metastasis
12967_2025_6351_MOESM4_ESM.pdf
Supplementary material 4: Figure S4. Interactions between neutrophil subpopulations, macrophage subpopulations, and exhausted T cells, related to Figure 6.N0-N3 neutrophils, Kupffer1-Kupffer3, and MAC1-MAC3 with exhausted T cells: receptor-ligand pairs elevated in MN.Expression levels of Cd274 and App in N0-N3 neutrophils.Expression levels of Cd274 and App in N1 neutrophils at different time points.N0-N3 neutrophils with Kupffer1-Kupffer3, and MAC1-MAC3: receptor-ligand pairs elevated in MN.Expression levels of Cd274 and Cd74 in Kupffer1-Kupffer3, and MAC1-MAC3.Expression levels of Cd274 and Cd74 in MAC2 at different time points. PMN refers to the pre-metastatic niche, while MN refers to the micro-metastatic niche.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
About this article
Cite this article
Jiang, Y., Long, G., Huang, X. et al. Single-cell transcriptomic analysis reveals dynamic changes in the liver microenvironment during colorectal cancer metastatic progression. J Transl Med 23, 336 (2025). https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12967-025-06351-3
Received:
Accepted:
Published:
DOI: https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12967-025-06351-3