The genes related with LYM were obtained from Human Gene Database (GeneCards, https://www.genecards.org/) and relevant review articles (19). A total of 466 LYM-related genes were included in this study. The detailed list of these genes can be found in Supplementary Table S1 .

The training dataset comprised transcriptomic profiles and corresponding clinical data for 412 GC patients and 36 control subjects, sourced from the TCGA- STAD database. For external validation, two independent cohorts were utilized: GSE84437 containing 433 GC patients and GSE84437 containing 357 GC patients, both retrieved from the Gene Expression Omnibus (GEO).

A total of 466 LYM-related genes were subjected to differential expression analysis comparing 412 GC samples with 36 normal tissue samples from the TCGA cohort. The “limma”R package was used to identify differentially expressed genes (DEGs) associated with LYM, and genes with an absolute log-fold change (|logFC|) >1 and an adjusted p-value < 0.05 were classified as statistically significant DEGs (20).

The R package “clusterProfiler” was utilized to conduct Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses to evaluate the potential biological pathways associated with LYM-related DEGs (21). The GO analysis comprised three categories: Biological Processes (BP), Molecular Functions (MF), and Cellular Components (CC), with statistical significance thresholds set at an adjusted p-value < 0.05 and q-value < 0.05.

Three independent cohorts comprising 385 GC patients from TCGA, 433 from GSE84437, and 357 from GSE84433 with available survival data were included in the analysis. Univariate Cox regression analysis was performed on the above three databases to identify genes with significant prognostic value (p < 0.05), and common genes were determined through intersection. The least absolute shrinkage and selection operator (LASSO) Cox regression method was employed to identify candidate genes, and the optimal signature was constructed using the “glmnet” package. Subsequently, the LYM score was calculated as follows: LYMS =∑ (βi Genei), where βi denotes the risk coefficient, and Genei represents the expression level of each gene. Patients were stratified into high- and low-LYMS groups based on median LYMS values. Kaplan Meier analysis was conducted using “survival” and “survminer” R packages to assess the association between LYMS and overall survival (OS). Finally, the LYMS was further validated in the GSE84437 and GSE84433 datasets.

An innovative prognostic nomogram was developed by integrating clinical characteristics, including age, TNM stages and LYMS through multivariate Cox and stepwise regression analyses. Calibration plots were used to evaluate the nomogram’s predictive accuracy. Additionally, time-dependent the receiver operating characteristic (ROC) curves, implemented via the R “timeroc” package, were utilized to assess the area under the curve (AUC) for the nomogram.

GSEA was conducted to identify differentially enriched pathways between GC patients with high-LYMS and low-LYMS. The “clusterProfiler” R package was utilized to conduct GSEA. In accordance with established GSEA guidelines, statistically significant results were defined as results meeting a significance threshold of p<0.05 and FDR<0.25 (22).

This study analyzed the association between the LYMS and immune cell infiltration utilizing various algorithms, including the CIBERSORT-ABS, EPIC, QUANTISEQ, TIMER, MCPCOUNTER, XCELL, and EPIC algorithms. Additionally, the stromal score, immune score, tumor purity, and ESTIMATE score were utilized to compare the tumor microenvironment (TME) between GC patients with high and low LYMS through the ESTIMATE algorithm (23).

Additionally, we employed the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu/) to predict immunotherapy responses between LYMS groups (24). Tumor mutation burden (TMB) between high- and low-LYMS groups was compared using somatic mutation data analyzed with the “maftools” package. Immunophenoscore (IPS), obtained from The Cancer Immunome Atlas (TCIA, https://tcia.at/home), was used to assess responses to immunotherapy across different risk groups. To evaluate drug sensitivity, the half maximal inhibitory concentration (IC50) values were calculated using data of GC obtained from Genomics of Drug Sensitivity in Cancer, with predictions generated via the “oncoPredict” package (25).

IHC analysis leverages antigen-antibody specificity binding to detect and localize target antigens in cellular and tissue contexts. We evaluated the expression of key genes in GC and normal tissues using IHC data from the Human Protein Atlas (HPA) database (http://www.proteinatlas.org/) (26).

The GES1, AGS, HGC27 and MKN1 cell lines were obtained from the General Surgery Laboratory of Chinese PLA General Hospital. The cell culture conditions were as follows: RPMI 1640 (G4535A, Servicebio, China) + 10% FBS (G8003, Servicebio, China) at 37 °C and 5% CO2. The steps of quantitative real-time PCR (qRTPCR) were as follows: Total RNA was isolated from the cells using FreeZol reagent (R711, Vazyme, China), and the extracted total RNA was reverse transcribed to cDNA using HiScript IV All-in-one Ultra RT SuperMix reagent (R433, Vazyme, China), and stained using SupRealQ Purple Universal SYBR qPCR Master mix (Q412, Vazyme, China) for staining, and the expression levels of the mRNAs of the 6 genes were detected according to the intensity of the fluorescent signals. β-actin mRNA expression level was used as an endogenous control. Three experiments were performed for each sample, and the results were used to calculate the expression values of the 4 genes according to Equation 2^-△△Ct. The primer sequences are shown in Supplementary Table S2 .

Statistical analyses were performed in R (v4.3.0) and Free Statistics software version 1.9.2. Group differences were analyzed using Student t-test or Wilcoxon test. Survival outcomes associated with categorical variables were assessed via Kaplan-Meier analysis. Statistically significant difference was defined as p<0.05.

In the TCGA-STAD cohort, 128 DEGs were identified between normal samples and gastric tumor samples using the limma analysis. This included 58 upregulated genes and 70 downregulated genes ( Figure 1A , Supplementary Table S3 ). The expression patterns of DEGs were illustrated in the heatmap presented in Figure 1B .

Then, we conducted GO and KEGG enrichment analysis based on the 128 DEGs. GO analysis results displayed the top 10 results in BP, MF and CC ( Figure 1C ). Based on the results of CC, DEGs related to LYM were predominantly enriched in the collagen-containing extracellular matrix (ECM). Interactions between the ECM and lymphatics, as well as the biophysical characteristics of the stroma, influence tumor formation, growth, and metastasis (27). ECM stiffness stimulates the expression of globin transcription factor (GATA) binding protein 2 and GATA2-dependent VEGFR-3, mediating the growth and migration of lymphatic endothelial cell in vivo ( 28). The results of the BP analysis indicated that the DEGs were closely associated to the chemotactic process, including cell chemotaxis, regulation of chemotaxis, and chemokine-mediated signaling pathway. The results of the MF analysis indicated that LYM was closely associated with cytokine activity. VEGFs are key cytokines involved in the LYM process, particularly VEGF-C and VEGF-D, which are known to be the main mediators of lymphatic endothelial cell proliferation and migration (29).

The results of KEGG analysis suggested that the LYM-related DEGs were significantly enriched in the following pathways: cytokine-cytokine receptor interaction, PI3K−AKT signaling pathway, proteoglycans in cancer, chemokine signaling pathway and focal adhesion ( Figure 1D ). These results provided novel insights in the mechanisms of LYM.

This study collected survival data from GC patients and performed further analyses. Univariate Cox regression analysis was employed to identify prognosis-associated genes across three gene databases. The analysis revealed significant associations between genes expression and OS in GC patients across three cohorts: 52 genes in TCGA (p < 0.05), 34 genes in GSE84437 (p < 0.05), and 28 genes in GSE84433 (p < 0.05). 16 genes were identified as common among TCGA, GSE84437, and GSE84433 cohorts ( Figure 2A , Supplementary Figure S1 ). Then, LASSO-Cox regression analysis was conducted on these genes to select the optimal penalty parameter (lambda value λ =0.025). The results of this analysis identified six genes: a disintegrin and metalloprotease with thrombospondin motifs 1 (ADAMTS1), sushi, von Willebrand factor type A, EGF and pentraxin domain containing 1 (SVEP1), caveolin-1 (CAV1), NADPH oxidase 4 (NOX4), neuronal pentraxin 1 (NPTX1), and secreted protein acidic and rich in cysteine (SPARC) ( Figures 2B, C ). Kaplan-Meier analysis demonstrated that each model gene exhibited a significant association with OS in GC patients (p<0.05, Supplementary Figure S2 ), indicating that these genes served as risk factors (HR > 1, p < 0.05). Finally, regression coefficients were calculated for each of the 6 genes. A patient’s risk score was determined utilizing the following formula: LYMS = 0.0797 * ADAMTS1 + 0.0078 * SVEP1 + 0.1338 * NOX4 + 0.0095 * CAV1 + 0.0336 * SPARC + 0.1071 * NPTX1.

According to the median LYMS (median=2.193), all GC patients in the TCGA database were divided into two groups: high-LYMS and low-LYMS. The distribution of risk scores, survival outcomes and differential expression of six signature genes between the two groups were visualized in Figure 2D . Kaplan-Meier survival analysis demonstrated a significantly poorer prognosis for high-LYMS patients compared to the low-LYMS group (p<0.001, Figure 2E ). Time-dependent ROC curve analysis showed that the LYMS predicted OS with AUC values of 0.627 (1-year), 0.660 (3-years), and 0.680 (5-years) ( Figure 2F ).

To validate the LYMS constructed from the TCGA cohort, patients in the GSE84437 and GSE84433 were also classified into high-LYMS and low-LYMS groups according to the median LYMS. Consistent with the training cohort, the high-LYMS group exhibited a significantly poorer prognosis than the low-LYMS group in both validation cohorts (p<0.001) ( Supplementary Figures S3A, B ). The distribution of risk scores and survival outcomes for samples in the validation cohort were illustrated in Supplementary Figures S3C, D , along with the differential expression of six genes in two groups. Additionally, the time-dependent ROC curves were presented in Supplementary Figures S3E, F .

VEGFs are key cytokines involved in the LYM process, particularly VEGF-C and VEGF-D, which are known to be the main mediators of lymphatic endothelial cell proliferation and migration (29). Further elucidating the relationship between LYMS and VEGFs helps to better understand LYMS. We analyzed the correlation between LYMS and VEGFs ( Supplementary Figures S4A–E ). The results revealed that LYMS exhibited positive correlations with VEGFB, VEGFC and VEGFD, but a negative correlation with VEGFA.

This study explored the association of the LYMS with clinical futures, including the T, M, TNM stage and OS in GC patients. The results were presented in Supplementary Figure S5 . Patients with poor prognosis, advanced T stage, and high TNM stage exhibited higher LYM scores.

We further investigated the prognostic significance of LYMS in GC patients. Univariate Cox regression analysis demonstrated that LYMS served as a significant risk factor (HR = 3.82, 95% CI 2.01-7.26, p<0.001, Figure 3A ). Multivariate Cox regression analysis further identified LYMS as an independent prognostic risk factor (HR = 4.62, 95% CI 2.30-9.19, p<0.001, Figure 3B ). Next, we incorporated factors such as LYMS, age, and TNM stage to develop a nomogram for the TCGA cohort. We employed multivariate Cox and stepwise regression analysis to predict the OS of GC patients ( Figure 3C ). The total points were negatively correlated with the survival rates of patients. Kaplan-Meier survival analysis indicated that GC patients with high total points on the nomogram had a worse prognosis compared to those with a low point (p<0.001, Figure 3D ). The calibration curves demonstrated the nomogram’s robust predictive performance for OS at 1, 3, and 5 years ( Figure 3E ). Additionally, when predicting 1, 3-, and 5-years OS in GC patients, the nomogram demonstrated strong ability with AUC values of 0.675, 0.743, and 0.815, respectively ( Figures 3F ). The nomogram showed superior predictive accuracy for OS compared to individual parameters, including age, TNM stage, and LYMS alone ( Supplementary Figure S6 ). Finally, external validation using the GSE84437 and GSE84433 databases confirmed the nomogram’s robust predictive performance ( Supplementary Figure S7 ). These findings suggest the nomogram holds significant potential for clinical prognostication in GC, particularly for long-term (5-year) survival outcomes. Clinically, LYMS can predict the prognosis of gastric cancer and serves as an independent risk factor.

To investigate the association between LYMS and tumor immune features, we examined the immune cell infiltration within the TME. Given the complexity of the TME, we conducted comprehensive analyses including tumor-associated stroma content (Stromalscore), immune cell infiltration levels (Immunescore), tumor purity and the overall characteristics of the tumor microenvironment (Estimatescore) to ensure rigorous and multidimensional evaluation. Compared with CIBERSORT ( Figure 4A ), the results of ( Figure 4B ), EPIC ( Figure 4C ), QUANTISEQ ( Figure 4D ), TIMER ( Figure 4E ), MCPCOUNTER ( Figure 4F ), and XCELL ( Figure 4G ) indicated that the high LYMS group exhibited a higher proportion of B cells, CD8+ T cells, natural killer (NK) cells, macrophages, and other immune cells. However, the results of CIBERSORT ( Figure 4A ) indicated no significant difference or even presented contrary results in T cells and NK cells between two groups. The observed phenomenon may be attributed to the elevated stromal cell abundance within in the TME of the high LYMS group, including fibroblasts, endothelial cells, and matrix components, as demonstrated in Figure 4G . An increased stromal cell proportion could diminish the relative representation of specific immune cells. TME analysis further demonstrated reduced tumor purity scores in the high-LYMS group, whereas stromal, immune, and estimate scores were significantly elevated ( Figures 4H–K ).

To further elucidate functional differences of LYMS, this article performed GSEA on patients with GC. The GSEA results indicated that cell adhesion molecules (CAMs), ECM receptor interactions, and focal adhesion were significantly activated in the high LYMS group ( Figure 5A ). In contrast, DNA replication, nitrogen metabolism, oxidative phosphorylation, and ribosome were inhibited in the high LYMS group ( Figure 5B ). These differential pathways may suggest the potential mechanisms underlying differences between two LYMS groups.

To further explore the role of LYMS in guiding GC immunotherapy, we performed correlation analyses between LYMS and three immunotherapeutic-related biomarkers: TIDE score, TMB, and IPS. TIDE results demonstrated that patients with elevated LYMS exhibited higher TIDE score ( Figures 6A, B ). High TIDE scores generally indicate that tumors possess robust immune evasion mechanisms, which suggests that patients with high LYMS may experience limited benefits from immunotherapy. To further elucidate the role of LYMS in predicting immunotherapy responsiveness, we conducted TIDE predictive analysis. The results demonstrated that GC patients who were responsive to therapy showed lower LYMS levels than non-responders ( Figure 6C ). Collectively, these findings indicated that GC patients with low-LYMS derived enhanced clinical benefits from immunotherapy relative to high-LYMS patients, which was in accordance with the TIDE results.

Subsequently, we examined the differences in TMB between two groups. Analysis indicated that the low LYMS group exhibited higher mutation frequency, including TTN, TP53, and LRP1B ( Figures 6D, E ). The TMB in the low LYMS group was significantly higher than that in the high LYMS group ( Figure 6F ). TMB exhibited a significant negative correlation with LYMS, with a correlation coefficient of -0.27 ( Figure 6G ). Furthermore, survival analysis further demonstrated that patients with high TMB exhibited a markedly improved prognosis compared to those with low TMB ( Figure 6H ). Finally, this study assessed the IPS of GC patients across different risk subgroups ( Figures 6I-L , Supplementary Figure S8 ). Higher IPS score correlate with stronger sample immunogenicity (30). The results indicated that in the cytotoxic T-lymphocyte-associated protein 4 (CTLA4) -positive and programmed death receptor 1 (PD1) -negative subgroup, the low-LYMS group exhibited significantly elevated IPS scores compared to the high-LYMS group. Clinically, these findings indicated that GC patients with low-LYMS might derive enhanced clinical benefits from ICIs.

To investigate the correlation between LYMS and sensitivity to antitumor drugs in GC, we conducted a drug sensitivity analysis. The results indicated that the IC50 values of common drugs including afatinib, gefitinib, dabrafenib, and lapatinib ( Figures 7A–D ) were positively correlated with LYMS. This indicated that patients with high LYMS exhibited insensitivity to these drugs. Conversely, dasatinib, JQ1, NU7441, JAK8517, OTX015 and alpelisib ( Figures 7E–J ) demonstrated a negative correlation with LYMS, indicating that patients with high LYMS were more sensitive to these antitumor drugs, providing references for clinical drug selection.

IHC was employed to validate the expression of LYMS model genes in GC tissues. Compared to normal tissues, ADAMTS1 and SVEP1 expression levels were significantly downregulated in GC tissues, whereas SPARC expression was markedly upregulated ( Figure 8 ). There were no significant differences for CAV1 and NPTX1. Additionally, HPA database lacks information on NOX4 expression in GC. Nevertheless, previous studies have reported significant upregulation of NOX4 in GC tissues compared to normal counterparts (31).

The qRT-PCR analysis revealed distinct expression profiles of six signature genes across gastric cancer cell lines ( Supplementary Figure S9 ). ADAMTS1 demonstrated significantly higher expression levels in HGC27 cells ( Supplementary Figure S9A ). Both NPTX1 and SPARC exhibited consistently low expression across all three cell lines (HGC27, AGS, and MKN1) ( Supplementary Figures S9B, C ). While CAV1 showed elevated expression in HGC27, its expression was markedly reduced in MKN1 ( Supplementary Figure S9D ). Notably, NOX4 displayed substantial upregulation in all tested cell lines (HGC27, AGS, and MKN1) ( Supplementary Figure S9E ). SVEP1 elevated expression was observed in HGC27 and MKN1 ( Supplementary Figure S9F ). The differential expression patterns of these six signature genes between gastric cancer cells and normal gastric mucosal cells further validate the feasibility of the LYMS.