Baseline characteristics

A total of 8035 and 1305 CRLM patients in the SEER cohort and Chinese cohort, respectively, who underwent surgical resection were included in this study. The baseline characteristics are shown in Supplementary Table 1. We conducted a comparative analysis of baseline characteristics between the LNM group and Non-LNM group, both before and after IPTW, in these two cohorts (Supplementary Table 2). After applying IPTW, the Non-LNM group and the LNM group were similar both in the SEER cohort and the Chinese cohort. This process ensured a well-balanced distribution of baseline characteristics between the two groups, including age at diagnosis, sex, race, primary tumor site, tumor histology type, postoperative radiotherapy and/or chemotherapy, surgical procedure, and tumor size (all p > 0.05).

LNM leads to worse CSS in patients with CRLM

To further explore the prognostic effect of LNM in patients with CRLM, we conducted Kaplan‒Meier survival analysis. Non-LNM remained associated with a better prognosis even after adjusting for other clinicopathological prognostic factors in the SEER cohort (HR = 1.72; 95% CI, 1.59–1.86; p < 0.001), making it an independent prognostic factor for stage IV CRC patients (Supplementary Table 3). The same conclusion was validated in the Chinese cohort (HR = 1.46; 95% CI, 1.24–1.73; p < 0.001). Both before and after IPTW, Kaplan‒Meier survival analysis revealed that patients with LNM had significantly worse CSS (all log-rank p < 0.001) than patients without LNM in the SEER cohort (Fig. 2a, b). In the Chinese cohort, patients with LNM also had significantly worse CSS (all log-rank p < 0.001) than patients without LNM (Fig. 2c, d).

Fig. 2: Survival comparisons between CRLM patients with and without LNM.

Kaplan‒Meier curve of stratified survival in patients in the SEER cohort: a before IPTW, b after IPTW. Kaplan‒Meier curve of stratified survival in patients in the Chinese cohort: c before IPTW, d after IPTW. CRLM colorectal cancer liver metastasis, LNM lymph node metastasis, SEER Surveillance, Epidemiology, and End Results, IPTW inverse probability of treatment weighting.

Gene expression characteristics of metastatic and nonmetastatic TDLNs

The gene expression levels of TDLNs were compared between the metastatic and nonmetastatic TDLN groups. The results revealed that 3116 genes were upregulated, while 1870 genes were downregulated in the metastatic TDLN group (Supplementary Fig. 2a). Among them, the epithelial cell differentiation-related KRT gene family, the cell cycle regulation gene FAM83H, and the tumor metastasis-related genes EPS8L1 and TNS4 were upregulated in metastatic TDLNs. Conversely, immune regulation-related genes such as CLEC4M and CD5L were downregulated in these lymph nodes.

To further characterize the biological attributes distinguishing the metastatic and nonmetastatic TDLN groups, we employed Metascape enrichment analysis. In the metastatic TDLN group, the pathways enriched in the upregulated genes were associated with stromal cell activation, including cell junction and adhesion, fibroblast-mediated extracellular matrix organization, endothelial cell-mediated blood vessel development and epithelial cell differentiation (Supplementary Fig. 2B). The downregulated genes were primarily enriched in pathways related to adaptive immunity, including the regulation of lymphocyte activation, differentiation, immune effector processes and antigen receptor-mediated signaling (Supplementary Fig. 2C).

Single-cell RNA sequencing reveals differences in the TDLN microenvironment

After preprocessing and quality control, our single-cell atlas contained high-quality transcriptomes from 84623 cells spanning both tumor and immune compartments. Unsupervised clustering analysis revealed 18 distinct clusters, including all major known epithelial, mesenchymal, and leukocyte lineages (Fig. 3a). Uniform manifold approximation and projection (UMAP) plots based on the sample of origin further revealed the impact of tumor metastases on the cellular composition of TDLNs (Supplementary Fig. 3A). Compared to nonmetastatic TDLNs, metastatic TDLNs showed a substantial presence of epithelial cells with similar characteristics to those of the primary tumor, consistent with the typical definition of LNM. For accurate cell type determination, we only used highly variable genes between cell types (Supplementary Fig. 3B). We then compared the microenvironmental composition at distinct sites of origin (Fig. 3b). In general, lymphocytes, including CD4+ T cells, CD8+ T cells, and B cells, represented the largest immune subpopulation within nonmetastatic TDLNs, while metastatic TDLNs contained a relatively greater proportion of fibroblasts, endothelial cells and macrophages.

Fig. 3: The immune infiltration profile of TDLNs.

a Uniform manifold approximation and projection (UMAP) analysis of the transcriptional profiles of TDLNs, colored by cell type. b Proportion plots of various cell types from different sample sources. c Box plots illustrating the differences in immune cell infiltration between negative and positive TDLNs using lymph node sc-RNA sequencing data as a reference based on CIBERSORTx. d Heatmap and box plots showing cell types highly expressed DEGs by the expression of upregulated DEGs (PosLN vs. NegLN) across cell types. NegLN negative lymph node, PosLN positive lymph node, TDLNs tumor-draining lymph nodes, sc-RNA single-cell RNA, DEGs differentially expressed genes.

To further validate the results of single-cell RNA sequencing, we conducted cell composition prediction using bulk RNA sequencing data from the aforementioned 173 TDLNs. Unlike traditional immune infiltration analyses, we conducted cell composition prediction using lymph node single-cell RNA sequencing data as a reference based on CIBERSORTx [27]. The numbers of lymphocytes, including CD4+ T cells, CD8+ T cells, and B cells, were increased in nonmetastatic TDLNs (Fig. 3c, p < 0.001, p < 0.001, and p < 0.001, respectively). In metastatic TDLNs, there was a significant increase in fibroblasts, epithelial cells, and macrophages (Supplementary Fig. 3C, p < 0.05, p < 0.001, and p < 0.001, respectively). Then, we identified DEGs that were upregulated or downregulated in each cell type by comparing samples from patients with and without LNM. The upregulated DEGs in patients with LNM were also highly expressed in epithelial cells, fibroblasts, macrophages, and endothelial cells (Supplementary Fig. 3D). The downregulated DEGs in patients with LNM were highly expressed in CD4+/CD8+ T cells and B cells (Fig. 3d).

In addition to changes in cell count and proportion, we further investigated the transcriptomic differences of B cells, CD4+ T cells, CD8+ T cells, and macrophages between metastatic and non-metastatic TDLNs (Supplementary Fig. 4A). First, we performed Gene Ontology (GO)/Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis on genes downregulated in each cell type in metastatic TDLNs compared to non-metastatic TDLNs. We found that in metastatic TDLNs, biological processes such as B cell activation and B cell-mediated immunity were downregulated, and molecular functions such as immune receptor activity were inhibited (Supplementary Fig. 4B). For both CD4+ T cells and CD8+ T cells, impairments in cell activation and differentiation were found (Supplementary Fig. 4C, D). For macrophages, the biological process of antigen processing and presentation was suppressed, leading to a reduced capacity for producing molecules that mediate immune responses and a diminished ability to regulate adaptive immunity (Supplementary Fig. 4E). Similarly, we also performed enrichment analysis on the upregulated genes, and the results are presented in Supplementary Fig. 5. In conclusion, after the metastasis of tumor cells, the transcriptome of TDLNs is altered to induce immune dysfunction and stromal cell proliferation.

The impact of LNM on the TDLNs of patients with CRLM

We first validated the above findings via proteomic sequencing. We found that the T-cell marker gene CD247 and the B-cell marker genes TCL1A and BANK1 were upregulated at the protein level, while the epithelial cell marker gene MACC1 and the macrophage marker genes FBP1, SPP1, and TFRC were downregulated (Fig. 4a and Supplementary Fig. 6A). A total of 309 proteins were differentially expressed, with 187 upregulated and 122 downregulated (Fig. 4b, c). To further explore the biological processes associated with LNM in the TDLNs of patients with CRLM, we conducted pathway enrichment analysis. KEGG enrichment analysis revealed that the downregulated proteins were mainly associated with immune-related pathways, including the T/B-cell receptor signaling pathway, helper T-cell differentiation pathway, and antigen processing and presentation pathway (Fig. 4d). GSEA also demonstrated the same results (Fig. 4e, g, Supplementary Fig. 6B, C). We also examined the expression of proteins related to these four downregulated pathways and found that they were all downregulated in metastatic TDLNs (Fig. 4f, h, Supplementary Fig. 6D). Additionally, we found that the upregulated proteins were mainly enriched in the extracellular matrix (ECM)-receptor interaction pathway, with ECM-related proteins such as COL1A1, FN1, and ITGB3 being upregulated in metastatic TDLNs (Supplementary Fig. 6E–G).

Fig. 4: Proteomic and pathological characteristics of positive and negative TDLNs.

a Violin plots verifying protein differences in immune cell-associated upregulated genes. b Heatmap of differentially expressed proteins in positive lymph nodes (PosLNs) and negative lymph nodes (NegLNs); c Volcano plot showing differentially expressed proteins in PosLNs and NegLNs; d Dot plot showing the KEGG pathways enriched in the proteins upregulated in NegLNs; e GSEA plot of the T-cell receptor signaling pathway gene set. f Violin plot showing the difference in protein expression between PosLNs and NegLNs in the T-cell receptor signaling pathway gene set. g GSEA plot of the B-cell receptor signaling pathway gene set. h Violin plot showing the difference in protein expression between PosLNs and NegLNs in the B-cell receptor signaling pathway gene set. i, j Immunohistochemical staining of CD3, CD8 and CD19 in PosLNs and NegLNs. i Pathological images. j Violin plot of the quantitative results. PosLN positive lymph node, NegLN negative lymph node, TDLNs tumor-draining lymph nodes, KEGG Kyoto Encyclopedia of Genes and Genomes, GSEA gene set enrichment analysis.

We further validated the impact of LNM on the immune function of TDLNs at the pathological level. When LNM occurred, the immune microenvironment of the TDLNs was comprehensively suppressed, leading to a decrease in the density of CD3+ T cells, CD8+ T cells, and CD19+ B cells (Fig. 4i, j, p < 0.001, p = 0.031, and p < 0.001, respectively). Morphologically, metastatic TDLNs were predominantly characterized by the lymphocyte depletion type (Supplementary Fig. 7, 36.4% vs. 9.8%, p < 0.05), while nonmetastatic TDLNs were predominantly characterized by the GC predominance type (41.2% vs. 14.5%, p < 0.05). No significant differences were observed between the LNM group and Non-LNM group in terms of lymphocyte predominance type (36.4% vs. 41.2%, p > 0.05). In summary, our proteomic and pathological analyses both revealed immune dysfunction in metastatic TDLNs.

The impact of LNM on primary tumors in patients with CRLM

Next, we validated the impact of LNM on CRLM at the level of the primary tumor. Only 52 proteins were differentially expressed in the primary tumor, with 33 upregulated and 19 downregulated DEGs (Supplementary Fig. 8A, B). GO enrichment analysis indicated that the downregulated proteins were mainly associated with immune response-related biological processes (Supplementary Fig. 8C, D). Unfortunately, we did not find any pathways enriched in differentially expressed proteins based on the KEGG enrichment analysis.

However, we did observe the impact of LNM on the primary tumor at the pathological level. Halo image analysis was applied to accurately assess the number of IF-positive cells according to Pearson correlation analysis (Fig. 5a, b, p < 0.001). Compared to that of those in the group without LNM, the density of CD3+ T cells and CD19+ B cells in the IF in the group with LNM was significantly lower (Fig. 5c–e, p = 0.018 and p = 0.017, respectively). Moreover, 30% of the primary tumors had a high KM grade in the non-LNM group, which was significantly greater than the proportion of tumors with a high KM grade in the LNM group (Fig. 5f, g, p < 0.05). The density and number of CLRs were greater in the group without LNM (Fig. 5h, i and Supplementary Fig. 8F, p = 0.006 and p = 0.038, respectively). The group without LNM had a greater proportion of patients with deficient mismatch repair (dMMR) status (Supplementary Fig. 8H, 10.7% vs. 4.1%, p < 0.05). We further investigated the prognostic impact of LNM under different MSI statuses in CRLM patients. Patients with pMMR without LNM had a better prognosis (Supplementary Fig. 9A, log-rank p < 0.0001), whereas among dMMR patients, there was no statistically significant difference in prognosis between those with and without LNM (Supplementary Fig. 9B, log-rank p = 0.64). However, the LNM group showed a trend toward poorer outcomes. No significant differences were observed between the two groups in the other subgroup analyses (Fig. 5e and Supplementary Fig. 8). In conclusion, we found that the TDLN status affected the immune status of the primary tumor, although this effect was spatially localized.

Fig. 5: Pathological characteristics of primary tumors with and without LNM.

a Halo image analysis automatically identified IHC-positive cells at the IF. b Pearson correlation analysis between manual counting and AI counting, showing the accuracy of Halo-image analysis. IHC staining of CD3 and CD19 at the IF of primary tumors with and without LNM: c CD3, d CD19. e Violin plot showing the quantitative results of IHC staining. Differences in the KM grade of primary tumors with and without LNM: f pathological images, g bar plot of quantitative results. Differences in CLR count between primary tumors with and without LNM: h pathological images, i violin plot of the quantitative results. AI artificial intelligence, LNM lymph node metastasis, IF invasive front, KM Klintrup–Makinen, CLR Crohn-like lymphoid reaction, IHC immunohistochemistry.

The impact of LNM on liver metastases in patients with CRLM

Finally, we further investigated the impact of LNM on liver metastasis. A total of 190 proteins were differentially expressed, with 168 upregulated and 22 downregulated DEGs (Fig. 6a, b). KEGG enrichment analysis revealed that the upregulated proteins were mainly enriched in metabolic pathways, including drug metabolism, protein metabolism, and oxidative phosphorylation, while the downregulated proteins were primarily associated with cell proliferation-related pathways, such as the TGF-β signaling pathway and the PI3K-A signaling pathway (Fig. 6c). The GO enrichment analysis also indicated that the downregulated proteins were mainly associated with immune response-related biological processes (Fig. 6d).

Fig. 6: Proteomic and pathological characteristics of liver metastasis with and without LNM.

a Heatmap of differentially expressed proteins in liver metastasis with and without LNM. b Volcano plot showing differentially expressed proteins in liver metastasis with and without LNM. c Lollipop diagram illustrating the KEGG pathways enriched in the proteins upregulated and downregulated in liver metastases. d Bar plot illustrating the GO terms enriched in proteins downregulated in samples from patients with liver metastasis with LNM; e, f Differences in the HGP of samples from patients with liver metastasis with and without LNM: e pathological images, f bar plot of quantitative results. BP biological process, CC cellular component, MF molecular function, LNM lymph node metastasis, HGP histological growth pattern, d desmoplastic, p pushing, r replacement, KEGG Kyoto Encyclopedia of Genes and Genomes, GO Gene Ontology.

Although we found no significant difference in the density of immune cells in liver metastatic lesions between the two groups (Supplementary Fig. 10), the morphological evaluation revealed that the HGP of liver metastatic lesions was predominantly desmoplastic HGP (dHGP) in the group without LNM, while it was predominantly a non-dHGP, including the pushing HGP (pHGP) and the replacement HGP (rHGP), in the group with LNM (Fig. 6e, f, p < 0.01).