RNA-sequencing transcriptomic information and associated clinical data for 15 patterns of PCD were obtained from The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov/; using lung adenocarcinoma as a key word). This dataset comprised 412 LUAD cases with 43 adjacent normal tissue samples. Clinical variables for the patients were collected, included sex, age, tumor stage and Tumor-Node-Metastasis (TNM) classification. In total, 1,670 PRGs were acquired via screening the GeneCards database (https://www.genecards.org/) and reviewing relevant literature. Expression data for these 1,670 PRGs were extracted from TCGA LUAD cohort for further analysis. To validate the findings obtained from TCGA data, four datasets from the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo) were employed: GSE8894 (24), GSE31201 (25), GSE30219 (26) and GSE42127 (27). A comprehensive list of the analyzed genes is presented in Table SI.

Wilcoxon test in R (version R 4.1.2; www.r-project.org) was used to determine differentially expressed PRGs between adjacent normal and LUAD tissues. The criteria for significance included a false discovery rate (FDR) <0.05 and an absolute log₂-fold change >1. The vioplot R package (https://rdocumentation.org/packages/vioplot/versions/0.4.0) was utilized to produce heatmaps and volcano plots, which visually depict the differential expression of PRGs in LUAD compared with that in adjacent normal tissue.

WGCNA is a comprehensive biological approach widely employed to analyze genetic association patterns across diverse samples, which is particularly effective for the identification of highly co-expressed genes. By examining the interrelatedness of genes and their associations with phenotypic traits, WGCNA facilitates the identification of potential candidate markers (28). Using the WGCNA R package (http://www.genetics.ucla.edu/labs/horvath/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/index.html), a gene co-expression network for PCD was developed. The selection of the soft threshold power was performed using the pickSoftThreshold function within the WGCNA package. Pearson's correlation analysis was then performed to evaluate the relationships between the LUAD index and module eigengenes. The maximum Pearson correlation coefficient was used to identify the key module most closely correlated with the LUAD index. Genes intersecting within this module were considered candidate hub genes relevant to PCD in LUAD. To further investigate the relationships among PRGs, Spearman's correlation analysis was conducted. Protein-protein interactions (PPIs) of PCD-related proteins were investigated using the STRING database (https://string-db.org/), resulting in the establishment of a PPI network. Critical modules and hub genes were identified using the cytohubba and MCODE plugins in Cytoscape 3.9.1 (https://cytoscape.org/). In addition, Gene Ontology (GO; http://geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.kegg.jp) pathway analyses were performed to functionally annotate the PRGs.

The associations between the PRGs and overall survival (OS) were assessed using univariate Cox regression analysis. To refine the selection of candidate genes and construct a prognostic model, the least absolute shrinkage and selection operator (LASSO) regression algorithm was employed via the R package glmnet (https://rdocumentation.org/packages/glmnet/versions/4.1–7). Ultimately, six genes, along with their associated coefficients, were retained, with the penalty parameter (λ) selected based on the minimum criteria. The risk score was calculated by multiplying the gene expression levels obtained from the LASSO Cox regression model by their corresponding coefficients. This risk score was then applied to classify the patients with LUAD into high- and low-risk groups. Kaplan-Meier survival analysis was performed to compare the OS probability between these groups, and receiver operating characteristic (ROC) curve analysis was conducted using the survminer (https://rdocumentation.org/packages/survminer/versions/0.4.9), timeROC (https://rdocumentation.org/packages/timeROC/versions/0.4) and survival (version 2.41-0; http://github.com/therneau/survival) R packages. In addition, relevant clinical variables across different risk groups were compared using the χ² test, and the results were visualized as a heatmap.

To validate the prognostic signature, the aforementioned GSE8894, GSE31201, GSE30219 and GSE42127 GEO datasets were used. Univariate and multivariate Cox regression analyses were conducted using the survivalROC and survival R packages to evaluate whether the risk score served as an independent prognosticative indicator. In addition, prognostic nomograms were developed using the rms (https://rdocumentation.org/packages/rms/versions/6.5–0) R package, incorporating risk scores and clinical stage to predict the OS of patients with LUAD. Calibration plots were constructed to evaluate the concordance between the predicted and observed 1-, 2-, 3- and 5-year OS rates.

GSEA was employed to identify potential mechanisms by analyzing the enriched pathways in the low- and high-risk groups. The reference gene sets included HALLMARK, C2:KEGG and C5:GO. The criteria for examination were a normalized enrichment score >1, q-value <0.25, and nominal P-value <0.05.

The profiles of immune components were assessed using several algorithms, including QUANTISEQ, TIMER, XCELL, EPIC, MCPcounter and CIBERSORT, with visualization facilitated by the limma R (https://rdocumentation.org/packages/limma/versions/3.28.14) package. Gene set variation analysis and single-sample GSEA (ssGSEA) were employed to evaluate immune function and immune cell subset infiltration. The ESTIMATE algorithm was utilized to calculate stromal, immune and tumor purity scores on the basis of the proportions of stromal and immune cells. Additionally, the expression of major histocompatibility complex (MHC) molecules was analyzed using specific signatures. Differentially expressed common immune checkpoints, including members of the tumor necrosis factor receptor superfamily, programmed cell death protein 1 (PD-1), cytotoxic T lymphocyte-associated protein 4 (CTLA4), programmed cell death ligand 1 (PD-L1) and T-cell immune receptor with immunoglobulin and immunoreceptor tyrosine-based inhibitory motif domains, were visualized using boxplots. The tumor immune dysfunction and exclusion (TIDE) database (tide.dfci.harvard.edu/) was used to calculate TIDE scores, which predict the response to immune checkpoint inhibitor therapies, for patients with LUAD in TCGA dataset. For single-cell analysis, t-distributed stochastic neighbor embedding was employed to visualize the expression of six signature DRGs within the TME. Single-cell sequencing aids in the identification of specific gene expression patterns across distinct cell populations, thereby enhancing our understanding of the functions of prognostic genes in different cell types. Subsequently, the correlations of the six DRGs with various immune cell interactions were analyzed, including B cells, T cells, natural killer (NK) cells, monocytes, macrophages, dendritic cells (DCs) and neutrophils, using the tidyverse R package. These immune cells were further categorized into naive, memory, activated and resting states. To investigate the impact of varying levels of BDNF gene expression on the response of patients with lung cancer to immune checkpoint inhibitor therapy, four cohorts of patients with NSCLC treated with PD-1/PD-L1 inhibitor were analyzed, including patients undergoing systemic therapy and neoadjuvant therapy: GSE207422 (29), GSE111414 (30), GSE126044 (31) and GSE135222 (32).

Stemness, epithelial-mesenchymal transition (EMT), angiogenic activity and tumorigenic cytokines serve as crucial indicators of malignant tumor properties. Scores for these characteristics in every tumor sample were calculated using the ssGSEA algorithm. Tumor stemness indices (TSIs) for patients diagnosed with LUAD were obtained from prior research (33), which identified that TSIs correlate with higher levels of tumor dedifferentiation and the presence of tumor stem cells. Somatic mutation data were obtained from TCGA, and gene mutation analysis was conducted employing the maftools R package. The tumor mutation burden (TMB) was computed for every patient and compared between the two risk groups. TMB score-based survival analysis was also performed. The cBioPortal for Cancer Genomics database (http://www.cbioportal.org/) was adopted to visualize somatic mutations in the key DRG signature.

The limma R package was used to detect differentially expressed genes between the two risk groups, which were subsequently used for small-molecule drug screening. The list of differentially expressed genes was submitted to the Connectivity Map database (CMap; http://clue.io/) to identify potential compounds relevant to the six-gene signature. The CMap database features gene expression profiles derived from nine cancer cell lines treated with 2,429 compounds, each with detailed annotations (34). Connectivity scores were calculated by matching CMap data with the expression levels of the six signature genes, in which these scores were negatively correlated with the treatment effects of the compounds. The pRRophetic (https://github.com/paulgeeleher/pRRophetic2) R package was also used to predict the half-maximal inhibitory concentration (IC₅₀) values for standard chemotherapeutics in different risk subgroups. Subsequently, the three-dimensional structures of potential therapeutic compounds were retrieved from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/).

Given the limited research on the role of honokiol in inducing PCD in lung adenocarcinoma and as honokiol was among the top 10 drugs identified by CMap with the highest enrichment scores, this compound was selected for further study. A549 cells were treated with honokiol, and sulforhodamine B (SRB) assays along with JC-1 staining were performed to validate the predictive accuracy of the CMap drug analysis. Honokiol and Taxol were obtained from Selleck Chemicals. For the SRB assay, A549 cells were seeded at a density of 8,000 cells per well in 96-well plates and incubated under standard conditions (37°C in 5% CO₂) for 24 h. Honokiol (0, 20, 40, 60 or 80 µM) and Taxol (0, 20, 40, 60 or 80 nM) were added to the treatment groups. Following 24 and 48 h of incubation, the SRB method was used to evaluate cell viability.

A549 cells were seeded into 6-well plates at a density of 1×10⁶ cells per well and subjected to the following treatments to form groups: Control (DMEM only), honokiol (50 µM), Taxol (50 nM), and a combination of honokiol (50 µM) + Taxol (50 nM). The cells were incubated at 37°C with 5% CO₂ for 24 h. Subsequently, staining was performed following the instructions provided with the JC-1 assay kit (IJ03009; Beijing Solarbio Science & Technology Co., Ltd.). Fluorescence images were captured using a fluorescence microscope (Axio Observer A1; Zeiss AG), and the fluorescence intensity ratio was quantified using ImageJ v1.8.0 (National Institutes of Health).

Immunohistochemistry data from the Human Protein Atlas (HPA; http://www.proteinatlas.org/) were used to validate the expression of the proteins corresponding to the prognostic genes in LUAD samples compared with normal lung samples.

Two-sample MR was employed to investigate the relationship between the risk of LUAD and hub genes, with single nucleotide polymorphisms (SNPs) defined as instrumental variables (IVs). Hub gene data were sourced from the Genome Wide Association Study database (https://gwas.mrcieu.ac.uk/). The core gene BDNF (id: prot-a-2122), and LUAD (id: ebi-a-GCST004744, containing 11,273 LUAD samples, 55,483 normal samples and 7,849,324 SNPs) as the representative disease, were selected for MR analysis. The analysis was conducted using the TwoSampleMR package, with inverse variance weighting (IVW) to evaluate the association between the hub gene expression value and the risk of LUAD. The MR-Egger method was applied for extra sensitivity analysis.

The lung cancer A549 (CCL-185^™), LLC (CRL-1642^™) and H1975 (CRL-5908^™) cell lines, and the normal human bronchial epithelial BEAS-2B (CRL-9609^™) cell line were purchased from the American Type Culture Collection. These cells were cultured in DMEM (C11875500BT; Gibco; Thermo Fisher Scientific, Inc.) into which 1% penicillin-streptomycin (15070063; Gibco; Thermo Fisher Scientific, Inc.) and 10% fetal bovine serum (10099-141; Gibco; Thermo Fisher Scientific, Inc.) were added. The cells were kept in an incubator at 37°C with 5% CO₂.

The small interfering RNA (siRNA) sequences and expression plasmids were designed and synthesized by Guangzhou RiboBio Co., Ltd. Specific siRNAs targeting BDNF (siBDNF-1 sense, 5′-TCCTTTTCCTTACTATGGTTATT-3′ and antisense, 5′-AAUAACCAUACUAAGGAAAAGGA-3′; and siBDNF-2 sense, 5′-TTCCTTACTATGGTTATTTCATA-3′ and antisense, 5′-UAUGAAAUAACCAUAGUAAGGAA-3′) and control siRNA (sense, 5′-TCCCAAATCGTCTGACCGATGATCCGTTCAAGAGACGGATCATCGGTCAGACGATTTT-3′ and antisense, 5′-AAAAUCGUCUGACCGAUGAUCCGUCUCUGAACGGAUCAUCGGUCAGACGAUUUGGGA-3′) were purchased from Santa Cruz Biotechnology, Inc. Prior to transfection, LLC cells were plated into 6-well culture plates and allowed to grow until reaching 70–80% confluence within 24 ho. According to the manufacturer's protocol, siRNA transfection was performed using Lipofectamine^® 3000 (Invitrogen; Thermo Fisher Scientific, Inc.) and siRNAs (50 nM). The transfection process was maintained in an incubator (37°C, 5% CO₂) for 8 h, after which the medium was replaced with fresh complete medium for an additional 48 h of incubation. Subsequently, the cells were harvested for further experiments.

Total RNA was isolated from the cell lines using TRIzol^® reagent (cat. no. 15596-026; Ambion; Thermo Fisher Scientific, Inc.). Total RNA was isolated from A549, LLC and BEAS-2B cells. The RNA purity and concentration were measured according to the ratio of absorbance at 260 and 280 nm, using a NanoDrop spectrophotometer (Thermo Fisher Scientific, Inc.). The RNA was then reverse-transcribed into cDNA using the PrimeScript^® RT reagent (Takara Bio, Inc.) according to the manufacturer's protocol. The thermocycling protocol was as follows: Amplification step, 94°C for 5 min; followed by 35 cycles of denaturation for 1 min at 94°C, 1 min of annealing at 55°C, elongation at 72°C for 1 min and a final extension step at 72°C for 1 min. The cDNA then served as the template for qPCR, with β-actin as the reference gene. Standard two-step qPCR amplification was conducted using SYBR Green I (10222ES60; Shanghai Yeasen Biotechnology Co., Ltd.) according to the manufacturer's protocol. A total of 20 µl qPCR reaction system contained 6 µl nuclease-free water, 10 µl SYBR Premix Ex Taq II (2X), 0.4 µl ROX Reference Dye II, 2 µl cDNA, 0.8 µl forward primer (10 µM) and 0.8 µl reverse primer (10 µM). The thermocycling conditions were as follows: 95°C for 30 sec; followed by 40 cycles of denaturation at 95°C for 5 sec, annealing at 60°C for 34 sec, elongation at 95°C for 15 sec, and extension at 60°C for 1 min. Gene expression levels were calculated using the 2^−ΔΔCq method (35). The primers used are listed in Table SII.

For the colony formation assay, transfected cells in the logarithmic growth phase were seeded at a density of 1,000 cells/well into a 6-well plate. After cultivation for 10 days (37°C in 5% CO₂), the cells were fixed with methanol for 20 min at room temperature, and then stained with 0.1% crystal violet for 15 min at room temperature. Images of the colonies formed in each well were captured and the number of colonies was manually counted. In this study, clusters of cells containing >50 cells were considered colonies. For the wound healing assay, transfected cells were seeded at a density of 2×10⁵ cells/well in a 6-well plate. After incubation for 24 h, a scratch tool was used to form a wound in the cell monolayer, and the medium was then replaced with serum-free medium. Images of the wound were captured at 0 h and 48 h using a light microscope. ImageJ v1.8.0 (National Institutes of Health) was used to evaluate the cell migration rate as follows: Wound closure surface area / wound total surface area ×100.

Total cellular proteins were extracted from transfected cells using RIPA lysis buffer, and their concentrations were measured with the Thermo BCA Protein Assay kit (Thermo Fisher Scientific, Inc.). The proteins were subsequently denatured at 95°C for 5 min in a water bath, separated on a 10% SDS-PAGE gel by electrophoresis, and transferred onto a 0.45-µm PVDF membrane. The membrane was then blocked with 5% skimmed milk at room temperature for 1 h, followed by the addition of primary antibody (pRIPK1, 1:1,000, cat. no. HY-p81539, MedChemExpress; RIPK1, 1:1,000, cat. no. 3493T, Cell Signaling Technology, Inc.; Cleaved CASP8, 1:1,000, cat. no. 8592T, Cell Signaling Technology, Inc.; CASP8, 1:1,000, cat. no. 4790T, Cell Signaling Technology, Inc.; Cleaved CASP1, 1:1,000, cat. no. HY-p80587, MedChemExpress; CASP1, 1:1,000, cat. no. HY-p81232, Cell Signaling Technology, Inc.; β-actin, 1:2,000, cat. no. 81115-1, Proteintech Group, Inc.) and incubation at 4°C overnight. After washing with TBST, the membrane was incubated with secondary antibody (anti-rabbit IgG, HRP-linked antibody, 1:3,000, cat. no. 7074P2, Cell Signaling Technology, Inc.) for 1 h at room temperature and washed again with TBST. Finally, a chemiluminescence reagent (Shanghai Biyuntian Biotechnology Co., Ltd.) was applied, and the membranes were developed using an automatic chemiluminescence imaging system (Odyssey; LI-COR Biotechnology).

All data were statistically analyzed using R software (version 4.2.1; www.r-project.org) with a two-tailed P-value of <0.05 considered to indicate a statistically significant difference. Kaplan-Meier analysis and log-rank test were conducted to compare the OS of the groups. Independent predictors were examined using Cox regression analysis. Hazard ratios (HRs) and 95% confidence intervals (CIs) were used to describe relative risk. Correlation analyses were performed using Spearman's method. Categorical variables of clinical factors (sex, stage, age and TNM status) were analyzed using the χ² test or Fisher's exact test. Student's t-test was used to compare two groups of normally distributed data, while ANOVA with Tukey's post hoc test was employed for multiple-group comparisons. Wilcoxon's test was employed for comparing ordinal and non-normally distributed data between subgroups. Data from cellular experiments were processed using GraphPad Prism 7.0 software (Dotmatics). The data are expressed as the mean ± SD from a minimum of three replicates.

Differential expression of the 1,670 PRGs identified from TCGA was analyzed between LUAD samples (n=412) and adjacent normal tissues (n=43). The volcano plot (Fig. 1A) and heatmap (Fig. 1B) demonstrated that 310 PRGs were significantly differentially expressed between the LUAD and normal tissues.

To determine the association between potential gene modules and LUAD, a WGCNA was performed on all candidate genes from the LUAD dataset from TCGA. The WGCNA R package was used to develop a system with the optimum efficacy value (β=3, R²=0.9; Fig. 1C and D). The 1,670 PRGs were grouped into seven modules based on their co-expression patterns, including the gray module. The clustering analysis results for all samples are shown in Fig. 1E, and the co-expression modules are visually represented in Fig. 1F using different colors. The results of the module correlation analysis indicate that the yellow module is the most strongly correlated with LUAD (r=−0.85; P<0.0001).

To determine the core differentially expressed PRGs in LUAD, the 310 differentially expressed genes from TCGA and the yellow module (150 genes) from the WGCNA analysis were imported into a Venn diagram. There were 84 overlapping genes that were identified as core genes for subsequent analysis (Fig. 2A). Subsequently, the differentially expressed genes were subjected to GO and KEGG pathway enrichment analyses. The GO analyses revealed that the genes were predominantly enriched in ‘cellular response to interleukin-1’, ‘cellular response to tumor necrosis factor’, ‘inflammatory response’, ‘apoptotic process’ and ‘positive regulation of apoptotic process’ (Fig. 2B). The KEGG analyses revealed that the genes were predominantly enriched in ‘PI3K-Akt signaling pathway’, ‘IL-17 signaling pathway’, ‘ECM-receptor interaction’, ‘Cytokine-cytokine receptor interaction’ and ‘PPAR signaling pathway’ (Fig. 2C). Correlations between the downregulated and upregulated genes are illustrated in Fig. 2D. A PPI network was then constructed utilizing Cytoscape software in conjunction with the STRING database (Fig. 2E). From this network, the top 10 hub genes were identified through degree analysis. These were PF4, IL6, CXCL12, BDNF, CAV1, FGF10, SPP1, GJA1, DCN and MMP13, as shown in Fig. 2F (colored nodes). Additionally, two modules were identified using the MCODE plugin (Fig. 2G and H).

Considering the association between PCD regulators and OS in patients with LUAD, univariate Cox regression analysis was conducted on the expression levels of 84 PRGs to assess their clinical relevance, and the top 20 PRGs with hazard ratios are displayed in a forest chart. The analysis revealed that adrenergic receptor b2 (ADRB2), microtubule associated protein 6 (MAP6) and recombinant integrin α8 functioned as protective factors, while BDNF, G protein-coupled estrogen receptor 1 (GPER1), angiopoietin-like protein 4 (ANGPTL4) and laminin b3 (LAMB3) were identified as risk factors (Fig. 3A). Additionally, correlation analyses showed that most of these protective and risk factor genes were interrelated (Fig. 3B).

To examine the clinical relevance of the PRGs further, patients with LUAD were grouped into subcategories according to their gene expression profiles. Optimal clustering was achieved with k=2, resulting in the separation of the patients into two distinct and non-overlapping groups (Figs. 3C-E). Significant variations in OS, age, stage and sex between these clusters were then assessed. The results indicated that patients in cluster 1 had an improved prognosis (P=0.003) compared with that of patients in cluster 2 (Fig. 3G). Additionally, cluster 1 exhibited lower N stages compared with cluster 2 (Fig. 3F). In summary, consensus clustering revealed a notable association between PRG expression patterns and different clinical parameters.

LASSO Cox regression analysis indicated that BNDF, GPER1, ANGPTL4 and LAMB3 have the highest predictive capabilities (Fig. 3H), as these were risk genes with HR ≥1. By contrast, ADRB2 and MAP6 were protective genes with HR <1. In addition, the LASSO Cox regression analysis identified six genes with the highest predictive power (Fig. 3I and J). By utilizing the coefficients derived from the LASSO algorithm, six optimal genes, namely BDNF, GPER1, ANGPTL4, LAMB3, ADRB2 and MAP6, were used to construct the risk model (Fig. 3K). Accordingly, the risk score derived from these coefficients may be calculated as follows: Risk score = (0.224877 × expression level of BDNF) + (0.060886 × expression level of GPER1)-(−0.08495 × expression level of ADRB2)-(−0.17324 × expression level of MAP6)-(0.005609 × expression level of ANGPTL4)-0.002494 × expression level of LAMB3). The associations between the risk score and the six genes are illustrated in Fig. 3L.

To assess the survival prediction value of these gene signature models, patients with LUAD were categorized into low- and high-risk groups according to the median risk score. A heatmap of clinically relevant characteristics was created to illustrate the differential expression of the six prognostic PRGs between the high- and low-risk groups (Fig. 4A). Significant variations in the clinical data were observed for sex (P<0.001), follow-up status (P<0.001), stage (P<0.05), lymph node metastasis (P<0.001) and tumor extent (P<0.01). Cox univariate analysis revealed that risk scores, tumor stage, tumor extent and lymph node metastasis were significantly associated with OS in patients with LUAD (P<0.001; Fig. 4B). Multivariate Cox regression analysis was then conducted to evaluate whether the risk score independently predicts OS in patients with LUAD, independent of other clinicopathological characteristics. The results demonstrated that the risk score (P<0.001), stage (P=0.008), and M stage (P=0.043) were independently associated with OS (P<0.001; Fig. 4C). In summary, the six-gene risk signature predicts the prognosis of patients with LUAD independently of other clinicopathological factors, including grade, histological age, sex and TNM stage.

Survival analyses were conducted using data from 371 patients with LUAD from TCGA, and the results showed that the OS of high-risk patients was significantly poorer compared with that of patients in the low-risk group (P<0.001; Fig. 4D). The 5-year OS rates were 48.3% for the low-risk group and 24.5% for the high-risk group. ROC curve analysis demonstrated that the area under the curve (AUC) for 1-, 3- and 5-year OS was 0.719, 0.712 and 0.628, respectively, indicating a strong predictive power for clinical outcomes. In additional, the distribution of risk scores among the patients with LUAD was determined, and survival status illustrated utilizing a dot matrix (Fig. S1).

The predictive capacity of the risk model was assessed in four GEO datasets for validation. Patients in the cohorts were classified into low- and high-risk groups based on the cutoff value of TCGA cohort. Survival analysis outcomes indicated that the low-risk group experienced significantly improved OS (GSE8894, P=0.033; GSE31201, P<0.001; GSE30219, P=0.011; and GSE42127, P=0.01) compared with that of the high-risk group (Figs. 4E-H). This finding is consistent with the results observed in TCGA cohort. The AUC for 1-year OS ranged from 0.566 to 0.719, for 3-year OS from 0.652 to 0.741 and for 5-year OS from 0.642 to 0.705. These findings validate the predictive accuracy of the risk model (Figs. 4E-H). Analysis of risk scores among various clinicopathological factors revealed that patients who died had a higher risk score than those who survived (P=0.044), and patients with more advanced disease, those with T2-4 tumors (P=0.0003), N1-3 status (P=0.0001) and stage II–IV (P=0.002), had higher risk scores. Consequently, higher risk scores were found to be associated with more advanced tumor stages (Fig. 4I).

To support the clinical implementation of the model, a comprehensive prognostic nomogram was developed incorporating the LUAD stage and risk score (Fig. 5A). This nomogram effectively predicted the 1-, 3- and 5-year OS of patients with LUAD. The AUC for the risk score was 0.742, indicating a strong association between the predictive effectiveness of the model and the tumor stage (Fig. 5B). These results confirm the good predictive value of the model for patients with LUAD. Calibration plots indicated a strong agreement between the actual and predicted survival rates at 1-, 2-, 3- and 5- years (Fig. 5C).

Enrichment analyses were performed using GO (Fig. 5D), KEGG (Fig. 5E) and HALLMARK (Fig. 5F) databases to identify PRGs with differential expression between the high- and low-risk groups and to explore subgroup survival benefits associated with biological processes and signal pathways. The analyses revealed that high-risk groups were associated with weaker immune response and inflammation-related functions, including ‘JAK-STAT signaling pathway’, ‘IL6-JAK-STAT signaling’, ‘cytokine receptor activity’, ‘inflammatory response’, ‘cytokine receptor activity’ and ‘interferon gamma response’. Conversely, cancer-promoting functions such as ‘DNA replication’, ‘MYC targets’ and ‘PI3K-AKT-MTOR signaling’ were upregulated in the high-risk group, indicating their association with worse prognosis. Selected results of the enrichment are depicted in Figs. 5G-J.

Functional enrichment analyses indicated that the primary roles of the PRGs encompass inflammation, immune response and PCD. Therefore, an immune infiltration analysis was performed to confirm these results. Using the CIBERSORT, MCPcounter, TIMER, QUANTISEQ, XCELL and EPIC algorithms, the levels of DCs (P<0.01), CD8⁺ T cells (CIBERSORT, P=0.006; and XCELL, P=0.003), T-helper 1 (Th1) cells (XCELL, P=0.017), NK (XCELL, P<0.001) and B cell plasma (XCELL, P=0.016) were observed to be lower in the high-risk group compared with those in the low-risk group, However, M2 macrophages (CIBERSORT, XCELL and QUANTISEQ, P<0.05) and cancer-associated fibroblast (CAF) (MCPCOUNTER and EPIC, P<0.05) levels were revealed to be higher in the high-risk group than in the low-risk group (Fig. 6A and B). The difference in expression of these factors suggests that immune infiltration may influence patient prognosis. Moreover, the quantification of enrichment fractions suggests that reduced immune function may impact the prognosis of patients in the high-risk groups, due to an association with cytolytic activity (P<0.05), C-C chemokine receptor (P<0.05), human leukocyte antigen (P<0.001), check-point (P<0.01), T cell co-stimulation and interferon response (P<0.05) (Fig. 6C). The ESTIMATE algorithm indicated higher tumor purity (P<0.001), as well as lower immune (P<0.001), ESTIMATE (P<0.001) and stromal (P<0.001) scores, in the high-risk group compared with the low-risk group (Fig. 6D). By contrast, the expression of MHC molecules was significantly higher in the low-risk group (Fig. 6E). Additionally, a distinct association was observed between the low-risk group and elevated expression of multiple immune checkpoints, such as TIGIT (P<0.01), CD276 (P<0.001), CTLA4 (P<0.01) and CD244 (P<0.001) (Fig. 6F), suggesting that these patients could be advantageously treated with immune checkpoint inhibitors. Notably, the low-risk group also exhibited higher tumor TIDE (P<0.001) and M2 tumor-associated macrophage scores (P=0.004), and lower levels of CAFs (P=0.004) and myeloid-derived suppressor cells (P<0.001) (Fig. 6G). Furthermore, lower TSIs, including differentially methylated probes-based stemness index (P=0.016), epiregulin-mDNA stemness index (P<0.001), enhancer-based stemness index (P=0.004) and mRNA stemness index (P<0.001), were observed in the low-risk group (Fig. 6H). The high-risk group had a lower stemness score (P=0.014) but a higher angiogenic activity score (P<0.001) (Fig. 6I). No significant differences in tumorigenic cytokines or EMT scores were identified between the high- and low-risk groups. Moreover, assessment of the correlation of the risk score with malignant features (Fig. 6J) indicated that the risk score had an association with angiogenic activity score (R=0.35, P<0.001) and stemness score (R=−0.15, P=0.004).

Single-cell analysis revealed limited expression of the six DRGs in immune cells (Fig. 7). For example, GPER1 exhibited low expression in CD8⁺ T cells, and ADRB2 was mainly expressed in macrophages, T cells and B cells, while LAMB3 exhibited slightly elevated expression in monocytes, NK cells and B cells. Combined with the previous immune infiltration analysis, these findings suggest that the six DRGs are primarily expressed in tumor cells and may serve as predictors for patient risk stratification. Further analysis of immune infiltration and single-gene associations (Fig. S2) revealed that BDNF expression was positively correlated with the mutual promotion of T cells CD4 memory activated, T cells CD8, macrophages and T cells CD4 memory activated. GPER1 expression was mainly positively correlated with the mutual promotion of monocytes and T cells CD4 memory resting, as well as mast cells resting and monocytes. ADRB2 was positively correlated with interactions between neutrophils and macrophages M2. MAP6 exhibited a positive correlation with the mutual promotion of monocytes, macrophages M2, DCs resting, DCs activated, mast cells resting and eosinophils. ANGPTL4 was positively correlated with the mutual promotion of mast cells activated, eosinophils and neutrophils. LAMB3 was positively correlated with interactions between T cells CD8, T cells CD4 memory activated and macrophages M1. These findings indicate that the expression of the six PRGs genes shows no significant variation in TME-associated immune cells, but is closely associated with immune cell function and interactions. Therefore, we hypothesize that PCD primarily occurs in tumor cells within the TME, which significantly promotes the activation and interactions of immune cells. This also suggests that the predictive model constructed in the present study not only estimates patient prognosis and risk stratification but also effectively reflects the immune functional status within the TME.

Four NSCLC immunotherapy cohorts were utilized to validate the impact of BDNF on the response to immune checkpoint inhibitors (Fig. S3). The results from the GSE135222 cohort indicated that patients with low BDNF expression treated with PD-1/PD-L1 inhibitors experienced improved progression-free survival compared with those patients with high BDNF expression (P<0.001; Fig. S3A). The findings from the GSE126044 cohort revealed that patients whose lung cancer was more sensitive to PD-1 inhibitors exhibited lower BDNF expression levels than non-responders (Fig. S3B). In the GSE111414 cohort, due to the limited sample size, the marked variability in BDNF expression within each group resulted in no significant difference being identified between the responder and non-responder groups (Fig. S3C). In the GSE207422 cohort, who received neoadjuvant therapy with PD-1 inhibitors combined with chemotherapy, patients with stage III lung cancer who received this treatment and exhibited a pathological partial response had lower BDNF expression than that of non-responders (Fig. S3D). These results indicate that the expression levels of BDNF can predict the response of patients with NSCLC to immune checkpoint inhibitor therapy.

Basic nucleotide variation data were retrieved from TCGA to investigate the differences in genomic mutations between the two risk groups. TTN (56%), TP53 (55%), CSMD3 (47%), MUC16 (41%) and RYR2 (40%) were the top five genes showing the highest mutation frequencies in the high-risk patients (Fig. 8A), and TP53 (39%), TTN (33%), MUC16 (35%), CSMD3 (29%) and RYR2 (32%) were the top five genes in the low-risk group (Fig. 8B). A higher mutation rate was detected in the high-risk group compared with that in the low-risk group (P<0.05). Analysis of somatic mutation interplay revealed widespread gene mutation co-occurrence across the majority of genes. However, KRAS and KEAP1 mutations were mutually exclusive in the high-risk group (P<0.05) (Fig. 8C), while EGFR mutations were mutually exclusive (do not co-occur with other genetic mutations within the same chromosomal region) in the low-risk group (P<0.05) (Fig. 8D). Additionally, the high-risk group had a higher TMB (P<0.001) (Fig. 8E). Among patients with LUAD, those in the high-TMB group had a higher OS time compared with those in the low-TMB group (P=0.009) (Fig. 8F). The 2- and 7-year survival rates were significantly lower for patients who were high-risk and had a low-TMB than for patients with other risk and TMB combinations (P<0.001) (Fig. 8G). These outcomes suggested that the current PRG-based model effectively optimizes TMB-based survival predictions. Finally, the mutation rates for the signature PRGs were observed to be as follows: BDNF (1.2%), GPER1 (5%), ADRB2 (1.4%), MAP6 (4%), ANGPTL4 (0.6%) and LAMB3 (7%) (Fig. 8H).

Sensitivity analysis of commonly used therapeutic drugs in clinical practice showed that bortezomib, BIBW2992, bleomycin, dasatinib, and bryostatin 1 were more efficacious in the low-risk groups, while camptothecin, CCT007093, ABT.888 and axitinib were more effective in the high-risk group (Fig. 9A). Next, the differentially expressed genes between the two risk groups were categorized into upregulated and downregulated subsets (Fig. 9B). These genes were then queried in the CMap database to identify compounds with potential for the treatment of LUAD. Applying an FDR <0.05 and standardized score for screening, 10 small-molecule compounds with potential therapeutic effects on LUAD were identified, where a negative enrichment score indicates an antitumor effect (Table I). The three-dimensional structures of the top five small-molecule drugs are displayed in the Fig. 9C. One of the top 10 ranked compounds, namely honokiol, was selected for cytotoxicity validation. An SRB assay (Fig. S4) revealed that at 24 and 48 h, honokiol inhibited the growth of A549 cells in a concentration-dependent manner within the concentration range of 2.5–80 µM, with an IC₅₀ of 50 µM. When used in combination with Taxol, honokiol was observed to potentiate the effects of Taxol, allowing for a marked reduction in the required effective dose while maintaining the same therapeutic efficacy. The results of a JC-1 assay indicate that honokiol (50 µM) disrupted the mitochondrial membrane of A549 cells, promoting early apoptosis. When combined with Taxol, this effect was significantly enhanced (Fig. S5). In summary, the PCD prognostic model was used to identify drugs with the potential to effectively treat lung cancer and a preliminary validation was conducted using cell experiments, which may serve as a reference for the development of novel lung cancer therapies.

The expression of certain proteins corresponding to the signature PRGs, namely LAMB3 and MAP6, was higher in normal lung tissue than in LUAD tissue (Fig. 10A). However, ANGPLT4 and GPER1 did not show clear differences between these tissue types. The expression of BDNF in LUAD tissue was higher than that in normal lung tissue. This analysis combined with the current predictive model indicates that the high expression of BDNF is negatively associated with the prognosis of patients with LUAD. Notably, it was not possible to evaluate ADRB2, as no immunohistochemistry data for this protein is available in the HPA database.

MR is a reliable method for inferring potential causal relationships, which uses SNPs as IVs to assess causal links between exposure factors and outcomes. MR leverages genetic variations strongly associated with exposure factors as IVs to infer causal effects. All SNPs used in the present study were robust IVs. The causal effects of each genetic variation on LUAD are illustrated in Fig. 10B and C.

A specific analysis of the causal association between BDNF levels and LUAD was performed. Using the IVW method, it was demonstrated that higher BDNF levels were associated with an increased risk of LUAD, with an odds ratio (OR) of 1.048 (95% CI, 1.002–1.096; P=0.04). However, the MR-Egger method, did not show statistical significance (OR, 1.0717; 95% CI, 1–1.1484; P=0.06). As depicted in Fig. 10D, a systematic MR analysis was conducted by iteratively removing each SNP, and the results remained consistent. This suggests that the causal inference was robust across all SNPs, indicating that no single SNP dominated the association between BDNF expression levels and LUAD, thereby validating the MR results. The funnel plot (Fig. 10E) shows a roughly symmetrical distribution, and the MR-Egger regression intercept did not reveal significant horizontal pleiotropy (P=0.416), further supporting the suggestion that that pleiotropy did not bias the causal effect.

The expression levels of BDNF in LLC (P<0.0001), H1975 (P<0.0001) and A549 (P<0.0001) cells were significantly higher than those in BEAS-2B cells (Fig. 11A). In addition, GPER1 was also more highly expressed in LLC (P<0.001), H1975 (P<0.0001) and A549 (P<0.0001) cells. However, ADRB2 as a protective gene, was more highly expressed in BEAS-2B cells than in the three lung cancer cell lines (P<0.0001, respectively). ANGPTL4 gene expression was higher in H1975 (P<0.0001) and A549 (P<0.01) cells than in BEAS-2B cells, and LAMB3 expression was elevated in H1975 (P<0.0001) cells compared with BEAS-2B cells. The expression levels of MAP6 were lower in H1975 cells compared with those in BEAS-2B (P<0.0001) and A549 (P<0.05) cells. These outcomes validate the differential expression of the six PRGs between normal and LUAD samples, highlighting their ability as predictive biomarkers.

Based on the results from PPI, LASSO regression and MR, BDNF was identified as the core gene in the six-gene model, as it exhibits a strong positive association with the occurrence and development of LUAD. Therefore, targeted knockdown of BDNF in LLC cells was performed to investigate its importance (Fig. 11B). A colony formation assay (Fig. 11C) revealed that the knockdown of BDNF suppressed the colony formation of LLC cells. Also, the wound healing time of transfected LLC cells was significantly slower compared with that of the si-NC group at 48 h after BDNF knockdown, indicating that BDNF is crucial for the migration of LLC cells (Fig. 11D).

Western blot analysis indicated that BDNF knockdown induces apoptosis in LLC cells, as evidenced by caspase-8 cleavage. Whether BDNF knockdown also induces necroptosis was also investigated by the measurement of phosphorylated receptor-interacting protein kinase 1 (p-RIPK1), a phosphorylated protein marker that regulates both apoptosis and necroptosis (36). Notably, LLC cells with BDNF-knockdown exhibited an increased p-RIPK1/RIPK1 ratio (P<0.0001). Similarly, BDNF-knockdown promoted the production of cleaved CASP1 (P<0.01) and cleaved CASP8 (P<0.0001), which are hallmarks of pyroptosis (37) and apoptosis (38) (Fig. S6). In summary, these data suggest that BDNF knockdown sensitizes LLC cells to multiple forms of PCD, including apoptosis, necroptosis and pyroptosis.