Integrated single-cell transcriptomic analysis of brain and lymph node metastases in non-small cell lung cancer

Abstract

BACKGROUND

Non-small cell lung cancer (NSCLC) is the leading cause of cancer-related mortality worldwide, and metastasis remains the primary driver of treatment failure and death. Despite recent advances in single-cell RNA sequencing (scRNA-seq), most studies have focused on isolated metastatic sites; there is still a lack of studies that simultaneously focus on the cellular and molecular characteristics of brain and lymph node metastases in NSCLC.

AIM

To explore features of brain and lymph node metastases in NSCLC and to nominate candidate pathways for subsequent functional interrogation.

METHODS

We integrated three public scRNA-seq datasets (GSE131907, GSE186344, and GSE256491) to construct a high-resolution transcriptomic atlas of NSCLC. These scRNA-seq cohorts were heterogeneous and not clinically matched primary-metastasis pairs, comprising a total of 45 samples (including 11 primary tumors, 7 lymph node metastases, 15 brain metastases, and 12 normal controls) and 208506 high-quality cells. We systematically dissected cellular heterogeneity, the distribution of tumor epithelial subpopulations in different metastatic sites, pathway activation, and immune microenvironment remodeling during metastasis. To assess the relationship between the expression of the tumor subclusters and overall survival of lung adenocarcinoma (LUAD) patients, we conducted a multivariate Cox regression analysis using The Cancer Genome Atlas (TCGA) data with ssGSEA gene-signature scores. Multiplex immunohistochemistry (mIHC) was performed on a cohort of 23 patients (Mets1: 12 cases of lymph node metastasis, including 2 with distant metastasis; Mets0: 11 cases of non-metastasis) to preliminarily validate macrophage infiltration and unfolded protein response (UPR)-related protein expression.

RESULTS

Metastatic lesions exhibited marked enrichment of T cells and myeloid cells. Immune landscape analysis revealed a dual immunosuppressive microenvironment comprising exhausted T cells overexpressing PDCD1, HAVCR2, CTLA4, and LAG3, and increased macrophage infiltration. The distribution of tumor epithelial subpopulations varied in different metastatic foci. The UPR was abnormally activated in both metastatic sites. TCGA data analysis revealed that high expression of the tumor epithelial subpopulation (C14_PBK_Tumor) was associated with inferior overall survival in LUAD patients (multivariate Cox regression analysis, hazard ratio = 2.39, 95%CI: 1.15-4.97, adjusted P = 0.02). Consistently, the mIHC results showed that in the metastasis group, the levels of CD68⁺CD163⁺ macrophages and cells triple-positive for UPR markers (XBP1, ATF6, and ATF4) were higher than those in the non-metastasis group. (Mann–Whitney U test; data are presented as medians; CD68⁺CD163⁺ cell rate: Met0: 1.913 vs Met1: 3.338, P = 0.0129; XBP1⁺/ATF6⁺/ATF4⁺ cell rate: Met0: 60.19 vs Met1: 73.26, P = 0.0028).

CONCLUSION

These findings suggest that in brain and lymph node metastatic sites of NSCLC, there is abnormal activation of the UPR, along with a dual immunosuppressive microenvironment and differential distribution of tumor epithelial subpopulations. However, further functional and prospective studies are needed to confirm these observations.

Key Words: Non-small cell lung cancer; scRNA-seq; T cells; Macrophages; Unfolded protein response; Brain metastasis; Lymph node metastasis

Core Tip: Our study offers a panoramic analysis of two distinct metastatic sites within the same disease context. Our main findings are: The distribution of tumor epithelial subpopulations differs across metastatic foci; abnormal activation of unfolded protein response (UPR) occurs in both metastatic sites; and a dual immunosuppressive microenvironment is characterized. Multivariate Cox regression analysis based on The Cancer Genome Atlas data revealed that, among the tumor epithelial subpopulations, the C14_PBK_Tumor subtype was associated with worse prognosis in lung adenocarcinoma patients. Importantly, multiplex immunohistochemistry was performed to assess macrophage infiltration and the expression of UPR-related proteins. These findings still require further verification through functional experiments and prospective studies.

INTRODUCTION

Lung cancer is the primary cause of cancer-related fatalities globally. Its pathogenesis is complex, particularly in the dynamic interaction between the immune system and tumors[1,2]. Non-small cell lung cancer (NSCLC) is the most common histological type of lung cancer, accounting for roughly 85% of all cases[3-5]. As a malignancy with high incidence and mortality worldwide, NSCLC causes nearly 1.8 million deaths annually, representing a major public health challenge[6,7]. Despite advances in surgery, chemoradiotherapy, targeted therapy, and immunotherapy, the prognosis of patients at advanced stages remains poor, with a 5-year survival rate below 20%[8-11]. Metastasis is the primary cause of treatment failure and death; however, the underlying molecular mechanisms remain incompletely understood[12-15].

The immune system plays a vital role in the progression of NSCLC. Tumor-infiltrating lymphocytes, especially T cells, are crucial mediators of the anti-tumor immune response[16]. However, tumor cells can reshape the tumor microenvironment (TME) to establish an immunosuppressive state, thereby evading immune surveillance[17]. The enrichment of regulatory T cells (Tregs) and M2-type tumor-associated macrophages (TAMs), along with the functional impairment of effector T cells, collectively forms an immunosuppressive network that promotes malignant progression[18,19].

Tumor cells adapt to distinct organ microenvironments to evade immune surveillance and form supportive niches. Powell et al[20] reviewed evidence showing that transcriptional regulatory adaptations are widely present during brain metastasis. During brain metastasis in breast cancer, subtype switching and TME interactions are key factors driving metastatic progression[21]. scRNA-seq has emerged as a powerful tool for exploring the mechanisms of tumor metastasis. Liu et al[22] comprehensively reviewed the advantages of scRNA-seq for elucidating the mechanisms of lymph node metastasis. Recent scRNA-seq studies have begun to unravel the heterogeneity of the tumor immune microenvironment (TIME) in NSCLC metastases. For example, Tagore et al[23] found a chromosomally unstable tumor subpopulation that drives brain metastasis, whereas Zhang et al[24] showed metabolic reprogramming in lymph node metastases. However, most studies focus on single metastatic sites; it remains unclear whether there are common molecular mechanisms between brain metastasis and lymph node metastasis.

We selected and integrated three public scRNA-seq datasets (GSE131907, GSE186344, and GSE256491) that collectively cover primary tumors, lymph node metastases, brain metastases, and normal tissues from NSCLC patients, with all samples derived from surgical resections or biopsies. To ensure comparability across cohorts, we applied uniform quality control and batch-effect correction using Harmony. Our main objective was to conduct single-cell transcriptomic analysis of brain metastasis and lymph node metastasis in NSCLC. We further sought to explore the immunosuppressive landscape across distinct microenvironments, identify potential activation pathways, and determine whether tumor subpopulation expression at different metastatic sites correlates with overall survival in lung adenocarcinoma (LUAD) patients. The research workflow is illustrated in Supplementary Figure 1.

MATERIALS AND METHODS

Single-cell data processing

The Seurat package (version 5.5.1)[25] was used for single-cell data processing and analysis. For three datasets (GSE131907, GSE186344, and GSE256491), quality control was performed individually. First, for GSE131907, a Seurat object was created using the parameters min.cell set to 3 and min.features set to 200. Quality control criteria included a mitochondrial gene ratio ≤ 20%, and only cells with RNA read counts from 1000 to 100000 and detected gene numbers ranging from 300 to 10000 were retained for further analysis. Second, for GSE186344, a Seurat object was similarly created with min.cells set to 3 and min.features set to 200. Quality control was performed for each sample. Third, for GSE256491, a Seurat object was created using min.cells set to 3 and min.features set to 50; quality control criteria included a mitochondrial gene ratio ≤ 25%, and only cells with RNA read counts from 1000 to 20000 and detected gene numbers ranging from 300 to 6000 were retained for further analysis. The three datasets were combined via the merge function. scRNA-seq dataset characteristics are listed in Table 1. Detailed information can be found in Supplementary Table 1.

Table 1 scRNA-seq cohorts.

Dataset	Sample	Tissue origin	Histology	Sex	Age (years)
GSE131907	BRONCHO_11	mLN	ADC	Female	81
GSE131907	EBUS_10	mLN	ADC	Male	68
GSE131907	EBUS_12	mLN	ADC	Female	60
GSE131907	EBUS_13	mLN	ADC	Male	53
GSE131907	EBUS_15	mLN	ADC	Male	64
GSE131907	EBUS_19	mLN	ADC	Female	68
GSE131907	EBUS_51	mLN	ADC	Male	43
GSE131907	LN_01	nLN	ADC	Male	61
GSE131907	LN_02	nLN	ADC	Female	66
GSE131907	LN_03	nLN	ADC	Male	67
GSE131907	LN_04	nLN	ADC	Male	67
GSE131907	LN_05	nLN	ADC	Female	56
GSE131907	LN_06	nLN	ADC	Female	63
GSE131907	LN_07	nLN	ADC	Male	68
GSE131907	LN_08	nLN	ADC	Female	77
GSE131907	LN_11	nLN	ADC	Male	74
GSE131907	LN_12	nLN	ADC	Male	62
GSE131907	LUNG_T06	tLung	ADC	Male	54
GSE131907	LUNG_T08	tLung	ADC	Female	60
GSE131907	LUNG_T09	tLung	ADC	Male	67
GSE131907	LUNG_T18	tLung	ADC	Male	62
GSE131907	LUNG_T19	tLung	ADC	Male	53
GSE131907	LUNG_T20	tLung	ADC	Male	55
GSE131907	LUNG_T25	tLung	ADC (Double)	Male	69
GSE131907	LUNG_T28	tLung	ADC (Double)	Male	73
GSE131907	LUNG_T30	tLung	ADC	Male	56
GSE131907	LUNG_T31	tLung	ADC	Male	65
GSE131907	LUNG_T34	tLung	ADC	Male	64
GSE131907	NS_02	mBrain	ADC	Female	46
GSE131907	NS_03	mBrain	ADC	Female	81
GSE131907	NS_04	mBrain	ADC	Male	39
GSE131907	NS_06	mBrain	ADC	Male	59
GSE131907	NS_07	mBrain	ADC	Female	64
GSE131907	NS_12	mBrain	ADC	Female	59
GSE131907	NS_13	mBrain	ADC	Female	50
GSE131907	NS_16	mBrain	ADC	Male	69
GSE131907	NS_17	mBrain	ADC	Female	56
GSE131907	NS_19	mBrain	ADC	Female	44
GSE186344	Lung-1	Frontal	Small cell neuroendocrine carcinoma	Male	76
GSE186344	Lung-2	Frontal	NSCLC	Female	57
GSE186344	Lung-3	Temporal	NSCLC	Female	62
GSE256491	TL_patient9	Right temporal lobe	NSCLC	Male	25
GSE256491	TL_patient10	Right temporal lobe	NSCLC	Male	24

NSCLC: Non-small cell lung cancer; ADC: Adenocarcinoma; nLN: Normal lymph nodes; mLN: Lymph node metastases; mBrain: Brain metastases..

Data were normalized by using NormalizeData with default parameters. We employed the variance-stabilizing transformation method to identify highly variable genes, with features set to 2000. The single-cell expression data were scaled and centered using the ScaleData function with default parameters. Dimensionality reduction was performed using RunPCA with npcs set to 50. Batch effects were corrected and data integration was carried out using IntegrateLayers by the HarmonyIntegration (v2.0.2)[26] method, applying the parameters theta = 2, sigma = 0.1, and lambda = 1. UMAP dimensionality reduction was applied using RunUMAP (dims = 1:20). The robustness of this integration was further supported by silhouette score analysis, which confirmed both the correction of batch effects and the retention of biological variability (Supplementary Figure 2A). A nearest neighbor graph was constructed using FindNeighbors (dims = 1:20), and cell clusters were identified using FindClusters with a resolution of 0.4.

To remove potential doublets, we applied DoubletFinder (v2.0.4)[27] to each sample independently with the following parameters: PN = 0.25, pK was automatically determined as the value corresponding to the peak of the BCmetric distribution using the find.pK function, sct = FALSE, reuse.pANN = FALSE, and nExp was set to 7.5% of the total cell number for each sample (i.e., total cells × 0.075). Doublets were then excluded from downstream analyses.

For biological annotation of the resulting cell clusters, the FindAllMarkers function in the Seurat package was used to identify specific marker genes for each cluster. To ensure the reliability of the identified markers, we applied the following criteria: The minimum detection proportion (min.pct) of a gene in the target cluster was at least 25%, and the average log fold change (logfc.threshold) of the gene between the target cluster and all other clusters was at least 0.5 (log² transformed fold change). Statistical significance of differential expression was evaluated via the Wilcoxon rank sum test, and P values from multiple hypothesis testing were corrected using the Bonferroni method to control the false positive rate.

Copy number variation inference

Copy number variation (CNV) inference was performed based on single-cell transcriptomic data using the R package InferCNV (v1.18.1). Gene position files derived from the human reference genome GRCh38 (hg38) were used for chromosomal ordering. Endothelial cells (ECs) and myofibroblast-like cells (Myo_fibro) were designated as the diploid reference baseline. The core analysis was executed via the infercnv:: Run() function with the following key parameter settings: Window_length = 101 (sliding window size, i.e., the number of genes included in each moving-average calculation), cutoff = 0.1 (a threshold appropriate for 10 × Genomics platform data, filtering genes with excessively low average expression in reference cells), denoise = TRUE (enabling denoising filtering based on the expression distribution of reference cells to suppress false-positive signals attributable to technical noise), cluster_by_groups = TRUE (performing hierarchical clustering separately according to cell annotation information), and num_threads = 10 (multi-threaded parallel computation). This method calculates relative log₂ fold-change values by comparing gene expression intensities within chromosomal sliding windows between tumor cells and reference cells. Following centering, smoothing, and noise reduction, an inverse-log-transformed visualization matrix was generated. The resulting chromosomal heatmap was used to identify malignant cells exhibiting genomic instability.

Cell subset analysis and annotation

Subset analyses were performed separately for epithelial cells, myeloid cells, and the T cell lineage. After data scaling with ScaleData, principal component analysis was carried out (npcs = 50), and the first 20 principal components were selected for UMAP dimensionality reduction and nearest neighbor graph construction. Clustering was conducted via the FindClusters algorithm, with resolutions determined by the lineage complexity (epithelial and myeloid cells: 0.4; T cells: 0.8). The optimal resolution for each cell type was determined by examining the cluster tree, and the robustness of this selection was supported by the clustering sensitivity analysis across a resolution range from 0.2 to 1.6 (Supplementary Figure 2B-D). Marker genes for each subcluster were identified using FindAllMarkers (min.pct = 0.25, logfc.threshold = 0.25), and cell subsets were annotated based on these markers combined with known cell type markers. A summary of cell type composition and counts is provided in Supplementary Table 2.

Analysis of cell type distribution and proportion differences across tissues

To systematically evaluate the distribution characteristics of each cell type across different tissues, we first calculated the relative proportion of each cell type in each sample using the prop.table function and visualized these proportions using boxplots (ggboxplot) from the ggpubr (v0.6.3) package. Intergroup comparisons were performed using the Wilcoxon rank-sum test, with P < 0.05 considered statistically significant.

Furthermore, to quantitatively assess the distribution preference of each cell type across tissues, we calculated the odds ratio (OR) for each cell type in each tissue using the Fisher’s exact test. The Benjamini-Hochberg method was applied to control the false discovery rate (FDR), with FDR < 0.05 considered statistically significant. An OR > 1 indicated enrichment of the cell type in the corresponding tissue, while an OR < 1 indicated depletion. Results were displayed as a heatmap.

Differential expression analysis and functional enrichment

To identify differentially expressed genes (DEGs) between the metastatic and primary groups, we employed two complementary strategies.

First, at the single-cell level, DEGs were identified using the FindMarkers function from the Seurat package using the Wilcoxon rank-sum test. The filtering criteria were: Adjusted P < 0.05, |avg_log2FC| > 0.5, and a detection rate (pct) > 0.25 in either group. Specifically, genes with avg_log2FC > 0.5, adjusted P < 0.05, and pct.1 > 0.25 were defined as upregulated; genes with avg_log2FC < -0.5, adjusted P < 0.05, and pct.2 > 0.25 were defined as downregulated. Results were visualized as volcano plots using the ggplot2 (v4.0.3) package.

Second, to validate the DEGs and mitigate the sparsity of single-cell data, we aggregated the single-cell count matrix into a pseudobulk expression matrix by sample ID using the AggregateExpression function from the Seurat package, yielding the sum of gene expression for all cells in each sample. DEGs were identified using the R package DESeq2 (v1.38.3). Briefly, the raw count matrix was extracted from single-cell or bulk RNA-seq data and used to construct a DESeqDataSet object. To ensure analytical reliability, genes with excessively low total counts across all samples were filtered out. Specifically, only genes with counts ≥ 2 in at least one sample within a given group were retained, thereby reducing the burden of multiple-testing correction imposed by low-expression genes.

The filtered count matrix was subsequently used to construct a DESeqDataSetFromMatrix object. The sample phenotypic information (colData) contained the grouping factor condition, with primary tumors designated as the reference level. The design formula was set to condition, indicating that intergroup differences were the primary focus. The DESeq() function was then called to execute the standard differential analysis pipeline, which sequentially performed the following steps: (1) Estimation of size factors for each sample using the “ratio” method to normalize for library size; (2) Estimation of gene-wise dispersion with parametric fitting (fitType = “parametric”), followed by dispersion shrinkage via the empirical Bayes approach to improve estimation stability under small sample sizes; and (3) Fitting of a negative binomial generalized linear model, with the Wald test (test = “Wald”) applied to assess the statistical significance of differential expression between the metastatic group and the primary tumor group. Throughout this process, the default betaPrior = FALSE was used to enable independent shrinkage of log₂ fold change (LFC) in subsequent steps. To mitigate false positives arising from low-expression and high-dispersion genes, the lfcShrink() function with the “apeglm” method was employed to obtain shrunken LFC estimates, yielding more robust effect sizes. Finally, the differential analysis result table was extracted using the results() function. Genes with P < 0.05 and |log2FoldChange| > 0.5 were defined as significantly differentially expressed.

Gene ontology functional enrichment analysis and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis were performed separately on the upregulated and downregulated genes identified by both strategies, using the clusterProfiler (v4.6.2) and org.Hs.eg.db (v3.16.0) packages. The top 10 significantly enriched pathways (adjusted P < 0.05) were selected for graphical representation.

Pathway analysis between tumor and normal groups

The human hallmark gene set (version h.all. v2022.1.Hs. symbols), comprising 50 representative pathways, was downloaded from the MSigDB database. After extracting the RNA expression matrix, the pathway enrichment scores for each sample were calculated using the ssGSEA algorithm, which is implemented in the GSVA (v1.50.5) method to assess pathway activity levels. Subsequently, differential pathway activity between the tumor and normal groups was compared using the empirical Bayes method from the limma (v3.54.2) package, with a threshold of P < 0.05 for selecting significantly different pathways. Finally, the identified differential pathways were visualized with the help of the ggplot2 (v4.0.3) package. Summary statistics for main comparisons are listed in Supplementary Table 3.

Prognostic and metastasis association analysis using The Cancer Genome Atlas data

First, marker genes for each cell subcluster were identified using the FindAllMarkers function in the Seurat package, with the parameter min.pct = 0.25 used to screen genes with relatively high expression frequency in the target cell cluster, and logfc.threshold = 0.25 as the default threshold. For each subcluster, the top 20 most significantly differentially expressed genes (ranked by descending avg_log2FC) were selected to construct a gene signature set for that subcluster. The RNA-seq expression matrix and clinical information (including survival time, survival status, and metastasis status; Mets0: Non metastatic, Mets1: Metastatic) for LUAD were downloaded from The Cancer Genome Atlas (TCGA) database. The ssGSEA algorithm was employed to calculate the enrichment score of each gene signature set in each TCGA sample, thereby evaluating the activity level of different cell subcluster signatures in individual tumor samples. Patients were divided into high score and low score groups based on the median ssGSEA score of each gene set. Kaplan-Meier survival curves (v3.5.3) were plotted, and survival probabilities were estimated using the survfit() function from the survival package. Differences in survival between the two groups were compared using the log rank test [survdiff() function], with P < 0.05 considered significant. Survival curves were visualized using the ggsurvplot() function from the survminer (v0.5.2) package.

For multivariate analysis, variables with P < 0.10 in the univariate Cox regression analysis (including sex, age, and clinical stage) were incorporated into a multivariate Cox proportional hazards regression model. The Wald test was used to evaluate the independent prognostic impact of each covariate, and results are presented as hazard ratios (HR) with 95% confidence intervals. P < 0.05 was considered statistically significant.

Additionally, samples were divided into metastatic and non-metastatic groups according to the patients’ metastasis status. Differences in the distribution of ssGSEA scores for each gene set between the two groups were compared using the Wilcoxon rank sum test [wilcox.test() function], with P < 0.05 as the significance threshold. Boxplots were created via the ggboxplot() function from the ggpubr (v0.6.3) package to display the differences in gene set activity scores between the metastatic and non-metastatic groups. Complete Cox regression results are shown in Supplementary Table 4. We conducted a verification using another dataset, GSE72094.

Expression characteristics of immune cell exhaustion-related genes and immune checkpoint molecules in different tissue microenvironments

To analyze the expression distribution of genes related to immune cell exhaustion and immune checkpoint molecules across various tissue microenvironments, we selected the exhaustion-related genes and immune checkpoint molecules as the detection indicators. Tissue type was used as the grouping variable, including mLN, nLN, mBrain, and nBrain. The average expression level and the percentage of cells expressing each gene in the indicated tissues were visualized using the DotPlot function. In the dot plot, each circle’s color indicates the average expression level of the gene in that tissue, and the size of the circle represents the percentage of cells expressing the gene in that tissue. This visualization enables a direct comparison of the distribution characteristics of immune exhaustion and inhibitory signals across different tissue microenvironments. We downloaded the new dataset of lymph node metastasis, GSE293360, and the original dataset, GSE131907, for integrated analysis.

mIHC staining of tissue sections

Tissue microarray chips used in this study were obtained from Hunan AiFang Biological Co., Ltd. All participants provided written informed consent, and the study protocol was approved by the Bioscience and Technology Ethics Review Committee of Hunan AiFang Biotechnology Co., Ltd. (approval No. HN20250401). mIHC cohort characteristics are listed in Table 2, and details of the sample composition are provided in Supplementary Table 5.

Table 2 Multiplex immunohistochemistry cohorts.

Group	Pathological number	Gender	Age (years)	Histology	Stage	TNM	Treatment regimens
Mets0	1037557	Male	58	ADC	I	T2aN0M0	Surgery
Mets0	1054412	Female	82	ADC	II	T2N0M0	Surgery
Mets0	1117778	Female	48	ADC	I	T2aN0M0	Surgery
Mets0	1117102	Female	61	ADC	I	T1N0M0	Surgery
Mets0	1123096	Male	58	ADC	I	T1N0M0	Surgery
Mets0	1121625	Male	51	ADC	II	T3N0M0	Chemotherapy
Mets0	1088996	Male	70	ADC	I	T1N0M0	Surgery
Mets0	1119756	Female	54	ADC	I	T1N0M0	Surgery
Mets0	1124448	Male	61	ADC	II	T2bN0M0	Surgery
Mets0	1136112	Male	68	ADC	II	T2bN0M0	Chemotherapy
Mets0	1135011	Male	60	ADC	II	T2bN0M0	Chemotherapy
Mets1	1090522	Female	48	ADC	III	T2aN2M0	Chemotherapy
Mets1	1092316	Female	53	ADC	II	T2aN1M0	Surgery
Mets1	641032	Male	50	ADC	IV	T1bN3M1b	Chemotherapy
Mets1	1099920	Female	50	ADC	II	T1bN1M0	Surgery
Mets1	1113168	Female	53	ADC	III	T2N2M0	Surgery
Mets1	1118269	Male	40	ADC	II	T2aN1M0	Surgery
Mets1	1120708	Female	58	ADC	III	T3N1M0	Surgery
Mets1	1117066	Female	79	ADC	IV	T2N1M1	Surgery
Mets1	1113386	Male	62	ADC	III	T3N2M0	Surgery
Mets1	1131710	Male	52	ADC	III	TIcN2M0	Surgery
Mets1	1128781	Male	63	ADC	III	T3N2M0	Surgery
Mets1	1145641	Male	43	ADC	II	T2N1M0	Surgery

ADC: Adenocarcinoma.

The collected tissue sections were dewaxed with standard xylene and hydrated with graded alcohol. Then, mIHC staining was done with a seven-color multiplex immunofluorescence staining kit (AFIHC027, AiFang Biological, China). Antigen retrieval was performed via microwave heating, and the sections were incubated with 3% hydrogen peroxide solution at room temperature for 15 minutes, and then incubated with 10% goat serum for 15 minutes. Primary antibody A was applied overnight at 4 °C, followed by PBST washes, and the sections were then incubated with Polymer-HRP anti-mouse/rabbit universal secondary antibody IgG (AFIHC001, AiFang Biological) for 30 minutes. After PBST washing, TYR fluorescent dye was added for 8 minutes, and the sections were washed three times with PBST. Antibody A was removed by microwave treatment, and the sections were washed three times with PBST. After blocking with goat serum, primary antibody B was added. The following antibodies were used for visualization: Anti-CD163 (AF20010, 1:200), anti-CD8 (AF20211, 1:200), anti-PD-1 (AF20083, 1:200), anti-ATF4 (AFRM0125, 1:2000), and anti-CD68 (AF20022, 1:4000) from AiFang Biological, as well as anti-ATF6 (Cat No. 24169-1-AP, 1:200) and anti-XBP1 (Cat No. 24868-1-AP, 1:2000) from Proteintech. After DAPI counterstaining at room temperature for 10 minutes in the dark and three PBST washes, the sections were mounted with an antifade mounting medium, and fluorescent images were captured using an eight-channel digital slide scanner (model AF-KL-20-8, AiFang Biological). Channel planning and exposure time are shown in Supplementary Table 6.

mIHC image acquisition and quantification

This study enrolled 23 patients, each serving as an independent biological replicate. Both image acquisition and data quantification were performed under blinded conditions. An eight-channel fluorescence digital slide scanner (AF-KL-20-8, AiFang Biological, China) was employed to capture images of mIHC slides, with the resulting images saved in .kfbf format. For image data analysis, we used the VISIOPHARM software (VISIOPHARM 2025.08.01 18881 x64, Hoersholm, Denmark) to perform quantitative assessments. To ensure consistency across analyses, fixed and uniform threshold parameters were established within the VISIOPHARM application for the identification and quantification of macrophages (CD68⁺CD163⁺), T cells (PD-1⁺/CD8⁺), and unfolded protein response (UPR) activation markers (XBP1⁺/ATF6⁺/ATF4⁺). Following the completion of the analyses, the resulting data were exported. All statistical evaluations were subsequently performed using GraphPad Prism 6.01. Patients were divided by TNM stage into Mets0 (no metastasis) and Mets1 (lymphatic/distant metastasis) groups. The Mann-Whitney U test was used to compare positive rate, positive density, and mean fluorescence intensity (MFI) between groups, with statistical significance set at P < 0.05. mIHC quantitative data are provided in Supplementary Table 6.

The following parameters were calculated: Total cell count (the total number of valid cells identified within the analysis region); positive cell count (the number of indicator-positive cells identified within the analysis area); positive rate = positive cell count/total cell count × 100%; positive cell density = positive cell count/area of the target tissue; and MFI I_i (signal intensity of the i^th valid fluorescent pixel within the target area), which is calculated using the following formula:

MFI=∑i=1nIin

where n is the total number of valid fluorescent pixels within the target area and ∑I_i is the sum of the signal intensities of all valid fluorescent pixels within the target area.

RESULTS

Single-cell transcriptomic landscapes of primary and metastatic NSCLC

To dissect the cellular heterogeneity and microenvironmental remodeling during NSCLC metastasis, we performed an integrated scRNA-seq analysis of adult human samples, including NSCLC primary tumors (Lung Cancer), lymph node metastases (mLN), normal lymph nodes (nLN), brain metastases (mBrain), and normal brain tissues (nBrain), using public datasets GSE131907, GSE186344, and GSE256491. After quality control and batch-effect correction, we obtained a high-resolution cellular atlas comprising all major cell types. Unsupervised clustering followed by canonical marker gene annotation identified T cells (CD3E, KLRB1), B cells (MS4A1, FCRLA), plasma cells (IGHG2, IGHG3), myeloid cells (FCN1, RETN, MARCO), mast cells (TPSAB1, CPA3), microglia (P2RY12, P2RY13), epithelial/tumor cells (FAM83A, GOLT1A), EC (CLEC14A, SOX18), myofibroblasts, neuroendocrine cells, oligodendrocytes (OPALIN, MOG, PLP1), and smooth muscle cells (Figure 1A-C).

Figure 1 Single-cell transcriptomic atlas of primary and metastatic non-small cell lung cancer. A: UMAP visualization of 174643 cells integrated from adult human non-small cell lung cancer primary tumors (lung cancer), lymph node metastases (mLN), normal lymph nodes (nLN), brain metastases (mBrain), and normal brain tissues (nBrain). Cells are color coded by annotated cell types: T cells, B cells, plasma cells, myeloid cells, mast cells, microglia, epithelial cells, endothelial cells, myofibroblasts (Myo_fibro), neuroendocrine cells, oligodendrocytes, and smooth muscle cells; B: Dot plot showing expression levels of selected canonical marker genes used for celltype identification; C: UMAP plot colored by tissue origin (nLN, nBrain, Lung Cancer, mBrain, and mLN); D: UMAP plot showing cell type distribution across the five tissue groups; E: Stacked bar plot showing cell proportion per tissue group; F and G: Bar plots highlighting changes in cell type composition between nBrain and mBrain (up) and between nLN and mLN (down). nLN: Normal lymph nodes; mLN: Lymph node metastases; mBrain: Brain metastases; nBrain: Normal brain tissues; ECs: Epithelial cells; SMC: Smooth muscle cells.

Comparison of cell-type proportions across different tissues revealed striking differences (Figure 1D). In nBrain, microglia were the predominant cell type, accounting for approximately 57% of all cells, while T cells and myeloid cells had almost disappeared. In mBrain, a profound shift in cellular composition was observed: Epithelial cells were the predominant cell type, accounting for approximately 38% of all cells, while the relative abundance of T cells increased to 18% and myeloid cells increased to approximately 6%. The proportions of oligodendrocytes, microglia, and SMC were also reduced in mBrain, reflecting the replacement of normal brain parenchyma by infiltrating tumor and immune cells (Figure 1E and F). Similarly, in mLN compared with nLN, we observed a significant increase in the proportions of T cells and myeloid cells, accompanied by a relative decrease in B cells (Figure 1E and G). As shown in Supplementary Figure 3, the UMAP projections of individual samples further confirmed the consistency of these cellular compositional shifts.

Taken together, the TIME is actively remodeled during metastatic progression. The striking increase in T cells and myeloid cells in metastatic sites prompted us to further characterize their functional states.

Distribution of tumor epithelial subpopulations in lymph nodes and brain metastases

To investigate whether specific tumor cell populations are associated with distinct metastatic sites, we extracted all epithelial/tumor cells from the integrated scRNA-seq dataset and performed unsupervised subclustering. This analysis identified two distinct epithelial subclusters and 14 distinct tumor subclusters, with the latter designated C1 through C14 and each defined by unique marker genes. The overall distribution of these subclusters across primary tumors and metastatic lesions is shown in UMAP plots (Figure 2A-C). In mLN, the tumor epithelial subclusters were predominantly C2-PSCA, C4-SPP1, and C7-TRPM8, indicating that these three subpopulations are selectively enriched in the mLN niche. In contrast, mBrain exhibited a broader and distinct composition. Notably, subclusters C1-CDH2 and C5-FABP7 were almost exclusively present in brain metastases (Figure 2D).

Figure 2 Different distribution of tumor epithelial subclusters in lymph node and brain metastases of non-small cell lung cancer. A: UMAP visualization of two epithelial subclusters and 14 tumor subclusters identified by subclustering analysis of 26587 epithelial/tumor cells derived from non-small cell lung cancer primary tumors (lung cancer), lymph node metastases (mLN), normal lymph nodes (nLN), brain metastases (mBrain), and normal brain tissues (nBrain). Each subcluster is annotated by its top marker genes; B: The UMAP shows the cell distribution of two epithelial subpopulations and 14 tumor subpopulations in the five tissue groups; C: Dot plot showing the expression levels of selected canonical marker genes used for cell type identification of the two epithelial subclusters and 14 tumor subclusters; D: Stacked bar plot showing the proportion of each tumor subcluster within different tissue types; E: UMAP plot showing the cell distribution of tumor and epithelial tissues in the five tissue groups; F: UMAP plot colored by tumor and epithelial tissues; G: UMAP plot colored by tissue origin (nLN, nBrain, Lung Cancer, mBrain, and mLN). nLN: Normal lymph nodes; mLN: Lymph node metastases; mBrain: Brain metastases; nBrain: Normal brain tissues.

UMAP analysis revealed that, in mLN and mBrain, the distribution of cell subpopulations varied (Figure 2E-G). The distribution of these subclusters across individual samples and the correlation analysis of representative subcluster signatures are presented in Supplementary Figure 4.

UPR is activated in both brain and lymph node metastases of NSCLC

To identify common molecular drivers of NSCLC metastasis independent of the metastatic site, we performed differential expression analysis comparing mBrain vs primary tumors (Figure 3A-C) and mLN vs primary tumors (Figure 3D-F), followed by functional enrichment analysis. Remarkably, despite the distinct cellular compositions and tumor subclusters observed in each metastatic site (Figure 2), both mBrain and mLN exhibited highly consistent enrichment of the UPR pathway. Specifically, genes and biological processes related to “cellular response to unfolded protein”, “cellular response to topologically incorrect protein”, and “chaperone-mediated protein folding” were among the upregulated categories in both metastatic sites compared with primary tumors (Figure 3B, C, E, and F). Differential expression and pathway analysis between metastatic lesions and primary lung cancer, including a volcano plot with key gene labels and a GSVA-based (v1.50.5) pathway comparison, are provided in Supplementary Figure 5A and B. These terms are core components of the UPR, which is activated upon endoplasmic reticulum stress to restore protein homeostasis[28,29]. To corroborate these scRNA-seq results, we applied pseudobulk differential expression analysis, which showed that the UPR was aberrantly activated in both brain and lymph node metastases of NSCLC (Supplementary Figure 5C and D). Analysis of an external dataset confirmed that UPR-related gene expression was elevated in metastatic tissues (Supplementary Figure 5E).

Figure 3 Upregulation of unfolded protein response and associated metabolic pathways in brain and lymph node metastases of non-small cell lung cancer. A: Volcano plot showing differentially expressed genes (DEGs) between brain metastases (mBrain) and primary lung tumors. Red dots (UP) indicate significantly upregulated genes, blue dots (DOWN) indicate downregulated genes, and gray dots (NOT) represent non-significant changes; B: Gene Ontology (GO) enrichment analysis of DEGs in mBrain compared to primary tumors; C: KEGG pathway enrichment analysis of DEGs in mBrain compared to primary tumors; D: Volcano plot showing DEGs between mLN and primary lung tumors. Red dots (UP) indicate significantly upregulated genes, blue dots (DOWN) indicate downregulated genes, and gray dots (NOT) represent non-significant changes; E: GO enrichment analysis of DEGs in mLN compared to primary tumors; F: KEGG pathway enrichment analysis of DEGs in mLN compared to primary tumors.

In addition to UPR activation, both types of metastases shared upregulation of several other pathways, including aerobic respiration, the electron transport chain, mitochondrial gene expression, and nucleoside triphosphate metabolic processes, indicating enhanced energy metabolism (Figure 3B-F). Notably, further analysis of additional gene sets revealed that glycolysis, mTOR signaling, fatty acid metabolism, and protein secretion pathways were also activated in both mBrain and mLN relative to primary tumors (Supplementary Figure 5). These pathways are known to support the high biosynthetic and energetic demands of metastatic tumor cells and may also contribute to the remodeling of the metastatic TME. Given the consistent enrichment of UPR-related genes in both metastatic sites, the UPR represents a biologically significant pathway meriting further investigation.

T cells in brain and lymph node metastases exhibit a profound exhausted phenotype

Given the marked increase in T cell abundance in both mBrain and mLN observed in Figure 1, we next characterized the functional states of tumor-infiltrating T lymphocytes. Unbiased subclustering of all T cells revealed multiple CD4⁺ and CD8⁺ T cell subsets, as well as their distribution across the different tissue types (Figure 4A and B). Strikingly, UMAP showed the distribution of all T cell subsets in lung cancer, mLN, nLN, mBrain, and nBrain, each defined by unique marker genes (Figure 4C-E). Compared with normal tissues, the proportion of CD4⁺/CD8⁺/γδT cells in primary and the two metastatic sites were altered (Figure 4F and G).

Figure 4 T cell landscape in lymph node and brain metastases of non-small cell lung cancer. A: UMAP visualization of 59450 T cells extracted from non-small cell lung cancer (NSCLC) primary tumors, normal lymph nodes (nLN), lymph node metastases (mLN), normal brain tissues, and brain metastases (mBrain). Major immune lineages are annotated, including CD4⁺ T cells, CD8⁺ T cells, and γδT cells; B: Distribution of CD4⁺ T cells, CD8⁺ T cells, and γδT cells across the five tissue groups; C: UMAP plot of T cell subclusters (generated by subclustering of all T cells), identifying multiple CD4⁺ and CD8⁺ subsets as well as γδT cells; D: UMAP plot showing the cell distribution of the 11 T cell subclusters across the five tissue groups; E: Dot plot showing the expression levels of selected canonical marker genes used for cell type identification of the 11 T cell subclusters; F: Stacked bar plot showing the proportion of each cell subset among the five tissue sample groups; G: Representative graph showing the proportional distribution of T cell subclusters across different tissue samples. Adjusted P values for between-group comparisons are shown at the top; H and I: Dot plot showing the elevated expression of exhaustion markers, including PD-1 (PDCD1), TIM-3 (HAVCR2), CTLA4, and LAG3, T cells of metastatic sites (mBrain and mLN) compared to normal tissues. nLN: Normal lymph nodes; mLN: Lymph node metastases; mBrain: Brain metastases; nBrain: Normal brain tissues.

Specifically, T cells in mLN and mBrain expressed significantly higher levels of exhaustion-associated inhibitory receptors, including PD-1 (PDCD1), TIM-3 (HAVCR2), CTLA4, and LAG3. In contrast, nLN and nBrain tissues exhibited low expression of these molecules. The elevated expression of multiple co-inhibitory receptors indicates that T cells infiltrating NSCLC metastases have entered a state of exhaustion, characterized by progressive loss of effector functions and impaired anti-tumor immunity. Importantly, this exhaustion phenotype was observed in two distinct metastatic organs (lymph node and brain) (Figure 4H and I), suggesting that a common immunosuppressive mechanism operates across different metastatic niches. The comprehensive subtyping, sample proportion distribution, and enrichment heatmap of T cells are shown in Supplementary Figure 6A-D. We then compared the expression of exhaustion-related molecules on T cells in another dataset, and the results were consistent with ours (Supplementary Figure 6E).

Taken together, these results demonstrate that NSCLC metastases establish a profoundly immunosuppressive T cell landscape that likely contributes to immune evasion and metastatic outgrowth.

Distinct myeloid profiles and immunosuppressive secretion in NSCLC metastatic sites

Myeloid cells are central regulators of the TIME and the cytokines that they secrete play a crucial role in tumor progression[30-32]. Similarly, we observed an increase in myeloid cell abundance in both mBrain and mLN as observed in Figure 1. To characterize their role in NSCLC metastasis, we performed subclustering analysis of all myeloid lineage cells, identifying major populations including macrophages, monocytes, dendritic cells (DCs), and neutrophils (Figure 5A and B). UMAP showed the distribution of myeloid cell subsets in lung cancer, mLN, nLN, mBrain, and nBrain, each defined by unique marker genes (Figure 5C-E). We next compared the composition of myeloid compartments in normal and metastatic tissues. nBrain was dominated by a single C8-SELENOP-Macro, whereas primary lung cancer exhibited diverse myeloid infiltration. Both mLN and mBrain displayed dramatic shifts in macrophage subpopulations, indicating remodeling of the myeloid landscape upon metastasis (Figure 5F and G).

Figure 5 Myeloid cell landscape in brain and lymph node metastases of non-small cell lung cancer. A: UMAP visualization of 15916 myeloid lineage cells, including macrophages, monocytes, dendritic cells (DC), and neutrophils, derived from non-small cell lung cancer (NSCLC) primary tumors, normal tissues, and metastatic lesions [lymph node metastases (mLN) and brain metastases (mBrain)]; B: UMAP visualization of macrophages, monocytes, DCs, and neutrophils across the five tissue groups; C: UMAP plot of myeloid cell subclusters (generated by subclustering of all myeloid cells), identifying multiple macrophage, monocyte, DC, and neutrophil subsets; D: UMAP plot showing the cell distribution of the 14 myeloid cell subclusters across the five tissue groups; E: Dot plot showing the expression levels of selected canonical marker genes used for cell type identification of the 14 myeloid cell subclusters; F: Stacked bar plot showing the proportion of each cell subset among the five tissue groups; G: Representative graph showing the proportional distribution of myeloid cell subclusters across different tissue samples. Adjusted P values for between-group comparisons are shown at the top; H-I: Dot plot showing the expression levels of key immunosuppressive in myeloid cells of metastatic sites (mBrain and mLN) compared to normal tissues. nLN: Normal lymph nodes; mLN: Lymph node metastases; mBrain: Brain metastases; nBrain: Normal brain tissues.

Next, we compared the expression levels of myeloid-derived inhibitory molecules between normal tissues and metastatic tumors. In mBrain, myeloid cells exhibited elevated expression of macrophage MIF and interleukin (IL)-10 (Figure 5H). MIF is a pro-inflammatory cytokine that also contributes to tumor progression and immune suppression; IL-10 is a well-established immunosuppressive and tumor-promoting cytokine that inhibits effector T cell function and promotes a protumorigenic environment. In mLN, myeloid cells exhibited elevated expression of transforming growth factor-beta (TGF-β), MIF, and IL-10 (Figure 5I). TGF-β and IL-10 are potent immunosuppressive cytokines. Within the TME, TGF-β, MIF, and IL-10 are predominantly secreted by immunosuppressive M2-polarized macrophages[33]. The combination of these factors creates a highly suppressive niche that effectively impairs anti-tumor immunity. The comprehensive subtyping, sample proportion distribution, and enrichment heatmap of myeloid cells are shown in Supplementary Figure 7A-C. We integrated the newly downloaded dataset GSE293360 with our original dataset GSE131907 for joint analysis. Comparison of immune-inhibitory molecules secreted by myeloid cells revealed consistent expression patterns across both datasets, corroborating our primary findings (Supplementary Figure 7D).

Taken together, myeloid cells derived from brain and lymph node metastases of NSCLC exhibit distinct molecular profiles accompanied by elevated secretion of immunosuppressive factors.

Prognostic value of tumor subpopulations and mIHC validation of key cellular and molecular features in the TME

To evaluate the clinical relevance of the cell subpopulations identified in our single-cell analysis, we queried bulk transcriptomic data from TCGA LUAD cohorts. In univariate Cox analysis, we found that several tumor epithelial subclusters were associated with prognosis; for example, C14_PBK_Tumor was associated with poor survival (HR = 1.68, 95%CI: 1.25-2.26, P < 0.001) (Supplementary Figure 8A and C-E). To evaluate the independent prognostic value of the identified tumor subpopulations, we performed multivariate Cox regression analysis incorporating the relative proportions of each tumor cluster, along with clinical covariates including gender, tumor stage, and age. Advanced tumor stage (III/IV) was strongly associated with poor overall survival (HR = 2.26, 95%CI: 1.65-3.10, P < 0.001), confirming its role as a well-established prognostic factor. Notably, among the various tumor subpopulations, C14_PBK_Tumor emerged as a significant independent adverse prognostic factor (HR = 2.39, 95%CI: 1.15-4.97, P = 0.020), indicating that patients with higher proportions of this specific cell cluster exhibited a substantially increased risk of death. No other tumor clusters, including C3_C1orf105_Tumor, C5_FABP7_Tumor, C10_FAM83D_Tumor, or C13_C6orf222_Tumor, reached statistical significance in this model. Likewise, gender and age did not show independent prognostic value (Figure 6A). We verified the impact of the C14_PBK_Tumor subtype on prognosis using another LUAD dataset, GSE72094. The results also showed that patients with tumors enriched for the C14_PBK_Tumor subtype had a worse prognosis (Supplementary Figure 8B). These findings suggest that the presence of C14_PBK_Tumor is associated with the prognosis of lung cancer, warranting further investigation into its biological role and implications.

Figure 6 Multiplex immunohistochemistry validation and prognostic value of tumor subpopulations. A: Multivariate Cox regression analysis of prognostic factors. Forest plot showing hazard ratio and P values for overall survival of lung adenocarcinoma patients from The Cancer Genome Atlas; B: Representative multiplex immunohistochemistry images showing the expression and spatial distribution of CD68⁺/CD163⁺, PD-1⁺/ CD8⁺, XBP1⁺/ATF6⁺/ATF4⁺ makers in Mets0 and Mets1 non-small cell lung cancer tissues. Scale bars: 50 μm; C: Comparison of positive cell density between Mets0 and Mets1 cancer regions; D: Comparison of positive cell rates between Mets0 and Mets1 cancer regions. Mets0 (n = 11) and Mets1 (n = 12); statistical significance was determined using Mann-Whitney U test; NS indicates not significant.

We performed mIHC on tissue sections from non-metastatic (Mets0) and metastatic (Mets1) NSCLC lesions, staining for markers of macrophages (CD68⁺/CD163⁺), T cell exhaustion (PD-1⁺/CD8⁺), and UPR activation (XBP1⁺/ATF6⁺/ ATF4⁺). Representative images are shown in Figure 6B. Quantification of positive cell rates revealed that, in Mets1 patients, the proportion of CD68⁺/CD163⁺ macrophages were higher than in Mets0 patients. Furthermore, the proportion of XBP1⁺/ATF6⁺/ATF4⁺ positive cells was higher in Mets1 patients (Figure 6C and D). Representative images of mIHC for multiple samples are shown in Supplementary Figure 9.

Combined with the enrichment of UPR-related gene signatures from scRNA-seq, this points to a role for the UPR in NSCLC brain and lymph node metastases that merits further investigation, although this requires subsequent large-scale validation.

DISCUSSION

In this study, we constructed a comprehensive scRNA-seq atlas of NSCLC primary tumors, mLN, and mBrain. Our main observations are threefold: (1) The tumor epithelial subpopulations are distributed differently in brain and lymph node metastatic lesions; (2) The UPR is abnormally activated across both metastatic sites; and (3) A dual immunosuppressive microenvironment. TCGA data analysis revealed that C14_PBK_Tumor was associated with worse overall survival. Notably, mIHC revealed that, in Mets1 patients, the proportion of CD68⁺/CD163⁺ macrophages and XBP1⁺/ATF6⁺/ATF4⁺ cells were higher than that in Mets0 patients.

Our observation of tumor epithelial subpopulations in the metastatic organs aligns with recent scRNA-seq studies. Tagore et al[23] reported a CIN-high neuronal-like subpopulation in NSCLC mBrain. We found that the C1-CDH2 and C5-FABP7 subgroups were enriched in the mBrain. However, our study additionally identified tumor epithelial subclusters within lymph node metastases (C2-PSCA, C4-SPP1, and C7-TRPM8) that were not characterized in their work. Zhang et al[24] described MOCCI-driven OXPHOS reprogramming in mLN. We found that the C7-TRPM8 subcluster shares a metabolic signature but also expresses distinct immune-related genes, suggesting further heterogeneity.

Bai et al[34] reported that brain metastases exhibit profound immunosuppressive reprogramming, with enrichment of HSP70-high stress-responsive T cells and PLTP⁺ TAMs. Our analysis revealed the presence of a dual immunosuppressive microenvironment in both brain metastases and lymph node metastases of NSCLC. Regarding T cell exhaustion, our finding of PD-1/TIM-3/CTLA4/LAG3 co-expression in both mBrain and mLN is consistent with previous reports[17,18]. Our scRNA-seq data uncovered the presence of exhausted CD8⁺ T cell subsets in the TIME of metastatic NSCLC. Although there was no significant difference in CD8⁺/PD-1⁺ expression among different groups as detected by mIHC, the results of scRNA-seq and mIHC were not completely consistent, likely owing to differences in the markers for progenitor-exhausted vs terminally-exhausted T cells. PD-1 expression is high in terminally-exhausted T cells but low in progenitor-exhausted T cells[35,36]. Therefore, CD8⁺/PD-1⁺ dual staining fails to sensitively identify early exhausted T cells. Zhou et al[37] reported SPP1⁺ macrophages in leptomeningeal metastases, while we found that in mBrain, MIF and IL-10 were upregulated, whereas in mLN, MIF, IL-10, and TGF-β all showed increased expression. These factors can be secreted by macrophages. The mIHC data further reveal that the CD68⁺CD163⁺ positive cell rate is higher in NSCLC metastatic lesions. Given that CD68 and CD163 are characteristic markers of macrophages and that previous studies have reported that myeloid cells in the metastatic niche tend to differentiate toward a macrophage phenotype, this finding suggests that macrophage infiltration is increased in NSCLC metastatic lesions.

While our current analysis primarily focused on T-cell exhaustion and macrophage-associated immunosuppression, the metastatic immune microenvironment encompasses a broader array of immune lineages whose functional states require further characterization. For example, Tregs are potent suppressors of antitumor immunity and have been shown to accumulate in metastatic niches, correlating with poor prognosis in NSCLC. Similarly, DCs, neutrophils, and natural killer (NK) cells all contribute critically to antitumor immune responses. Moreover, ligand–receptor interactions between tumor cells and infiltrating immune cells are essential for deciphering the molecular crosstalk that drives immune evasion and metastatic progression. Although our current dataset and analytical scope do not permit a comprehensive dissection of all these components, this is an important direction for our future research.

To our knowledge, our study offers a panoramic analysis of two distinct metastatic sites within the same disease context. We found that the UPR was abnormally activated in both brain metastases and lymph node metastases of lung cancer. While UPR activation has been implicated in individual cancer types[38,39], its consistent upregulation across two distinct metastatic organs in the same disease has not been shown. To corroborate these scRNA-seq results, we conducted pseudobulk differential expression analysis and observed a high degree of consistency between the two approaches. The mIHC results showed that the XBP1⁺/ATF6⁺/ATF4⁺ cells were higher in the Mets1 group, further supporting the biological significance of this pathway. This observation alone does not demonstrate a causal role in metastasis. UPR-related genes were consistently enriched in both metastatic sites; however, this association may be influenced by tissue and cohort composition. Despite this, we propose UPR activation as a candidate pathway that merits functional testing. Finally, in multivariate Cox regression analysis, the C14_PBK_Tumor subcluster was identified as an independent predictor of poor overall survival in LUAD patients. Nevertheless, further validation in larger independent patient cohorts is warranted to assess its translational potential and clinical value.

Several limitations should be acknowledged. First, the scRNA-seq data were derived from public cohorts with limited clinical annotations; prospective validation with larger sample sizes is needed. Second, UPR-related genes were consistently enriched in both metastatic sites; this association may be influenced by tissue and cohort composition. Third, our mIHC validation, while supportive, was performed on a small sample cohort and should be considered preliminary. Fourth, several confounding factors need to be considered. Different treatment exposure statuses may affect gene expression and immune infiltration patterns. Although all samples were adenocarcinomas, there was inconsistency in the histological subtypes. Additionally, there were differences in tumor purity among the samples, which may affect the detection of intrinsic signals of tumor cells. Furthermore, the presence of brain-resident microglia in brain metastases complicates the interpretation of the myeloid cell population. Finally, although we integrated multiple datasets and experimental platforms to enhance the robustness of our findings, cross-study batch effects remain a non-negligible concern. In addition, incomplete clinical matching between the scRNA-seq public cohorts and our mIHC cohort limits the generalizability of our integrated observations. Despite these limitations, our findings establish a foundation for future large-scale validation studies.

Looking forward, we propose a validation pathway and work plan. First, our key findings should be validated in independent, clinically annotated cohorts of mBrain and mLN from treatment-naïve patients, with sufficient sample sizes for histology-stratified analyses. Second, spatial validation should be performed via multiplex immunofluorescence or spatial transcriptomics on tissue sections. Third, in vivo and in vitro experiments should be performed to verify the function of the UPR pathway.

In summary, we performed single-cell transcriptomic analysis of brain and lymph node metastases in NSCLC. We further explored the prognostic relevance of tumor subclusters using TCGA data and performed preliminary mIHC validation of macrophage infiltration and UPR-related protein expression, demonstrating the feasibility of detecting these markers at the protein level. From a translational perspective, as a pivotal stress-responsive pathway, the UPR merits in-depth exploration in the context of tumor progression.

CONCLUSION

Our study integrates multiple public scRNA-seq datasets to construct a transcriptomic atlas of NSCLC primary tumors and metastatic sites, identifying tumor epithelial subpopulations in different metastatic organs, common activation of the UPR pathway, and immunosuppressive remodeling characterized by exhausted T cells and macrophage enrichment. Among the tumor epithelial subpopulations, the C14_PBK_Tumor subtype was associated with worse prognosis in LUAD patients. However, further functional and prospective studies are needed to confirm these observations.