We integrated publicly available single-cell and bulk RNA sequencing datasets. A total of seven single-cell RNA-seq (scRNA-seq) datasets were obtained from the Gene Expression, including: E-MTAB-13526, GSE164789, GSE207422 (Hu et al., 2023), GSE117570 (Wang et al., 2021), GSE241934 (Chen M. et al., 2025), GSE131907 (Kim et al., 2020) and GSE148071 (Wu et al., 2021). All datasets were downloaded from the GEO repository and preprocessed using standardized workflows. Three bulk RNA-seq datasets were incorporated for integrative analysis. The TCGA cohort (n = 133), with gene expression profiles and clinical annotations obtained from the National Cancer Institute’s Genomic Data Commons (GDC) Data Portal. Data regarding clinical outcomes and survival probabilities for the GSE115002 (Liao et al., 2025) and GSE13213 (Tomida et al., 2009) cohorts were retrieved from the GEO database to ensure the consistency and comparability of clinical variables. Data on immune checkpoint inhibitors in lung cancer were obtained from su2c cohort (Ravi et al., 2023). All dataset information were collected in Supplementary Table S1.

For RNA-seq datasets, gene expression values were transformed to log2 (TPM + 1), whereas for microarray datasets, normalized expression matrices provided by the original studies were used. The batch effects across datasets were eliminated by means of the implementation of the Combat algorithm from the SVA R package (Zhang et al., 2020). Following correction, the three bulk cohorts were merged into a unified meta-cohort for downstream integrative analyses and feature extraction.

For single-cell transcriptomic analysis, initial quality control was performed in accordance with the following criteria: (1) removal of ribosomal and kinase genes; (2) mitochondrial gene content <15%. The detail information of scRNA-seq data were collected in Supplementary Table S2. The raw expression matrix was subsequently normalized to eliminate technical bias. To address batch effects arising from multi-sample integration, Harmony algorithm was applied to the principal component analysis (PCA) results for batch correction. The top 30 principal components were selected for downstream analysis. Unsupervised clustering was conducted to identify cell populations, and UMAP was used for dimensionality reduction and visualization. To ensure the accuracy of cell-type annotation, clusters were annotated at a resolution of 0.6 using lineage-specific marker genes. The annotation of each cluster to specific cell types was based on established canonical markers. Myeloid subpopulations were then subjected to further annotation, after which differential gene expression analysis was performed using the COSG algorithm. Pathway enrichment analysis was conducted for each subpopulation using the KEGG database via the clusterProfiler R package (v4.12.1), with Benjamini–Hochberg false discovery rate correction.

Samples from the TCGA cohort were stratified based on pathological lymph node status (PN), with N0M0 samples designated as the control group and N > 0M0 samples designated as the treatment group. The differential expression analysis between the treatment and control groups was performed using the limma package (v3.56.2). The ranking of genes was based on adjusted p-values and log₂ fold changes. The top 100 up- and downregulated genes were selected to construct up- and downregulated gene sets, respectively. Subsequently, gene set enrichment scores were calculated in single-cell RNA sequencing data using the irGSEA framework, based on these two gene sets. Based on the up- and downregulated gene sets, cells with enrichment scores in the top 10% were defined as LNM cells and non-LNM cells, respectively.

To evaluate the functional contributions of distinct cell types within bulk RNA-seq datasets, deconvolution analysis was performed on nine annotated myeloid subpopulations using InstaPrism (version 0.1.6) to assess their infiltration levels. Based on the estimated infiltration abundance of each cell type, samples were stratified for prognostic evaluation, and survival analyses were conducted within the bulk cohorts. Specifically, infiltration scores were compared between different nodal statuses (N > 0 vs. N = 0) and clinical stages (Stage I-II vs. Stage III) using the Wilcoxon rank-sum test.

Regulatory network inference and transcription factor (TF) activity analysis were conducted using the pySCENIC framework (v0.11.2) (Van de Sande et al., 2020) with default parameters. Initially, GRNBoost2 was applied to construct gene regulatory networks by identifying co-expression modules and putative TF-target gene relationships. RcisTarget was then used to perform motif enrichment analysis to confirm TF binding site enrichment. Regulon activity scores were subsequently computed using the AUCell algorithm. For downstream analysis, only the top ten most specific regulons were retained for each cell subtype. To identify RELB target genes, the top 50 genes ranked by importance score were selected. The irGSEA was performed at the single-cell level to evaluate enrichment patterns across individual cells, with these RELB target genes being the basis for the study.

The limma package (Ritchie et al., 2015) was used to identify genes differentially expressed between the high- and low-RELB score groups. KEGG pathway enrichment analysis was conducted via GSEA using the gseKEGG function in clusterProfiler (v4.12.1) (Demirel et al., 2020), with multiple testing correction performed using the Benjamini–Hochberg false discovery rate (FDR) method.

Human neutrophil cells (Cat# JN-CC6058) was obtained from Jining Bio Co., LTD. The human non-small cell lung cancer (NSCLC) cell lines A549 (Cat# CBP60053) and H1299 (Cat# CBP60053) were also purchased from Cobioer Biosciences Co., LTD. Upon thawing, all cell lines were centrifuged at 1,000 rpm for 5 min to remove the upper layer, followed by resuspension in complete DMEM medium supplemented with 10% fetal bovine serum (FBS) and 1% penicillin–streptomycin. Cells were maintained in a humidified incubator at 37 °C with 5% CO₂. The culture medium was replaced every 2 days, and subculturing was performed every 3 days.

For gene silencing experiments, neutrophils cells were transfected with two independent small interfering RNAs targeting RELB (siRELB#1 and siRELB#2), as well as a negative control siRNA (siCtrl), using Lipofectamine™ RNAiMAX transfection reagent according to the manufacturer’s instructions. Cells were harvested 24 h post-transfection for downstream analyses. The siRNA sequence, used were as follows: siCtrl (UUC​UCC​GAA​CGU​GUC​ACG​U), siRELB#1 (GAAGAAGGAAAUUGAAGCU), siRELB#2 (GCU​ACG​GUG​UGG​ACA​AGA​A).

Total RNA was extracted from cultured cells using TRIzol reagent according to the manufacturer’s protocol. RNA concentration and purity were assessed spectrophotometrically. Reverse transcription was performed using the ReverTra Ace qPCR RT Kit, following standard procedures. Quantitative real-time PCR (qPCR) was conducted using the 7500 Fast DX Real-Time PCR System with gene-specific primers. The primer sequences were in Supplementary Table S3. Gene expression levels were quantified using the 2^−ΔΔCt method, with β-actin serving as the internal control. All qPCR reactions were performed in triplicate, and expression differences were statistically analyzed.

Human neutrophil cells were seeded into 96-well plates. Cell viability was assessed at 0, 24, 48, and 72 h post-transfection using the Cell Counting Kit-8 (CCK-8). At each time point, 10 μL of CCK-8 solution was added to each well, followed by incubation at 37 °C for 2 h. Absorbance was measured at 450 nm using a microplate reader to evaluate cell viability. All measurements were performed in triplicate.

Human neutrophil cells transfected with siRELB#1 or siRELB#2 were co-cultured with A549 or H1299 cells for 14 days, or until most colonies contained more than 50 cells. During the incubation period, the culture medium was replaced every 3 days. At the end of the experiment, cells were washed once with phosphate-buffered saline (PBS), fixed with 4% paraformaldehyde, and stained with crystal violet. Colony images were captured using a digital camera, and colony formation was quantified for analysis.

A549 and H1299 cells were enzymatically dissociated, rinsed, and resuspended in serum-free medium at 4 × 10⁶ cells/mL. Subsequently, 500 μL of human neutrophils were placed in the lower chambers of 24-well plates, upon which Matrigel-coated Transwell inserts (8.0-μm pores) were mounted to facilitate direct contact between the 2 cell populations. Subsequently, 200 μL of A549 or H1299 cell suspension (8 × 10⁵ cells) was added to the upper chamber. After 24 h of incubation, non-migrated and non-invaded cells on the upper surface of the membrane were carefully removed. Cells on the lower surface were fixed with methanol, stained with crystal violet, and quantified under a light microscope at 100× magnification. Three representative fields were selected per insert for statistical analysis. Graphs were generated and statistical significance was determined using GraphPad Prism.

All statistical analyses were performed using R software (v 4.3.1). Survival analyses were conducted using the survival (v3.5–7) and survminer (v0.4.9) packages. Kaplan–Meier survival curves were generated to compare survival outcomes, and Cox proportional hazards regression models were used to estimate hazard ratios. Statistical significance was assessed using the Wilcoxon rank-sum test, with multiple comparisons adjusted via the Holm method. Differences in qPCR expression levels and cell counts were evaluated using unpaired t-tests. The infiltration scores difference between different nodal statuses (N > 0 vs. N = 0) and clinical stages (Stage I-II vs. Stage III) were accessed using the Wilcoxon rank-sum test. Drug-target pathway enrichment was assessed using hypergeometric testing, with p-values corrected for multiple comparisons using the Benjamini–Hochberg (BH) method.

To decipher the cellular basis of lymph node metastasis (LNM) in EGFR-wildtype lung adenocarcinoma (LUAD), we integrated multiple single-cell RNA sequencing datasets from this patient subset. The batch correction was performed using Harmony algorithm (Supplementary Figure S1A). We employed the UCell algorithm to project the LNM status, defined from bulk RNA-seq data of the TCGA-LUAD (EGFR-wildtype) cohort, onto the integrated single-cell transcriptomes. This approach successfully defined LNM, non-LNM, and NS cells (Figure 1A). The UCell scores for LNM-upregulated and LNM-downregulated gene sets were visualized on UMAP plots for each cell (Supplementary Figure S1B,D) and the distribution of these scores across cell types were displayed using violin plots (Supplementary Figure S1C,E). Notably, the myeloid subpopulation exhibited the highest enrichment for the LNM-upregulated gene set (Supplementary Figure S1C). UMAP visualization revealed the comprehensive cellular composition of the EGFR-wildtype tumor microenvironment (Figure 1B), validated by the expression of canonical cell-type marker genes (Figure 1C). Proportional analysis identified myeloid cells as significantly enriched in the LNM group, a shift illustrated by sankey and bar plots (Figures 1D,E). This indicates a specific association between myeloid cell abundance and the LNM status in this subtype. Deconvolution analysis of bulk TCGA data further confirmed a significantly higher infiltration score for myeloid cells in the LNM group (Figure 1F).

Additionally, the analysis of distribution distance on the UMAP space suggested that myeloid cells stayed in a relatively stable state between LNM and non-LNM conditions (Figure 1G). This suggests that the observed enrichment is not due to major transcriptional state shifts but rather an expansion of a conserved myeloid population. Collectively, these outcomes designate myeloid cells as a significant cell population related to LNM in EGFR-wildtype LUAD, encouraging their examination in upcoming analyses.

By performing sub-clustering, we aimed to analyze the heterogeneity of myeloid cells in the EGFR-wildtype context and identify specific subsets that drive LNM, leading to the identification of nine different myeloid subpopulations (Figure 2A). Proportional analysis indicated a distinct increase in the ELANE + Neu subpopulation within the LNM group, whereas other subsets like FCGR3B + Neu and LAD1+DC exhibited relative decreases (Figures 2B,C). The data suggests a targeted growth of ELANE + Neu within metastatic environments. Pseudo-bulk analysis of differentially expressed genes in LNM myeloid cells underscored their participation in inflammatory and immunosuppressive pathways (Figure 2D). This functional bias indicates a role in establishing a supportive microenvironment for metastasis. The analysis of individual samples showed that ELANE + Neu was significantly higher in LNM samples, with notable alterations in other subpopulations, including an increase in SPP1+TAM and OLFML3+TAM and a decrease in CD1C + DC, FCGR3B_Neu, and LAD1+DC (Figure 2E). This confirms the sample-level consistency of the ELANE + Neu expansion. Pathway analyses, including KEGG and Hallmark, provided additional insights into the specific functional roles of each myeloid subpopulation (Figures 2F,G). Specifically, ELANE + Neu was notably enriched in oxidative phosphorylation, phagocytosis, and ferroptosis, indicating a metabolic and inflammatory profile that supports metastasis. This unique functional profile potentially underlies its pro-metastatic capacity. To evaluate the clinical relevance of these subpopulations, we deconvoluted their infiltration scores into multiple independent EGFR-wildtype cohorts and performed survival analysis. A high ELANE + Neu infiltration score was significantly associated with shorter overall survival (Figure 2H). This establishes a direct link between ELANE + Neu abundance and poor patient outcome. Clinicopathological correlation further confirmed that this score was elevated in patients with advanced N stage and higher clinical stage (Figure 2I; Supplementary Figure S2).

Together, these results identify ELANE + Neu not merely as an associated cell type but as a pivotal myeloid subpopulation driving metastasis and poor prognosis in EGFR-wildtype LUAD. They also highlight the unique functional and metabolic profile of this subset within the tumor microenvironment.

To elucidate the intrinsic regulatory mechanism governing the pro-metastatic function of the ELANE + Neu subpopulation, we performed transcription factor (TF) regulon analysis using pySCENIC. The results identified RELB as a top TF specifically enriched not only in the LNM group (Figure 3A) but also within the ELANE + Neu subpopulation itself (Figure 3B). This dual enrichment strongly suggests RELB is a key regulator specific to this pro-metastatic context and cell type. The activity of the RELB (+) regulon was significantly elevated in ELANE + Neu cells (Figure 3C) and across the LNM myeloid compartment (Figure 3D), with its activity score significantly higher in the LNM group (Figure 3E). This confirms heightened RELB-driven transcriptional activity in metastasis-associated myeloid cells. Consistent with its established role in promoting metastasis and immunosuppression in other cancer types such as hepatocellular carcinoma and diffuse large B-cell lymphoma (Tan et al., 2015; Nuan-Aliman et al., 2022), our findings highlight RELB as a key transcriptional regulator in EGFR-wildtype LUAD. Enrichment analysis of RELB target genes revealed significant involvement in NF-κB signaling, inflammatory response, and immune regulatory pathways (Figure 3F). This indicates that RELB exerts its function through these well-known pro-tumorigenic pathways.

To functionally validate the role of RELB, we knocked down RELB in neutrophils and co-cultured them with two tumor cell lines (A549 and H1299). This caused a reduction in the expression of NF-κB-associated and inflammatory genes IDO-1 and IL-6 within the tumor cells (Figure 3G). This suggests that RELB within neutrophils may modulate the tumor cell phenotype to be more inflammatory.

Furthermore, silencing RELB led to a significant reduction in tumor cell proliferation (Figure 3H), colony formation (Figure 3I), migration (Figure 3J), and invasion (Figure 3K) in both A549 and H1299 cell lines. These functional validations suggest that RELB activity in neutrophils may play an important role in enhancing the aggressive behaviors of tumor cells.

Collectively, these results suggest that RELB may function as a key transcription factor influencing the pro-metastatic behavior of the ELANE + Neu subpopulation, in part through activation of NF-κB and inflammatory pathways. This regulatory program is likely to contribute to lymph node metastasis in EGFR-wildtype LUAD.

To evaluate the transcriptional output of the RELB pathway and its clinical relevance beyond the single-cell level, we constructed a RELB signature score based on its target gene set using the ssGSEA method and applied it to bulk transcriptomic data. In three independent EGFR-wildtype LUAD cohorts, a high RELB score consistently predicted shorter overall survival (Figure 4A). The score was significantly elevated in tumors with lymph node metastasis and advanced clinical stage (Figure 4B; Supplementary Figure S3), underscoring its close association with disease progression. We then stratified the TCGA cohort into high- and low-RELB score groups and performed differential expression analysis. The volcano plot illustrated considerable transcriptomic alterations in the high-score group as opposed to the low-score group, with 2011 genes upregulated and 3620 downregulated (Figure 4C), reflecting the broad regulatory impact of RELB on downstream gene expression.

The extensive influence emphasizes RELB as a primary regulator of tumor traits. GSEA of the upregulated genes showed significant activation across various pathways. The extensive influence emphasizes RELB as a primary regulator of tumor phenotype. Gene set enrichment analysis (GSEA) of the genes that were upregulated identified notable activation of various pathways, including the IL-17 signaling pathway, TNF signaling pathway, PI3K-AKT signaling pathway, pathways in Cancer, MAPK signaling pathway, and NF-κB signaling pathway (Figure 4D), indicating that RELB plays an important role in promoting inflammatory responses and cell survival pathways. The observed pro-tumorigenic functions align with this pathway activation profile. GSVA additionally confirmed the activation of hypoxia, EMT and other malignancy-related hallmark pathways in tumors with high scores (Figure 4E). Together, these analyses establish that high RELB activity is associated with a more aggressive tumor phenotype characterized by specific transcriptional reprogramming. Immune infiltration analysis showed that a high RELB score was correlated with features of an immunosuppressive microenvironment, characterized by enhanced myeloid-mediated immunosuppression, enrichment of pro-tumor cytokines, increased infiltration of cancer-associated fibroblasts, and reduced enrichment of immune-activation features such as anti-tumor cytokines and effector cells (Figure 4F). This suggests that the RELB signature captures an immunosuppressive TME state, potentially explaining its association with poor prognosis and metastasis.

In conclusion, this study demonstrates that the RELB signature score may serves as a prognostic biomarker and delineates a malignant tumor phenotype characterized by transcriptomic reprogramming, enhanced invasiveness, and immunosuppression in EGFR-wildtype LUAD.

We next assessed the translational potential of the RELB signature in predicting therapy response. Distinct genetic modifications were identified in the EGFR-wildtype TCGA cohort’s mutational landscape analysis (Figures 5A,B), along with a significantly lower tumor mutational burden in the high RELB score group (Figure 5C). This inverse correlation indicates distinct biological differences between tumors with high and low RELB levels. By integrating the RELB score with TMB, patients were effectively categorized, revealing a subgroup with the worst prognosis (Figure 5D). This demonstrates the complementary prognostic value of RELB score and TMB. In the SU2C immunotherapy cohort (predominantly non-EGFR mutant), a low RELB score was associated with better therapeutic response, as confirmed by ROC analysis (Figures 5E,F), and longer overall survival (Figure 5G). This provides direct evidence that the RELB signature can predict inferior response to immune checkpoint inhibitors.

A drug sensitivity analysis was performed using the calcPhenotype function from the oncoPredict R package. It was found that there were several targeted anti-tumour agents that may be effective in patients with low RELB scores, such as NU7441, Doramapimod, AZD1208, etc. (Figure 5H). Patients with elevated RELB scores were sensitized to several chemotherapeutic agents (Docetaxel, 5-Fluorouracil, Irinotecan, Topotecan, etc.) and targeted anti-tumor agents (Crizotinib, Luminespib, AZD1332, MIM1, etc.) (Figure 5H). The differential pattern of drug sensitivity suggests potential treatment strategies that can be tailored according to the RELB activity status. Notably, several ADC drug targets were expressed at significantly higher levels in the high RELB score group (Figure 5I), suggesting a potential therapeutic vulnerability for this aggressive tumor subtype. This discovery suggests that ADCs could be a promising treatment option for tumors with high RELB expression that resist immunotherapy.

In conclusion, these findings indicate that the RELB signature may be a strong biomarker for forecasting prognosis, response to immunotherapy, and unique drug sensitivity patterns, emphasizing its potential in guiding personalized treatment strategies in EGFR-wildtype LUAD.