Dissecting origin factors of lymph node metastasis in non-small cell lung cancer via multimodal omics

Domain: biology
Published: 2026-04-13
Editorial outcome: major_revision

Abstract
Cancer metastasis accounts for the majority of cancer-related deaths and lymph node metastasis is critical in cancer staging and prognosis. However, the mechanisms involved in lymph node metastatic non-small cell lung cancer (NSCLC) remain unclear. Here, we delineate the cellular and spatial landscape of metastatic lymph nodes and their matched primary tumors in NSCLC using multimodal omics. In results, tumor core and peripheral regions respectively exhibit high stemness and elevated expression of immunoglobulins. Increased expressions of CD274 (PD-L1) are observed in the immunoglobulin-high malignant cells and SPP1+IFI30+ macrophages. Notably, residual and tertiary lymphoid structures are identified. However, their anti-tumor immunity is compromised by CXCL13+ T cell exhaustion and immune exclusion. Furthermore, we identify the ferroptosis signaling, with its key factor FTH1 showing strong associations with cancer cell stemness and metastasis. This study elucidates the mechanisms of metastasis and immune evasion in lymph node metastatic NSCLC, providing insights for targeted therapies.

Introduction
Lung cancer is the most common cause of cancer-related death around the world, and non-small cell lung cancer (NSCLC) is the most dominating histologic subtype 1. Approximately 80–90% of lung cancer-related deaths are caused by metastasis 2. Lymph nodes (LNs) play critical roles in the arrest and proliferation of metastatic cells, facilitate secondary dissemination to distant organs, and are integral to immune modulation in cancer 3, 4, 5. It has been validated that LN metastasis can predict poor prognostic outcomes in various types of cancer 6, 7, and is also a very common factor in NSCLC 1. Lymph node dissection is regarded as a standard surgery for NSCLC, contributing to accurate tumor staging and improved prognosis 8.

The potential molecular basis of LN metastasis is complicated and involves multiple processes and factors 3, 9, and current researches are mainly focused on lymph-angiogenesis 10, 11, malignant cell migration towards LNs 12, formation of premetastatic niche 13, immune evasion 14 and metabolic adaptation 15. Nonetheless, the mechanism underlying LN metastasis is still limited understood, more precise investigations are needed to obtain comprehensive insights into the LN and long distance metastasis.

Recent advances in single cell and spatial omics provide effective tools to explore the heterogeneous molecular functions and spatial organization of tumor microenvironment (TME) at a high resolution, promoting the understanding of cellular landscape and key mechanisms in cancer 16, 17, 18, 19, 20. A collection of single cell RNA-sequencing (scRNA-seq) studies have been conducted on NSCLC 21, 22, 23, 24, 25, 26, 27, 28, 29. These studies have presented the cellular and transcriptional landscape of cancer cells, immune cells and stromal cells in lung adenocarcinoma (LUAD) and lung squamous-cell carcinoma (LUSC) 21, 27, early-stage LUAD harboring EGFR mutations 22, metastatic LUAD 23 and treatment relevant lung cancers 25, 30, 31. Differences between LUAD and LUSC have also been investigated 21, 28. Meanwhile, spatial transcriptomics (ST) studies have revealed some spatial features of NSCLC 20, 21, 24. Several scRNA-seq studies also explored the cellular features of metastatic lymph nodes (mLNs) in breast cancer 32 and head and neck cancer 33. However, the cellular and spatial characteristics of mLNs in NSCLC remain unclear, with only limited study about the cellular characteristics of mLNs in LUAD 23. There is currently a lack of systematic collection and integrative analysis of patient-matched mLNs and primary tumor tissues (mLN+ PTs) of NSCLC, which impedes the ability to trace malignant cells and TME between mLNs and their corresponding primary tumors, thereby limiting the depth and scope of such analyses.

Here, we aim to investigate the TME characteristics in both mLNs and primary tumor tissues that may promote metastasis to LNs, using multimodal omics including scRNA-seq, ST and bulk RNA-seq. The cellular features of mLNs are described by integrative and comparative analysis of the scRNA-seq data of patient-matched mLN, mLN+ PT and mLN-positive adjacent normal lung (mLN+ N) samples, as well as the other contrastive samples including normal LN (nLN), mLN negative primary lung tumor (mLN- PT) and adjacent normal lung (mLN- N) samples. Meanwhile, the spatially specific characteristics, particularly the lymphoid structure like (LS-like) and immune exclusion spatial organizations, were also described by integrative analyses of the ST data.

Our study reveals that the LN metastasis of NSCLC involves a complicated TME where cancer cell stemness, malignant-cell-derived-immunoglobulins, cancer associated fibroblast (CAF) cells, epithelial–mesenchymal transition (EMT), immune evasion mechanisms as well as cellular and spatial interactions are involved. Additionally, it serves as a comprehensive data resource for understanding the mechanisms of metastatic and immune evasion involved in LN metastatic NSCLC.

Results
The TME landscape of matched mLNs and mLN+ PTs in NSCLC
To characterize the TME characteristics involved in LN metastasis of NSCLC, we collected mLN (n= 22) and matched mLN+ PT samples (n=20) from 19 NSCLC patients (P0-P18) for scRNA-seq analysis, where matched mLN+ N samples (n=4) were also collected (Fig. 1A, created by BioGDP.com 34). In addition, to explore relevant spatial characteristics, we collected mLNs (n=4) and matched mLN+ PTs (n=4) from 4 of the 19 NSCLC patients for ST analysis (Fig. 1A, Methods). Meanwhile, public scRNA-seq datasets of mLN- samples including mLN- PT, mLN- N and nLN samples, and ST data of mLN- PT samples (De Zuani et al. 21, Zhang et al., 28, Kim et al. 23) were collected for integration and comparison (Fig. 1A). The clinicopathological information of the collected LN metastatic NSCLC patients was provided in Supplementary Data 1.

For the scRNA-seq data, 848337 single cells (334956 cells from the self-generated dataset, 513381 cells from the public datasets) passed quality controls in total. We confirmed the batch effects were corrected (Supplementary Fig. 1A, the bias of the nLN cells were due to the absence of epithelial cells in nLNs). As results, 9 major cell types were identified based on whether each cluster show specific expressions of canonical cell markers: B (CD79A, MS4A1), T (CD3D, CD3E, TRBC2), endothelial (CLDN5, PECAM1), epithelial (EPCAM, KRT17, KRT19), fibroblast (ACTA2, COL1A1, COL1A2), mast (TPSAB1, TPSB2), myeloid (LYZ, CD68, APOE), natural killer (NK) (NKG7, KLRB1, KLRD1), and plasma (MZB1, JCHAIN) 23, 27 (Fig. 1B and Supplementary Fig. 1B, 1C). Differences in cellular compositions across different samples, particularly between mLN+ PTs and mLN- PTs, as well as mLNs and nLNs, were observed. Notably, mLN+ PTs exhibited a higher prevalence of epithelial and plasma cells, while mLNs showed an increased proportion of B cells. In contrast, mLN- PTs had a greater abundance of myeloid cells, and nLNs had a higher proportion of T cells (Fig. 1C). These observations suggest that the TME undergoes alterations at both the primary and metastatic sites during LN metastasis.

For the ST data, there were 48535 spatial spots across the 4 pairs of mLNs (15401 spots) and mLN+ PTs (14082 spots) and the other 8 mLN- PTs (19052 spots). BayesSpace-based 35 spatial clustering analysis identified 21 different clusters (Supplementary Fig. 1D), and the clusters were annotated based on their marker genes (Fig. 1D, 1E and Supplementary Fig. 1D, 1E, Methods). For example, spatial cluster SC1 showed elevated expressions of SFTPC, SCIGB3A1, SCGB3A2, and SFTPB - all classical marker genes of lung epithelial cells - and was therefore annotated as epithelial-dominant (Supplementary Fig. 1D). For spatial clusters SC18 and SC19, which showed co-expressions of marker genes for both myeloid (LYZ, C1QB) and plasma cells (IGHG1, JCHAIN), the dominant cell types were defined as Myeloid+Plasma (Supplementary Fig. 1D), implying the co-localization of myeloid cells and plasma cells. Notably, two of the spatial clusters (SC3 and SC4) showed simultaneous high expressions of MS4A1, TRBC2 (Fig. 1D and Supplementary Fig. 1D), indicating the aggregation of T and B cells 36. The two clusters also showed high expressions of CXCL13, CCL19 and CCL21 (Supplementary Fig. 1D), all of which contribute to the formation of secondary and tertiary LS (TLS) 37, 38, 39. These two clusters were annotated as LS-like spots, representing the residual LSs in the mLNs, and potential TLSs in the mLN+ PTs or mLN- PTs. The spatial organizations of LS-like spots showed differences, with LS-like spots being more spatially aggregated in the mLN+ PTs than the mLN- PTs (Fig. 1F and Supplementary Fig. 1F). Additionally, epithelial cells typically exhibited spatial aggregation and with close proximity with fibroblasts, while sparse immune cell infiltrations into epithelial cells were observed, especially for the mLNs and mLN+ PTs (Fig. 1D, 1F).

Epithelial cells in mLNs and mLN+ PTs exhibit enhanced cancer stemness
We next re-clustering the epithelial cells and identified 7 types of epithelial cell sub-populations based on the canonical markers of lung epithelial cells 40: alveolar type I (AT1) (AGER, CLIC5), alveolar type II (AT2) (SFTPA1, SFTPA2, SFTPB, SCGB3A2), basal (KRT5, KRT15), ciliated (CAPS, TPPP3), club (SCGB1A1, SCGB3A1), and the other uncharacterized epithelial cells which were mainly patient specific (Fig. 2A and Supplementary Fig. 2A-2C, Methods). LUAD and LUSC were respectively dominant by AT2 and basal cells, as previously reported 28 (Fig. 2A, 2B). As expected, the epithelial cell compositions of mLNs were similar with those from the matched mLN+ PTs (Fig. 2B).

Cancer stem cells (CSCs) are important for tumor initiation, proliferation and metastasis 41. We evaluated cell stemness of the epithelial cells using CytoTRACE which can predict differentiation states of cells from scRNA-seq data, with higher scores representing higher stemness (Fig. 2C, 2D, Methods). Our findings revealed that epithelial cells in mLNs exhibited the highest CytoTRACE scores, and mLN+ PTs also demonstrated elevated CytoTRACE scores compared to other samples (Fig. 2C). This indicates a higher degree of stemness in the epithelial cells within both mLNs and mLN+ PTs.

CD44, EpCAM and CD166 (gene names CD44, EPCAM and ALCAM) are widely recognized markers of CSCs in NSCLC 42, 43. Here, epithelial cells exhibiting high CytoTRACE scores also showed a significant higher expressions of CD44, EPCAM and ALCAM (Fig. 2E). Additionally, a series of genes, such as FTH1, FTL, ENO1 and PKM, displayed strong correlations with the CytoTRACE based cell stemness (Fig. 2D, 2F). FTL and FTH1, which encode ferritin light chain and heavy chain respectively, are responsible in iron storage and involved in ferroptosis 44. Meanwhile, FTL and FTH1 along with multiple other ferroptosis relevant genes (PCBP1, PCBP2, SLC7A11, SLC3A2 and GPX4 45) were up-regulated in the epithelial cells of mLN and mLN+ PT samples compared to the mLN+ N samples (Supplementary Fig. 2D). The interplay between iron metabolism and CSC stemness or metastasis has been reported 46, 47. ENO1 has been demonstrated to improve the stemness of CSCs by promoting glycolysis in gastric cancer 48. Pyruvate kinase M2 isoform (PKM2) is capable of modulating the activity of key stemness regulators like c-Myc, and β-catenin, which are crucial for maintaining CSC phenotype 49, 50.

Furthermore, we utilized gene set enrichment analysis (GSEA) to identify CSC relevant pathways (Methods). Our analysis revealed that the cell stemness in both LUAD and LUSC was significantly associated with MYC target signatures and a variety of metabolic pathways, including oxidative phosphorylation (OXPHOS), fatty acid metabolism, glycolysis and reactive oxygen species (ROS) (Fig. 2G). C-Myc has been reported to play an importance role in cell stemness 50, 51. CSCs exhibit high metabolic plasticity and can switch between OXPHOS and glycolysis in response to TME conditions 52, 53. ROS and fatty acid metabolism are also crucial for CSC survival and are linked to iron metabolism 46, 54, 55. We also established a CSC score for the bulk RNA-seq data based on the CSC relevant genes (Methods), and observed that a high CSC score was related with un-favorable overall survival in multiple cancer types (Fig. 2H and Supplementary Fig. 2E, 2F).

Elevated expression of immunoglobulin genes in epithelial cells of mLNs and mLN+ PTs
To further explore the specific properties of epithelial cells within the matched mLNs and mLN+ PTs, we compared the gene expressions of epithelial cells from the matched mLNs and mLN+ PTs to the other samples including mLN- PTs, mLN- Ns and mLN+ Ns. As expected, some of the CSC-relevant genes like FTL, HSPB1 and PABPC1 also showed higher expressions in the matched mLNs and mLN+ PTs than the other samples (Fig. 3A). Interestingly, we also observed elevated expression of multiple immunoglobulin relevant genes including IGHA1, IGHG1, IGHG4 and IGLC2, considering both AT2 and basal cells (Fig. 3A, 3B). The ST data also demonstrated that the elevated expression of immunoglobulin genes in the epithelial spots of the mLNs and mLN+ PTs when compared to mLN- PTs (Supplementary Fig. 3A). The immunoglobulin genes are typically expressed by B/plasma cells, which are known for their ability to produce and secrete antibodies. However, some studies have reported the presence of immunoglobulin expressions in tumor cells, challenging the notion that B/plasma cells are the sole source of immunoglobulins 56. To validate the elevated immunoglobulin gene expressions, we performed comparative analysis using the two public scRNA-seq datasets (De Zuani et al. [21], Zhang et al., [28], Methods). The results confirmed that the elevated expressions of immunoglobulin genes in the epithelial cells of mLN+ PTs than mLN- PTs (Supplementary Fig. 3B, 3C). These results underscore the potential involvement of immunoglobulins in LN metastatic of NSCLC.

Additionally, epithelial cells in the matched mLNs and mLN+ PTs were further divided into malignant and non-malignant based on copy number variations (CNVs) (Supplementary Fig. 3D, Methods). The CNV profiles also showed that malignant cells from one mLN were similar with the mLN+ PT from the same individual, while different patients exhibited divergences (Supplementary Fig. 3D). Malignant cells also showed higher expressions of immunoglobulin genes than the non-malignant epithelial cells (Fig. 3C). Meanwhile, we conducted parallel clustering analyses of malignant cells - one preserving patient-specific batch effects and one with batch effect correction, to ensure robust results (Methods). Both results showed that some of the clusters presented elevated expressions of immunoglobulin genes than the other clusters (Fig. 3D-3F and Supplementary Fig. 3E-3G). Notably, the malignant cells expressing high levels of immunoglobulin genes displayed significantly lower CytoTRACE scores than those with low immunoglobulin expressions (Fig. 3G and Supplementary Fig. 3H). This suggests that immunoglobulin genes are more likely to be expressed by cells with lower stemness, indicating the independent roles of CSCs and immunoglobulin-high malignant cells. Meanwhile, GSEA results showed these immunoglobulin-high malignant cells showed down-regulations of P53 pathway and apoptosis pathways (Fig. 3H), indicating its high proliferation ability. Various CSC-relevant metabolism pathways like fatty acid metabolism, ROS were down-regulated in the immunoglobulin-high malignant cells (Fig. 3H).

Importantly, these immunoglobulin-high malignant cells also exhibited significantly elevated expressions of CD274 (Fig. 3I and Supplementary Fig. 3I) which encodes the key immune checkpoint ligand PD-L1. The up-regulations of PD-L1 can facilitate T cell suppression. In addition, an integrative analysis incorporating scRNA-seq data of metastasis brain (mBrain) from LUAD manifested that the immunoglobulin genes were also highly expressed by mBrain sites (Fig. 3J and Supplementary Fig. 3J, Methods), indicating the potential role of malignant-cell-derived-immunoglobulins in secondary dissemination to distant organs.

We also investigated spatially dependent characteristics of epithelial cells in the mLNs and mLN+ PTs. The epithelial spots were further categorized into core and peripheral tumor regions based on their proximity to non-epithelial spots (Fig. 3K and Supplementary Fig. 3K, Methods). Notably, the peripheral region showed significant up-regulations of immunoglobulin genes, collagen genes (COL1A1, COL1A2) and major histocompatibility class I subunit genes (B2M and HLA-A) (Fig. 3L, 3M). In contrast, the tumor core regions had higher CytoTRACE scores compared to the peripheral regions (Supplementary Fig. 3L). Pathway enrichment analysis revealed that the up-regulated genes in the peripheral regions were enriched by EMT, and immune relevant pathways like interferon gamma response, while the core regions were enriched by CSC-relevant pathways like MYC targets, fatty acid metabolism and OXPHOS (Fig. 3N). NicheNet analyses identified several ligands including TGFB1, TGFB2 and BMP4 which were likely driving the target gene expressions of the EMT-dominated tumor peripheral regions (Fig. 3O).

Fibroblasts form a cell transition niche and immune barrier in the leading edge of tumors in mLNs and mLN+ PTs
Stromal cells were categorized into sub-populations based on specific marker genes of fibroblasts 57 and endothelial cells (ECs) 58, including matrix CAF (mCAF) (MMP11, COL1A1), inflammatory CAF (iCAF) (CFD, C3), vascular CAF (vCAF) (COL18A1, NOTCH3), tumor-like CAF (tCAF) (GAPDH, ENO1), smooth muscle cell (SMC) (TAGLN, ACTA2), pericyte (RGS5, NOTCH3), aerocytes (HPGD), arterial EC (HEY1, IGFBP3), capilary EC (CA4), EC progenitor (CD34, KDR), lymphatic (LYVE1, CCL21), venous (ACKR1) and proliferation stromal (MKI67). Additionally, CST1+ fibroblasts showed specific high expressions of CST1, WIF1, COL8A1 and COL10A1 (Fig. 4A and Supplementary Fig. 4A, 4B, Methods). CST1 (Cystatin 1) is a protease inhibitor. CST1 overexpression has been observed in various cancer types 59, 60, 61. Its role in inhibiting cysteine proteases suggests a potential function in modulating the extracellular matrix. The CST1+ fibroblasts have also been identified in tumor tissues of esophageal squamous cell carcinoma, and show prognostic values 61. Here, the CST1+ fibroblasts were specific to mLNs and mLN+ PTs in LUAD. Moreover, the mCAFs predominated in primary tumors, iCAFs were more enriched by normal lung tissues, vCAFs were most abundant in the mLNs (Fig. 4B). Endothelial cell populations in the mLNs closely resembled those in the mLN+ PTs, with venous ECs being dominant, whereas LECs were most common in the nLNs (Fig. 4C). Additionally, cell-cell communication analysis revealed that the mCAFs/vCAFs interacted with the malignant cells through various collages in both mLNs (Fig. 4D) and mLN+ PTs (Supplementary Fig. 4C). The main receptors of the malignant cells included SDC1, SDC4, CD44 and integrin genes. Previous studies have reported that SDC1 plays as role in the activation of CSCs 62. Meanwhile, the collagen-integrin interactions can help malignant cells to invade blood vessels and develop a metastatic niche in distant organs 63. Among the fibroblast sub-populations, the mCAFs and CST1+ fibroblasts had higher cell transition scores than the others (Fig. 4E, Methods).

Notably, the mLNs or mLN+ PTs exhibited elevated expressions of the vCAF marker genes (COL8A1, COL18A1 and IGFBP7) and several cell transition relevant genes like TGFBI, AMBRA1 and PLEKHA7 compared to the other sample types in the stromal cells (Fig. 4F and Supplementary Fig. 4D). AMBRA1 plays a key role in tumorigenesis and progression through regulating the autophagy process, cell transition and a variety of processes 64. PLEKHA7 regulates the assembly of cytoplasmic apical zonula adherens junctions, thereby maintaining cell-cell adhesion and tight-junction 65. TGFBI was up-regulated in the mLN+ PTs than the mLN- PTs (Fig. 4F). TGF-β as a cytokine can secret and is critical for epithelia cell transition. Meanwhile, its role in the exclusion of CD8+ T cells from tumor core areas has been demonstrated 66. ST-based analysis revealed that fibroblast spots were with significantly higher CCNs with epithelial spots (p=0.00012, W = 249, 95 % confident interval = [-366, -135], by unpaired two-sided Wilcoxon test) than other spots, suggesting that fibroblasts are with closer proximity to the epithelial cells than the other cells, particularly in the mLN+ PTs (Fig. 4G). Immune cells were relatively far from the tumor cores, being separated by the fibroblasts, rather than infiltrated into the tumor cells in mLNs and mLN+ PTs (Fig. 1D). Meanwhile, the proportions of fibroblasts, especially mCAF and tCAF cells showed negative correlations with both T and B cells (Fig. 4H). These findings suggest the fibroblasts may form an immune exclusion barrier in the leading edge of tumor cells during LN metastasis in NSCLC.

Considering the specificity of the CST1+ fibroblasts in the mLNs and mLN+ PTs, we examined the spatial expressions of CST1 in different samples. The mLNs and mLN+ PTs also exhibited elevated expression of CST1 compared to the mLN- PTs based on the ST data (Fig. 4I, 4J and Supplementary Fig. 4E). Pathway enrichment analysis revealed that expression of CST1 was significantly associated with cell transition process (Supplementary Fig. 4FL). Additionally, several other cell transition-relevant genes like CTGF 67, ISLR 68 and SMOC2 69 also showed elevated expressions in the mLN+ PTs based on the ST data (Fig. 4I). These up-regulated genes were significantly enriched by pathways like EMT and interferon response pathways (Fig. 4K). NicheNet analysis identified the potential ligands for CST1 and the other up-regulated genes of mLN+ PTs included JAG1, INHBA, IL25 and IL4 (Fig. 4L). Meanwhile, fibroblasts in close proximity to epithelial cells showed significantly higher CST1 expressions than those located far away (Fig. 4M). These findings suggest the role of fibroblasts, especially the CST1+ fibroblasts, in promoting cell transition in LN metastatic NSCLC.

Role of myeloid cells in immune suppression and LS formation in mLNs and mLN+ PTs
Re-clustering analysis of the myeloid cells identified 11 sub-populations based on the marker genes of macrophage (MP), dendritic cell (DC), monocyte 70 and some specific genes including SPP1, FABP4, GPNMB and LAMP3 (Fig. 5A-5C and Supplementary Fig. 5A, Methods). Both SPP1 71 and GPNMB 72 have been reported to facilitate the polarization of M2 macrophages which supports tumor progression by producing immunosuppressive cytokines in MPs, while FABP4 is mainly produced by alveolar macrophages which play a key role in innate immune response against infection 73. Here, FABP4+ MPs were abundant in normal lungs, while SPP1+ MPs were more enriched by primary tumors and mLNs, and GPNMB+ MPs were more enriched by the mLNs and mLN+ PTs of LUAD (Fig. 5C, 5D). Further comparing mLNs to nLNs and mLN+ PTs to mLN- PTs, it was observed that both mLNs and mLN+ PTs exhibited larger proportions of SPP1+ MP and neutrophil (Neu) cells compared to their respective controls (Fig. 5D). Meanwhile, the high expression of SPP1 was associated with unfavorable prognosis in LUAD (Fig. 5E, Methods). Notably, the proportion of SPP1+ MP cells, rather than the other sub-populations of MP cells, showed a strong correlation with the proportion of malignant cells (Fig. 5F).

Further comparisons of the gene expressions of the SPP1+ MPs across different sample types discovered that certain genes including IFI30, FOS, HSPA1A, IGHG4, CCL3, CCL4, etc., were specifically up-regulated in the mLNs and mLN+ PTs (Fig. 5G). Meanwhile, clustering analysis of the SPP1+ MPs identified a cluster predominantly enriched by mLNs and mLN+ PTs, with a particular prevalence in the mLNs (Supplementary Fig. 5B-D). This cluster also showed significantly up-regulated expressions of IFI30 along with the other aforementioned genes (Supplementary Fig. 5D). IFI30, also known as gamma interferon-inducible lysosomal thiolredoxin, is mainly expressed in antigen-presenting cells, including MP, DC and B cells 74. Recent studies have revealed that IFI30 can promote tumor metastasis and immune evasion in cancer 75, 76. Remarkably, overexpression of IFI30 has been shown to promote the M2 polarization of MPs and increase the expression of PD-L1 in breast cancer 75. In our findings, we also observed that CD274 was up-regulated in the SPP1+ IFI30+ MPs (Supplementary Fig. 5E). Additionally, the SPP1+IFI30+ MPs also showed up-regulated expressions of GABARAP, MIF, TNFRSF13, etc. (Fig. 5H). These genes have been demonstrated to promote the M2 polarization of MPs 77, 78, 79. Among them, MIF is a pro-inflammatory cytokine that plays a crucial role in various immune processes 78, and TNFSF13 (also known as APRIL), a proliferation-inducing ligand of the tumor necrosis factor (TNF) superfamily, has emerged as a promising target for countering immunosuppression 79. Moreover, cell-cell communication analysis revealed that the SPP1+ MPs interacted with the malignant cells through MIF (Supplementary Fig. 5F, 5G) and interacted with the T cells, through CXCL9-CXCR3 and CXCL16-CXCR6 axes (Fig. 5I). Spatial examination of cellular distributions indicated a non-co-localization pattern between SPP1+ MP and malignant cells, suggesting potential long-range interactions mediated by paracrine or endocrine signaling 80 (Supplementary Fig. 5H, 5I). In contrast, SPP1+ MP and T cell sub-populations exhibited co-localization, indicating short-range cellular crosstalk (Supplementary Fig. 5H, 5I).

Based on the ST data, the myeloid-relevant spots (myeloid/plasma+myeloid spots) were categorized into far and near groups, based on their distance from epithelial cells (see Methods). The analysis showed that myeloid spots in close proximity to epithelial cells exhibited higher expression levels of CXCL9, PGK1 and SPP1, while those farther away showed higher expressions of IFI30, APOE and GPNMB (Fig. 5J), indicating the spatial distinct characteristics of myeloid cells. Moreover, the myeloid cell spots were divided into 7 spatial clusters (SCs) (Methods), and the spatial cluster 7 (SC7), which was in close proximity to LS-like spots, showed remarkably high expression of CCL21 and CCL18 (Fig. 5K, 5L and Supplementary Fig. 6A). Meanwhile, the SC7 was more enriched by the mLNs (Supplementary Fig. 6B). Most of the marker genes were primarily expressed by the GPNMB+ MPs or SPP1+ MPs, while CCL21 and CPVL were also highly- expressed by DCs (Supplementary Fig. 6C). The chemokines CCL21 and CCL18 play important roles in immune cell trafficking and are both involved in the formation and function of TLSs 81. CCL18 can also contribute to immune suppression when it attracts T regulatory (Treg) cells 82 and M2-like MPs 83.

Reduction of anti-tumor immunity through T cell exhaustion and immune exclusion
T and NK cells were further classified into 12 sub-populations (Fig. 6A, Methods). In addition to some well-defined sub-populations including CD4 effector memory T, CD4 naïve T, CD4 resident effector memory T, Treg, CD8 cytotoxic T, CD8 exhausted T as well as NK cells, we also identified several other sub-populations (Fig. 6A and Supplementary Fig. 7A, Methods). Sub-populations showing specific high expressions of ALDOA, ARHGAP15 and CXCL13 were defined by the gene names (Fig 6A and Supplementary Fig. 7A, 7B). A small sub-population which also showed specific high expressions of HAVCR2, TIGIT, CTLA4, classical markers of T cell exhaustion, but with differential expression profiles from the CD8 exhausted T cells, was defined as “other exhausted T cells” (Supplementary Fig. 7A). This sub-population was specific in the mLNs (Supplementary Fig. 7C). Some CXCL13+ T cells are characterized as follicular helper (Tfh) 84 or exhausted T cells 85 by recent studies. Here, the CXCL13+ T cells showed specific high expressions of CXCL13, ICA1, IL6ST, CTLA4 and PDCD1, however, BCL6 and CXCR5, two classical marker genes of Tfh cells 86, were not highly expressed (Supplementary Fig. 7A, 7B). The CXCL13+ T cells also exhibited the highest exhaustion scores among the T cells (Fig. 6C, Methods), and showed high expressions of exhaustion marker genes like CTLA4 87 and PDCD1 88 (Fig. 6D and Supplementary Fig. 7B), highlighting its exhausted state. Both mLNs and mLN+ PTs contained more Tregs and CXCL13+ T cells than other sample types (Fig. 6B and Supplementary Fig. 7C). ALDOA, or Aldolase A, is an enzyme crucial for glycolysis and gluconeogenesis which are important for maintaining cellular energy balance, and ALDOA expression may influence T cell behavior 89, potentially affecting their ability to effectively target and eliminate cancer cells 90. The ALDOA+ T cells were predominantly enriched by nLNs, but significantly reduced in mLNs (Fig. 6B), strongly suggests a state of metabolic impairment or dysfunction, contributing to the ineffective anti-tumor response in mLNs. ARHGAP15, a gene encoding a Rho GTPase-activating protein, can inhibit Rac1 activity 91. Rac1 is a small GTPase involved in various cellular processes, including T cell activation and migration 92. ARHGAP15+ T cells were more enriched in mLNs and mLN+ PTs (Fig. 6B), indicating suppression of T cell activation and motility, offering a plausible mechanism for their trapped and dysfunctional state. The mLNs were mainly composed of CD4 T cells, especially CD4 naïve T cells, while the mLN+ PTs contained more CD8 exhausted T cells (Fig. 6B and Supplementary Fig. 7C). Moreover, the LS-like spots exhibited higher expressions of the CXCL13+ T relevant genes including CTLA4, CXCL13, PDCD1 and TRBC2 than other spatial spots (Fig. 6E). The CXCL13+ T cells were primarily enriched by the LS-like spots based on the deconvolution analysis of ST data (Fig. 6F and Supplementary Fig. 7D). Notably, the LS-like spots in the mLNs and mLN+ PTs showed significantly higher exhaustion scores than the mLN- PTs (Fig. 6G). Taken together, the anti-tumor cytotoxicity of T cells in both mLNs and mLN+ PTs appears to be reduced by Tregs and exhausted T cells.

CXCL13 is strongly associated with the aggregation of T and B cells, and it plays a pivotal role in the formation of secondary LSs in LNs and TLSs in tumors 93. The gene ontology (GO) based pathway enrichment analysis revealed the CXCL13+ T cells were associated with B cell receptor signaling pathway, indicting its ability to recruit B cells. Cell-cell communication analysis revealed that the CXCL13+ T cells can interact with the B cells via the CXCL13-CXCR5 axis (Fig. 6H), while the CAF and LAMP3+ DC cells can interact with the CXCL13+ T cells through the CCL19-CCR7 axis (Fig. 6I). CARD-based deconvolution analysis of the ST data also showed CXCL13+ T cells were significantly co-localized with key interacting cells, including LAMP3+ DCs, GC B cells, and activated B cells, within LS-like spots (Supplementary Fig. 7E, 7F). While spatial distributions of iCAF and mCAF cells were broader (Supplementary Fig. 7E), and did not show significant enrichment in LS-like spots (Supplementary Fig. 7F), suggesting that their interaction with CXCL13+ T cells might be less constrained to the core LS-like architecture. Additionally, we also found the CXCL13+ T cells and Tregs were more enriched by the mLN+ PTs from NSCLC patients without smoking history comparing to those with smoking history (Supplementary Fig. 7G).

The B and plasma cells were further characterized into sub-populations based on the re-clustering analysis (Fig. 6J and Supplementary Fig. 8A, Methods). The mLNs contained more activated B cells, while the mLN+ Ns and mLN+ PTs contained a large fraction of plasma cells (Fig. 6K), indicating the activation of B cell mediated adaptive immune responses during LN metastasis in NSCLC. The plasma cells in the mLN+ PTs showed frequent interactions with malignant cells and mCAFs (Fig. 6L). However, the excess activated B cells and plasma cells in mLNs and mLN+ PTs indicated prolonged activation of B cells which can lead to immune dysfunction or exhausted B cell states 94. Plasma cells can produce antibodies to recognize and attack tumor cells. Further clustering analysis of the plasma cells revealed that the mLN+ PTs and mLN- PTs were respectively enriched by the Plasma_C2 and Plasma_C1 clusters (Fig 6M, 6N). The two clusters showed distinct immunoglobin expression patterns, with Plasma_C2 showing significantly higher expressions of IGHG4 and IGHV4-34, whereas Plasma_C1 showing higher expressions of IGHA2, IGHG2 and IGKC (Fig. 6O and Supplementary Fig. 8B), indicating remarkable differences between the antibody isotypes of plasma cells from mLN+ PTs and mLN- PTs. Moreover, since the malignant cells in the mLNs and mLN+ PTs also exhibited high expressions of IGHG4 (Fig. 3B), it is possible that the malignant-cell-derived IgG neutralizes the plasma-cell-derived IgG to inhibit complement-dependent cytotoxicity, thereby preventing host immune system and promoting tumor growth, as revealed previously 95. Deconvolution analysis of B cell sub-populations revealed that the LS-like spots were enriched with germinal center (GC) B cells, a hallmark of mature TLS 96, while few plasma cells in the LS-like spots (Supplementary Fig. 8C). The proportions of GC B cells in the mLNs were significantly higher than the mLN+ PTs (Supplementary Fig. 8D), indicating the TLSs in the mLN+ PTs may be immature.

To further characterize the TME surrounding the LS-like spots, we extracted the two-step neighbors of each LS-like spot (Fig. 6P, Methods). In the mLNs and mLN+ PTs, the neighbors of the LS-like spots were primarily other LS-like spots (Fig. 6P, 6Q), indicating the spatial aggregation of the LS-like spots. Meanwhile, the LS-like regions in the mLN- PTs were more likely to be infiltrated into epithelial cells (Fig. 6P), however, fewer LS-like regions were infiltrated into the tumor areas in the mLNs and the mLN+ PTs, with fibroblast or myeloid cells forming an immune-excluded 97 interface between the tumor and immune cells (Fig. 1E, 4G, 6P). These findings together described a complex immunosuppressive TME that caused by the immune exhaustion states and spatial organizations.

Distinct cellular mechanisms involved in lymph node metastasis between LUAD and LUSC
To further investigate the differences between LUAD and LUSC regarding LN metastasis, we quantified cell composition changes of mLN+ PT comparing to either normal lung or mLN- PT tissues in LUAD and LUSC samples respectively. Both conserved and specific cellular features were identified (Supplementary Fig. 9). Notably, mLN+ PTs of LUAD displayed significant accumulation of ARHGAP15+ T, CXCL13+ T, plasma and Treg cells. In contrast, mLN+ PTs of LUSC displayed depletion of various effector immune cells, particularly CD4 effector T cells and NK cells (Supplementary Fig. 9). CAFs, such as vCAF and tCAF, were consistently enriched by the mLN+ PT samples in both LUAD and LUSC, suggesting a similar stromal response. These divergent immune landscapes suggest fundamentally distinct metastatic mechanisms - where LUAD is more likely to develop an immunosuppressive niche through Treg and exhausted lymphocytes, while LUSC appears to evade immunity more by immune cell exclusion.

Validations of key TME characteristics associated with LN metastasis in NSCLC
First, we examined the CSC relevant gene signatures in multiple independent scRNA-seq datasets of NSCLC. In the GSE127465 dataset 98, the epithelial cells from the primary tumors at the N1 stage exhibited significantly higher CSC scores and up-regulated expressions of the CSC relevant genes compared to those at the N0 stage (Fig. 7A). Meanwhile, analysis of datasets GSE123904 99 and GSE131907 23 demonstrated that epithelial cells derived from distant metastatic sites and mLNs consistently showed higher CSC scores and increased expression of CSC-associated genes than both primary tumor cells and normal lung tissues (Fig. 7B, 7C and Supplementary Fig. 10A). These results indicate the CSCs play important roles in both LN and distant metastasis.

Furthermore, to investigate the general expressional patterns and prognostic effects of the various signature genes identified in this study, we performed clustering analysis of the TCGA-LUAD tumor samples based on the bulk RNA-seq data of the signature genes (Methods). Four bulk RNA-seq data based clusters (BCs) were identified, where the cluster BC3 showed the most un-favorable prognosis (Fig. 7D, 7E and Supplementary Fig. 10B). Genes up-regulated in BC3 included several CSC relevant genes like FTH1, GPX2, AKR1C1, AKR1C2 and AKR1C3. Up-regulation of these genes like GPX2, AKR1C1 have been demonstrated to inhibit ferroptosis 100, 101. FTH1, a key regulator in iron storage, was identified as one top-ranked CSC-relevant gene in both LUAD and LUSC (Fig. 2D, 2E). Moreover, mBrain and mLN samples showed higher expressions of FTH1 than normal lung tissues, and primary tumors in the epithelial cells (Fig. 7F). The prognostic significance of BC3 was validated by the other two independent bulk RNA-seq datasets, which respectively included NSCLC and LUAD tumor samples (Supplementary Fig. 10C, 10D). Meanwhile, the BC3 displayed a higher prevalence of metastatic LNs than the other BCs in the two datasets (Supplementary Fig. 10E).

To validate FTH1 and immunoglobulin gene relevant functions both in vitro and in vivo, we conducted loss-of function analyses of FTH1 and IGHG4 in lung cancer cell lines (Fig. 7G and Supplementary Fig. 10F). Notably, FTH1 and IGHG4 knockdown reduced the several EMT marker genes expression of Vimentin, N-cadherin, MMP2/9, ZEB1/2 and Twist in A549 and H460 cells, while increasing E-cadherin expression increased in FTH1 knockdown cells (Fig. 7G and Supplementary Fig. 10F). Furthermore, two independent FTH1 short hairpin RNAs (shRNAs) significantly decreased the proliferation of A549 and H460 cells (Supplementary Fig. 10G). However, IGHG4 knockdown did not influenced A549 and H460 cell proliferation (Supplementary Fig. 10H). Interestingly, treatment with the ferroptosis inducer Erastin did not further enhance ferroptosis-dependent growth inhibition in FTH1 knock-downed lung cancer cells (Supplementary Fig.10I, 10J). Additionally, FTH1 and IGHG4 knockdown decreased cell invasion in both A549 and H460 cells (Fig. 7H and Supplementary Fig. 10K, 10L). In vivo, FTH1 knockdown cells formed fewer metastatic lung nodules compared to control cells, with no significant changes in body weight observed during the metastatic colonization assay (Fig. 7I, 7J and Supplementary Fig. 10M). Collectively, these findings prove that FTH1 and IGHG4 play crucial regulatory roles in lung cancer metastasis.

We further analyzed the expression profiles of SPP1+IFI30+ MP and CXCL13+ T cells across the validation datasets. Notably, the primary tumors at the N1 stage demonstrated significantly higher scores considering the SPP1+IFI30+ MP and CXCL13+ T relevant gene signatures in the myeloid and T cells compared to those at the N0 stage (Fig. 7K). This pattern was consistently observed in distant metastatic lesions and mLN samples, which exhibited higher scores for both immune cell types (Fig. 7L, 7M). The expressional signatures of these two immune cell types, especially the SPP1+IFI30+ MP, showed strong positive correlations with the CSC signature (Supplementary Fig. 10N). Supporting these findings, both parental and sixth-generation LN samples derived from LN metastasis mouse models 14 maintained elevated expression scores for these immune cell populations (Supplementary Fig. 10O). These results highlight the potential significance of SPP1+IFI30+ MP and CXCL13+ T cells in facilitating NSCLC metastasis. Moreover, analysis of an independent NSCLC cohort 102 revealed that elevated signature scores of SPP1+ IFI30+ macrophages and CXCL13+ T cells in pre-treatment non-responders to immuno-chemotherapy (Supplementary Fig. 11A, 11B), suggesting their associations with resistance to immuno-chemotherapy.

We further examined our findings in the mLNs of breast cancer based on ST data 32. Similar to NSCLC, mLNs of breast cancer also exhibited residual LSs (LS-like spots defined based on high expressions of CXCL13, CCL19, CCL21 and MS4A1, Supplementary Fig. 12A-12F), with significantly elevated CXCL13+ T cell signature scores (Supplementary Fig. 12G). In addition, the mLNs of breast cancer also displayed significant higher expressions of SPP1 and IFI30, along with elevated SPP1+IFI30+ MP signature scores compared to nLNs (Supplementary Fig. 13A-13E). These conserved features suggest a potential universal role for SPP1+IFI30+ MP and CXCL13+ T cells in LN metastasis across multiple cancer types.

Methods
Patient samples
All procedures performed in our study were in accordance with the ethical standards of our institutional research committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. This study was approved by the Institutional Ethics Committee of Liaoning Cancer Hospital & Institute before the initiation of the study, and all patients signed informed consent. Matched metastatic lymph node (mLN) and primary tumor tissue samples (mLN+ PT) were collected from 19 patients diagnosed with NSCLC with lymph node metastasis from Jan-2023 to Jan-2024 at Liaoning Cancer Hospital & Institute in China. The normal lung tissues from 4 of the 19 patients were also collected and termed as mLN+ normal (mLN+ N) lung tissues. All hematoxylin and eosin slides of the samples were subjected to evaluation for histopathological morphology and tumor components by expert pathologists. The samples enrolled in this study should meet the following criteria: (1) all patients have underwent surgery,(2) All tumor and lymph node specimens were definitively identified as non-small cell lung carcinoma and metastatic lymph nodes, respectively, following intraoperative frozen section pathological examination and (3) All patients had not been subjected to any preoperative antitumor therapy. The clinical information of the enrolled NSCLC patients were also recorded, including age, TNM stage and smoking history. Tumor and non-involved lung samples were obtained from surgical specimens of patients undergoing resection at Liaoning Cancer Hospital & Institute. Patients were consecutively enrolled from a single institution based on the above criteria. While this single-center surgical cohort may introduce geographic and treatment-selection biases, the core findings were validated using independent multi-center public datasets to enhance generalizability. All available clinicopathological patient information were provided in Supplementary Data 1.

Cell culture
Human lung cancer cell lines A549 and H460 were obtained from American Type Culture Collection (ATCC, Manassas, VA, USA; A549: CRM-CCL-185, H460: HTB-177) and cultured in RPMI 1640 medium (GIBCO, USA). Cell lines were maintained in culture supplemented with 10% (v/v) FBS (GIBCO, USA) and 1% penicillin/streptomycin (Thermo) at 37°C with humidified 95% air and 5% CO2 in a humidified incubator.

In vivo animal studies
For tumor lung metastasis assay, male BALB/c nude mice ( 6-8 weeks old) were obtained from China Medical University. All mice were maintained under specific pathogen-free conditions with a 12-hour light/dark cycle, at an ambient temperature of 20-24°C and a relative humidity of 50-60%. A549 cells transfected with control (shControl) or shRNA against FTH (shFTH1_1) ( 2×106 per mice) in 150 μL phosphate-buffered saline (PBS) were injected into the tail vein. After 10 weeks, the lungs were removed and photographed followed by nodules counting. All experiments were carried out according to the regulations set by the Ethics Committee of Liaoning Cancer Hospital & Institute (Protocol number: LNCH20210104). Randomization was performed using a simple random number generator to assign mice to cages and treatment groups before the experiment. Investigators were blinded to the group information during the assessment of lung metastatic nodules. The maximal tumor burden permitted by the ethics committee is no more than 1500 mm3. When the tumor burden reached 1500 mm3, mice were euthanized, and the tumors were dissected for further analysis. All the tumors were removed, photographed and weighed.

Transfection and lentivirus production
HEK293T cells were used for transfection and lentivirus packaging. HEK293T cells were transfected using PEI (Polyethylenimine, Polysciences 24765). The transfection mass ratio of plasmids to PEI was 1:3 in Opti-MEM medium (Gibco). After 48 hours transfection, protein expression was tested.
For lentivirus packaging, HEK293T cells were infected with shRNA expressing lentivirus vectors plus the packaging plasmids psPAX2 and pVSVG using PEI. Viral supernatant was harvested after 48 hours, filtered through a 0.45 μm membrane to remove cell debris, then the virus containing medium was used to infect target cells with 5 μg/ml polybrene for 8-24 hours. Puromycin or blasticidine were used to obtain positive cells, based on the resistance gene encoded by the lentiviral vector. Knockdown of each shRNA target was confirmed by immunoblotting. The sequences of shRNAs used were as follows:
Human FTH shRNA1: CCTGTCCATGTCTTACTACTT;
Human FTH shRNA2: GCCGAATCTTCCTTCAGGATA;
Human IGHG4 shRNA1: GTCCAGTTCAACTGGTACGTG;
Human IGHG4 shRNA2: GCTGTCACTCACATGCTCTGT;

Immunoblotting assay
Cells were harvested from plates and lysed using RIPA lysis buffer [50 mM Tris-Cl PH7.4, 150 mM NaCl, 1 mM EDTA PH8.0, 0.25% DOC (deoxycholic acid), 10% glycerol, 1% NP-40 and 1% Triton-X100], containing protease inhibitors (Bimake) and phosphatase inhibitors (Bimake). The protein concentration was estimated by EASY Ⅱ Protein Quantitative Kit (Transgene), after centrifuging at 13,000 r.p.m. for 15 minutes at 4 ℃ to remove debris. Equal amounts of protein for each sample was electrophoresed on SDS-PAGE gel and transferred to PVDF membranes (Millipore, Germany), followed by blocking with 5% non-fat milk for 1 hour. Then the proteins on the PVDF membrane were immunoblotted with indicated primary antibodies overnight at 4°C. Primary antibodies used include FTH1 (Cell Signaling Technology, 3998; 1:1000), E-cad (BD Biosciences, 610182; 1:1000), Vimentin (Proteintch, 60330-1-Ig; 1:1000), MMP2 (Proteintch, 10373-2-AP; 1:1000), MMP9 (Proteintch, 10375-2-AP; 1:1000), ZEB1 (Proteintch, 21544-1-AP; 1:1000), Twist (Proteintch, 25465-1-AP; 1:1000), GPX4 (Cell Signaling Technology, 52455; 1:1000), IGHG4 (Proteintch, 66408-1-Ig; 1:5000), N-cadherin (Proteintch, 22018-1-AP; 1:2000)，ZEB2 (Proteintch, 14026-1-AP; 1:4000)，Vinculin (Santa Cruz, sc-73614; 1:1000), Actin (Proteintch, 66009-1-Ig; 1:20000). Secondary antibodies used include goat anti-rabbit (Jackson, 115-005-003, 1:5000) or goat anti-mouse (Jackson, 115-035-003, 1:5000).

Cell proliferation assay
Cell proliferation assays were performed using cell counting kit-8 (CCK-8) assay according to the manufacturer’s requirement. Briefly, 0.8×103 or 1.5×103 cells were suspended in 100 μL of RPMI 1640 medium seeded in 96-well plates. Then cells were cultured under 37℃ and 5% CO condition overnight. 2
At different time points (0, 24, 48, 72, and 96 hours), cells were incubated with a 100 μL relative culture medium containing CCK-8 reagent (10 μL/well) incubated for 4 hours at 37 °C. Plates were measured at OD 450 nm using a plate reader (Cytation5, Biotech).

Transwell migration assay
Cell migration ability was assessed using a transwell assay. The cell transfer chamber was obtained from Corning Company (Costar 3422). 5×104 cells were seeded in the upper chamber in a 150 μL medium without serum, while 600 μL of medium containing 15% fetal bovine serum was added to the lower chamber as a chemoattractant. After incubation for 24 hours, cells were fixed with paraformaldehyde, stained with 0.1% crystal violet (Sorabio, Beijing, China) for 10 min. Non-migrated cells were removed with cotton swabs, and the rest cells were pictured and counted.

H&E Staining
Lung samples were fixed in 4% paraformaldehyde, embedded in paraffin and then sectioned (5 μm) to expose the main nodules. Hematoxylin and eosin (H&E) staining was performed using a standard protocol to evaluate inflammation. The images were captured with a microscope.

Single cell tissue dissociation, preparation and sequencing
The fresh tissues were stored in the sCelLiveTM Tissue Preservation Solution (Singleron) on ice after the surgery within 30 mins. The specimens were washed with Hanks Balanced Salt Solution (HBSS) for three times, minced into small pieces, and then digested with 3 mL sCelLiveTM Tissue Dissociation Solution (Singleron) by Singleron PythoN™ Tissue Dissociation System at 37 °C for 15 min. The cell suspension was collected and filtered through a 40-micron sterile strainer. Afterwards, the GEXSCOPE® red blood cell lysis buffer (RCLB, Singleron) was added, and the mixture [Cell: RCLB=1:2 (volume ratio)] was incubated at room temperature for 5-8 min to remove red blood cells. The mixture was then centrifuged at 300 × g 4 ℃ for 5 mins to remove supernatant and suspended softly with PBS. Finally, the samples were stained with Trypan Blue and the cell viability was evaluated microscopically. Single-cell suspensions (2×105 cells/mL) with PBS (HyClone) were loaded onto microwell chip using the Singleron Matrix® Single Cell Processing System. Barcoding Beads are subsequently collected from the microwell chip, followed by reverse transcription of the mRNA captured by the Barcoding Beads and to obtain cDNA, and PCR amplification. The amplified cDNA is then fragmented and ligated with sequencing adapters. The scRNA-seq libraries were constructed according to the protocol of the GEXSCOPE® Single Cell RNA Library Kits (Singleron) 118 . Individual libraries were diluted to 4 nM, pooled, and sequenced on Illumina novaseq 6000 with 150 bp paired end reads.

Spatial transcriptomics sample preparation, tissue optimization and sequencing
Fresh tissues were concurrently frozen and embedded in optical cutting tissue (OCT) compound in liquid nitrogen. The RNA quality of OCT embedded block was assessed by Agilent 2100. RNA integrity number (RIN) of tissues greater than 7 were used for visium spatial gene expression experiments. Cryosections were performed on a Leica CM3050 and brightfield images were taken on a Leica Aperio Versa8 whole-slide scanner at 20× resolution.

The Visium Spatial Tissue Optimization Slide & Reagent kit (10X Genomics) was used to optimize the tissue permeabilization conditions according to Visium Spatial Tissue Optimization User Guide (CG000238, 10X Genomics). Briefly, the workflow includes placing tissue sections on 7 capture areas on a Visium Tissue Optimization slide. The sections are fixed, stained, and then permeabilized for different times. mRNA released during permeabilization binds to oligonucleotides on the capture areas. Fluorescent cDNA is synthesized on the slide and imaged. The permeabilization time that results in maximum fluorescence signal with the lowest signal diffusion is optimal. If the signal is the same at two time points, the longer permeabilization time is considered optimal.

The Visium Spatial Gene Expression Slide & Reagent kit (10X Genomics) was used to construct sequencing libraries according to the Visium Spatial Gene Expression User Guide (CG000239, 10X Genomics). A 10um frozen tissue section was placed on one of the Visium gene expression slide capture areas in a slide. After tissue hematoxylin and eosin (H&E) staining, bright-field images were acquired. Reverse transcription experiment was conducted and sequencing libraries were prepared following the manufacturer’s protocol. Sequencing was performed with a Novaseq PE150 platform according to the manufacturer’s instructions (Illumina) at an average depth of 300 million read-pairs per sample.

Single-cell RNA pre-processing
Raw reads from scRNA-seq were processed to generate gene expression matrixes using CeleScope (https://github.com/singleron-RD/CeleScope) pipeline. Raw reads were first processed with CeleScope to remove low quality reads to trim poly-A tail and adapter sequences. Cell barcode and UMI were extracted. After that, we used STAR v2.6.1a 119 to map reads to the reference genome GRCh38 (ensembl version 92 annotation). UMI counts and gene counts of each cell were acquired with featureCounts v2.0.1 120, and used to generate expression matrix. Doublet scores of the cells in each sample were caculated using R package scDblFinder (v1.12.0) 121 with default parameters.

Integrative clustering analysis of scRNA-seq data
In addition to the scRNA-seq data of the above mLN relevant samples (mLN+ N, mLN+ PT, mLN), we also collected three additional scRNA-seq datasets to do integration analysis. The first one was from Zhang’s integrative study about primary LUAD and LUSC 28 (https://lungcancer.chenlulab.com/#/download), where the normal and tumor tissue samples of the mLN-negative (N stage was recorded as N0) patients were taken as the mLN- N and mLN- PT samples for integration. The second dataset was from Kim’ study (GSE131907 in GEO) 23, where the scRNA-seq count matrix of normal lymph node (nLN) tissue samples of LUAD patients were collected for the integration. The third one was from DCe Zuani’s study 21 which was also about primary LUAD and LUSC (https://www.ebi.ac.uk/biostudies/, with accession number E-MTAB-13526), where the tumor samples from mLN-negative patients were also taken as the mLN- PT samples for integration. The count matrix from different studies were merged as one. For quality control (QC), we excluded cells either expressing fewer than 1000 nCount_RNA or 300 nFeature_RNA or with a mitochondrial proportion higher than 25%. Meanwhile, poor QC cells (cell clusters with low nFeature_RNA numbers, high doublet scores, or enriched by abnormal marker genes) identified during the clustering and sub-population analysis were also removed during the analysis.

Then, the integrative analysis was mainly performed using R package Monocle3 (v1.3.4) 122 which showed a relatively well overall performance among several scRNA-seq analysis tools 123. Firstly, preprocess_cds function was applied for data normalization and principle component analysis (PCA). The batch effects among different patients from different studies were removed by matching mutual nearest neighbors 124 based on the align_cds function. Then, reduce_dimension function was applied for UMAP reduction where 20 principle components (PCs) were utilized. Clustering of the cells was based on a community detection technique as part of the phenoGraph algorithm using the cluster_cells function. The marker genes, i.e., genes most specifically expressed in each cluster were identified using the top_markers function, where the specificity was based on the Jensen-Shannon divergence, and P-value was based on likelihood ratio test, adjusted by Bonferroni correction method. Boxplots comparing the normalized expressions of genes in different sample types or conditions were generated using the plot_genes_violin function.

Identification of major cell types
The major cell types were determine based on classical marker genes of epithelial (EPCAM, KRT19, KRT17), fibroblast (COL1A1, COL1A2, ACTA2), endothelial (PECAM1, CLDN5, VMF), myeloid (LYZ, CD68), T (CD3D, CD3E), NK (NKG7, KLRD1, KLRB1), B (CD79A, MS4A1) and plasma (MZB1, JCHAIN) cells. Each cluster was assigned with one of the major cell types based on whether the clusters showed specific high expressions of the classical marker genes (p.adjust < 0.01, likelihood ratio test, two-sided, adjusted by Bonferroni). The cell types were also confirmed and validated based on the prediction results by SingleR (V 1.10.0) 125 and the previously annotated cell types of the mLN- samples 28.

Identification of epithelial cell sub-populations
Sub-populations of the epithelial cells were identified by extracting the count matrix of the epithelial cells, and re-clustering the epithelial cells following the pre-processing, PCA, batch effect correction, uniform manifold approximation and projection (UMAP) reduction, clustering, identifying cluster markers and annotation as the above integrative analysis work-flow. The annotation was based on the defined marker genes of AT1 (AGER, CLIC5), AT2 (SFTPB, SFTPA1, SFTPA2, SFTPC), basal (KRT5, KRT15), ciliated (TPPP3, CAPS), club (SCGB1A1, MUC5AC) and the clusters without the above marker genes but showing patient-specific enrichment were annotated as Patient Specific clusters. The other clusters which showed poor QC data (low nFeature_RNA or high doublet score) were removed.

CytoTRACE analysis about the epithelial cells
The stemness of the epithelial cells were inferenced using R package CytoTRACE (V 0.3.3) with default parameters. The CytoTRACE method can predict differentiation states of cells from scRNA-seq data by a two-step procedure of smoothing gene counts signature on the basis of transcriptional covariance among single cells 126. The CytoTRACE method assigns a score (range from 0 to 1) to each individual cell, with higher scores representing higher stemness or lower differentiation. The CytoTRACE scores of the epithelial cells from LUAD and LUSC patients were respectively calculated. The the Pearson correlations between gene expressions and the CytoTRACE scores were also calculated, and genes with correlation coefficient > 0.5 for both LUAD and LUSC were taken as the cancer-stem-cell (CSC) relevant gene signatures. The ranked gene list based on the correlations with CytoTRACE scores were also used as the input of gene set enrichment analysis (GSEA), and the significantly enriched pathways were identified.

CNV inference
The copy number variation (CNV) inference was performed using R package inferCNV (V 1.12.0) (https://github.com/broadinstitute/InferCNV). For each individual patient, if there was one mLN+ N samples, we created one inferCNV object with the raw count matrix of epithelial cells of all the patient’s samples as the input matrix, and the normal epithelial cells from the mLN+ N samples were utilized as a reference of normal cells; if there was no mLN+ N samples, the raw count matrix of epithelial and main immune cells (T, B and myeloid cells) of the patient’s samples were utilized as the input matrix, with the immune cells as the reference of normal cells. Then, the CNV of the epithelial cells in the mLN and mLN+ PT samples were estimated with default parameters (cutoff=0.1, cluster_by_groups=TRUE, denoise=TRUE, HMM = FALSE) using the run function.

Identification and clustering analysis of malignant cells
For each individual patient, the epithelial cells of the mLN samples were malignant cells metastasis from the primary tumors, showing significant CNVs comparing to the normal cell reference, so the epithelial cells of the mLN samples were all marked as malignant. For the epithelial cells of each mLN+ PT scRNA-seq sample, they were clustered based on the inferred CNV scores by hierarchical clustering (Ward’s method, cluster number = 10), and the CNV-based clusters with significantly fewer CNVs than the epithelial cells of mLN samples were considered as non-malignant, while the remained clusters were taken as malignant.

The malignant cells were also re-clustered following the pre-processing, PCA, UMAP reduction, clustering, and identification of marker genes work-flow. We initially maintained patient-specific biological variation by omitting batch correction, as malignant cells consistently clustered by patient origin. Based on this result, we defined Ig-high malignant cell clusters based on significant overexpression of immunoglobulin genes including IGHG1, IGHG4, IGHG2 or IGHA1 ( examined by likelihood ratio test, two-sided, Bonferroni adjusted P < 0.01 and mean_expression > 10). To systematically evaluate potential batch effects, we performed parallel analyses using harmony based batch correction on the PCA reduction and re-clustered these malignant cells, subsequently evaluating the persistence of both Ig-high/low cluster distinctions and their association with CytoTRACE-based stemness scores.

Validation of the high expressions of immunoglobulin-relevant genes in mLN+ PTs by different cohorts
The scRNA-seq data of mLN+ PT (N0 stage recorded as N1 or N2) and mLN- PT (N stage recorded as N0) samples in either De Zuani’s or Zhang’s study were respectively compared for validation. Notably, the epithelial cells were extracted based on the cell type annotations provided by the corresponding studies, and the mLN+ PTs in these studies were not included in the initial integrative analysis. Consequently, the comparisons in Supplementary Fig. 2B and S2C were independent of the comparison results in Fig. 3A. For each cohort, we checked and compared the normalized expressions of IGHA1, IGHG1, IGHG4 and IGLC2 in mLN+ PT and mLN- PT by function plot_genes_violin of Monocle3.

Integrative analysis about the metastatic Brain and mLN+ epithelial cells
The scRNA-seq data of metastatic Brain tumor tissues (10 samples, 15463 single cells) of LUAD were also obtained from GEO131907 dataset. The cell type annotations were obtained from the meta data. We extracted and merged the count matrix of the epithelial cells of the mBrain dataset and our self-generated mLN and mLN+ PT dataset. Then, we checked and compared the normalized expressions of IGHA1, IGHG1, IGHG4 and IGLC2 in mLN, mLN+ PT and mBrain by function plot_genes_violin of Monocle3.

Identification of stromal cell sub-populations
Sub-populations of the stromal cells were also identified by re-clustering the fibroblast cells following the pre-processing, PCA, batch effect correction, UMAP reduction, clustering, identifying cluster markers and annotation procedures as the above integrative analysis work-flow. The annotation was based on the defined marker genes of sub-populations of cancer associated fibroblast (CAF) cells 57 and endothelial cells (ECs) 58, including matrix CAF (MMP11, COL1A1), inflammatory CAF (iCAF) (CFD, C3), vascular CAF (vCAF) (COL18A1, NOTCH3), tumor-like CAF (tCAF) (GAPDH, ENO1), SMC (TAGLN, ACTA2), Pericyte (RGS5, NOTCH3, COX4I2), aerocytes (HPGD), arterial EC (HEY1, IGFBP3), capilary EC (CA4), EC progenitor (CD34, KDR), lymphatic (LYVE1, CCL21), venous (ACKR1) and proliferation stromal (MKI67) , as well as one cluster with specific high expression of CST1. The other clusters which showed aberrant high expressions of marker genes of the other non-stromal major cell types were removed.

Identification of myeloid cell sub-populations
Sub-populations of the myeloid cells were identified by re-clustering the myeloid cells following the pre-processing, PCA, batch effect correction, UMAP reduction, clustering, identifying cluster markers and annotation procedures as the above integrative analysis work-flow. The annotation was based on the defined marker genes 70 of macrophage (MP) (CD68, MARCO), cDC1 (CLEC9A, AIF1), cDC2 (CD1C, FCER1A, CD1B, CCL17), CD1C-DC (LST1, AIF1, FCGR3A), pDC(TCF4, IRF8, IRF7), neutrophil (Neu) (S100A8, S100A9), granulocyte (Gran) (S100A12) 127, proliferation myeloid (MKI67, STMN1) as well as the specific marker genes of the clusters including SPP1, FABP4, GPNMB and LAMP3. The other clusters which showed aberrant high expressions of marker genes of the other non-myeloid major cell types were removed.

Identification and analysis of T/NK cell sub-populations
Sub-populations of T/NK cells were identified by re-clustering the initially defined T and NK cells following the pre-processing, PCA, batch effect correction, UMAP reduction, clustering, identifying cluster markers and annotation procedures as the above integrative analysis work-flow. The cluster annotation was based on the defined marker genes 70 of CD4 naive (CCR7, SELL, TCF7, LEF1), CD4 resident effector memory (LTB, S100A4, KLRB1, IL32), CD4 effector memory (RORA, IFNG), Treg (FOXP3, CTLA4), CD8 cytotoxic (CD8A, CD8B, GZMK, CCL5, CCL4, GZMA), CD8 exhausted (CD8A, CD8B, LAG3, TIGIT, PDCD1) 70, γδT (TRDC), NK (NKG7, XCL1, FCGR3A) 128, other exhausted T (HAVCR2, TIGIT) as well as the specific marker genes including ARHGAP15, CXCL13 and ALDOA. The other clusters which showed aberrant high expressions of marker genes of the other major cell types were removed.

Identification of B cell sub-populations
Sub-populations of B cells were identified by re-clustering the initially defined B and plasma cells following the pre-processing, PCA, batch effect correction, UMAP reduction, clustering, identifying cluster markers and annotation procedures as the above integrative analysis work-flow. The cluster annotation was based on the defined marker genes 129 of native B (TCL1A, SELL, IGHD), memory B (BANK1, LTB, BLK) , activated B (CD69, CD83), germinal center (GC) B (AICDA), plasma (JCHAIN, MZB1, IGHG1, IGHG5, IGLC2), as well as the specific marker genes ZBTB20. The other clusters which also showed aberrant expressions of marker genes of the other major cell types were removed.

Spatial transcriptomics pre-processing
Raw FASTQ files and histology images were processed by sample with the Space Ranger (version spaceranger-1.2.0, 10X Genomics) software with default parameters. The filtered gene-spots matrix and the fiducial-aligned low-resolution image was used for down-streaming analyses.

Integrative analysis of ST data
In addition to the self-generated ST data about mLN and mLN+ PT samples, we also collected the public ST data of mLN- PT samples (BioStudies, https://www.ebi.ac.uk/biostudies/, accession number E-MTAB-13530) for integration. The ST count matrix as well as the spot location information were read by R package Seurat (v4.3.0) 130, and merged as one. Next, data normalization, PCA, batch effect correction and UMAP reduction were performed using the SCTransform, RunPCA, RunHarmony, RunUMAP functions of Seurat. Then, clustering of the spatial spots was performed using R package BayesSpace (v1.6.0) 35, with the harmony-corrected PCA reduction results as the input data of spatialCluster function in BayesSpace, and the number of clusters was set as q=21 based on the qTune function.

Spatial annotation of ST data
Over expressed marker genes of each spatial cluster were identified using FindAllMarkers function in Seurat, with the following thresholds: avg_log2FC> 0.1, P.adj<0.01, pct.1 > 0.3, and ranked top-8 based on log2FC. The spatial clusters were annotated based on whether the cluster marker genes included the classical marker genes of the major cell types including epithelial (SCGB3A1, KRT5, SFTPB, SFTPA1,S100A2, AKR1C1, CES1 ), fibroblast (COL1A1, COL1A2, ACTA2, FN1), myeloid (LYZ, C1QB, TYROBP), plasma (IGHG1, IGHG3, MZB1, JCHAIN) and LS-like structures (CXCL13, CCL19, CCL21, MS4A1 and TRBC2) 36, 37, 38, 39. One spatial cluster can be assigned with more than one cell type owing to the co-localization of different cell types spatially. Meanwhile, spatial clusters showing high expressions of mitochondria genes were annotated as MT-high.

Cell neighborhood analysis
Each spot was assigned with its dominant cell types according to the spatial annotation results. For the spots on each ST sample, we checked on every un-duplicated spatially adjacent spot-spot pair, and counted the number of cell-cell neighbors (CCNs) based on the dominant cell type of the two spots. For instance, if one epithelial spot is adjacent with two epithelial spots and four myeloid spots, then the epithelial-epithelial CCN adds 2, and the epitheilal-myeloid CCN adds 4. Meanwhile, there is no difference in the order of cell types, i.e., epitheilal-myeloid or myeloid-epithelial means the same CCN, and only counts once.

Cell type deconvolution analysis of the ST data
The deconvolution of the cell types on the mLN and mLN+PT ST data were performed using R package CARD (V 1.0) 131. For each ST sample, we created a CARD object using createCARDObject function, with the count matrix as well as the cell type annotations from the matched scRNA-seq sample (from the same patient and the same tissue type with the ST data) as the reference, and the ST data as the deconvolution target. Next, the deconvolution analysis was done using CARD_deconvolution function with default parameters, and each spot was assigned with the proportions of different cell types.

Determination of tumor core and peripheral partitions for epithelial spots
Based on the assumption that non-epithelial spots primarily represent non-tumor components of the microenvironment, for each epithelial spot on every ST slice, we calculated the minimum Euclidean distance between this spot and all the other non-epithelial cells using the spots’ coordinates. Here, a small distance indicates closer proximity to non-tumor regions and, by extension, far from the tumor core. Then, epithelial spots with minimum distances below the median value of all epithelial spots on the same ST slice were classified as belonging to the tumor periphery; the remaining epithelial spots were assigned to the tumor core region.

Definition of epithelial-near/far based on the distances with epithelial spots
On each ST slice, for each fibroblast/myeloid spot, we computed its minimum Euclidean distance (based on spots’ coordinates, recorded as [0,0], [0,2], [1,1] etc. ) to all epithelial spots within the same ST slice (Fig. 4M and 5J). According to the spatial coordinates of Visium, the distance between two adjacent spots in the same line equals 2, and between two adjacent spots from two different lines equals 1.414, and two-step neighbors (i.e., neighbors’ neighbors) were with min-distance <= 4. Consequently, fibroblast/myeloid spots with minimum distance <= 4 to epithelial spots were designated as epithelial-near, while the others were epithelial-far.

Cell-cell communication analysis
The ligand-receptor based cell-cell communications among malignant cells, non-malignant epithelial cells and all the sub-populations of the immune cells in the the mLN/mLN+ PT scRNA-seq data were deduced using R package CellChat (v1.6.1) 132. The list of known ligand-receptor pairs was derived from CellChatDB, a repository for ligands, receptors, and their interactions. Over-expressed signaling genes and ligand-receptor interactions were calculated using the identifyOverExpressedGenes and identifyOverExpressedInteractions functions, communication probability and cell communication network was deduced by computeCommunProb and computeComunProbPathway functions. We also utilized the NicheNet algorithm to identify potential ligands that contribute to the signature gene expression of tumor peripheral regions based on the ST data. Using the signature genes showing specific high expressions in the tumor peripheral regions than the core regions, potential ligands from the non-epithelial spots were predicted using predict_ligand_activities function.

Further clustering analysis of the myeloid/plasma+myeloid spots in mLN and mLN+ PT ST samples
We extracted the count matrix of spots annotated as myeloid or plasma+myeloid, and create a new cell_data_set object for Monocle3 analysis. Then, pre-processing and PCA, batch effect correction, UMAP reduction, clustering, and identification of cluster markers were performed using preprocess_cds, align_cds, reduce_dimension, cluster_cells and top_markers functions of Monocle3. The resulted spatial clusters were names as SCi, with i from 1 to the total cluster number.

Determination of the two-step neighbors of the lymphoid structure-like spatial spots
For each lymphoid structure (LS)-like spatial spot of each ST sample, we collected its direct neighbors, i.e., the adjacent spots, as the one-step neighbors. Typically, one spot from the Visium ST data is surrounded by six neighbors. Then, these one-step neighbors together with their adjacent spots (i.e., neighbors’ neighbors) were all taken as the two-step neighbors of the LS-like spots.

Computation of the signature scores of CSC, SPP1+IFI30+ MP and CXCL13+ T cells in the validation datasets
We defined the signature genes of CSC, SPP1+IFI30+ MP and CXCL13+ T according to the identified cellular features. The CSC relevant genes were the top-ranked genes based on the correlation coefficients with the CytoTRACE scores for LUAD and LUSC. The SPP1+IFI30+ MP gene signatures included SPP1, IFI30 as well as the other top-ranked genes showing high expressions in the IFI30-high MP cells. The CXCL13+ T cell signature genes included the top-ranked marker genes of the cell type among all sub-types (identified by top_marker function of Monocle3, top-20 ranked by specificity, p<0.001, likelihood test). These signature genes were listed in Supplementary Data 2. We computed the signature scores for the scRNA-seq/ST data from GSE127465 (scRNA-seq data of primary tumors in N0 and N1 stages of NSCLC), GSE123904 (scRNA-seq data of normal lung, primary and metastasis tumor tissues of LUAD), GSE131907 (scRNA-seq data of normal lung, primary tumor, mbrain and mLNs of LUAD ), GSE117529 (scRNA-seq data of LNs from tumor naive, parental (non-metastatic) tumor, and LN6 [LN-metastatic] tumor mice ), Zenodo ( https://zenodo.org/records/8227624, scRNA-seq data of pre/post-treatment NSCLC patients with different responses to immuno-chemotherapy treatments) and GSE190811 (ST data of mLNs in breast cancer) using the AddModuleScore function in Seurat. For the validation data, the cell type/spot annotations for each validation dataset were performed following above scRNA-seq/ST integrative analysis work-flow. When cell type information was provided alongside the matrix data, we directly utilized the pre-annotated cell types.

Clustering analysis of the TCGA-LUAD tumor samples based on bulk RNA-seq data
We performed clustering analysis of the TCGA-LUAD tumor samples based on the bulk RNA-seq data. The signature genes, including the marker genes of different cell sub-populations (Supplementary Data 3, identified by top_marker function of Monocle3, top-50 ranked and p<0.001 for each sub-population) as well as the CSC-relevant genes were considered here. Consensus clustering analysis of the bulk RNA-seq data of the signature genes were performed using ConsensusClusterPlus (v1.60.0) 133, where the expressional matrix was log2-transformed and scaled, the number of clusters was set as 4, and clustering algorithm was set as ‘pam’.

Validation of the prognostic significance of BC3 by other datasets
We trained a sub-type predictor using glmnet (v 4.1-7) 134 to classify BC3 and the other clusters based on the expression matrix for clustering analysis, and the labels of the clustering results. The parameters for glmnet were optimized through 10-fold crossover-validation. This predictor was then applied on the log2-transformed and scaled mRNA expression matrix from the other datasets to predict cluster labels. The prognostic difference between BC3 and the other clusters was examined by log-rank test and KM survival curves.

Pathway enrichment analysis
Pathway enrichment analysis was performed using R package clusterProfiler (V4.4.4). Two kinds of pathway enrichment analysis including pathway over-presentation and GSEA algorithm were utilized. The pathway over-representation algorithm was performed using enricher function, with the input gene list setting as the top-ranked up-regulated genes obtained by specific comparisons. The GSEA algorithm was performed using GSEA function, with the order ranked gene list based on the Log2-fold change of genes between two different conditions as the input. HallMark pathway and Gene ontology information was downloaded from MSigDB website (https://www.gsea-msigdb.org/gsea/msigdb/human/collections.jsp#H), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway information was extracted from the KGML files from KEGG pathway website (https://www.kegg.jp/kegg/pathway.html).

Calculation of cell transition and exhaustion scores
The cell transition and exhaustion scores were calculated using AddModuleScore function in Seurat based on the gene expression of the corresponding signature genes. The cell transition signature genes were obtained from the HallMark-EMT pathway. The exhaustion markers included CXCL13, HAVCR2, PDCD1, TIGIT, LAG3, CTLA4, LAYN, RBPJ, VCAM1, GZMB, TOX and MYO7A 135.

TCGA and GEO data collection
The normalized gene expression data in the form of RNA-Seq by Expectation Maximization(RSEM) and survival information of TCGA-LUAD, TCGA-LUSC, TCGA-LIHC, TCGA-KIRC samples were obtained by R package TCGAbiolinks (v 2.14.1) 136 using GDCquery and GDCdownload functions. GEO datasets (GSE30219 and GSE68465 ) including both the bulk transcriptomics and the patient clinical information were obtained by GEOquery (v 2.6 4.2) 137.

Survival analysis
In each cancer type, the patients were stratified into high and low groups based on the median value of CSC-score or SPP1 expression, where CSC-score was calculated by the mean normalized expressions of the CSC relevant gene signatures determined by CytoTRACE analysis. Kaplan-Meier curves were plotted to visualize survival curves and compare the prognosis differences between the high and low groups using R package survival (v3.2-13) and survminer (v0.4.9), where log-rank test was utilized to examine the significance of difference.

Statistical analysis
Additional analysis or visualization was performed using R Version 4.2.1, SPSS 19.0 and GraphPad Prism 5.0. Group comparisons were performed using Student’s t test, Wilcoxon test or likelihood ratio test. Specificity was assessed using Jensen-Shannon divergence. Correlations were assessed using the Pearson or Spearman correlation coefficient.

Data availability
The raw scRNA-seq and ST data generated in this study have been deposited at the Genome Sequence Archive at the National Genomics Data Center (Beijing, China) 138 under accession ID HRA009872 (https://ngdc.cncb.ac.cn/gsa-human/browse/HRA009872). The raw sequencing data are available under restricted access for research purposes only, access can be obtained by the DAC (Data Access Committees) of the GSA-human database. According to the guidelines of GSA-human (https://ngdc.cncb.ac.cn/gsa-human/policy), all non-profit researchers can obtain access to the data, and the Principle Investigator of any research group is allowed to apply the data. The access authority can be obtained for Research Use Only. The user can also contact the corresponding author directly. Access requests are typically processed within two weeks. Once access has been approved, the data will be available to download for 2 months. The pre-processed scRNA-seq and ST data were deposited on Science Data Bank (https://doi.org/10.57760/sciencedb.19197; https://doi.org/10.57760/sciencedb.19202). The public single-cell RNA sequencing datasets used in this study can be found in Gene Expression Omnibus (GEO) with the accession number GSE127465 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE127465), GSE123904 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE123904), GSE131907 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE131907) and GSE117529 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE117529), in BioStudies with accession number E-MTAB-13526 (https://www.ebi.ac.uk/biostudies/arrayexpress/studies/E-MTAB-13526), from Zenodo (https://zenodo.org/records/8227624) or in the LungCancer website (https://lungcancer.chenlulab.com/#/download). The public ST data of NSCLC can be found in BioStudies with accession number E-MTAB-13530 (https://www.ebi.ac.uk/biostudies/arrayexpress/studies/E-MTAB-13530). The public ST data of mLNs in breast cancer can be found in GEO with accession number GSE190811 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE190811). The ST data of nLN was from 10X Genomics website (https://www.10xgenomics.com/datasets/human-lymph-node-1-standard-1-1-0). Bulk RNA-seq data was from the cancer genome atlas (TCGA) and GEO (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE30219; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE68465). Source Data are provided with this paper.

Code availability
Code is available from github (https://github.com/diChen310/mLN_NSCLC/, DOI: http://doi.org/10.5281/zenodo.18690124) 139.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This study by Chen et al. aims to understand lymph node metastasis in non-small cell lung cancer (NSCLC) by integrating single-cell RNA-seq (scRNA-seq) and spatial transcriptomics (ST) data from matched primary tumors and metastatic lymph nodes (mLNs). Main findings include the identification of two distinct mechanisms underlying lymph node metastasis and immune evasion: one driven by cancer stemness and ferroptosis resistance, and the other mediated by malignant-cell-derived immunoglobulins at the tumor periphery. The study also reveals a key role of fibroblasts in forming a physical barrier that separates immune cells from tumor epithelial cells, while also contributing to cell transitions. Furthermore, the authors find that tumor microenvironment (TME) in both mLNs and metastatic-primary tumors (mLN⁺ PTs) are enriched with immunosuppressive cell types, including SPP1⁺IFI30⁺ macrophages, CXCL13⁺ exhausted T cells, and hyperactivated B cells. The authors use independent single-cell datasets, TCGA-LUAD bulk RNA-seq data, and functional experiments to validate the findings. This work not only provides a comprehensive view of the cellular and molecular landscape of lymph node metastasis but also offers valuable datasets to study cancer biology, immunology, and therapeutic strategy for NSCLC and potentially other tumors. The manuscript is written clearly overall. However, we have several specific comments that the authors should address to enhance the clarity of interpretation and the strength of the conclusions.\nMajor concerns:\n1. For “stemness” analysis, the authors use CytoTRACE to infer stemness scores, where a higher score reflects greater stemness or lower differentiation. However, given that Monocle3 was used for pre-processing the data and generating trajectories. If their results are consistent, they seem redundant: why not just using Monocle3 to infer pseudotime for inferring stemness? If the results are different, how to interpret the discrepancies? In other words, how robust are the CytoTRACE-derived stemness scores? Has any form of validation (e.g., correlation with known stemness markers other than CD44 and ENO1) been conducted?\n2. Cell-cell interaction analysis results need more clarification and/or validation. For example, in the cell neighborhood analysis, the number of cell-cell neighbors (CCNs) are counted based on the dominant cell type of two spots. How are dominant cell types within each spot defined? For example, if a spot contains a 50/50 mix of epithelial and myeloid cells, how is this ambiguity resolved?\n3. The authors identified interactions between CXCL14⁺ T cells and B cells, and between CAFs and LAMP3⁺ dendritic cells based on CellChat results. Given the availability of spatial transcriptomic data, is there any quantitative assessment of spatial proximity to support these inferred interactions? Similarly, on page 21, the interaction between SPP1⁺ macrophages and both malignant cells and exhausted T cells is noted. Can these interactions be spatially validated?\n4. Spatial transcriptomics data analysis requires clearer and more rigorous presentation. For example, on page 18, the authors refer to fibroblasts being “in close proximity” to epithelial cells. How is this spatial proximity quantitatively defined? Likewise, on page 21, how are spots categorized as “near” or “far” from myeloid cells? What is the spatial threshold or metric used to make this classification?\n5. Biological relevance of identified marker genes needs better explanation and discussion. For example, the analysis of fibroblasts highlights CST1⁺ fibroblasts as being specific to mLNs and mLN⁺ PTs. However, the functional relevance of CST1 is not discussed. Providing biological context or prior evidence for its significance would strengthen the rationale for focusing on this subpopulation. Similar clarification is warranted for ALDOA⁺, ARHGAP15⁺, and CXCL13⁺ T cells. Why were these genes selected for emphasis, and what is their biological or clinical importance?\nMinor concerns:\n6. On page 24, the term “other exhausted T cells” is used, and these cells appear to be specific to mLNs. Given their apparent specificity, is there anything interesting about them?\n7. Figure 1E: Spots with higher feature values are displayed with lighter colors, which is visually counterintuitive.\n8. Figure 2D: It is unclear whether the color gradient reflects CytoTRACE scores or relative gene expression. Additionally, the term \"relative expression\" is ambiguous. Does it refer to normalized expressions or something else? The legends should be clarified, since it currently refers to \"correlated gene expression.\"\n9. Figure 4H: It would be helpful to clarify what is being correlated. Is it the abundance of the fibroblast subpopulation versus the abundance of other cell types within each mLN⁺ sample, or something else?\n10. Figure 4G: The y-axis label should be \"CCN\" (cell-cell neighbors), not \"CNN.\"\n11. Figure 5D: The x-axis label should be consistently presented as \"mLN vs nLN\" instead of \"mln\".", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_3": "In this manuscript, Chen et al. present a comprehensive investigation into the mechanisms of lymph node (LN) metastasis in NSCLC using a multimodal omics approach, including single-cell RNA sequencing (scRNA-seq), spatial transcriptomics (ST), and bulk RNA sequencing. The study characterizes the tumor microenvironment (TME) of metastatic LNs (mLNs) and matched primary tumors (mLN+ PTs), providing novel insights into cancer stemness, immune evasion, and spatial organization driving metastasis. The findings offer significant translational potential for understanding NSCLC progression and developing novel therapies.\nThe work provides valuable contributions to our understanding of LN metastasis by providing a comprehensive multimodal omics resource for dissecting TME heterogeneity. This study holds potential implications for immunotherapy. Most notably, the identification of FTH1-driven ferroptosis resistance as a key driver of cancer stemness and metastasis builds upon emerging ferroptosis research while offering novel insights through its specific focus on LN metastasis and innovative integration with spatial transcriptomics data. The discovery of immunoglobulin expression by malignant cells represents another intriguing finding with potential implications for immune evasion. Furthermore, the detailed characterization of the immunosuppressive TME, particularly the mechanisms of immune exclusion, significantly advances the field and highlights promising therapeutic targets such as FTH1 and CXCL13. Overall, the study is comprehensive and well-structured by multimodal data, validation across independent datasets, as well as in vitro and in vivo experiments.\nConcerns and suggestions for improvement:\n- The potential link between ferroptosis/FTH1 and LN metastasis is compelling but remains speculative, if feasible, experimental validation would strengthen the claim.\n-The proposed functions of tumor-derived immunoglobulins also remain speculative, if feasible, functional assays could clarify their impact.\n- Given histological distinctions in NSCLC, comparing LUAD and LUSC could reveal subtype-specific metastatic mechanisms.\n-Please define key abbreviations (e.g., mLN, mLN+ PT, LN+N) in figure legends when first introduced.\n-Please clarify criteria for defining tumor core vs. periphery in ST data\n-Please Justify thresholds for gene expression cutoffs (eg. how \"immunoglobulin-high\" malignant cells were classified)\n-If feasible, may consider analyzing whether SPP1+IFI30+ macrophages or CXCL13+ T cells correlate with immunotherapy resistance in clinical cohorts.\n-In the discussion section, please consider further addressing clinical translation and therapeutic implications, as well as comparing NSCLC metastasis to other solid tumors.", "reviewer_4": "The manuscript profiles matched metastatic lymph nodes and primary tumors from 19 NSCLC patients using scRNA-seq (848,337 cells) and Visium ST (4 matched pairs), integrating public datasets to dissect cellular and spatial programs underlying LN metastasis. The authors report high stemness in mLN/mLN+ PT epithelial cells, discovery of immunoglobulin-high malignant cells expressing PD-L1, fibroblast-driven immune exclusion at tumor edges, enrichment of SPP1+IFI30+ macrophages and CXCL13+ exhausted T cells in LS-like regions, and implicate ferroptosis/FTH1 in metastasis with functional validation. The work is timely and offers an integrative cellular–spatial view of LN metastasis. The dataset is valuable, but several analytical choices, and reporting gaps limit reproducibility and interpretability.\nMajor concerns:\n-Only 4 mLN/mLN+ PT pairs underwent ST; no normal LN or additional metastatic sites were profiled spatially, limiting generalizability and statistical power for spatial claims in the paper.\n-LS-like spots are defined by co-expression of MS4A1, TRBC2, CXCL13, CCL19, CCL21 without any orthogonal validation (imaging, IHC). The methodology for “two-step neighbors” is computationally defined and may overestimate aggregation. The authors should provide histologic confirmation (such as CD20, CD3, BCL6 staining) or cite prior validation of this marker-based LS/TLS call.\n-Batch correction is explicitly skipped for malignant cells “to retain patient specificity,” which risks inflating inter-patient differences in downstream clustering and DE. The authors should show that their key findings (Ig-high clusters, stemness gradients) persist under an alternative integration (eg Harmony, etc) using additional sensitivity analysis.\n-The methodology for spatial analysis need to be clarified. For example, the edge/core classification uses median distance thresholds, which could be sample-size dependent. The authors need to provide justification.\n-Also, LS like spot identification and two-step neighbor definition are described, but criteria for LS like designation are not detailed (which gene sets? morphology?)\n-Numerous Wilcoxon tests and GSEA are reported, FDR control is stated but not consistently shown for all comparisons. Clarify all thresholds used.\nMinor points:\n-Some typos need to be corrected: “matirx” ; “contral” ; “princeple” ; “probalility” ; “up-regualted” . Please run a thorough copy edit.\n-There are issues with acronyms, please define all at first use (e.g., “mCAF”, “tCAF”, “SC7”, “LS-like”).\n-Provide exact gene lists for CSC, SPP1⁺IFI30⁺ MP, and CXCL13⁺ T signatures used in validations for reproducibility."}, "author_response": "Here, according to the Editor and Reviewers’ important suggestions and comments, we have revised carefully and have made corrections which we hope to meet with approval. We attached the response letter and revised manuscript and the corrections are colored in red to make the revisions notable. The “Calibri” fonts cited the Editor and Reviewers’ comments, and the “Times new roman” and “Blue” fonts are our replies. The “Red” fonts are cited in the manuscript for helping the editor and reviewers easily read. If any question arises, please do not hesitate to let us know.\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author): bioinformatic analysis of omics data including ST\nThis study by Chen et al. aims to understand lymph node metastasis in non-small cell lung cancer (NSCLC) by integrating single-cell RNA-seq (scRNA-seq) and spatial transcriptomics (ST) data from matched primary tumors and metastatic lymph nodes (mLNs). Main findings include the identification of two distinct mechanisms underlying lymph node metastasis and immune evasion: one driven by cancer stemness and ferroptosis resistance, and the other mediated by malignant-cell-derived immunoglobulins at the tumor periphery. The study also reveals a key role of fibroblasts in forming a physical barrier that separates immune cells from tumor epithelial cells, while also contributing to cell transitions. Furthermore, the authors find that tumor microenvironment (TME) in both mLNs and metastatic-primary tumors (mLN⁺ PTs) are enriched with immunosuppressive cell types, including SPP1⁺IFI30⁺ macrophages, CXCL13⁺ exhausted T cells, and hyperactivated B cells. The authors use independent single-cell datasets, TCGA-LUAD bulk RNA-seq data, and functional experiments to validate the findings. This work not only provides a comprehensive view of the cellular and molecular landscape of lymph node metastasis but also offers valuable datasets to study cancer biology, immunology, and therapeutic strategy for NSCLC and potentially other tumors. The manuscript is written clearly overall. However, we have several specific comments that the authors should address to enhance the clarity of interpretation and the strength of the conclusions.\nRe: Thanks a lot for the affirmation of our work and the valuable suggestions. We have revised the manuscript according to the comments and suggestions.\nMajor concerns:\n1. For “stemness” analysis, the authors use CytoTRACE to infer stemness scores, where a higher score reflects greater stemness or lower differentiation. However, given that Monocle3 was used for pre-processing the data and generating trajectories. If their results are consistent, they seem redundant: why not just using Monocle3 to infer pseudotime for inferring stemness? If the results are different, how to interpret the discrepancies? In other words, how robust are the CytoTRACE-derived stemness scores? Has any form of validation (e.g., correlation with known stemness markers other than CD44 and ENO1) been conducted?\nRE: We appreciate the reviewer’s insightful question regarding the relationship between CytoTRACE-derived stemness scores and Monocle3 pseudotime. Below, we clarify our rationale and present evidence supporting the robustness of our approach:\nWhile Monocle3 is powerful for trajectory inference, its pseudotime calculation requires user-specified root nodes (as noted in the Monocle3 tutorial, https://cole-trapnell-lab.github.io/monocle3/docs/trajectories/#learn-graph). This introduces subjectivity, as the choice of root nodes (often based on prior knowledge or assumptions) directly influences pseudotime values. In contrast, CytoTRACE objectively infers stemness solely from gene expression signatures, without requiring user input on root nodes.\nTo address the reviewer’s concern, we explicitly tested the correlation between CytoTRACE scores and Monocle3 pseudotime, where the cells with high CytoTRACE scores were selected as the root nodes. The results were displayed as below:\nFigure R1_1. Umap visualization of the CytoTRACE and Pseudotime values of the epithelial cells in LUAD and LUSC respectively.\nWe calculated the correlations between CytoTRACE scores and the pseudotime values:\nFor LUAD, Spearman’s ρ = -0.40 (p < 2.2e-16).\nFor LUSC, Spearman’s ρ = -0.38 (p < 2.2e-16).\nThe significant negative correlations confirm that cells with high CytoTRACE scores (greater stemness) align with early pseudotime states (root nodes), as expected biologically.\nThough consistent, CytoTRACE was chosen as the primary metric of cell stemness because the following reasons:\n1) CytoTRACE does not require manually selected root nodes\n2) CytoTRACE specializes in stemness prediction (unlike Monocle3’s broader trajectory focus),\n3) Monocle3 here was mainly applied for the pre-processing, reduction and clustering analysis of the scRNA-seq data, since a benchmark study for scRNA-seq analysis tools demonstrated that Monocle3 also displayed good performance for the general clustering analysis of the scRNA-seq data [1].\n4) We further validated CytoTRACE scores by confirming their agreement with known stemness markers (e.g., CD44, EpCAM, CD166), see Figure 2E and Page 8-9, Line 155-157, as below:\nCD44, EpCAM and CD166 (gene names CD44, EPCAM and ALCAM) are widely recognized markers of CSCs in NSCLC [2, 3]. Here, epithelial cells exhibiting high CytoTRACE scores also showed a significant higher expression of CD44, EPCAM and ALCAM (Figure 2E).\nFigure 2E. Box plot showing the differential expressions of stemness marker genes CD44, EPCAM and ALCAM in epithelial cells with high or low CytoTRACE scores in LUAD and LUSC respectively. High and low is separated by the median value of CytoTRACE scores. ****p < 0.0001 by Wilcox-test, unpaired, two-sided.\nReferences:\n1. Yu, L., et al., Benchmarking clustering algorithms on estimating the number of cell types from single-cell RNA-sequencing data. Genome Biol, 2022. 23(1): p. 49.\n2. Zakaria, N., et al., Human non-small cell lung cancer expresses putative cancer stem cell markers and exhibits the transcriptomic profile of multipotent cells. BMC Cancer, 2015. 15: p. 84.\n3. Zhang, W.C., et al., Glycine decarboxylase activity drives non-small cell lung cancer tumor-initiating cells and tumorigenesis. Cell, 2012. 148(1-2): p. 259-72.\n2. Cell-cell interaction analysis results need more clarification and/or validation. For example, in the cell neighborhood analysis, the number of cell-cell neighbors (CCNs) are counted based on the dominant cell type of two spots. How are dominant cell types within each spot defined? For example, if a spot contains a 50/50 mix of epithelial and myeloid cells, how is this ambiguity resolved?\nRE: We appreciate the reviewer's insightful question regarding the definition of dominant cell in spatial spots. Below we provide a detailed clarification of our methodology and the relevant results.\nThe definition of dominant cell types for the spots were based on the spatial clustering results. Here, the dominant cell types of each spatial cluster was annotated based on whether the cluster showed significantly higher expressions of classical marker genes of the major cell types. A spatial cluster could be assigned more than one cell type (e.g., Plasma+Myeloid) if markers from multiple cell types met the thresholds, reflecting the spatial co-localization of different cell types. The spatial annotation results are shown in Figure 1D and S1D.\nWe have clarify the relevant definition and analysis methods in the Methods sections, in Page 51, Line 1043-1052 , as below:\nSpatial annotation of SD data\nOver expressed marker genes of each spatial cluster were identified using FindAllMarkers function in Seurat, with the following thresholds: avg_log2FC> 0.1, P.adj<0.01, pct.1 > 0.3, and ranked top-8 based on log2FC. The spatial clusters were assigned to major cell types based whether the cluster over-expressed classical marker genes of the major cell types or LS-like structures: epithelial (SCGB3A1, KRT5, SFTPB, SFTPA1,S100A2, AKR1C1, CES1 ), fibroblast (COL1A1, COL1A2, ACTA2, FN1), myeloid (LYZ, C1QB, TYROBP), plasma (IGHG1, IGHG3, MZB1, JCHAIN) and LS-like structures (CXCL13, CCL19, CCL21, MS4A1 and TRBC2). One spatial cluster can be assigned with more than one cell type owing to the co-localization of different cell types spatially.\nWe also added explanations of the dominant cell types in the results section, Page 6, Line 108-113 and Figure S1D, as below:\nFor example, spatial cluster SC1 showed elevated expressions of SFTPC, SCGB3A1, SCGB3A2, and SFTPB - all classical marker genes of lung epithelial cells - and was therefore annotated as epithelial-dominant (Figure S1D). For spatial clusters SC18 and SC19, which showed co-expressions of marker genes for both myeloid (LYZ, C1QB) and plasma cells (IGHG1, JCHAIN), the dominant cell types were defined as Myeloid+Plasma (Figure S1D), implying the co-localization of myeloid cells and plasma cells.\n3. The authors identified interactions between CXCL13⁺ T cells and B cells, and between CAFs and LAMP3⁺ dendritic cells based on CellChat results. Given the availability of spatial transcriptomic data, is there any quantitative assessment of spatial proximity to support these inferred interactions? Similarly, on page 21, the interaction between SPP1⁺ macrophages and both malignant cells and exhausted T cells is noted. Can these interactions be spatially validated?\nRE: Thank you for this important question regarding spatial validation of cell-cell interactions. Since the resolution of ST data here is not single-cell level, we utilized CARD method to estimate the abundance of different cell sub-populations within each spatial spot, and then investigated the spatial relationships.\nFor CXCL13+ T, we calculated the summed abundance scores of CXCL13+ T cells and each interacting cell type (including LAMP3+ DC, GC B,activated B, iCAF and mCAF cells) to quantitatively assess their co-localization potential. As shown in Figures S7E and S7F, these interacting partners showed significantly higher co-localization scores with CXCL13+ T cells within LS-like spots compared to other regions, supporting the CellChat-predicted interactions in a spatially resolved context.\nFigure S7(E-F)\n(E) Spatial distribution of CXCL13+ T cells and their interacting cells in P6T (mLN+ PT) and P6L (mLN) samples. The cell abundances were deconvoluted using CARD.\n(F) Violin plot showing the summed scores (derived from CARD deconvolution) between CXCL13+ T cells and their interacting cell types based on the ST data of mLN+ PT and mLN samples. Wilcox-test, un-paired, two-sided. ****P<0.0001.\nExplanations have been added in Page 25, Line 479-482, as below:\nCARD-based deconvolution analysis of the ST data also showed the co-localization of the CXCL13+ T cells with the interacting cells like LAMP3+ DC, GC B, activated B, iCAF and mCAF cells in the LS-like spots (Figures S7E and S7F).\nSimilarly, to evaluate spatial interactions involving SPP1+ MP cells and their interaction cells identified using scRNA-seq data, like malignant cells and different T cell sub-populations, we visualized their spatial distributions (Figure S5H) and computed the Spearman correlation coefficients between the cell abundance scores (Figure S5I), as below:\nFigure S5(H-I)\n(H) Spatial distribution of SPP1+ MP, malignant epithelial cell and CD8 Exhausted T cells in mLN+ PT (P3T, P5T) and mLN (P3L, P5L) samples. The cell abundances were deconvoluted using CARD.\n(I) Heatmap showing the Spearman correlations of the CARD-based cell abundance among different cell sub-populations.\nWe observed a clear spatial dissociation between SPP1+ MP and malignant cells (Figure S5H), suggesting that their interaction may occur through long-range signaling. In contrast, SPP1+ MP showed notable co-localization with T cell subsets, indicative of short-range interactions or direct cell-cell contact (Figure S5H and S5I).\nThe results have been described in Page 20, Line 395-399, as below:\nSpatial examination of cellular distributions indicated a non-co-localization pattern between SPP1+ MP and malignant cells, suggesting potential long-range interactions mediated by paracrine or endocrine signaling (Figures S5H and S5I). In contrast, SPP1+ MP and T cell sub-populations exhibited co-localization, indicating direct short-range cellular crosstalk (Figures S5H and S5I).\n4. Spatial transcriptomics data analysis requires clearer and more rigorous presentation. For example, on page 18, the authors refer to fibroblasts being “in close proximity” to epithelial cells. How is this spatial proximity quantitatively defined? Likewise, on page 21, how are spots categorized as “near” or “far” from myeloid cells? What is the spatial threshold or metric used to make this classification?\nRE: We appreciate the reviewer's request for clarification regarding our spatial proximity analyses. We have revised the manuscript to provide more rigorous descriptions of the methods used to define cellular proximity relationships under different conditions.\nFor fibroblast-epithelial cell proximity (Figure 4G),\nWe compared the average number of cell-cell neighbors (CCNs) between epithelial spots and different kinds of spatial spots. The significantly higher CCNs between fibroblast and epithelial spots (p<0.05 by Wilcoxon test) quantitatively supports their spatial co-localization. We have revised the relevant text to explicitly state in Page 17, Line 318, as below:\nST-based analysis revealed that there were significantly higher CCNs between fibroblast and epithelial spots (p<0.05 by Wilcoxon test), suggesting that fibroblasts are with closer proximity to the epithelial cells than the other cells, particularly in the mLN+ PTs (Figure 4G)\nFor myeloid cell distance classification (Figure 5J),\nThe \"near\" versus \"far\" categorization of myeloid spots relative to epithelial spots was determined as:\nFor each myeloid spot, we computed its minimum distance (based on spots’ relative coordinates, recorded as [0,0],[0,2], [1,1] etc. ) to all epithelial spots within the same tissue section. According to the spatial arrangement principle of Visium, the distance between two adjacent spots in the same line equals 2, and between two adjacent spots from two different lines equals 1.414, and two-step neighbors were with min-distance <= 4, as shown in the figure below:\nFigure R1_2. Schematic of the relative coordinates and distances of spatial spots in Visium ST data.\nConsequently, myeloid spots with minimum distance <= 4 from epithelial spots on the same ST slice were designated as \"epithelial-near\", while the other myeloid spots were \"epithelial-far\". The Methods section now includes full computational details in Page 53, Line 1076-1083, as below:\nDefinition of epithelial-near/far based on the distances with epithelial spots\nOn each ST slice, for each fibroblast/myeloid spot, we computed its minimum Euclidean distance (based on spots’ coordinates, recorded as [0,0], [0,2], [1,1] etc.) to all epithelial spots within the same ST slice (Figures 4M and 5J). According to the spatial coordinates of Visium, the distance between two adjacent spots in the same line equals 2, and between two adjacent spots from two different lines equals 1.414, and two-step neighbors (i.e., neighbors’ neighbors) were with min-distance <= 4. Consequently, fibroblast/myeloid spots with minimum distance <= 4 to epithelial spots were designated as epithelial-near, while the others were epithelial-far.\n5. Biological relevance of identified marker genes needs better explanation and discussion. For example, the analysis of fibroblasts highlights CST1⁺ fibroblasts as being specific to mLNs and mLN⁺ PTs. However, the functional relevance of CST1 is not discussed. Providing biological context or prior evidence for its significance would strengthen the rationale for focusing on this sub-population. Similar clarification is warranted for ALDOA⁺, ARHGAP15⁺, and CXCL13⁺ T cells. Why were these genes selected for emphasis, and what is their biological or clinical importance?\nRE: We agree that providing biological context is crucial for interpreting our findings. The marker genes (CST1, ALDOA, ARHGAP15, CXCL13) were selected not only based on their specific expression in distinct population (see the figure below, see also Figure S4A for CST1, Figure S7A for ALDOA, ARHGAP15 and CXCL13), but, more importantly, because their known biological functions provide direct mechanistic hypotheses for the stromal modulation and immunosuppressive microenvironment in mLNs and mLN+ PTs. We have now expanded the discussion to include this rationale.\nFigure R1_3. Violin plot showing the differential expressions of genes across different cell sub-populations.\nCST1 (Cystatin 1) is a protease inhibitor. CST1 overexpression has been observed in various cancer types [4-6]. Its role in inhibiting cysteine proteases suggests a potential function in modulating the extracellular matrix. The CST1+ fibroblasts have also been identified in tumor tissues of esophageal squamous cell carcinoma, and show prognostic values [6]. Here, the CST1+ fibroblasts were specific to mLNs and mLN+ PTs in LUAD. (Page 16, Line 296-300)\nALDOA, or Aldolase A, is an enzyme crucial for glycolysis and gluconeogenesis which are important for maintaining cellular energy balance, and ALDOA expression may influence T cell behavior [7], potentially affecting their ability to effectively target and eliminate cancer cells [8]. The ALDOA+ T cells were predominantly enriched by nLNs, but significantly reduced in mLNs (Figure 6B), strongly suggests a state of metabolic impairment or dysfunction, contributing to the ineffective anti-tumor response in mLNs. (Page 24, Line 457-461)\nARHGAP15, a gene encoding a Rho GTPase-activating protein, can inhibit Rac1 activity [9]. Rac1 is a small GTPase involved in various cellular processes, including T cell activation and migration [10]. ARHGAP15+ T cells were more enriched in mLNs and mLN+ PTs (Figure 6B), indicating suppression of T cell activation and motility, offering a plausible mechanism for their trapped and dysfunctional state. (Page 24, Line 462-465)\nReferences:\n4. Dai, D.N., et al., Elevated expression of CST1 promotes breast cancer progression and predicts a poor prognosis. J Mol Med (Berl), 2017. 95(8): p. 873-886.\n5. Zhang, L., et al., CST1 promoted gastric cancer development by activating the AKT pathway. Clinics (Sao Paulo), 2025. 80: p. 100561.\n6. Dinh, H.Q., et al., Integrated single-cell transcriptome analysis reveals heterogeneity of esophageal squamous cell carcinoma microenvironment. Nature Communications, 2021. 12(1).\n7. Chang, Y.C., et al., Roles of Aldolase Family Genes in Human Cancers and Diseases. Trends Endocrinol Metab, 2018. 29(8): p. 549-559.\n8. Guo, M., et al., Molecular, metabolic, and functional CD4 T cell paralysis in the lymph node impedes tumor control. Cell Rep, 2023. 42(9): p. 113047.\n9. Costa, C., et al., The RacGAP ArhGAP15 is a master negative regulator of neutrophil functions. Blood, 2011. 118(4): p. 1099-1108.\n10. Faroudi, M., et al., Critical roles for Rac GTPases in T-cell migration to and within lymph nodes. Blood, 2010. 116(25): p. 5536-5547.\nMinor concerns:\n6. On page 24, the term “other exhausted T cells” is used, and these cells appear to be specific to mLNs. Given their apparent specificity, is there anything interesting about them?\nRE: The number of this sub-population of T cells is small, showing up-regulations of two exhaustion relevant marker genes, HAVCR2, TIGIT, however, low expression of CXCL13 and CD8A/CD8B, making it different from CXCL13+ T and CD8+ exhausted T cells (Figure R1_4, see also Figure S6A), consequently, we termed this sub-population as other exhausted T. They may represent a rare unconventional T cell subset that has acquired an exhausted phenotype within the immunosuppressive mLN microenvironment. The small cell number precluded extensive validation, we did not do other in-further analysis.\nFigure R1_4. Violin plot showing the differential expressions of some representative genes across different T cell sub-populations.\n7. Figure 1E: Spots with higher feature values are displayed with lighter colors, which is visually counterintuitive.\nRE: We have revised the colors,\nFigure 1F (original Figure 1E). Heatmap showing the average number of cell-cell neighbors (avg_CCNs) in different ST samples. One CCN refers to a pair of spatially adjacent spots, and the CCN type is determined based on the annotations of the two spots.\n8. Figure 2D: It is unclear whether the color gradient reflects CytoTRACE scores or relative gene expression. Additionally, the term \"relative expression\" is ambiguous. Does it refer to normalized expressions or something else? The legends should be clarified, since it currently refers to \"correlated gene expression.\"\nRE: To be more concise, we have revised the “relative expression” into “scaled gene expression”, and annotated in Figure 2D, as below\nFigure 2D. Heatmap showing CytoTRACE scores (first line) and the correlated gene expressions of epithelial cells in LUAD and LUSC samples (top-20 genes based on Pearson correlations). Only 500 randomly-sampled cells are visualized for each tumor type. Each column represents one cell, and the cells from left to right are ranked based on the CytoTRACE scores. The heatmap displays min-max scaled gene expression values: raw read counts were first log2-transformed and then scaled to a [0,1] range using min-max normalization.\n9. Figure 4H: It would be helpful to clarify what is being correlated. Is it the abundance of the fibroblast sub-population versus the abundance of other cell types within each mLN⁺ sample, or something else?\nRE: We have revised the figure legend:\nFigure 4H. Heatmap showing the Spearman correlation coefficients between the proportions of different fibroblast sub-populations and other major cell types in each mLN+ sample based on the scRNA-seq data.\n10. Figure 4G: The y-axis label should be \"CCN\" (cell-cell neighbors), not \"CNN.\"\nRE: We have revised this.\n11. Figure 5D: The x-axis label should be consistently presented as \"mLN vs nLN\" instead of \"mln\".\nRE: We have revised this.\nReviewer #2 (Remarks to the Author): an ECR co-reviewer\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nReviewer #3 (Remarks to the Author): lung cancer metastasis and immune evasion\nIn this manuscript, Chen et al. present a comprehensive investigation into the mechanisms of lymph node (LN) metastasis in NSCLC using a multimodal omics approach, including single-cell RNA sequencing (scRNA-seq), spatial transcriptomics (ST), and bulk RNA sequencing. The study characterizes the tumor microenvironment (TME) of metastatic LNs (mLNs) and matched primary tumors (mLN+ PTs), providing novel insights into cancer stemness, immune evasion, and spatial organization driving metastasis. The findings offer significant translational potential for understanding NSCLC progression and developing novel therapies.\nThe work provides valuable contributions to our understanding of LN metastasis by providing a comprehensive multimodal omics resource for dissecting TME heterogeneity. This study holds potential implications for immunotherapy. Most notably, the identification of FTH1-driven ferroptosis resistance as a key driver of cancer stemness and metastasis builds upon emerging ferroptosis research while offering novel insights through its specific focus on LN metastasis and innovative integration with spatial transcriptomics data. The discovery of immunoglobulin expression by malignant cells represents another intriguing finding with potential implications for immune evasion. Furthermore, the detailed characterization of the immunosuppressive TME, particularly the mechanisms of immune exclusion, significantly advances the field and highlights promising therapeutic targets such as FTH1 and CXCL13. Overall, the study is comprehensive and well-structured by multimodal data, validation across independent datasets, as well as in vitro and in vivo experiments.\nRE: We sincerely thank the reviewer for their positive assessment of our work and their insightful comments regarding the significance and implications of our findings.\nWe have revised the manuscript according to the concerns and suggestions.\nConcerns and suggestions for improvement:\n- The potential link between ferroptosis/FTH1 and LN metastasis is compelling but remains speculative, if feasible, experimental validation would strengthen the claim.\nRE: Thank you for your suggestion. During the revision period, we attempted to examine FTH1 and ferroptosis-related protein expression, as well as perform pathological analyses in human lymph node (LN) metastatic samples. However, obtaining these samples required approval from the institutional research ethics committee, and the ethical review process is expected to take longer than the revision timeframe. Moreover, direct observation of ferroptosis or FTH1 regulation in murine LN metastatic samples proved impractical. According to assessment of our pathologists, it is difficult to accurately delineate metastatic regions within mouse lung tissues.\nIn the initial submission, we conducted an in vivo metastasis assay, where FTH1 knockdown cells formed markedly fewer metastatic lung nodules compared to control cells, without significant changes in body weight (original Figure 7I and Figure S8J). Furthermore, histological examination of metastatic lung samples by H&E staining revealed extensive metastatic lesions in the control group (new Figure 7J in the revised manuscript). Although we were unable to directly visualize ferroptosis/FTH1 regulatory events in LN metastasis, our collective data support that FTH1 and ferroptosis may play a significant role in promoting both lung and LN metastasis.\nFigure 7J.\nCells transfected with control (shControl) or shRNA against FTH1(shFTH1_1) were injected into the null mice with tail vein for 10 weeks. H&E staining of the lungs with tail vein injected shControl or shFTH1_1 cells. Scale bar, 600 μM. n=3 biologically independent samples.\n-The proposed functions of tumor-derived immunoglobulins also remain speculative, if feasible, functional assays could clarify their impact.\nRE: Thank you for your suggestion. We have generated IGHG4 knockdown cell lines by two independent IGHG4 short hairpin RNAs (shRNAs) in A549 and H460 cells. And we found that IGHG4 knockdown reduced the several EMT marker genes expression of Vimentin, N-cadherin, MMP9 and ZEB2 in A549 and H460 cells (new Figure S10F in revised manuscript). However, IGHG4 knockdown did not influenced A549 and H460 cell proliferation (New Figure S10H in revised manuscript). Additionally, IGHG4 knockdown decreased cell invasion in both A549 and H460 cells (New Figure S10L in revised manuscript). Taken together, these findings suggest that IGHG4 may play a crucial regulatory role in cell invasion and metastasis.\nFigure S10(F, H, L).\n(F) Cells were transfected with control or shRNAs against IGHG4. Immunoblotting analysis of IGHG4, N-cad, Vimentin, MMP9, and ZEB2 protein levels in A549 and H460 cells. ACTIN was used as a loading control.\n(H) Cells transfected with control or shRNAs against IGHG4 were suspended in 100 μL of RPMI 1640 medium seeded in 96-well plates at a density of 0.6×103 cells overnight. Relative cell viability was assessed by the CCK8 assay in A549 and H460 cells (n=3). Error bars indicate mean ± SD. One-way ANOVA test.\n(L) Cells were transfected with control or shRNA against IGHG4. Representative images (left panel) and data quantification (right panel) of transwell invasion assay in A549 and H460 cells (n=3). Unpaired two-sided Student’s t-test. Scale bar, 100 μm.\n- Given histological distinctions in NSCLC, comparing LUAD and LUSC could reveal subtype-specific metastatic mechanisms.\nRE: Previous analyses have identified differences between LUAD and LUSC in terms of cell types such as epithelial cells (Figure 2B, LUAD and LUSC were respectively dominated by AT2 and Basal cells) and CST1+ fibroblast (Figure 4B, CST1+ fibroblast cells were specifically enriched in LUAD). We have supplemented the relevant comparative analyses between LUAD and LUSC to address this important question. Our comparative analysis of LUAD and LUSC has uncovered both conserved and subtype-specific cellular features during lymph node metastasis. The results have been shown in Figure S9, as below:\nFigure S9. Barplot showing changes in cell composition of mLN+ PTs compared to normal lung or mLN- PTs in LUAD and LUSC respectively. Cell compositions were determined based on proportions of all cell sub-populations identified in scRNA-seq data. For each cell sub-population, Log2FC values were calculated for the comparisons between mLN+ PTs and normal lung (mLN+ PT vs normal lung) and between mLN+ PTs and mLN- PTs (mLN+ PT vs mLN- PT). Only statistically significant changes (Wilcoxon test, unpaired, two-sided; P < 0.01) are shown as bars.\nThe results have been described in Page 26-27, Line 516-516, as below:\nDistinct cellular mechanisms involved in lymph node metastasis between LUAD and LUSC\nTo further investigate the differences between LUAD and LUSC regarding LN metastasis, we quantified cell composition changes of mLN+ PT comparing to either normal lung or mLN- PT tissues in LUAD and LUSC samples respectively. Both conserved and specific cellular features were identified (Figure S9). Notably, mLN+ PTs of LUAD displayed significant accumulation of ARHGAP15+ T, CXCL13+ T, plasma and Treg cells. In contrast, mLN+ PTs of LUSC displayed depletion of various effector immune cells, particularly CD4 effector T cells and NK cells (Figure S9). CAFs, such as vCAF and tCAF, were consistently enriched by the mLN+ PT samples in both LUAD and LUSC, suggesting a similar stromal response. These divergent immune landscapes suggest fundamentally distinct metastatic mechanisms - where LUAD is more likely to develop an immunosuppressive niche through Treg and exhausted lymphcytes, while LUSC appears to evade immunity more by immune cell exclusion.\n-Please define key abbreviations (e.g., mLN, mLN+ PT, LN+N) in figure legends when first introduced.\nRE: We have added the full names in the figure legend of Figure 1:\nmLN: metastatic lymph node; mLN+ N: mLN-positive adjacent normal lung; mLN+ PT: mLN-positive primary tumor; mLN- N: mLN-negative adjacent normal lung; mLN- PT: mLN-negative primary tumor; nLN: normal lymph node.\n-Please clarify criteria for defining tumor core vs. periphery in ST data\nRE: We have revised the method section in Line ** ,as below:\nDetermination of tumor core and peripheral partitions for epithelial spots\nBased on the assumption that non-epithelial spots primarily represent non-tumor components of the microenvironment, for each epithelial spot on every ST slice, we calculated the minimum Euclidean distance between this spot and all the other non-epithelial cells using the spots’ coordinates. Here, a small distance indicates closer proximity to non-tumor regions and, by extension, far from the tumor core. Then, epithelial spots with minimum distances below the median value of all epithelial spots on the same ST slice were classified as belonging to the tumor periphery; the remaining epithelial spots were assigned to the tumor core region.\n-Please Justify thresholds for gene expression cutoffs (eg. how \"immunoglobulin-high\" malignant cells were classified)\nRE: we have revised the figure legend of Figure 3F:\nFigure 3(F) Dot plot showing the mean expression of classical marker genes of AT2/basal and immunoglobulin relevant genes in different clusters of malignant cells. Clusters demonstrating significantly elevated expression of immunoglobulin genes (IGHG4, IGHG1, IGHA1, or IGHG2; likelihood ratio test, Benjamini–Hochberg adjusted P < 0.01, mean expression > 10) are highlighted in red and designated as Ig-high clusters; all others are classified as Ig-low clusters.\nIg-high: immunoglobulin-high; Ig-low: immunoglobulin-low.\n-If feasible, may consider analyzing whether SPP1+IFI30+ macrophages or CXCL13+ T cells correlate with immunotherapy resistance in clinical cohorts.\nRE: This is a very valuable suggestion. To preliminarily explore the potential role of SPP1+IFI30+ macrophages and CXCL13+ T cells in immunotherapy resistance, we analyzed a publicly available scRNA-seq dataset from a previous immuno-chemotherapy study in NSCLC [1]. In this dataset, both post and pre treatment samples, with different responses to treatment were collected. We calculated the IFI30+SPP1+ MP and CXCL13+ T cell signature scores in the macrophage and T cells respectively. As the results show, both signatures showed elevated scores in pre-treatment samples and in non-responders, suggesting a possible association with treatment resistance (Figure S11, see below). We have now included this exploratory analysis as a new supplementary figure (Figure. S11) and added a brief description in the Results section (Page 31, Line 603-606) to highlight the potential clinical relevance of our findings, as below:\nMoreover, analysis of an independent NSCLC cohort [1] revealed that elevated gene signature scores of SPP1+ IFI30+ macrophages and CXCL13+ T cells in pre-treatment non-responders to immuno-chemotherapy (Figures S11A and S11B), suggesting their associations with resistance to immuno-chemotherapy.\nFigure S11. Investigation of IFI30+ SPP1+ MP and CXCL13+ T cell signatures in an immuno-chemotherapy NSCLC cohort.\n(A-B) Distributions of IFI30+ SPP1+ macrophage (A) and CXCL13+ T cell (B) signature scores in the macrophage/T cells from the pre and post- treatment samples with different responses to immuno-chemotherapy treatments. PostMPR: Post-treatment, major pathologic response; PostNMPR: Post-treatment, non-major pathologic response; PreMPR:Pre-treatment, major pathologic response; PreNMPR:Pre-treatment, non-major pathologic response. The preprocessed ST data were obtained from Zenodo (https://zenodo.org/records/8227624).\nReference:\n1. Yan, Y., et al., Multi-omic profiling highlights factors associated with resistance to immuno-chemotherapy in non-small-cell lung cancer. Nat Genet, 2025. 57(1): p. 126-139.\n-In the discussion section, please consider further addressing clinical translation and therapeutic implications, as well as comparing NSCLC metastasis to other solid tumors.\nRE: we have added further discussion:\nPage 34, Line 669-672:\nThese findings carry important clinical implications, as FTH1 represents a promising therapeutic target. Its inhibition could disrupt iron metabolism and sensitize tumor cells to ferroptosis inducers, particularly against metastatic CSCs.\nPage 36, Line 703-705:\nThe prominent roles of SPP1+IFI30+ macrophages and CXCL13+ T cells in mediating immunosuppression also highlight their potential as novel immunotherapy targets for countering immune evasion in metastatic NSCLC.\nPage 36, Line 706-713:\nInterestingly, the immune evasion mechanisms we identified in NSCLC extend beyond this malignancy, reflecting broader relevance across diverse cancer types. For example, SPP1+ MP have been implicated in immune suppression and metastasis in melanoma [2], breast [3], and colorectal cancer [4]. Similarly, CXCL13+ T cells—often associated with tertiary lymphoid structures and immune regulation—have been observed in breast cancer [5, 6] and gastric cancer [7], suggesting a potential conserved mechanism of immune modulation across tumors. Moreover, ferroptosis resistance mediated by FTH1 or lipid metabolism has been reported in liver and ovarian cancers [8, 9], suggesting the broad relevance of metabolic adaptation in metastasis.\nReferences:\n1. Wei, C.Y., et al., Delineating the early dissemination mechanisms of acral melanoma by integrating single-cell and spatial transcriptomic analyses. Nature Communications, 2023. 14(1).\n3. Cha, S.M., et al., SPP1+ macrophages in HR+ breast cancer are associated with tumor-infiltrating lymphocytes. NPJ Breast Cancer, 2024. 10(1): p. 83.\n4. Wu, Y.C., et al., Spatiotemporal Immune Landscape of Colorectal Cancer Liver Metastasis at Single-Cell Level. Cancer Discovery, 2022. 12(1): p. 134-153.\n5. Liu, T., et al., Single cell profiling of primary and paired metastatic lymph node tumors in breast cancer patients. Nat Commun, 2022. 13(1): p. 6823.\n6. Zhang, Y., et al., Single-cell analyses reveal key immune cell subsets associated with response to PD-L1 blockade in triple-negative breast cancer. Cancer Cell, 2021. 39(12): p. 1578-1593 e8.\n7. Chen, J., et al., Single-Cell Profiling of Tumor Immune Microenvironment Reveals Immune Irresponsiveness in Gastric Signet-Ring Cell Carcinoma. Gastroenterology, 2023. 165(1): p. 88-103.\n8. Zhang, Q., et al., Iron promotes ovarian cancer malignancy and advances platinum resistance by enhancing DNA repair via FTH1/FTL/POLQ/RAD51 axis. Cell Death Dis, 2024. 15(5): p. 329.\n9. Zhang, X.Y., et al., CircPIAS1 promotes hepatocellular carcinoma progression by inhibiting ferroptosis via the miR-455-3p/NUPR1/FTH1 axis. Mol Cancer, 2024. 23(1): p. 113.\nReviewer #4 (Remarks to the Author): ST, TME, TIME, computational biology/bioinformatics\nThe manuscript profiles matched metastatic lymph nodes and primary tumors from 19 NSCLC patients using scRNA-seq (848,337 cells) and Visium ST (4 matched pairs), integrating public datasets to dissect cellular and spatial programs underlying LN metastasis. The authors report high stemness in mLN/mLN+ PT epithelial cells, discovery of immunoglobulin-high malignant cells expressing PD-L1, fibroblast-driven immune exclusion at tumor edges, enrichment of SPP1+IFI30+ macrophages and CXCL13+ exhausted T cells in LS-like regions, and implicate ferroptosis/FTH1 in metastasis with functional validation. The work is timely and offers an integrative cellular–spatial view of LN metastasis. The dataset is valuable, but several analytical choices, and reporting gaps limit reproducibility and interpretability.\nRE: We sincerely thank the reviewer for their thorough evaluation of our work and their constructive feedback. We have revised the manuscript according to the suggestions.\nMajor concerns:\n-Only 4 mLN/mLN+ PT pairs underwent ST; no normal LN or additional metastatic sites were profiled spatially, limiting generalizability and statistical power for spatial claims in the paper.\nRE: We sincerely thank the reviewer for raising this important point regarding the sample size and spatial profiling scope in our study. We fully acknowledge that the limited number of spatial samples (n=4 mLN/mLN+ PT pairs) will affect the generalizability of our spatial findings.\nTo address this concern and enhance the robustness of our conclusions, we performed an independent validation using publicly available spatial transcriptomics dataset of breast cancer lymph node metastases and normal lymph node (nLN). Notably, we observed consistent spatial features in breast cancer samples, including: recurrent presence of LS-like structures characterized by high expression of CXCL13, CCL19, and MS4A1, elevated CXCL13+ T cell signature scores in the LS-like spots compared to other spots of mLNs, and elevated SPP1+ IFI30+ macrophage signature scores in mLNs compared to nLNs.\nThe results have been added in Figures S12-S13, as below:\nFigure S12. Evaluation of residual LSs and CXCL13+ T cell signatures in mLNs of breast cancer using ST data.\n(A-F) H&E images (A), annotations of LS-like spots (B) and spatial distribution of CCL21 (C), CXCL13 (D) and MS4A1 (E) expressions as well as CXCL13+ T cell signature scores (F) in 4 mLNs of breast cancer and one normal lymph node (nLN). Spatial clustering was performed on the four mLN samples using BayesSpace, and clusters exhibiting elevated expression of LS-related markers (MS4A1, TRBC2, CXCL13, CCL19, CCL21) were annotated as LS-like spots. In the nLN sample, all spots were defined as LS spots, given its origin from histologically normal lymph node tissue. The LS or LS-like spots (LS/LS-like) are represented as red nodes and the others as blue (B).\n(G) Violin plot showing the difference of CXCL13+ T cell signature scores between LS-like spots and the other spots in the 4 mLNs of breast cancer. ****p < 0.0001 by Wilcox-test, unpaired, two-sided.\nFigure S13. Evaluation of SPP1+ IFI30+ MP cell signature in mLNs of breast cancer using ST data.\n(A-C) H&E images (A) and spatial distributions of IFI30 (B) and SPP1 (C) expressions in 4 mLNs of breast cancer and one nLN.\n(D-E) Violin plot showing the differences of IFI30 and SPP1 expressions (D) and SPP1+ IFI30+ MP signature scores (E) between mLNs of breast cancer and one nLN.\n****p < 0.0001 by Wilcox-test, unpaired, two-sided.\nThe results have been described in Page 31, Line 607-613, as below:\nWe further examined our findings in the mLNs of breast cancer based on ST data [32]. Similar to NSCLC, mLNs of breast cancer also exhibited residual LSs (LS-like spots defined based on high expressions of CXCL13, CCL19, CCL21 and MS4A1, Figures S12A-S12F), with significantly elevated CXCL13+ T cell signature scores (Figure S12G). In addition, the mLNs of breast cancer also displayed significant higher expressions of SPP1 and IFI30, along with elevated SPP1+IFI30+ MP signature scores compared to nLNs (Figures S13A-S13E). These conserved features suggest a potential universal role for SPP1+IFI30+ MP and CXCL13+ T cells in LN metastasis across multiple cancer types.\n-LS-like spots are defined by co-expression of MS4A1, TRBC2, CXCL13, CCL19, CCL21 without any orthogonal validation (imaging, IHC). The methodology for “two-step neighbors” is computationally defined and may overestimate aggregation. The authors should provide histologic confirmation (such as CD20, CD3, BCL6 staining) or cite prior validation of this marker-based LS/TLS call.\nRE: We thank the reviewer for raising these important points regarding the definition and validation of LS-like structures. Below we provide a detailed response addressing both the marker selection and histological correlation.\nBiological basis for LS/TLS marker selection:\nMS4A1 (CD20) is a cell surface protein primarily found on B cells and is a common marker used to identify B cell follicles within TLS [1]. TRBC2 encodes a component of the T cell receptor, and its expression is indicative of T cell presence within the TLS[2]. The chemokines CXCL13, CCL19, and CCL21 are vital for recruiting immune cells, particularly B cells and T cells, into the tumor microenvironment to form TLS. CXCL13 is produced by T cells and follicular dendritic cells within TLS, attracting CXCR5+ B and T cells [2]. CCL19 and CCL21, produced by fibroblasts and other cell types, also contribute to cell recruitment [3]. Consequently, this combination of the markers (MS4A1, TRBC2, CXCL13, CCL19, CCL21) are supported by the literature for identifying TLS and lymphoid aggregates. We have cited these references to validate this marker-based LS/TLS annotations.\nReferences:\n1. Fridman, W.H., et al., Tertiary lymphoid structures and B cells: An intratumoral immunity cycle. Immunity, 2023. 56(10): p. 2254-2269.\n2. Sautes-Fridman, C., et al., Tertiary lymphoid structures in the era of cancer immunotherapy. Nature Reviews Cancer, 2019. 19(6): p. 307-325.\n3. Zhang, Y., et al., CCL19-producing fibroblasts promote tertiary lymphoid structure formation enhancing anti-tumor IgG response in colorectal cancer liver metastasis. Cancer Cell, 2024. 42(8): p. 1370-1385.e9.\nMethodological clarification:\nThe LS-like spots identified in our analysis encompass both residual lymphoid structures and tertiary lymphoid structures (TLS). The LS-like spots were annotated according to the spatial clustering results, rather than “two-step neighbors” method, we have revised the method part to more explicitly state about spatial annotation method, in Page 51, Line 1043-1051, as below:\nSpatial annotation of ST data\nOver expressed marker genes of each spatial cluster were identified using FindAllMarkers function in Seurat, with the following thresholds: avg_log2FC> 0.1, P.adj<0.01, pct.1 > 0.3, and ranked top-8 based on log2FC. The spatial clusters were annotated based on whether the cluster marker genes included the classical marker genes of the major cell types including epithelial (SCGB3A1, KRT5, SFTPB, SFTPA1,S100A2, AKR1C1, CES1 ), fibroblast (COL1A1, COL1A2, ACTA2, FN1), myeloid (LYZ, C1QB, TYROBP), plasma (IGHG1, IGHG3, MZB1, JCHAIN) and LS-like structures (CXCL13, CCL19, CCL21, MS4A1 and TRBC2) [20-23]. One spatial cluster can be assigned with more than one cell type owing to the co-localization of different cell types spatially. Meanwhile, spatial clusters showing high expressions of mitochondria genes were annotated as MT-high.\nThe description about the LS-like spots were in Page 4, Line 113-118, as below:\nNotably, two of the spatial clusters (SC3 and SC4) showed simultaneous high expressions of MS4A1, TRBC2 (Figures 1D and S1D), indicating the aggregation of T and B cells [20]. The two clusters also showed high expressions of CXCL13, CCL19 and CCL21 (Figures 1D and S1D), all of which contribute to the formation of secondary and tertiary LS (TLS) [36-38]. These two clusters were annotated as LS-like spots, representing the residual LSs in the mLNs, and potential TLSs in the mLN+ PTs or mLN- PTs.\nHistopathological validation:\nWith guidance from clinical pathologists, we marked LS-like regions in the spatial transcriptomics data based on corresponding H&E-stained tissue sections.\nFurthermore, we have included zoomed-in views of two representative regions (new Figures 1D and 1E) to demonstrate strong concordance between the computationally annotated spatial clusters and histopathological assessment.\nFigure 1. Overview of the TME landscape of LN metastatic NSCLC\n(D) Hematoxylin and eosin (H&E) images and ST based spot-annotation results of the paired mLN+ PTs and mLNs from four patients. Each ST spot is colored according to the clustering and annotation results.\n(E) Enlarged regions of interest (ROIs) containing LS-like spots, including TLS in mLN+ PT and residual LS in mLN. Corresponding ROIs are outlined in red boxes in panel (D).\n-Batch correction is explicitly skipped for malignant cells “to retain patient specificity,” which risks inflating inter-patient differences in downstream clustering and DE. The authors should show that their key findings (Ig-high clusters, stemness gradients) persist under an alternative integration (eg Harmony, etc) using additional sensitivity analysis.\nRE: We thank the reviewer for raising this important point. In response, we have performed Harmony-based batch correction for the analysis of malignant cells to evaluate the robustness of our key findings. These results found that Ig-high clusters and CytoTRACE scores showed similar trends as identified by the previous analysis.\nWe have added these results in Figures S3F-S3I, as below:\nFigure S3. Elevated expression of immunoglobulin genes in epithelial cells of mLNs and mLN+ PTs\n(F) Umap plot showing the harmony-corrected clustering results of malignant cells. Harmony correction was applied on the PCA reductions to remove batch effects among different patients.\n(G) Dot plot showing the mean expression of classical marker genes of AT2/basal and immunoglobulin genes across the harmon-corrected malignant cell clusters. Clusters exhibiting significantly elevated expressions of Ig-genes (IGHG1, IGHG4, IGHA1 or IGHG2; examined by likelihood ratio test, Benjamini Hochberg adjusted P < 0.01 and mean_expression > 10) are highlighted by red fonts and termed as Ig-high clusters, and the others as Ig-low clusters.\n(H) Box plot of the CytoTRACE scores in the harmon-corrected Ig-high and Ig-low malignant cell clusters.\n(I) Violin plot comparing the normalized expression of CD274 in the harmon-corrected Ig-high and Ig-low malignant cell clusters. The point in the middle marks the median value.\nand we have added description in the results section, Page 12, Line 224-228, as below:\nMeanwhile, we conducted parallel clustering analyses of malignant cells - one preserving patient-specific batch effects and one with batch effect correction, to ensure robust results (Methods). Both results showed that some of the clusters presented elevated expressions of immunoglobulin genes than the other clusters (Figures 3D-3F and S3E-S3G). Notably, the malignant cells expressing high levels of immunoglobulin genes displayed significantly lower CytoTRACE scores than those with low immunoglobulin expressions (Figures 3G and S3H).\n-The methodology for spatial analysis need to be clarified. For example, the edge/core classification uses median distance thresholds, which could be sample-size dependent. The authors need to provide justification.\nRE: We thank the reviewer for this insightful comment regarding the spatial partitioning method. We have revised the Methods section to provide a more detailed description of the tumor core/periphery partition method, Page 52, Line 1068-1075 as follows:\nDetermination of tumor core and peripheral partitions for epithelial spots\nBased on the assumption that non-epithelial spots primarily represent non-tumor components of the microenvironment, for each epithelial spot on every ST slice, we calculated the minimum Euclidean distance between this spot and all the other non-epithelial cells using the spots’ coordinates. Here, a small distance indicates closer proximity to non-tumor regions and, by extension, far from the tumor core. Then, epithelial spots with minimum distances below the median value of all epithelial spots on the same ST slice were classified as belonging to the tumor periphery; the remaining epithelial spots were assigned to the tumor core region.\nHere, the median distance was computed independently for each ST slice, ensuring that the partition is based on a sample-specific threshold, not influenced by sample-size.\n-Also, LS like spot identification and two-step neighbor definition are described, but criteria for LS like designation are not detailed (which gene sets? morphology?)\nRE: We appreciate the reviewer’s request for further clarification regarding the criteria used to define LS-like spots. As we have explained in the previous question about LS-like spots, LS-like spots were annotated based on the high expressions of marker genes of LS/TLS including MS4A1, CXCL13, CCL19, CCL21 and TRBC2. This marker set is consistent with previously published definitions of tertiary lymphoid structures (TLS) and lymphoid aggregates [1-3].\nTo ensure biological relevance and structural correspondence, all computationally identified LS-like spots were cross-referenced with H&E morphology by an experienced pathologist. Representative regions are provided in Figure 1E.\nTwo-step neighbors of the LS-like spots were used to investigate the neighborhood surrounding LS-like regions, but not for definition of LS-like spots. we have revised the method part to more explicitly state about spatial annotation method, in Page 51, Line 1043-1051.\nReferences:\n1. Fridman, W.H., et al., Tertiary lymphoid structures and B cells: An intratumoral immunity cycle. Immunity, 2023. 56(10): p. 2254-2269.\n2. Sautes-Fridman, C., et al., Tertiary lymphoid structures in the era of cancer immunotherapy. Nature Reviews Cancer, 2019. 19(6): p. 307-325.\n3. Zhang, Y., et al., CCL19-producing fibroblasts promote tertiary lymphoid structure formation enhancing anti-tumor IgG response in colorectal cancer liver metastasis. Cancer Cell, 2024. 42(8): p. 1370-1385.e9.\n-Numerous Wilcoxon tests and GSEA are reported, FDR control is stated but not consistently shown for all comparisons. Clarify all thresholds used.\nRE: We thank the reviewer for raising this important point regarding multiple testing correction. We have carefully reviewed and clarified all statistical thresholds applied throughout our study, as detailed below:\nValidation-oriented single gene tests (e.g., marker gene validation or hypothesis-driven comparisons of pre-specified genes) were assessed using uncorrected Wilcox-test, as is conventional for focused, confirmatory analyses.\nExploratory multi-gene comparisons (e.g., differential expression across cell types or conditions involving the whole gene list), p-values were subjected to Benjamini–Hochberg correction.\nGene Set/pathway enrichment analysis was performed in an exploratory manner, and all reported pathway enrichment results include Benjamini–Hochberg adjusted P-values.\nAll figure legends have been updated to unambiguously indicate the specific statistical thresholds and correction methods used in each analysis.\nMinor points:\n-Some typos need to be corrected: “matirx” ; “contral” ; “princeple” ; “probalility” ; “up-regualted” . Please run a thorough copy edit.\nRE: We have revised them.\n-There are issues with acronyms, please define all at first use (e.g., “mCAF”, “tCAF”, “SC7”, “LS-like”).\nRE: We have checked acronyms and make sure we define all at first use.\n-Provide exact gene lists for CSC, SPP1⁺IFI30⁺ MP, and CXCL13⁺ T signatures used in validations for reproducibility.\nRE: The gene lists were provided in Table S2.\nHere, according to the Editor and Reviewers’ important suggestions and comments, we have revised carefully and have made corrections which we hope to meet with approval. We attached the response letter and revised manuscript and the corrections are colored in red to make the revisions notable. The “Calibri” fonts cited the Editor and Reviewers’ comments, and the “Times new roman” and “Blue” fonts are our replies. The “Red” fonts are cited in the manuscript for helping the editor and reviewers easily read. If any question arises, please do not hesitate to let us know.\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nMost concerns have been sufficiently addressed, and the revised manuscript has been substantially improved. However, there are several minor issues with the additional cell type co-localization analysis.\n1. In Figure S7E, the spatial distribution of CXCL13⁺T cells and iCAF appears largely non-overlapping, with minimal spatial co-localization. This observation seems inconsistent with the main text claiming their co-localization within LS-like spots. Although cell-cell interactions do not necessarily require stringent co-localization, as exemplified by the proposed long-distance interaction between SPP1⁺macrophages and malignant epithelial cells in Figure S5H, the discrepancy between visual spatial pattern and the claim of colocalization requires clarification or quantitative spatial proximity measurement to substantiate the interpretation.\nRE: We sincerely thank the reviewer for this comment and the opportunity to clarify. The reviewer is absolutely correct in their visual assessment of Figure S7E, the distribution of iCAFs is indeed broader and less perfectly overlapping with the dense cores of CXCL13+ T cells.\nIn our previous analysis, we calculated the co-localization scores across all spots, and compared the results in spot level, which did not account for sample-level variability. To address this, we have now performed a revised analysis at the sample level. Specifically, for each biological sample, we calculated the average co-localization scores separately for LS-like spots and all other spots within the sample.\nThe revised figure and figure legend are as below:\nFigure S7F. Violin plots showing the sample-averaged co-localization scores between CXCL13+ T cells and their interacting cells, calculated specifically within LS-like spots or all other spots across mLN and mLN+PT samples. The co-localization score between CXCL13+ T cell and a partner cell type for one spot was calculated as the sum of the CARD-derived abundance scores of CXCL13+ T cell and the partner cell type. T-test, un-paired, one-sided.\nThe corresponding results have been revised as below (Page 25, Line 481-485):\nCARD-based deconvolution analysis of the ST data also showed CXCL13+ T cells were significantly co-localized with key interacting cells, including LAMP3+ DCs, GC B cells, and activated B cells, within LS-like spots (Figures S7E and S7F). While spatial distributions of iCAF and mCAF cells were broader (Figure S7E), and did not show significant enrichment in LS-like spots (Figure S7F), suggesting that their interaction with CXCL13+ T cells might be less constrained to the core LS-like architecture.\n2. In Figure S7F, the y-axis appears to be mislabeled. It should be “Sum score with CXCL13⁺T” instead of “CSCL”.\nRe: We apologize for this error. We have revised this.\n3. The method for calculating the summed abundance score requires clarification. Is it defined as the sum of the abundance scores of CXCL13⁺ T cells and the interacting cell type, or based on another formulation?\nRe: We thank the reviewer for prompting us to clarify this. The method is defined as the sum of the abundance scores of CXCL13+ T cells and the interacting cell type. To make this clear, we have updated the figure legend for S7F, as below:\nThe co-localization score between CXCL13+ T cell and a partner cell type for one spot was calculated as the sum of the CARD-derived abundance scores of CXCL13+ T cell and the partner cell type.\n4. Since the hypothesis specifically predicts higher co-localization between these cell types and CXCL13⁺ T cells within LS-like regions, a one-sided Wilcoxon test would be more appropriate than the current two-sided test.\nRe: We thank the reviewer for this valuable suggestion. We fully agree that a one-sided test is more appropriate. Accordingly, we have re-analyzed the data using a one-sided, unpaired T-test. We selected the T-test because the new comparative variable (the sample-averaged co-localization score) is continuous and, despite the small sample size, the data met the assumptions for parametric testing in these specific comparisons. This choice offers greater statistical power than a non-parametric test. For completeness, we also performed a one-sided Wilcox test, which yielded consistent significance conclusions. The results from the T-test are now reported in Figure S7F.\nReviewer #2 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nReviewer #3 (Remarks to the Author):\nThe revised manuscript is significantly strengthened and addresses most prior concerns. The integration of new analyses and experiments alongside extensive methodological clarifications, has improved both the scientific rigor and clarity of the study. The manuscript now provides a more robust and reproducible multimodal dissection of LN metastasis in NSCLC with clear translational angles.\nRe: Thanks a lot for the affirmation of our revision and the valuable suggestions. We have revised the manuscript according to the comments and suggestions.\n- The claim that FTH1 upregulation promotes ferroptosis resistance and metastatic stemness is intriguing and largely supported by the data. However, it remains somewhat inferred that FTH1’s effect is indeed via ferroptosis resistance. It would strengthen the manuscript if the Discussion could explicitly note this remaining gap – i.e. that the link between FTH1 and LN metastasis, though plausible, is not yet causally shown in situ in LN tissues. If possible, add a short analysis of lipid peroxidation/ferroptosis gene programs in LN lesions to further buttress the mechanism. Nonetheless, given the practical constraints, the authors’ approach is reasonable.\nRe: We agreed that it remains somewhat inferred that FTH1’s effect is indeed via ferroptosis resistance. We have added discussion about this issue, in Page 34, Line 668-675, as below:\nWhile our data strongly associate FTH1 up-regulation with enhanced cancer stemness and lymph node metastasis in NSCLC, the causal role of FTH1 in conferring ferroptosis resistance within the lymph node microenvironment remains to be fully elucidated. Recent intravital imaging of melanoma LN metastasis has revealed that lymph-derived oleic acid protects disseminated tumour cells from ferroptosis [104] and fatty acid metabolite of 7-DHC promotes metastasis [105]，supporting the plausibility of our hypothesis in the LN metastasis. Future investigations using LN-specific ferroptosis models could help to better clarify whether FTH1 acts through ferroptosis evasion to facilitate LN metastasis.\nMeanwhile, we checked on the ferroptosis relevant gene expressions in the epithelial cells from metastasis lymph node (mLN), mLN+ primary tumor (mLN+ PT), and the paracancerous normal lung tissue samples (mLN+ N). Multiple ferroptosis-inhibition relevant genes (FTL, FTH1, PCBP1, PCBP2, SLC7A11, SLC3A2 and GPX4) were up-regulated in mLN and mLN+ PT samples compared to the mLN+ N samples.\nThe figure has been added as Figure S2D, as below:\nFigure S2D. Violin plot showing the up-regulations of some ferroptosis-inhibition relevant genes in the epithelial cells of mLN and mLN+ PT samples compared to normal lung samples (mLN+ N). P: Wilcox-test, unpaired, one-sided.\nWe also added text explanations in Page 9, Line 159-162, as below:\nMeanwhile, FTL and FTH1 along with multiple other ferroptosis relevant genes (PCBP1, PCBP2, SLC7A11, SLC3A2 and GPX4 [45]) were up-regulated in the epithelial cells of mLN and mLN+ PT samples compared to the mLN+ N samples (Figure S2D).\n- The discovery that a subset of malignant epithelial cells express immunoglobulin genes (IgG heavy chains in particular) is a novel and intriguing finding with implications for tumor biology and immune evasion. The authors have improved the evidence for this claim by performing IGHG4 knockdown experiments in NSCLC cell lines, which showed reduced expression of EMT/invasion markers and decreased invasiveness. The authors might consider briefly discussing potential mechanisms (for example, autocrine signaling etc) by which cancer-cell-produced immunoglobulins facilitate metastasis.\nRe: We have added some discussion, in Page 35, Line 685-688:\nExperimental analysis validated the regulatory role of IGHG4 in lung cancer metastasis, however, the exact mechanisms remain to be fully elucidated. Tumor-derived immunoglobins may be a target to prevent and treat tumor-associated thrombosis because it interacted directly with FcγRIIa and activated platelets [107]. Meanwhile, tumor-derived-immunoglobins can neutralize plasma-cell-derived immunoglobins to inhibit complement dependent cytotoxicity [95].\n-The spatial evidence for CAF-associated exclusion at tumor edges is convincing. If available, one orthogonal stain (αSMA/collagen or CD3) in a representative case would further anchor this conclusion; otherwise, acknowledge as a limitation.\nRe: We acknowledged this as a limitation in the discussion, Page 37, Line 727-729, :\nThe definitions of tumor core and peripheral using ST data, based on a limited cohort (n=4), should be considered preliminary and require further validation.\nReviewer #4 (Remarks to the Author):\nThe authors have sufficiently addressed my concerns and have made significant improvements to the manuscript. I do not have any further comments.\nRe: Thanks a lot for the affirmation of our revision\nREVIEWERS' COMMENTS\nReviewer #1 (Remarks to the Author):\nThe authors have fully addressed all the concerns we raised. We do not have any further comments.\nRe: We sincerely thank the reviewer for the positive feedback and for taking time to review our manuscript.\nReviewer #2 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nRe: We thank the reviewer for their collaborative effort in evaluating our manuscript.\nReviewer #3 (Remarks to the Author):\nThank you for the careful and substantive revision. The manuscript is notably stronger. I do not have any further comments.\nRe: We sincerely thank the reviewer for the positive feedback and for taking time to review our manuscript."}, {"round": 2, "reviewer_comments": {"reviewer_1": "Most concerns have been sufficiently addressed, and the revised manuscript has been substantially improved. However, there are several minor issues with the additional cell type co-localization analysis.\n1. In Figure S7E, the spatial distribution of CXCL13⁺ T cells and iCAF appears largely non-overlapping, with minimal spatial co-localization. This observation seems inconsistent with the main text claiming their co-localization within LS-like spots. Although cell-cell interactions do not necessarily require stringent co-localization, as exemplified by the proposed long-distance interaction between SPP1⁺ macrophages and malignant epithelial cells in Figure S5H, the discrepancy between visual spatial pattern and the claim of colocalization requires clarification or quantitative spatial proximity measurement to substantiate the interpretation.\n2. In Figure S7F, the y-axis appears to be mislabeled. It should be “Sum score with CXCL13⁺ T” instead of “CSCL”.\n3. The method for calculating the summed abundance score requires clarification. Is it defined as the sum of the abundance scores of CXCL13⁺ T cells and the interacting cell type, or based on another formulation?\n4. Since the hypothesis specifically predicts higher co-localization between these cell types and CXCL13⁺ T cells within LS-like regions, a one-sided Wilcoxon test would be more appropriate than the current two-sided test.", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_3": "The revised manuscript is significantly strengthened and addresses most prior concerns. The integration of new analyses and experiments alongside extensive methodological clarifications, has improved both the scientific rigor and clarity of the study. The manuscript now provides a more robust and reproducible multimodal dissection of LN metastasis in NSCLC with clear translational angles.\n- The claim that FTH1 upregulation promotes ferroptosis resistance and metastatic stemness is intriguing and largely supported by the data. However, it remains somewhat inferred that FTH1’s effect is indeed via ferroptosis resistance. It would strengthen the manuscript if the Discussion could explicitly note this remaining gap – i.e. that the link between FTH1 and LN metastasis, though plausible, is not yet causally shown in situ in LN tissues. If possible, add a short analysis of lipid peroxidation/ferroptosis gene programs in LN lesions to further buttress the mechanism. Nonetheless, given the practical constraints, the authors’ approach is reasonable.\n- The discovery that a subset of malignant epithelial cells express immunoglobulin genes (IgG heavy chains in particular) is a novel and intriguing finding with implications for tumor biology and immune evasion. The authors have improved the evidence for this claim by performing IGHG4 knockdown experiments in NSCLC cell lines, which showed reduced expression of EMT/invasion markers and decreased invasiveness. The authors might consider briefly discussing potential mechanisms (for example, autocrine signaling etc) by which cancer-cell-produced immunoglobulins facilitate metastasis.\n-The spatial evidence for CAF-associated exclusion at tumor edges is convincing. If available, one orthogonal stain (αSMA/collagen or CD3) in a representative case would further anchor this conclusion; otherwise, acknowledge as a limitation.", "reviewer_4": "The authors have sufficiently addressed my concerns and have made significant improvements to the manuscript. I do not have any further comments."}, "author_response": "Re: Thanks a lot for the affirmation of our revision and the valuable suggestions. We have revised the manuscript according to the comments and suggestions.\n- The claim that FTH1 upregulation promotes ferroptosis resistance and metastatic stemness is intriguing and largely supported by the data. However, it remains somewhat inferred that FTH1’s effect is indeed via ferroptosis resistance. It would strengthen the manuscript if the Discussion could explicitly note this remaining gap – i.e. that the link between FTH1 and LN metastasis, though plausible, is not yet causally shown in situ in LN tissues. If possible, add a short analysis of lipid peroxidation/ferroptosis gene programs in LN lesions to further buttress the mechanism. Nonetheless, given the practical constraints, the authors’ approach is reasonable.\nRe: We agreed that it remains somewhat inferred that FTH1’s effect is indeed via ferroptosis resistance. We have added discussion about this issue, in Page 34, Line 668-675, as below:\nWhile our data strongly associate FTH1 up-regulation with enhanced cancer stemness and lymph node metastasis in NSCLC, the causal role of FTH1 in conferring ferroptosis resistance within the lymph node microenvironment remains to be fully elucidated. Recent intravital imaging of melanoma LN metastasis has revealed that lymph-derived oleic acid protects disseminated tumour cells from ferroptosis [104] and fatty acid metabolite of 7-DHC promotes metastasis [105]，supporting the plausibility of our hypothesis in the LN metastasis. Future investigations using LN-specific ferroptosis models could help to better clarify whether FTH1 acts through ferroptosis evasion to facilitate LN metastasis.\nMeanwhile, we checked on the ferroptosis relevant gene expressions in the epithelial cells from metastasis lymph node (mLN), mLN+ primary tumor (mLN+ PT), and the paracancerous normal lung tissue samples (mLN+ N). Multiple ferroptosis-inhibition relevant genes (FTL, FTH1, PCBP1, PCBP2, SLC7A11, SLC3A2 and GPX4) were up-regulated in mLN and mLN+ PT samples compared to the mLN+ N samples.\nThe figure has been added as Figure S2D, as below:\nFigure S2D. Violin plot showing the up-regulations of some ferroptosis-inhibition relevant genes in the epithelial cells of mLN and mLN+ PT samples compared to normal lung samples (mLN+ N). P: Wilcox-test, unpaired, one-sided.\nWe also added text explanations in Page 9, Line 159-162, as below:\nMeanwhile, FTL and FTH1 along with multiple other ferroptosis relevant genes (PCBP1, PCBP2, SLC7A11, SLC3A2 and GPX4 [45]) were up-regulated in the epithelial cells of mLN and mLN+ PT samples compared to the mLN+ N samples (Figure S2D).\n- The discovery that a subset of malignant epithelial cells express immunoglobulin genes (IgG heavy chains in particular) is a novel and intriguing finding with implications for tumor biology and immune evasion. The authors have improved the evidence for this claim by performing IGHG4 knockdown experiments in NSCLC cell lines, which showed reduced expression of EMT/invasion markers and decreased invasiveness. The authors might consider briefly discussing potential mechanisms (for example, autocrine signaling etc) by which cancer-cell-produced immunoglobulins facilitate metastasis.\nRe: We have added some discussion, in Page 35, Line 685-688:\nExperimental analysis validated the regulatory role of IGHG4 in lung cancer metastasis, however, the exact mechanisms remain to be fully elucidated. Tumor-derived immunoglobins may be a target to prevent and treat tumor-associated thrombosis because it interacted directly with FcγRIIa and activated platelets [107]. Meanwhile, tumor-derived-immunoglobins can neutralize plasma-cell-derived immunoglobins to inhibit complement dependent cytotoxicity [95].\n-The spatial evidence for CAF-associated exclusion at tumor edges is convincing. If available, one orthogonal stain (αSMA/collagen or CD3) in a representative case would further anchor this conclusion; otherwise, acknowledge as a limitation.\nRe: We acknowledged this as a limitation in the discussion, Page 37, Line 727-729, :\nThe definitions of tumor core and peripheral using ST data, based on a limited cohort (n=4), should be considered preliminary and require further validation."}, {"round": 3, "reviewer_comments": {"reviewer_1": "The authors have fully addressed all the concerns we raised. We do not have any further comments.", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_3": "Thank you for the careful and substantive revision. The manuscript is notably stronger. I do not have any further comments."}, "author_response": "Re: We sincerely thank the reviewer for the positive feedback and for taking time to review our manuscript."}]

Decision letter