Single-cell and spatial transcriptome analyses reveal tumor immunometabolism in lymph node metastasis of lung cancer

Cell lines and cell culture

The human lung cancer cell lines A549 (RRID: CVCL_0023) and PC9 (RRID: CVCL_B260) were purchased from American Type Culture Collection (ATCC) (Shanghai, China). A549 and PC9 lung cancer cells were cultured in Roswell Park Memorial Institute (RPMI) 1640 medium supplemented with 5% fetal bovine serum (FBS, Gibco) and 1% penicillin–streptomycin antibiotic (Beyotime). All cell lines were tested and authenticated, and were cultured at 37°C in a humidified atmosphere containing 5% CO₂. The cell lines were Mycoplasma-free and authenticated through PCR analysis monthly.

Tumor sample handling and dissociation into a single-cell suspension for scRNA-seq

Tumor and LN samples were resected by an experienced physician, to ensure that the tumor samples contained neoplastic, stromal, and para-cancerous tissue, and the LN samples consisted of neoplastic and lymphoid tissue with complete cross-sections. Fresh samples were washed twice with precooled 0.04% bovine serum albumin (BSA, MACS) in RPMI 1640 medium and evenly divided into 2 parts for scRNA-seq and spatial transcriptomic (ST) analyses. Each sample for scRNA-seq was cut into approximately 1 mm³ pieces and enzymatically digested with 10 mL digestion medium containing 0.2% collagenase type IV (Gibco), 0.05% hyaluronidase type I-S (Sigma-Aldrich), and 0.002% DNase I (Applichem). After 45 min of digestion at 37°C in a shaking water bath, the enzymatic hydrolysate was filtered through a 40-μm pore size nylon mesh and centrifuged at 1,200 rpm for 5 min to obtain a single-cell suspension. Next, erythrocytes were lysed with red blood cell lysis solution (MACS) for 5 min. The cell suspension was centrifuged at 1,500 rpm for 5 min. After removal of the supernatant, the cell pellet was washed twice with RPMI 1640 with 0.04% BSA, then re-suspended in sorting buffer (PBS supplemented with 0.04% FBS). Cell concentration and viability were assessed with a Luna cell counter.

This study was approved by the Ethics Committee and Institutional Review Boards of Fudan University Shanghai Cancer Center (FUSCC) (IRB No. 2008223-9).

scRNA-seq and read processing

The concentrations of fresh cell suspensions for each sample were adjusted to 700–1,200 cells/μL. The cell suspensions were processed with a Chromium Next GEM Single Cell 3′ Reagent Kit v3.1 (10x Genomics, Pleasanton, CA) for library preparation according to standard protocols. The single cell libraries were sequenced on an Illumina NovaSeq 6000 system with paired-end sequencing (150 bp). The Cell Ranger software pipeline (version 3.1.0) from 10x Genomics was used to demultiplex cellular barcodes; map reads to the genome and transcriptome with the STAR aligner; and down-sample reads as required to generate normalized aggregate data across samples, thus producing a matrix of gene counts vs. cells.

Quality control and batch effect correction of scRNA-seq data

The Unique Molecular Identifier (UMI) count matrix was processed with the R package Seurat (version 4.4.1)¹⁶. To remove low-quality cells and likely multiplet captures, a major concern in microdroplet-based experiments, we applied the following criteria to filter cells with UMI/gene numbers out of the limit of mean ± 2 standard deviations, by assuming a Gaussian distribution of each cell’s UMI/gene number. After visual inspection of the distribution of cells according to the fraction of mitochondrial genes expressed, we discarded low-quality cells in which >10% of the counts belonged to mitochondrial genes. After application of these quality control criteria, 308,985 single cells were included in the downstream analyses. Library size normalization was performed with the NormalizeData function in Seurat to obtain the normalized count. To remove the batch effects in single cell RNA sequencing data, we used the RunHarmony function from the R package Harmony. The top 2,000 highly variable genes in the dataset were identified with the FindVariableFeatures function. Finally, the scaled and batch effect-corrected expression profiles of all samples were obtained for downstream analyses.

Unsupervised clustering and dimensionality reduction

The top principal components (PCs) were computed with the RunPCA function, according to the gene expression profiles of the top 2,000 highly variable genes. The ElbowPlot function in Seurat was used to select the optimal number of PCs for further analysis, as recommended by Seurat. The FindNeighbors and FindClusters functions in Seurat were both applied for cell clustering. The RunTSNE and RunUMAP functions were both performed for visualization, when appropriate. Clustree was applied to visualize changes in clustering patterns from a resolution of 0.1 to 1.5 (Figure S1E). The emergence of new clusters at lower resolutions and the formation of multiple parent clusters suggested instability in these regions. On the basis of prior literature¹⁷, we concluded that these new clusters were unlikely to represent genuine groups within the dataset and constitute over-clustering. Consequently, the resolution was ultimately set to 0.1, and 12 clusters were then identified from 308,985 single cells. To resolve the high heterogeneity in tumor tissue, we conducted CSV analysis to identify finer clusters within epithelial cells.

Identification of marker genes for cell clusters and cell type annotation

The marker genes in each cluster were identified with the FindAllMarkers function in Seurat. For a given cluster, FindAllMarkers identified markers positive with respect to all other cells. Subsequently, the markers were used to infer the cell of origin of each of the single cells independently and identify cell types. In the “low-resolution” clustering, epithelial cells (EPCAM, CDH1, KRT8, and KRT18), endothelial cells (VWF, PECAM1, CDH5, and CD34), plasma cells (MZB1, IGHG1, and SDC1), myeloid cells (CD14, CD68, AIF1, ENPP, and CSF1R), mast cells (MS4A2, KIT, and ENPP3), tubulin (PPBP, TUBB1, and PF4), B cells (CD19, CD79B, and MS4A1), T cells (CD3E, CD3D, CD8A, and CD4), and fibroblasts (DCN, COL1A1, COL1A2, and LUM) were identified.

Single-cell copy-number variation (CNV) evaluation

CNV evaluation for each chromosomal region in each cell was conducted with the inferCNV R package (version 1.20.0) (https://github.com/broadinstitute/inferCNV)¹⁸, according to gene expression in the single-cell transcriptome data. The CNV levels of epithelial cells were calculated with other cells serving as the reference. The genes were sorted by chromosomal location, and the average gene expression was calculated with 101 genes as the sliding window size. Normal cell expression was used as a control to generate the final CNV result file after denoising. The parameters of inferCNV analysis included “denoise”, default hidden Markov model settings, and a value of 0.1 for the “cutoff”. To decrease the false-positive CNV cells, we implemented the default Bayesian latent mixture model to identify the posterior probabilities of the CNV alterations in each cell, with the default value of 0.5 as the threshold.

Identification of signature genes among cell clusters

Differentially expressed genes (DEGs) were identified with the FindMarkers function of the Seurat package. The significance levels of these signature genes were determined with the Wilcoxon rank-sum test along with Bonferroni correction. P-value <0.05 and |log₂foldchange| >0.25 were set as thresholds for significant differential expression.

Inference of cell fate through trajectory analysis

Trajectory analysis was performed with the monocle R package (version 2.32.0)¹⁹ to reveal the cell-state transitions. The newCellDataSet function converted the Seurat object to a CellDataSet object. The estimateSizeFactors and the estimateDispersions functions were used to estimate size factors and dispersions. The ordering genes in the trajectory analysis were determined with the function FindVariableFeatures. The setOrderingFilter function was then applied to mark the genes to be used. The reduceDimension function in Monocle2 was applied to reduce the dimensionality. All functions used default settings. Because the normal ductal epithelial cells were separated and clearly defined in this study, those cells were set as the root state, and the orderCells function in Monocle2 was performed.

Tissue processing and Visium data generation

Lung cancer tissues were gently washed with cold 1× PBS (Gibco) and evenly divided to prepare a single-cell suspension and ST. Fresh frozen tumor tissues for ST were stored at −80°C and embedded in OCT (Sakura, Alphen aan den Rijn, Netherlands) before cryosectioning for RNA extraction. Each tumor cryosection was cut to a thickness of 10 μm with a cryostat (Leica CM1950, Germany) and mounted onto ST microarrays. For processing, the tissue was dehydrated with isopropanol for 1 min and stained with H&E. Bright-field images were obtained with a whole slide scanner (Panoramic MIDI FL, 3DHISTECH) at 20× resolution. H&E-stained sections of each sample were carefully reviewed by 2 experienced pathologists to confirm the pathology, and were then manually annotated by a trained pathologist to identify tumor, stromal, and normal ductal regions.

Responsiveness analysis

Augur (v1.0.3) is a method to prioritize cell types involved in the response to biological perturbations within scRNA-seq data. By training a machine-learning model for each cell type and evaluating their separability in the space of molecular measurements, Augur identified the responsiveness of each cell type.

ST barcoded microarray slide information

Library preparation slides were purchased from the ST team (https://www.spatialtranscriptomics.com). Each spot printed onto the array was 55 μm in diameter and 100 μm from center to center, covering an area of 6.5 × 6.5 mm². Each slide included 4 capture zones. Each capture zone had approximately 5,000 unique gene expression spots.

Tissue permeabilization, cDNA synthesis, library construction, and sequencing

Tissue sections were processed with Visium Spatial Gene Expression Reagent Kits (10x Genomics, Pleasanton, CA) for library preparation according to standard protocols. First, fixed and stained tissue sections were permeabilized with 70 μL permeabilization enzyme at 37°C for 18 min. Subsequently, each well was washed with 100 μL SSC, and 75 μL reverse transcription master mix containing reverse transcription reagents was added for cDNA synthesis. After incubation with reagents, full-length cDNAs with spatial barcodes generated from polyadenylated mRNAs were obtained. At the end of the first-strand synthesis, RT master mix was removed from the wells, and 75 μL 0.08 M KOH was added and incubated for 5 min at room temperature. The KOH was then removed from the wells, which were then washed with 100 μL EB buffer. A 75 μL volume of second strand mix was added to each well for second-strand synthesis. Finally, spatial barcodes were encoded, and full-length cDNA was amplified with PCR to obtain sufficient amounts for library construction. Libraries were prepared according to the Visium Spatial Gene Expression User Guide (CG000239, Rev B, 10x Genomics). The constructed library was sequenced on an Illumina NovaSeq 6000 system with a sequencing depth of at least 100,000 reads per spot with a paired-end 150 bp reading strategy.

SMI CosMx sequencing

SMI CosMx Sequencing technology is a powerful approach for in situ transcriptomics and proteomics analyses at subcellular resolution. This enzyme-free, amplification-free, hybridization-based single-molecule barcoding method enables direct interrogation of intact FFPE samples on standard pathology glass slides. The SMI approach uses a two-component probe design consisting of a target-binding domain and a readout domain. The target-binding domain is a 35- to 50-nucleotide sequence that hybridizes to the target RNA of interest, whereas the readout domain contains 4 consecutive 10- to 20-nucleotide reporter-landing sequences. Subsequent hybridization of these landing sequences with fluorescently labeled branched reporters enables a 64-bit barcoding scheme through combinatorial signaling across 4 fluorescence channels. The SMI CosMx sequencing workflow involves standard tissue preparation, which is followed by probe hybridization and cyclic reporter readout within a flow cell integrated with the SMI instrument. Notably, our use of a gene panel encompassing 6,000 targets enabled comprehensive transcriptomic profiling. After each round of reporter hybridization and imaging, a UV cleavage step efficiently quenched the fluorescence signals and allowed for subsequent incubation with the next set of reporters.

Spatial neighborhood analysis

A neighborhood graph was built with the generic spatial coordinates of the SMI CosMx data with the function gr.spatial_neighbors in squidpy (v1.6.1). Distances between 2 cells were compared with a given radius; if the distance was short, the 2 cells were considered neighbors. Through analysis of the graph, cells that surrounded specific cell clusters or types within a given radius were visualized. The results were normalized and are presented as proportions.

Animal experiments

An orthotopic lung cancer model was established through intrapulmonary injection in 4-week-old male nude mice. The mice were anesthetized with pentobarbital (50 mg/kg sodium) and positioned in right lateral decubitus. After an 8- to 10-mm incision was made in the left thoracic wall, the pleural cavity was exposed to visualize the left lung. A cell suspension comprising 1 × 10⁶ PC9 cells with >90% viability was prepared in 40 μL PBS supplemented with 40 mg growth factor-reduced Matrigel. The mice were randomly assigned to receive an injection of either control or shMOCCI preparations, according to a computer-generated random number sequence. This randomization ensured that any inherent variability among the mice was evenly distributed across all experimental groups, thereby minimizing its potential impact on the outcomes. The suspension was directly inoculated into the left lung parenchyma via the pleura. After injection, the surgical site was stapled, and the mice were positioned in left lateral decubitus for monitoring until they had fully recovered. At 4 weeks after orthotopic lung tumor implantation, the mice underwent bioluminescence imaging to evaluate primary tumor establishment and metastatic dissemination. The surgeons and analysts were blinded to all subsequent outcome assessments, including the measurement of primary tumors, and the dissection and weighing of LNs for metastasis, until all data were collected and analyzed.

Animal maintenance and all experiments were performed in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals and the Institutional Animal Care and Use Committee (IACUC), and were approved by the Shanghai Experimental Animal Center (Shanghai, China) (approval No. FUSCC-IACUC-S2025-1006).

Determination of cell type enrichment/depletion with multimodal intersection analysis (MIA)

MIA is an analytical method integrating single-cell expression profiles with spatial transcriptome data¹⁶. We first delineated sets of cell-type-specific and tissue-region-specific genes, then determined whether their overlap was higher (enrichment) or lower (depletion) than expected by chance. After obtaining 2 types of gene sets, we queried the significance of the overlap between ST genes and cell type marker genes with the hypergeometric cumulative distribution, with all genes as the background, to compute the P-value. In parallel, we conducted a test for cell type depletion by computing −log₁₀ (1–P).

Cell–cell communication analysis

To infer the potential intercellular communication network between clusters, we used the CellChat R package (version 2.1.2)²⁰ to quantitatively measure networks through the law of mass action, according to the average expression values of a ligand by one cell group and that of a receptor by another cell group, as well as their cofactors. The normalized expression matrix was imported, and a cellchat object was created with the createCellChat function. After data preprocessing with identifyOverExpressedGenes, identifyOverExpressedInteraction, and smoothData, we used the computeCommunProb, filterCommunication (min.cells = 10), and computeCommunProbPathway functions to calculate potential ligand–receptor interactions. Finally, with the aggregateNet function in CellChat, the aggregated cell–cell communication network was calculated.

KEGG enrichment analysis, hallmark pathway enrichment analysis, GSVA, and ssGSEA of DEGs

To determine the biological and molecular functional processes of the prioritized gene list as well as their significantly enriched pathways, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis and Hallmark pathway enrichment analysis with clusterProfiler (version 4.12.6), with P < 0.05. Single sample gene set enrichment analysis (ssGSEA) was used to define an enrichment score to represent the degree of absolute enrichment of a gene set in each sample, including KEGG, Gene Ontology Biological Process, and Reactome pathways. Normalized enrichment scores were calculated for the OXPHOS or glycolysis pathway. The ssGSEA analysis was performed in the R package GSVA (version 1.52.3)²¹. A P-value <0.05 was considered statistically significant. The gene sets used in ssGSEA were obtained with the function keggGet from the R package KEGGREST (v1.44.1). In detail, the OXPHOS and glycolysis gene sets were attached with hsa00190 and hsa00010, respectively. The database used to perform enrichment analyses was attached with the R package org.Hs.eg.db (v3.19.1). The gene signatures used to investigate the cellular states (activation, cytotoxicity, stress, and inhibition) for CD8+ T cells, CD4+ T cells, and NK cells were collected and curated from previous research²². Epithelial–mesenchymal transition (EMT) gene signatures for EMT score assessment were attached from previous research²³.

Plasmid/short hairpin RNA construction and viral infection

Human COX6C and DHRS2 were amplified with reverse-transcribed cDNA from the A549 and PC-9 cell lines and cloned into the lentiviral vector pSIN-Flag (puromycin-resistant; Addgene). Authenticity was verified by sequencing. The short hairpin RNA sequences of COX6C and DHRS2 were purchased from Sigma-Aldrich and cloned into the lentiviral vector pLKO.1 (Addgene). A highly efficient lentiviral system was used to generate the viruses. The cell lines were infected with the lentiviruses, and stable cell lines were established. The lentiviral transfection efficiency was >90% in all cell lines.

Western blot analysis

Total protein was extracted with RIPA buffer (Thermo Fisher Scientific, Waltham, MA, USA) and detected with BCA assays (Thermo Fisher Scientific, Waltham, MA, USA). Protein (30 μg per sample) was separated with SDS-PAGE and transferred onto a PVDF membrane (Gene Molecular Biotech, Inc., Shanghai, China). After being blocked with 5% milk for 2 h at room temperature, the membrane was incubated overnight at 4°C with primary antibodies to the following: C15orf48 (2 μg mL⁻¹, Abcam, ab110267), β-catenin (1:5,000, Abcam, ab32572), vimentin (1:1,000, Abcam, ab92547), E-cadherin (1:10,000, Abcam, ab40772), Snail (1:1,000, Abcam, ab216347), and β-actin (1:1,000, Abcam, ab179511). The membrane was subsequently incubated with HRP-conjugated goat anti-rabbit IgG secondary antibodies (1:10,000, Abcam, ab205718) or HRP-conjugated goat anti-mouse IgG secondary antibodies (1:10,000, Abcam, ab205719) for 1 h at room temperature, and expression levels were measured with an ECL kit (Roche Diagnostics, Basel, Switzerland) and Western blot imaging system.

Transwell assays

Transwell migration and invasion assays were conducted with Corning Transwell inserts (8.0 μm). For Transwell migration assays, 1.5 × 10⁴ transfected cells suspended in 50 μL serum-free medium were placed in the upper chamber, and 600 μL medium (10% FBS) was placed in the lower compartment. The cells were incubated at 37°C for 24 h. Successfully translocated cells were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet, then counted in 4 randomly chosen fields (200×) under a microscope. For Transwell invasion assays, 1.2 × 10⁵ cells were seeded on Transwell coated with 50 μL Matrigel (dilution of 1:4 with 0.2% BSA). Matrigel was used to coat the membranes for 12 h at 37°C before the invasion assays. The culture conditions were the same as those in the Transwell migration assays. The cells on the lower surface were fixed, stained, and photographed under a microscope after 48 h.

Immunohistochemical staining of formalin-fixed, paraffin-embedded tissue

Formalin-fixed, paraffin-embedded tumor tissues of 2 additional cases were separately cut into 5-μm sections and mounted on glass slides. The slides were baked at 65°C overnight. After deparaffinization and hydration, the slides were boiled in citrate buffer at 100°C for 15 min. Subsequently, 3% H₂O₂ solution was used to block endogenous peroxidase activity for 20 min. To prevent nonspecific antibody binding, we incubated the slides with 5% normal goat serum for 1 h at room temperature. The slides were subsequently incubated at 4°C overnight with anti-C15orf48 primary antibody (Abcam, ab110267) at 1:100, which were validated for immunohistochemistry (IHC) by the manufacturer in HeLa cells and PC9 cells, respectively. After being washed with TBST 3 times, the slides were incubated with HRP-conjugated goat anti-rabbit/mouse secondary antibodies (GeneTech, GK500705) for 1 h at room temperature. Sections were stained with DAB, then counterstained with hematoxylin according to the manufacturer’s instructions.

Quantification and statistical analysis

All statistical analyses were performed in R (version 4.4.1) and GraphPad Prism 8.0. Assays were repeated at least 3 times, and the data are presented as mean ± standard deviation. Statistical tests were dependent on data distribution: normality was assessed with the Shapiro–Wilk or Kolmogorov–Smirnov test, depending on sample size (< 50, Shapiro–Wilk test or otherwise Kolmogorov–Smirnov test). Homogeneity of variance was validated with analysis of variance. The t-test and Wilcoxon rank-sum test were applied accordingly, as indicated in each figure legend. Spearman rank correlation coefficients, which are less sensitive to outliers and more robust, were calculated with the cor.test function in the base R package. This non-parametric measure of rank correlation assesses the strength and direction of the association between 2 variables. Differentially expressed genes were identified with two-sided Wilcoxon rank-sum test with Bonferroni correction in Seurat. In pathway enrichment analyses, the Benjamini–Hochberg method was used for multiple comparison correction. For CellChat cell-cell interaction analyses, receptor–ligand interactions between cell types were identified according to the specific expression of a receptor by one cell type and a ligand by another cell type. The significance of specific ligand–receptor interactions between 2 cell subsets was calculated with permutation tests. Wilcoxon rank-sum tests were used to identify genes significantly overexpressed (P < 0.05) in specific cell clusters compared with all other cells.

scRNA-seq identifies a primary metastatic cancer cell cluster in paired primary lung tumors and metastatic lymph nodes

To comprehensively identify the cellular composition and architecture of lung cancer primary tumors and metastatic LNs, we performed single-cell RNA sequencing (scRNA-seq, 10x Genomics) on a cohort comprising 7 matched primary tumor, MetLN, uiLN, and PBMC samples, as well as 3 pairs of primary tumor, uiLN, and PBMC samples collected from 10 patients in the FUSCC cohort during surgical resection (Figure 1A, B). Among these patients, 7 had lung squamous cell carcinoma, and 3 had lung adenocarcinoma (LUAD). Two pathologists confirmed all LNs, to ensure metastasis. After filtering the scRNA-seq data to exclude damaged or dead cells and putative cell doublets, we constructed an atlas consisting of 308,985 single cells, with a median of 1,732 expressed genes that passed the stringent quality filtering.

• Download figure
• Open in new tab
• Download powerpoint

• Download figure
• Open in new tab
• Download powerpoint

Figure 1

Expression profiling of 308,985 single cells in metastatic LNs, uninvolved LNs, primary lung tumors, and peripheral blood. (A) Graphical overview of the workflow. (B) Heatmap showing the clinical features of scRNA-seq data at the patient level. (C) UMAP plot of 308,985 cells, showing 9 major cell types (left), LN+ patients and LN− patients (middle), and 2 tissue types (right). Each cell type is shown in a unique color. (D) Bubble plots of marker gene expression and expressed percentages in major cell types. Dot size reflects the percent of a cell type expressing a specific gene, whereas dot color represents expression level. (E) Bar plots demonstrating the relative proportions of cell types across 10 patients in 4 tissue types. Each cell type is colored differently. LN, lymph node; LN+, lymph-node–positive; LN−, lymph-node–negative; UMAP, Uniform Manifold Approximation and Projection.

After integrating the transcriptional profiles from all captured cells, we first applied low-resolution Uniform Manifold Approximation and Projection (UMAP) clustering and generated 2D graphs with 12 clusters across 37 samples from 10 patients (Figure S1A). To delineate the principal cellular constituents within this atlas, we annotated each cluster according to the expression of canonical marker genes (Figures 1C and S1B). All cells were divided into 9 cell types: epithelial cells (EPCAM, KRT19, and PROM1), plasma cells (IGHG1, MZB1, and SDC1), myeloid cells (CD68, CD14, AIF1, and CSF1R), mast cells (MS4A2, KIT, and ENPP3), tubulin (PPBP, TUBB1, and PF4), B cells (CD19, CD79B, and MS4A1), T cells (CD3D, CD3E, CD8A, and CD4), fibroblasts (COL1A1, COL1A2, DCN, and LUM), and endothelial cells (PECAM1 and VWF) (Figures 1D and S1C). Clustree¹⁷ was applied to avoid over-clustering (see Materials and Methods) (Figure S1D). In MetLNs, the proportions of B and T cells were higher than observed in primary tumors but remained lower than observed in uiLNs. Myeloid and mast cells were enriched in PBMCs and primary tumors, but were comparatively low in MetLNs and uiLNs (Figures 1E and S1E). In summary, we identified multiple cell populations with distinct distribution patterns in different lung cancer tissues.

Given that lung cancer cells originate from epithelial cells, we applied CNV analysis to delineate malignant epithelial cells from their non-malignant counterparts. CNV analysis is a widely used approach in scRNA-seq to explore tumor evolution and disease progression. We evaluated the CNV level of each epithelial cell cluster (Figures 2A and S2A). Among all epithelial cell clusters, C2 exhibited markedly lower CNV levels than the others (Figure 2A) and therefore appeared to correspond to a non-malignant epithelial population, whereas the remaining clusters probably represented malignant phenotypes. To characterize the trajectory of epithelial transformation, we performed pseudo-time analysis across the 5 epithelial clusters, which revealed a bifurcated trajectory originating from C2 and progressing toward 2 distinct malignant endpoints (Figures 2B and S2B). In agreement with its low CNV level, C2 was located at the bottom left corner of the trajectory curve, representing the inferred origin of epithelial tumor evolution. From this origin, the trajectory diverged into 2 branches, referred to as Cell Fate 1 and Cell Fate 2. C4, which constituted the highest proportion of LN metastatic cells and a smaller fraction of primary tumor cells, was therefore classified as PLMCs (Figure S2C). Pseudo-time analysis showed that tumor cells from MetLNs tended to be present in the middle of the trajectory and also appeared at the ends of both branches (Figure S2D, E). Although clusters C1, C3, and C4 co-occupied the terminal trajectory space, the limited abundance of C1 and C3 in LNs suggested that these clusters did not represent the phenotype of LN metastatic tumors.

• Download figure
• Open in new tab
• Download powerpoint

• Download figure
• Open in new tab
• Download powerpoint

Figure 2

Characteristics of metastatic epithelial cells and MOCCI. (A) Violin plot of the CNV levels of all epithelial cells in 5 clusters. (B) Pseudo-time trajectory plot of all epithelial cells. Two distinct cell fates are identified by trajectory. Arrows show the potential evolutionary orientation of the trajectory. (C) Bubble plot of enriched KEGG pathways between PLMC and NPLMC clusters. Dot color shows the adjusted P-value of each pathway, and dot size represents the gene counts of each pathway. The difference was assessed with Wilcoxon signed-rank test. (D) Heatmap of the selected genes between PLMC and NPLMC clusters along 2 cell fates, which clustered into 3 patterns. The color bar, from blue to red, indicates expression levels from low to high. (E) Dot plot of the dynamic expression of marker genes in OXPHOS and glycolysis along the pseudotime trajectory. (F, G) Kaplan–Meier curves demonstrating overall survival of patients with LUAD from FUSCC cohorts (F) and TCGA datasets (G). Patients were divided into 2 groups according to MOCCI expression. CNV, copy number variation; DEG, differentially expressed genes; FUSCC, Fudan University Shanghai Cancer Center; KEGG, Kyoto Encyclopedia of Genes and Genomes; LUAD, lung adenocarcinoma; MOCCI, Modulator of Cytochrome C Oxidase during Inflammation; NPLMC, non-primary lymph node metastatic cell; OXPHOS, oxidative phosphorylation; PLMC, primary lymph node metastatic cell.

A metabolic switch between OXPHOS and glycolysis is a key regulator in PLMC clusters

To delineate the transcriptional alterations in PLMCs, we first performed DEG analysis between PLMCs and the remaining tumor cell populations (NPLMCs). Significantly modulated genes were subjected to KEGG pathway enrichment analysis, which revealed key pathways altered during LN metastasis (Figure 2C). The top enriched pathways, including OXPHOS, ATP biosynthetic process, and ATP synthesis coupled electron transport, were predominantly associated with energy metabolism, thus collectively suggesting the potential metabolic reprogramming process accompanying the transition from primary to metastatic tumor states. Notably, the significant OXPHOS pathway enrichment observed in PLMCs suggested novel metabolic traits of LN metastasis in lung cancer.

To further explore this transition, we analyzed the expression dynamics of pathway-associated genes along the inferred pseudo-time trajectory. A heatmap depicted changes in gene expression from the initial state to Cell Fate 1 or 2, and identified 2 different transformation patterns (Figure 2D). Although few lung cancer studies have focused on OXPHOS pathway switching, similar transitions have been observed in breast cancer. Liu et al.²⁴ have recently described a switch toward OXPHOS during early metastatic seeding, in agreement with our findings. Representative OXPHOS markers [ATP synthase membrane subunit e (ATP5ME), dehydrogenase/reductase 2 (DHRS2), and NADH ubiquinone oxidoreductase subunit A3 (NDUFA3), NDUFA4] exhibited an up-then-down expression pattern, whereas marker genes of glycolysis pathways (LDHA) followed an inverse, down-then-up pattern (Figures 2E and S2F, G). In parallel, the expression of genes associated with cell adhesion and migration (CLDN3 and F11R) also exhibited an up-then-down pattern indicating dynamic changes in cellular behavior during dissemination (Figure S2H). These findings suggested that this metabolic switch in lung cancer cells might trigger changes in EMT during the metastatic process. Furthermore, the reversal of metabolic preference might help meet the energy demands of cell proliferation after metastatic seeding. Of note, among the leading pathway-associated genes identified along the pseudo-time trajectory (Figure 2D), MOCCI (also known as C15orf48) was the sole gene that both correlated with OXPHOS and co-varied precisely with the OXPHOS metabolic shift (Figures 2E and S2I). Recent work has shown that MOCCI functions as a component of mitochondrial complex IV, in which it displaces the NDUFA4 protein under inflammatory conditions, and its expression is positively associated with the severity of various inflammatory diseases^25,26. We validated the associations among MOCCI, the OXPHOS pathway, and mediastinal LN metastasis (known as the N descriptor in the TNM staging system of lung cancer) through ssGSEA in our FUSCC cohort, comprising bulk RNA-sequencing and whole-genome sequencing data from 875 lung cancer cases²⁷. Compared with patients without LN metastasis, those with LN metastasis demonstrated significantly higher expression of MOCCI, higher OXPHOS pathway scores, and lower glycolysis pathway scores (Figure S2J–L). Notably, A significant correlation was observed between MOCCI expression and higher nodal (N) status (Figure S2J). Furthermore, analysis of gene expression and clinical outcome data from The Cancer Genome Atlas (TCGA) revealed that MOCCI upregulation was significantly associated with poorer survival in LUAD (Figure 2F, G). These results suggested that MOCCI might play a critical role in regulating OXPHOS reprogramming. In summary, our findings uncovered a potential metabolic transition between OXPHOS and glycolysis in the evolution of lung cancer metastases, thereby implicating OXPHOS as a critical metabolic feature in this process.

T-cell landscape and PLMC signaling reveal T-cell exhaustion during metastasis and distinct intercellular communication profiles

To investigate the potential interactions between PLMCs and immune cells within the TME, we first extracted all T cells, NK/NKT cells, and B cells for re-clustering, aiming to characterize immune composition and functional transitions within MetLNs, in comparison with both primary tumors and uiLNs. T cells from all samples were integrated and re-clustered into 13 subsets according to marker genes (Figure S3A, B). T cells were initially categorized according to the conventional CD4⁺/CD8⁺ paradigm. Although no statistical significance was found, uiLNs appeared to have the highest proportion of CD4⁺ T cells and the lowest proportion of CD8⁺ T cells (P = 0.078), whereas MetLNs showed a similar trend to primary tumors (P = 0.11) (Figure 3A). To further characterize functional states, we applied 4 commonly used gene signatures (activation, cytotoxicity, stress, and inhibition) to CD8⁺ T cells (Table S1), CD4⁺ T cells (Table S2), and NK/NKT cells (Table S3). Intriguingly, CD8⁺ T cells from primary tumors exhibited the highest scores for activation, cytotoxicity, and inhibitory markers among the 4 tissues, but had stress scores comparable to those in MetLNs (Figure 3B). A similar trend was observed in CD4⁺ T cells (Figure S3C): primary tumors showed the highest scores across all 4 dimensions, and were followed by MetLNs, uiLNs, and PBMCs. Interestingly, the cytotoxicity scores for CD4⁺ T cells were comparable between uiLNs and MetLNs. In NK/NKT cells, primary tumors also exhibited lower inhibitory scores and higher cytotoxicity scores than MetLNs and uiLNs (Figure S3D). These findings provided a rationale for primary tumors having slightly more stimulated and functional T cells than MetLNs, and also confirmed the immune reserve ability from uiLNs.

• Download figure
• Open in new tab
• Download powerpoint

• Download figure
• Open in new tab
• Download powerpoint

Figure 3

Evolution of the tumor microenvironment in the transformation process and its interaction with malignant T cells. (A) Bar plots (left) showing the proportions of each tissue type in CD4/CD8 T cells. Line plots (right) representing the difference in MetLN-uiLN (top) and MetLN-tumor (bottom) proportions in CD4+/CD8+ T cells. (B) Box plots comparing the activation (top left), cytotoxicity (top right), inhibitory (bottom left), and stress (bottom right) scores of CD8⁺ T cells among 4 tissue types. The boxes represent medians (indicated by the horizontal line with a numerical value), and encompass the second to third quartiles (forming the body of the box); Tukey-style whiskers extend beyond the box to depict the data range. The Wilcoxon rank-sum test was applied to assess the difference. (C) Bar plots (left) showing the proportion of each tissue type assigned to a T cell type. Line plots (right) representing the differences in MetLN-uiLN (top) and MetLN-tumor (bottom) proportions assigned to a T cell type. (D) Heatmap showing the mean expression of 7 gene sets (co-stimulatory, cytotoxicity, inhibitory, MHC-II, naïve, pre-dysfunction, and resident) from each T cell type among 3 tissue types. Color bar coding from blue to red indicates expression levels from low to high. (E) Bar plots of the number and strength of inferred interactions of T cells, PLMCs, and NPLMCs in both MetLNs and tumors. (F, G) Heatmaps of differential interaction strength (F) and differential number of interactions (G) of T cells, PLMCs, and NPLMCs between MetLNs and tumors. (H) Circle plots showing cell-cell interactions among PLMC clusters, NPLMC clusters, and T cell types in MetLNs. Arrows and colors indicate the direction of interaction. Dot size indicates the proportion of cells in each cell group. Edge thickness indicates the number (left) and strength (right) of interactions. (I) Bubble plots showing the significantly differentially expressed ligand-receptor pairs of T cells, PLMCs, and NPLMCs in MetLNs. Dot size reflects P-values, and dot color represents communication probabilities. MetLN, metastatic lymph node; NPLMC, non-primary lymph node metastatic cell; PLMC, primary lymph node metastatic cell; uiLN, uninvolved lymph node.

On the basis of the expression of selected marker genes, we annotated the 15 T cell clusters into distinct functional categories: naïve/memory-like (clusters 2, 3, 7, 8, 10, and 11, including 1 with a more naïve phenotype, 3 central memory, and 2 expressing IFN genes; T-CM), exhausted/regulatory-like (clusters 5, 6, and 14, 2 comprising CD4⁺ T cells and 1 comprising CD8⁺ T cells; T-EXH), effector (clusters 1, 9, 13, and 15; T-EFF), and follicular helper T cells (cluster 12; Tfh) (Figure S3E). At single-cell resolution, the evolution of T cell composition was delineated across 4 distinct paired tissues (Figure 3C). As expected, we found a consistent increase in T-EFF and T-EXH and a decrease in T-CM and Tfh cells from MetLNs to tumors. We also observed a significant increase in T-EFF cell clusters during the transformation from uiLNs to MetLNs, whereas other clusters, including T-CM, T-EXH, and Tfh cells, appeared to be invariant, thus suggesting an immune reserve ability from uiLNs.

We next investigated the transcriptional signatures of T cell subsets across paired tissue types (Figure 3D). T-CM cells in tumors and uiLNs showed a more naïve signature, characterized by elevated expression of CCR7, TCF7, LEF1, SELL, and IL7R. In contrast, MetLN-derived T-CM cells appeared to exhibit diminished reserve ability. Tumor-derived T-EFF cells had the highest levels of cytotoxicity markers (GZMB, PRF1, IFNG, and GNLY) and tissue-residency markers (CD69 and ITGAE), and were followed by MetLNs, thus implying their robust cytotoxic potential and high propensity to differentiate into tissue-resident memory T (Trm) cells, probably driven by the microenvironment. Although uiLNs demonstrated a relatively quiescent state, T-EFF cells from uiLNs retained considerable expression of GZMA and CST7, thereby suggesting anti-tumor capability. Both uiLN- and MetLN-derived T-EFF cells exhibited low inhibitory markers, whereas tumor-derived T-EFF cells displayed markedly high expression of LAG3, PDCD1, and HAVCR2. The T-EXH cells from tumors also showed elevated inhibitory and co-stimulatory markers, such as CTLA4, TIGIT, TNFRSF9, and CD28, thus further corroborating their functionally exhausted phenotype. Interestingly, MetLN-derived Tfh cells displayed broadly lower expression of MHC class II (MHC-II) genes than observed in primary tumors and uiLNs; these findings indicated relative suppression of antigen presentation in MetLNs. Differential gene analysis of immune cells between MetLNs and uiLNs also demonstrated that T-EFF cells in MetLNs, compared with uiLNs, significantly overexpressed the cytotoxic gene GZMB, thus indicating potent cytotoxic potential (Figure S3F).

To explore potential intercellular communication between PLMCs and immune cells within the TME, we performed CellChat analysis to identify cell-cell signaling links. Although the strength of interaction among T cell subtypes was higher in primary tumors, the number of interactions between the 2 groups was very close (Figure 3E). PLMC clusters showed fewer interactions with T cells in MetLNs, except for Tfh cells, whereas the strength of interaction was similar to that in primary tumors (Figure 3F–H). An investigation of the specific ligand–receptor interactions between tumor subtypes (PLMCs and NPLMCs) identified a markedly different pattern between primary tumors and MetLNs (Figure 3I). In MetLNs, the macrophage migration inhibitory factor pathway was highly active between T cells and either PLMCs or NPLMCs, thus indicating its potential role in overall tumor progression. However, PLMCs showed more active ligand-receptor interactions than NPLMCs in MetLNs, such as AREG-(EGFR+ERBB2), AREG-EGFR, and BMP4-(BMPR1B+BMPR2); PLMCs therefore might more strongly induce EMT^28,29, angiogenesis³⁰, and metastases³¹. Similarly, PLMCs demonstrated more ligand receptors that tended to induce metastases in primary tumors, such as ANGPTL4-(ITGA5+ITGB1)³² and CCL20-CCR6³³. CellChat analysis of potential intercellular communication from immune cells to PLMCs within the TME also revealed that PLMCs exhibited more active BMP4-(BMPR1B+BMPR2) ligand-receptor pair interactions compared to NPLMCs with T cells in MetLNs, thus suggesting enhanced EMT-inducing ability (Figure S3G)²⁹. Additionally, in tumors, T-EFF cells exhibited active ligand-receptor pair interactions with PLMCs via CCL4-CCR5, thereby suggesting PLMC participation in tumorigenesis and promotion of tumor progression (Figure S3H)³⁴. Remarkably, unlike other ligand-receptor interactions, GDF15-TGFBR2 demonstrated a strong interaction in MetLNs specifically not in primary tumors, thus suggesting its potential role in LN-specific tumor-T cell interactions. Overall, the CellChat analysis demonstrated that PLMCs had an intercellular communication profile distinct from that of NPLMCs and characterized by enhanced metastatic signaling activity.

Integrative analysis of the B-cell landscape and PLMC signaling reveals enhanced metastatic potential in lymph node metastasis

We extracted and re-clustered all B cells into 12 subsets according to canonical marker genes. Further annotation based on the expression of selected marker genes allowed us to categorize the 12 B cell clusters as memory-like (clusters 1, 3, 7, and 9; memory-B), naïve-like (clusters 2, 6, 10, and 11; naïve-B), germinal-center-like (cluster 8; GC-B), and plasma-like (cluster 4, 12; plasma cell) groups (Figure S4A–C). Compared with tumors, MetLNs exhibited a higher proportion of memory-B, GC-B, and naïve-B cells, but a lower proportion of plasma cells (P = 0.16) (Figure 4A). However, compared with uiLNs, MetLNs contained slightly fewer memory-B (P = 0.3) and naïve-B cells (P = 0.078), and a greater proportion of plasma cells, thus suggesting a more mature B cell state in MetLNs than uiLNs. A stress signature score (Table S4) was used to evaluate the microenvironment across tissues, and little difference was found (Figure S4D).

• Download figure
• Open in new tab
• Download powerpoint

• Download figure
• Open in new tab
• Download powerpoint

Figure 4

Characteristics of malignant B cells and their interactions with other cells. (A) Bar plots (left) showing the proportion of each tissue type assigned to a B cell type. Line plots (right) representing the differences in MetLN-uiLN (top) and MetLN-tumor (bottom) proportions assigned to a B cell type. (B) Bar plots of the number and strength of inferred interactions of B cells, PLMCs, and NPLMCs in both MetLNs and tumors. (C) Heatmaps of differential interaction strength (left) and differential number of interactions (right) of B cells, PLMCs, and NPLMCs between MetLNs and tumors. (D) Circle plots showing cell-cell interactions among PLMCs, NPLMCs, and B cells in MetLNs. Arrow and color indicate the direction of interaction. Dot size illustrates the proportion of cells in each cell group. Edge thickness indicates the number (left) and strength (right) of interactions. (E) Bubble plots showing the significantly differentially expressed ligand-receptor pairs of B cells, PLMCs, and NPLMCs in the MetLNs. Dot size reflects P-values, and dot color represents communication probabilities. (F) Lollipop chart (left) illustrating the Augur scores of cell types across all cells. The actual value of the Augur score for each cell type is indicated. Heatmap (right) demonstrating correlations between 2 different cell types. Color from blue to red indicates the intensity of correlations from low to high. (G) Correlations between the proportions of selected cell clusters. Redder hues indicate stronger correlations between the 2 cell subpopulations. Spearman correlation coefficient was used. MetLN, metastatic lymph node; NPLMC, non-primary lymph node metastatic cell; PLMC, primary lymph node metastatic cell; uiLN, uninvolved lymph node.

To comprehensively characterize B cell receptor (BCR) repertoire diversity, we analyzed paired scBCR-seq data from 7 metastatic patients, with 28 matched samples (from 4 tissue types) that were also included in the scRNA-seq cohort, and identified 68,793 cells showing productive VDJ rearrangements. Compared with primary tumors, the diversity in BCR repertoires was markedly greater in MetLNs and uiLNs (Figure S4E), including a higher proportion of expanded clonotypes present at elevated clonal frequencies (Figure S4F). In both MetLNs and uiLNs, memory-B cells contained the most clone types, whereas plasma cells were predominantly clonally expanded in tumors (Figure S4G, H). We also measured how the BCR repertoire in a sample from a given geospatial location overlapped with that from another location, according to the Jaccard index (Figure S4I). The Jaccard index of shared clones (0.014) between primary tumors and MetLNs was similar to that between primary tumors and uiLNs (0.012).

To investigate the potential interactions between PLMCs and B cells, we performed CellChat analysis. Whereas the strength of interactions among B cell subtypes was higher in primary tumors, the number of interactions between the 2 groups was very close, similarly to observations in T cells (Figure 4B). In MetLNs, PLMC clusters showed significantly more interactions with GC-B, naïve-B, and memory-B cells compared to those in primary tumors, but similar strength in interactions with B cells (Figure 4C, D). We then examined ligand–receptor interactions between PLMC and NPLMC tumor subtypes (Figure 4E). In MetLNs, NAMPT-(ITGA5+ITGB1), which is known to activate the PI3K-AKT pathway³⁵, was highly active between B cells and either PLMCs or NPLMCs, thus indicating its potential role in overall tumor progression. However, PLMCs showed more active ligand-receptor interactions, such as CXCL8-CXCR2 and FASL-FAS, than NPLMCs in MetLNs. Therefore, as previously suggested, PLMCs appeared to more strongly induce B cell transition³⁶ and Breg-regulated immune tolerance³⁷. Similarly, PLMCs demonstrated more ligand receptors that tended to induce metastases, such as ANGPTL 4-(ITGA5+ITGB1) and CCL20-CCR6. Interestingly, strong activation of GDF15-TGFBR2 was observed exclusively in MetLNs, but not in primary tumors, similarly to the findings in T cells. Overall, CellChat analysis demonstrated that PLMCs exhibited a distinct B cell interaction profile suggesting greater metastatic potential than that of NPLMCs.

To investigate the underlying causes of differences in the TME between MetLNs and primary tumors, we applied the Augur algorithm²² to compare the transcriptional responses of various cell subtypes in these regions (Figure 4F). In MetLNs, B cells and T cells displayed higher response scores than myeloid cells or fibroblasts; therefore, B cells and T cells might play a more important role in LN metastasis. Although mast cells exhibited an exceptionally high response score, their very low abundance in LNs might have inflated their rank. Among all subsets of immune cells, because naïve-B cells scored highest (and were followed by T-CM and Tfh cells), their activation and differentiation appeared to be key features of the MetLN microenvironment. Furthermore, we conducted Spearman correlation analysis to explore the relationships among various cell clusters. Although no significant positive correlation was observed between immune cells and NPLMCs, various subtypes of B cells and T cells, including T-EFF and plasma cells, showed significant positive correlations with PLMCs (Figure 4G).

Combined spatial transcriptome and scRNA-seq analyses reveal spatial features of PLMC clusters

To further assess the spatial organization of epithelial cells, we analyzed 27 paired cryosections from 10 patients, comprising 10 tumor, 7 MetLN, and 10 uiLN samples, all of which were also included in the scRNA-seq cohort. CosMx SMI³⁸ analysis was performed to characterize spatial interactions between tumor cells and the surrounding TME at single-cell resolution (see Materials and methods). A total of 362,428 cells were obtained at a median depth of 645 genes/cell. We identified 17 distinct cell types: epithelial cells, endocrine cells, B cells, plasma cells, pDC cells, macrophages, mast cells, monocytes, cDC cells, endothelial cells, fibroblasts, lymphatic vessels endothelial cells, pericytes, CD4⁺ T cells, CD8⁺ T cells, NK cells, and Treg cells (Figures 5A and S5A, B). The cellular composition notably varied across tissue compartments (Figure S6A, B).

• Download figure
• Open in new tab
• Download powerpoint

• Download figure
• Open in new tab
• Download powerpoint

Figure 5

Spatial features of the PLMC cluster. (A) Unsupervised clustering of ST cells (left), H&E staining of tissue sections (middle left), PLMC and NPLMC distribution of ST cells (middle right), and OXPHOS scores of ST cells (right) in lymph nodes and tumors. (B) UMAP plot of 40,053 epithelial cells in ST data, showing 16 clusters. Each cluster is shown in a different color. (C) Heatmap showing the mean expression of 7 gene sets (co-stimulatory, cytotoxicity, inhibitory, MHC-II, naïve, pre-dysfunction, and resident) from T cells around PLMCs and NPLMCs at a radius of 40 μm in ST data. Color bar coding from blue to red indicates expression levels from low to high. (D) Representative images of Fixable Viability Dye 640/770, CD4, and CXCR5 staining to identify Tfh cells, and PD-1 expression to assess exhaustion in MetLNs from mice bearing orthotopic MOCCI-positive or control lung tumors. (E) Bubble plots showing the significantly differentially expressed ligand-receptor pairs of B cells, T cells, mast cells, PLMCs, and NPLMCs in the lymph nodes of ST. Dot size reflects P-values, and dot color represents communication probabilities. (F) Schematic illustration. (G) Bar plots of the proportions of GDF15+/− cells in PLMCs, NPLMCs, epithelial cells, and other ST cells. (H) Heatmap showing mean expression of 7 gene sets (co-stimulatory, cytotoxicity, inhibitory, MHC-II, naïve, pre-dysfunction, and resident) from T cells around GDF15+/− PLMCs at a radius of 40 μm in ST data. Color bar coding from blue to red indicates expression levels from low to high. (I) Representative example of a LUAD section stained by mIF with anti-CXCR5 (blue), anti-PanCK (light blue), anti-TGFBR2 (yellow), anti-GDF15 (red), anti-CD4 (purple), and anti-MOCCI (green). Arrows indicate specific cell types. MOCCI+ cancer cells (MOCCI+/PanCK+) with GDF15 secretion (orange arrow); Tfh cells (CD4+/CXCR5+) (red arrow) with GDF15-TGFBR2 binding (white arrow). MetLN, metastatic lymph node; MIA, multimodal intersection analysis; NPLMC, non-primary lymph node metastatic cell; OXPHOS, oxidative phosphorylation; PLMC, primary lymph node metastatic cell; ST, spatial transcriptomes; Tfh, follicular helper T cell; uiLN, uninvolved lymph node; UMAP, Uniform Manifold Approximation and Projection.

To explore the spatial features of PLMC clusters in detail, we applied UMAP clustering to the epithelial cells from the ST data in MetLNs, and we identified 16 epithelial clusters (Figure 5B). We performed MIA to map single-cell-defined clusters onto the spatially resolved data (Figure S6C)¹⁶. Similarly, UMAP clustering and MIA were also performed on epithelial cells from primary tumors (Figure S6D, E). In line with our previous findings, PLMC clusters (I, II, and VI) exhibited greater OXPHOS and EMT signature scores than NPLMCs in both tumor and MetLN tissues (Figures 5A, S6F, G, and Table S5), a finding highlighting their distinct metabolic and phenotypic states.

To better dissect the TME surrounding PLMCs, we used Squidpy³⁹ to delineate the cellular neighborhood of each malignant cell in MetLNs at 3 spatial scales (radii of 110, 340, and 648 pixels), which yielded approximately 3 (directly touching), 30 (local niche), and 100 (distal niche) neighboring cells, respectively⁴⁰. This multiscale analysis enabled us to distinguish cells in direct physical contact with tumor clones from those occupying the local niche and those situated at more distal locations. Notably, the most common neighboring cell type was other tumor cells, and a substantial proportion of cells neighboring PLMCs were also PLMCs (59.8%, 46.6%, and 41.1% at the respective distances). Similarly, for NPLMCs, a considerable percentage of neighbors (75.2%, 70.2%, and 66.5% at the respective distances) were also NPLMCs (Figures S6H and S7A, B). As expected, PLMCs had higher proportions of cell neighboring types, such as CD4⁺ T cells, CD8⁺ T cells, B cells, and fibroblasts, at all distances, thus suggesting a globally intensified inflammatory microenvironment.

To further elucidate how PLMCs adapt to and potentially reshape this immunologically enriched niche and facilitate metastatic colonization, we further stratified the B cells and T cells in LNs, and investigated the changes in PLMC neighboring cells (Figure S8A–F). Neighboring cells were selected at a 40 μm radius to capture only local niche interactions, and overlapping neighbors shared by both PLMCs and NPLMCs were excluded. DEG and KEGG pathway enrichment analyses were performed between the neighboring B cells of PLMCs and NPLMCs. As expected, the B cells around PLMCs, compared with NPLMCs, showed markedly lower enrichment in the ‘regulation of protein-containing complex assembly’ pathway, which was indicative of a less active phenotype (Figure S8G). The T-EFF cells around PLMCs also showed a less cytotoxic phenotype with lower expression of GZMA and CST7 compared to the T-EFF cells around NPLMCs (Figure 5C). Of note, when comparing the phenotypes of immune cells within vs. outside a 40 μm radius around PLMCs, we observed that T cells located beyond this 40 μm radius exhibited a significantly higher cytotoxic signature than those located within it (Figure S8H). Likewise, B cells beyond this 40 μm boundary showed significantly greater enrichment in the ‘B cell activation’ and ‘leukocyte proliferation’ pathways, thus suggesting a suppressive phenotype within a 40 μm radius (Figure S8I). Together, these findings suggested a potential mechanism through which PLMCs might establish an immunosuppressive niche in the local microenvironment.

We observed markedly greater PD-1 expression on CD4⁺CXCR5⁺ Tfh cells within metastatic LNs and orthotopic lung tumor tissue from mice bearing MOCCI-positive orthotopic lung tumors than control tumors, according to flow cytometry (Figures 5D and S8J). Therefore, the presence of MOCCI+ tumors promoted Tfh cell exhaustion within the draining LN microenvironment as well as in orthotopic lung tumor tissue. These findings suggested that MOCCI-derived cytokines intrinsically reshape the immune landscape of MetLNs and orthotopic tumor tissue, and drive Tfh exhaustion during early metastatic dissemination. To identify key signaling axes involved in TME remodeling, we systematically screened all activated ligand–receptor pairs predicted by CellChat in the ST dataset (Figure 5E). Among these interactions, GDF15 emerged as the most differentially expressed ligand between PLMC and NPLMC populations, and showed notably high expression in PLMCs. Interestingly, this interaction was also observed in our scRNA-seq data (Figures 3I and 4H). Given the previously reported immunosuppressive role of GDF15 in tumor progression^41,42, we hypothesized that PLMCs might secrete GDF15, thereby establishing an immunosuppressive niche in the local microenvironment (Figure 5F). As expected, PLMCs exhibited a higher frequency of GDF15+ cells (30.6%) than did NPLMCs (15.9%), normal epithelial cells (17.4%), and other cells (7.9%); PLMCs therefore might have a critical role as the primary source of GDF15 (Figure 5G). Survival analysis based on TCGA data indicated that the elevated GDF15 expression significantly correlated with poorer survival in LUAD and lung squamous cell carcinoma (Figure S8K, L). Notably, GDF15+ PLMC neighboring cells had a significantly lower proportion of plasma cells, and a higher proportion of naïve-B cells and T-EXH cells, than GDF15− PLMC neighboring cells, thus suggesting a locally immunosuppressive environment (Figure S8M). To demonstrate how GDF15+ PLMCs modulate neighboring immune cells through direct interaction, we selected TGFBR2+ neighboring cells (the predicted activated receptor from CellChat, also reported previously⁴³) within a 40 μm radius. We investigated the functional genes of TGFBR2+ T cell subsets within GDF15+ and GDF15− neighboring cells (Figure 5H). Within the GDF15+ subset, Tfh cells exhibited the greatest gene-expression variation, characterized by marked up-regulation of inhibitory signatures (CTLA4, LAG3, and PDCD1) and down-regulation of MHC-II signatures, indicative of an exhausted phenotype of Tfh. Given the well-established role of Tfh cells in facilitating B cell activation ^44,45, we inferred that PLMCs, probably via secreting GDF15 acting through TGFBR2, fostered Tfh cell exhaustion and consequently suppressed the activation of neighboring B cells (Figure 5F). Multiplex immunofluorescence histochemistry also indicated GDF15-TGFBR2 ligand-receptor binding, and close contact between GDF15-expressing MOCCI cancer cells and Tfh cells (Figure 5I). Therefore, PLMCs rebuild a locally immunosuppressive niche via the GDF15-TGFBR2 axis and ultimately contribute to a metastasis-permissive microenvironment.

MOCCI regulates OXPHOS reprogramming and drives lymph node metastasis through intracellular ROS levels

To evaluate the role of MOCCI in lung cancer metastasis, we used short hairpin RNA to generate PC9 and A549 lines with low MOCCI expression. In Transwell assays, knockdown of MOCCI significantly inhibited both migration and invasion of lung cancer cells, with reductions exceeding 50% compared to control cells (Figure 6A). These findings indicated that MOCCI promoted lung cancer migration and invasion. Similarly, wound healing assays showed that shMOCCI cells exhibited significantly slower gap closure compared to control cells, thus demonstrating a lower speed of collective motion of cells (Figure 6B). IHC staining was next performed on LNs and primary tumors from patients with or without LN metastasis. MOCCI protein expression was frequently positive in metastatic LNs. Notably, primary tumors from patients with LN metastasis also showed a higher proportion of positive staining than those from patients without metastasis (Figure 6C). Moreover, intrapulmonary delivery of PC9 cells overexpressing MOCCI (PC9-ovMOCCI) significantly enhanced lymphatic dissemination in vivo: serial bioluminescence imaging indicated markedly elevated photon flux in the LN metastatic mice. Analysis of dissected lungs and LNs confirmed a greater number of metastatic LNs in the ovMOCCI group than the control group (Figure 6D, E). H&E staining was also used to verify LN metastasis (Figure 6F).

• Download figure
• Open in new tab
• Download powerpoint

• Download figure
• Open in new tab
• Download powerpoint

Figure 6

Upregulation of MOCCI promotes lung cancer cell invasion, migration, OXPHOS, and LN metastasis. (A) Bar plots showing that downregulation of MOCCI significantly inhibits the cell invasion and migration ability of PC9 and A549 cells in Transwell assays (right). Representative images randomly selected from PC9 and A549 cells (left). (B) Bar plots showing that downregulation of MOCCI significantly inhibits the cell migration ability of PC9 and A549 cells in wound healing assays (right). Representative images randomly selected from PC9 and A549 cells (left). (C) IHC images of MOCCI expression in primary tumors with or without lymph node metastasis. Three independent experiments were performed and generated similar results. Scale bar = 100 μm. D, E, F) Workflow of the LN metastasis mouse model (D, left). Bioluminescence imaging indicated markedly elevated photon flux in an LN metastatic mouse (D, right). The dissected LNs confirmed a greater number of metastatic lymph nodes than observed in control mice (E), as confirmed by H&E staining (F). (G) Line chart for Seahorse XF mitochondrial stress testing in PC9 cells (left). Bar plot showing significant differences in basal respiration, maximal respiration, and ATP production in the shMOCCI, ovMOCCI, and control groups (right). (H) Representative fluorescence images of DCFDA (total ROS) and mitoSOX™ Red (mitochondrial ROS) with or without NAC rescue of ROS levels (left). Box plots showing mean fluorescence intensity quantification of DCFDA and mitoSOX™ Red among 5 groups (right). Boxes represent the median (indicated by the horizontal line), and encompass the second to third quartiles (forming the body of the box); Tukey-style whiskers extend beyond the box to depict the data range. The Wilcoxon rank-sum test was applied to assess the difference. DCFDA, 2′,7′-dichlorodihydrofluorescein diacetate; IHC, immunohistochemistry; LN, lymph node; MOCCI, Modulator of Cytochrome C Oxidase during Inflammation; NAC, N-acetyl-L-cysteine; OXPHOS, oxidative phosphorylation; ROS, reactive oxygen species.

We then tested how MOCCI affected metabolic reprogramming. Seahorse XF mitochondrial stress testing revealed significantly greater basal oxygen-consumption rate, ATP-linked oxygen-consumption rate, and maximal respiratory ability in PC9-ovMOCCI cells than PC9 controls (n = 3, P < 0.01), thus resulting in significantly greater spare-respiratory ability (Figure 6G). Moreover, A549-shMOCCI and PC9-shMOCCI showed significantly down-regulated mitochondrial OXPHOS ability (Figures 6G and S9A). Therefore, a switch to a more bioenergetically flexible state produced by MOCCI met higher ATP demand and endowed the cells with a greater ability to cope with the metastatic process. To further explore the regulation mechanism of MOCCI, we performed JC-1 staining, which showed extensive depolarization of mitochondria in PC9 and A549 cells in the shMOCCI group, as evidenced by a shift from red aggregates to green monomers (Figure S9B). Intracellular ROS was tested with DCFH-DA staining. Mean DCF fluorescence was significantly elevated in PC9-ovMOCCI cells and significantly diminished in PC9-shMOCCI cells. In agreement with these findings, the significantly elevated mitochondrial ROS in ovMOCCI cells, as revealed by MitoSOX staining, indicated that MOCCI augments ROS generation in both cytosolic and mitochondrial compartments (Figure 6H). The finding that pre-treatment with 40 μM N-acetyl-L-cysteine (NAC) completely normalized the DCF and MitoSOX signals confirmed that the observed changes were ROS dependent. Finally, functional rescue experiments demonstrated that 40 μM NAC fully restored the migration and invasion abilities of PC9-ovMOCCI in 24 or 48 h (Figure S9C, D). We therefore concluded that MOCCI induces LN metastasis through increasing mitochondrial ROS levels and up-regulating OXPHOS. ELISA revealed significantly higher levels of GDF15 in the ovMOCCI group of both PC9 and A549 cells compared to control group, thus suggesting that OXPHOS reprogramming correlated with enhanced GDF15 secretion (Figure S9E). We then treated the ovMOCCI groups of both PC9 and A549 cells with OXPHOS antagonists and measured GDF15 secretion via ELISA. GDF15 levels markedly decreased after OXPHOS inhibition. This key result demonstrated that the upregulation of GDF15 was dependent on enhanced OXPHOS activity. Therefore, we concluded that GDF15 induction was not an independent, parallel event but a direct consequence of the MOCCI-driven boost in mitochondrial OXPHOS (Figure S9E). However, further validation of this possibility is necessary.