Single-cell RNA-seq data (GSE127465) from GEO database (7 LUAD samples, 31,179 cells) were analyzed alongside TCGA data (515 tumors, 59 normals). External validation employed GEO datasets GSE13213 and GSE11969.

Single-cell data were processed using MAESTRO v1.5.1, including quality control, batch correction, clustering, and annotation (23). Cells with <1,000 reads or <500 genes were excluded. PCA (top 2,000 genes) and clustering (KNN/Louvain; 30 PCs, resolution=1) were performed, followed by t-SNE visualization (25 clusters). DEGs were identified (Wilcoxon; |logFC|≥0.25, FDR<1×10^-5) and annotated using published markers.

TCGA LUAD samples were stratified into lymph node-positive (N1+N2+N3) and -negative (N0) groups. DEGs were analyzed using “limma”, with comparative assessment of MUC5B versus other MUC family genes. Survival analysis was performed after dichotomizing patients by median MUC5B expression.

We performed WGCNA on lung adenocarcinoma samples to identify modules of genes associated with lymph node metastasis (24). The expression matrix was filtered using “goodSamplesGenes” (soft threshold=10), identifying 12 gene modules. Genes with module membership>0.7 and gene significance>0.1 were retained for lymph node metastasis analysis. Pearson correlation identified candidate genes associated with both MUC5B and metastasis.

A PPI network was constructed using GeneMANIA (25), with hub genes identified (degree>30) and further screened by MNC/Betweenness centrality. NetworkAnalyst analyzed miRNA-MUC5B (TarBase v8.0) and TF-MUC5B (JASPAR) interactions, integrated into a miRNA-MUC5B-TF network via Cytoscape. Functional enrichment was performed using GO/KEGG.

Core genes were analyzed using LASSO (“glmnet”) and SVM-RFE (“e1071”) (20) for feature selection. SHAP visualized model interpretations, while “pROC” calculated AUC values. Nomograms (“rms”) and decision curves (“rmda”) were generated for clinical utility assessment.

Immune-related pathways were identified through literature review. “GSVA” calculated immune scores, with significant changes assessed (Pearson’s test, P<0.05). Correlations between immune scores/diagnostic genes and core genes/immune infiltration were analyzed. For immune pathway correlation analyses, statistical significance was determined using the p.adjust function in R to correct for multiple comparisons (26).

This retrospective study analyzed 65 LUAD patients diagnosed at the Cancer Hospital, Chinese Academy of Medical Sciences (January 2017-December 2018). Inclusion criteria: 1) surgical treatment; 2) pathological LUAD confirmation; 3) complete clinical/follow-up data. Exclusion criteria: 1) other malignancies; 2) non-LUAD pathology; 3) lost to follow-up. The study was approved by the Institutional Ethics Committee (NCC2019C-167). All patients underwent 5-year postoperative follow-up (or until death) with informed consent.

For IHC analysis, LUAD specimens were deparaffinized in xylene and rehydrated through graded alcohols. After endogenous peroxidase blocking (3% H₂O₂) and antigen retrieval (citrate buffer), sections were incubated with anti-MUC5B antibody (1:100, HUABIO) at 4°C overnight, followed by HRP-conjugated secondary antibody (30 min, RT). DAB (Servicebio) was used for visualization. Staining was quantified using the H-score formula: H-Score = Σ (pi × i) = (% weak cells × 1) + (% moderate cells × 2) + (% strong cells × 3), where i represents staining intensity (0=negative, 1=weak, 2=moderate, 3=strong) and pi is the percentage of cells at each intensity level (27, 28). The IHC scoring was independently performed by two senior pathologists with extensive experience in pathological diagnosis of thoracic tumors.

H1975 and A549 cells were maintained in RPMI-1640 and DMEM (10% FBS) at 37°C/5% CO₂, respectively. Cells were transfected with either negative control (NC) or MUC5B-targeting siRNA (siMUC5B; GenePharma) using the following sequences: Sense: 5′-GCAGCUACGUUCUGUCCAATT-3′; Antisense: 5′-UUGGACAGAACGUAGCUGCTT-3′. Transfected cells were cultured under standard conditions until analysis.

Total RNA was extracted from tissues or cells with TRIzol (Invitrogen, USA), according to the manufacturer’s instructions. Total RNA was reverse transcribed using a reverse transcription kit (Vazyme Biotech Co., Ltd, China), according to the manufacturer’s protocol. β-actin was used to normalize the RNA levels using the 2-ΔΔCt method. The sequences of the primers for qRT-PCR were: MUC5B-F: 5′- CCCGTGTTGTCATCAAGGC -3′; MUC5B-R: 5′- CAGGTCTGGTTGGCGTATTTG -3′; β-actin-F: 5′- TGGCACCCAGCACAATGAA-3′; β-actin-R: 5′- CTAAGTCATAGTCCGCCTAGAAGCA-3′; GINS1-F: 5′-ACGAGGATGGACTCAGACAAG-3′; GINS1-R: 5′-TGCAGCGTCGATTTCTTAACA-3′; GINS2-F: 5′-CCCTGGTTTACCCGTGGAAG-3′; GINS2-R: 5′-GGGAGCAGGCGACATTTCT-3′; GINS3-F: 5′-ACTTTTATCGGACGTTTTCGCC-3′; GINS3-R: 5′-TCTCCATCTCGTCTAGCCTGG-3′; GINS4-F: 5′-AGTTGGCCTTTGCCAGAGAG-3′; GINS4-R: 5′-GAACTGCCCGAAAGAGGTCC-3′.

Cell growth was measured using a CCK-8 kit (Dojindo, Kumamoto, Japan) according to the instructions. Briefly, approximately 1×10³ cells per well were plated into 96-well plates and cultured in the indicated medium. Cell proliferation, measured at 450 nm, was examined every day for four days according to the manufacturer’s protocol according to the manufacturer’s protocol.

1×10³ cells per well were seeded in 6-well plates, and the culture was terminated when the cells formed obvious clones under the microscope. The colonies were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet for 30 min at room temperature. Finally, the colony number was quantified.

H1975 and A549 cells were seeded into 6-well plates. Once the cells reached full confluence, a scratch was created using a 200 µl pipette tip. The dislodged cells were removed with serum-free medium. Images of the wounded area were captured at 0 hours and 24 hours. The cell migration rate was calculated using the following formula: %migration rate = ((0 h wound area - 24 h wound area) × 100)/0 h wound area.

Invasion assay:

H1975 (2.5×10⁴) and A549 (5×10⁴) cells in serum-free medium were seeded into Matrigel-coated Transwell chambers (8 µm; Corning). Lower chambers contained 10% FBS medium. After 48h incubation (37°C/5% CO₂), non-invading cells were removed, membranes were fixed (4% PFA) and stained (0.1% crystal violet, 30min). Invading cells were counted from three random fields

Migration assay:

Performed similarly without Matrigel coating.

RNA sequencing was conducted using the Illumina NovaSeq 6000 system (Shanghai Personal Biotechnology Co., Ltd). The transcriptomic data analysis workflow included several steps: initial quality filtering with fastp (v0.22.0) to obtain high-quality FASTQ files; read mapping to the reference genome using HISAT2 (v2.1.0); gene expression quantification with HTSeq (v0.9.1), applying FPKM for normalization; identification of differentially expressed genes through DESeq2 (criteria: fold change > 2, p < 0.05); and functional annotation using topGO (v2.50.0) and ClusterProfiler (v4.6.0) for GO and KEGG enrichment analyses. All bioinformatics analyses were performed on the GenesCloud platform provided by Personalbio (https://www.genescloud.cn).

Statistical analyses were performed using R (v4.2.3) and GraphPad Prism (v10.1.2). Diagnostic accuracy was evaluated via ROC curves (“survivalROC” package) with AUC calculation. Nomograms and calibration curves were generated using “rms”. Functional enrichment (GO/KEGG) was conducted with “clusterProfiler”. Continuous variables were compared using Wilcoxon rank-sum (two groups) or Kruskal-Wallis (≥3 groups) tests, while categorical variables used Fisher’s exact/Chi-square tests. Unpaired t-tests analyzed experimental comparisons. P<0.05 was considered statistically significant.

After quality control, a total of 31,797 genes and 31,179 cells were obtained. We applied the MAESTRO pipeline for single-cell data analysis, identifying the top 2,000 variable features and using PCA for dimensionality reduction, KNN, and Louvain algorithms to identify clusters. To capture cell differences across varying datasets, we set the number of principal components to 30 and used a clustering resolution of 1. We employed t-SNE method to further reduce dimensions and visualize clustering results. We first determined the malignancy of these clusters, identifying 26,055 immune cells, 3,995 malignant cells, and 1,129 stromal cells, as shown in Figure 1A. The results of the cell clustering are displayed in Figure 1B, where 25 clusters were identified. The number of cells in each cluster is presented in Supplementary Table 1. Using the marker-based annotation method in MAESTRO, cell clusters were annotated based on differentially expressed (DE) genes, with marker genes collected from published resources. The annotated clustering results are shown in Figure 1C. Finally, we analyzed MUC family gene expression across different immune cells, finding the highest expression in malignant cells (Figure 1D). Due to dataset limitations (lacking normal/metastatic lymph node data), we analyzed peripheral blood mononuclear cells (PBMCs) and tumor cells to characterize MUC family gene expression patterns. Subgroup analyses by N-stage revealed metastatic-associated expression differences. Notably, MUC family genes showed significantly elevated expression in tumor tissues, particularly malignant cells, compared to normal cells (Supplementary Figure 1).

We performed differential gene expression analysis on lung cancer samples, identifying 158 DEGs (adj.P.Val < 0.05, |log1.3FC| > 1) between lymph node metastasis (N1+N2+N3) and non-metastasis (N0) samples. Among these, 141 genes were upregulated in the lymph node metastasis group, while 17 genes were downregulated (Figure 2A). A heatmap of differential gene expression is shown in Figure 2B. We also analyzed MUC family gene expression differences between metastasis and non-metastasis groups, finding significant differences for MUC5B and MUC21 (Figure 2C).

Further analysis revealed that only MUC5B could stratify lung adenocarcinoma patients into high- and low-risk groups. The expression in tumor versus normal tissues is shown in Figure 3A, and the prognostic Kaplan-Meier curve is shown in Figure 3B. MUC21 was not significantly associated with prognosis. We extended progression-free survival (PFS) analysis using TCGA survival data, revealing significant association between MUC5B and PFS in lung adenocarcinoma (p<0.05), reinforcing its lymph node metastasis relevance (Supplementary Figure 2).

To identify genes associated with lymph node metastasis, we performed weighted gene co-expression network analysis (WGCNA) on 16,033 genes, with MAD > 0.15. Hierarchical clustering of samples is shown in Figure 4A. We used the R package WGCNA to construct a weighted co-expression network, with a soft threshold of 10 for module selection. The network showed characteristics of a scale-free network, with a negative correlation between the log(k) of node connectivity and the log(P(k)) of node occurrence. Using average-linkage hierarchical clustering, we grouped genes into modules, with a minimum of 50 genes per module. We identified 12 modules (Figure 4C), excluding the gray module (unclustered genes). Correlation clustering of these modules is shown in Figure 4D. We analyzed the correlation between each module and immune scores, finding strong correlations for the blue (1,217 genes), green yellow (114 genes), and pink (271 genes) modules with lymph node metastasis. Gene significance (GS) and module membership (MM) relationships are shown in Figure 4F. We retained genes with MM > 0.7 and GS > 0.1, identifying 276 co-expressed genes related to lymph node metastasis.

Pearson correlation analysis with MUC5B revealed 3,532 genes significantly associated with MUC5B, from which 25 candidate genes were identified based on overlap with lymph node metastasis-associated genes (Figure 5A).

Next, we used the GeneMANIA database to construct the protein-protein interaction (PPI) network of 25 candidate genes (Figure 5B). Using the cytoHubba plugin, we identified core genes within this network, retaining 23 genes with a degree greater than 30 as key genes (Figure 5C). To enhance the PPI analysis, we included Maximal Clique Centrality (MNC) and Betweenness centrality for gene screening. Key genes consistently ranked highly across all parameters (Degree, MNC, Betweenness), demonstrating their network significance (Supplementary Figure 3). Finally, we performed functional enrichment analysis of these key genes using GO and KEGG databases (Figures 5D, E). Our GO analysis reveals that the model’s key genes are significantly enriched in mitochondrial membrane organization, DNA replication/repair, and chromosomal structure maintenance pathways, suggesting two clinically relevant mechanisms underlying lung adenocarcinoma metastasis. The strong association with inner mitochondrial membrane organization and intermembrane space components indicates potential metabolic reprogramming in metastatic cells, consistent with clinical observations of enhanced oxidative phosphorylation and altered mitochondrial morphology in circulating tumor cells. The DNA-related processes, including break-induced repair and replication preinitiation, reveal molecular vulnerabilities that may explain the high mutational burden in lymph node metastases and poor chemotherapy response observed in clinical cohorts. Notably, the enrichment of protein-DNA complexes highlights potential therapeutic targets, such as PARP inhibitors for tumors with DNA repair defects and MCM complex inhibitors currently in clinical trials. These findings collectively suggest that mitochondrial dysfunction and genomic instability jointly contribute to metastatic progression, with DNA replication/repair mechanisms serving as both biomarkers for metastatic risk and potential targets for therapeutic intervention. KEGG analysis revealed key pathways including cell cycle, DNA replication, and ubiquitin-mediated proteolysis for proliferation control, along with galactose and nucleotide sugar metabolism for cellular energetics. Notably, Human T-cell leukemia virus 1 infection pathway enrichment suggests potential immune evasion mechanisms. These findings collectively indicate our model genes coordinate tumor growth through proliferation, metabolic reprogramming, and immune modulation - hallmarks of cancer progression.

In addition to WGCNA, two machine learning methods were used to identify predictors of lymph node metastasis in lung adenocarcinoma. LASSO regression selected features from 23 core genes, identifying 23 optimal variables (Figure 6A). SVM-RFE further refined the selection, yielding 13 genes (Figure 6B). Their intersection produced 13 key model genes (Figure 6C), including TIMM10, EXOSC5, SF3B5, CDC20, PAK1IP1, CHCHD3, GINS4, GALK1, ANAPC11, GINS2, TIMM8A, RPL39L, and GINS1. SHAP analysis assessed the contribution of each gene (Figures 6D, E).

Immune scores were calculated using GSVA, and Pearson correlation was performed between the 14 key genes (13 model genes plus MUC5B) and immune components (Figure 6F). Most genes showed positive correlations with Th2 cells and negative correlations with mast cells and type II interferons. Th2 cells secrete cytokines such as IL-4, IL-5, and IL-13, which not only induce the polarization of M2-type macrophages but also directly suppress the activation and proliferation of CD8⁺ cytotoxic T cells. As core anti-tumor immune cells, the impaired function of CD8⁺ T cells significantly weakens the body’s ability to kill tumor cells, thereby clearing immune obstacles for tumor metastasis. On the other hand, the negative correlation between high model gene expression and mast cells reduces the release of chemokines (e.g., CXCL10, CCL3, CCL5) by mast cells. These chemokines are key signaling molecules that recruit CD4⁺ T cells and CD8⁺ T cells to the tumor immune microenvironment; reduced expression of these chemokines leads to insufficient infiltration of anti-tumor immune cells, further weakening immune surveillance. Additionally, the negative correlation between model genes and type II interferon reduces the inhibitory effect of type II interferon on Th2 cell differentiation, while also weakening its ability to promote the maturation of antigen-presenting cells. Impaired maturation of APCs prevents effective presentation of tumor antigens to T cells, ultimately blocking the initiation of anti-tumor immune responses. These findings suggest that high expression of model genes may indicate an immunosuppressive microenvironment conducive to immune evasion.

Finally, NetworkAnalyst was used to explore miRNAs and transcription factors (TFs) regulating MUC5B (Figure 6G). miRNAs such as hsa-miR-484, hsa-miR-335-5p, and hsa-miR-182-5p may modulate MUC5B through interactions with TFs. These miRNAs are potentially involved in LUAD metastasis and immune evasion by affecting tumor cell adhesion and migration. Key TFs regulating MUC5B included ELF5, SREBF2, STAT3, PPARG, and CREB1, suggesting MUC5B may drive lymph node metastasis via these transcriptional networks.

To further evaluate the diagnostic accuracy of the identified genes for predicting lymph node metastasis in lung adenocarcinoma, we performed ROC analysis and validated the model using the GSE13213 and GSE11969 dataset. ROC analysis demonstrated an AUC of 0.931 in the training set (Figure 7A). The prognostic Kaplan-Meier curve (Figure 7B) showed that the model successfully stratified patients into high- and low-risk groups. In the GSE13213 validation set, the model achieved an AUC of 0.925 (Figure 7C), and the KM curve (Figure 7D) confirmed the model’s ability to separate patients into distinct risk groups. AUC of ROC in GSE11969 was 0.894 (0.809–0.947). These results suggest the 13-gene model has high diagnostic and prognostic value.

We collected tumor tissues from 65 patients with lung adenocarcinoma and performed immunohistochemical (IHC) staining to assess MUC5B expression. Among these patients, 30 had no lymph node metastasis, while 35 had lymph node metastasis. The IHC results showed that MUC5B expression was significantly higher in lung adenocarcinoma tissues with lymph node metastasis compared to those without metastasis (Figures 8A, B). The positive signal of MUC5B was mainly localized in the cytoplasm of lung adenocarcinoma cells, showing a granular or diffuse brownish-yellow distribution. Furthermore, we calculated the H-score for each sample, and the H-score for MUC5B expression in the group without lymph node metastasis was significantly lower than in the metastasis-positive group (P = 0.0125, Figure 8C). Based on the H-scores, patients were divided into two groups: MUC5B_high (top 50%) and MUC5B_low (bottom 50%). Kaplan-Meier survival analysis demonstrated that patients with higher MUC5B expression had a significantly poorer prognosis compared to those with lower expression (HR: 2.08, 95% CI: 1.01-4.3, P = 0.0427, Figure 8D). To avoid the impact of basic factors on survival analysis, we collected data on patients’ gender, age, TNM stage, Grade, chemotherapy, and radiotherapy, and supplemented our analysis by performing a multivariate Cox regression analysis comparing patients with high and low H-score. The results showed that only the H-score of MUC5B was a significant prognostic factor (HR: 2.42, 95%CI: 1.11–5.25, P = 0.026). These findings suggest that MUC5B may play an important role in promoting lymph node metastasis and is associated with worse overall survival in lung adenocarcinoma patients.

We used siRNA targeting MUC5B (siMUC5B) and a negative control siRNA (siNC) to transfect both A549 and H1975 lung adenocarcinoma cells. Cells were co-transfected with FAM-labeled siNC, and the transfection efficiency was observed to be over 80% via fluorescence microscopy at 48 hours post-transfection. qRT-PCR analysis confirmed a substantial reduction in MUC5B mRNA expression in both cell lines following siMUC5B transfection, compared to the control groups (Figures 9C, D). The CCK-8 viability assay demonstrated a significant decrease in the proliferation rate of both A549 and H1975 cells after MUC5B knockdown (Figures 9A, B). The colony formation assay revealed that MUC5B knockdown drastically inhibited the colony-forming ability of both cell lines, highlighting its key role in promoting LUAD cell proliferation (Figures 9E, F). In the wound healing assay, MUC5B knockdown markedly reduced the migration capacity of both.

cell lines, evidenced by decreased wound closure at 24 hours post-scratch (Figures 9G, H). Additionally, transwell migration and invasion assays further confirmed that MUC5B depletion led to a significant reduction in both the migratory and invasive abilities of A549 and H1975 cells (Figures 9I, J). Overall, these findings collectively suggest that MUC5B plays a pivotal role in facilitating lung adenocarcinoma cell proliferation, migration, invasion, and colony formation, underscoring its potential as a therapeutic target in LUAD.

To explore the downstream molecular effects of MUC5B, we performed RNA-seq analysis comparing lung adenocarcinoma cells with and without MUC5B knockout. A total of 1,431 genes were upregulated and 970 genes were downregulated following MUC5B deletion (Figure 10A). Gene Ontology (GO) and KEGG enrichment analyses of differentially expressed genes showed significant enrichment in pathways related to DNA replication and cell cycle regulation (Figures 10B, C). Notably, three genes from our 13-gene predictive model—GINS1, GINS2, and GINS4—were significantly downregulated upon MUC5B knockout (Figure 10D). These genes are all core members of the GINS family and together with GINS3, form a heterotetrameric complex essential for DNA replication initiation and elongation (29). This observation prompted further investigation. We subsequently validated the expression of GINS1, GINS2, GINS3 and GINS4 using quantitative PCR (qPCR) in A549 cells. Consistently, all four genes exhibited markedly decreased expression in MUC5B-deficient cells. We hypothesize that MUC5B may promote lymph node metastasis in lung adenocarcinoma by modulating the expression of the GINS gene family, thereby enhancing DNA replication and cell cycle progression.