Abstract Background Lung adenocarcinoma (LUAD) is a high-risk malignancy. Telomeres- (TRGs) and aging-related genes (ARGs) play an important role in cancer progression and prognosis. This study aimed to develop a novel prognostic model combined TRGs and ARGs signatures to predict the prognosis of patients with LUAD. Methods LUAD patient’s sample data and clinical data were obtained from public databases. The prognostic model was constructed and evaluated using the least absolute shrinkage and selection operator (LASSO), multivariate Cox analysis, time-dependent receiver operating characteristic (ROC), and Kaplan-Meier (K-M) analysis. Immune cell infiltration levels were assessed using single-sample gene set enrichment analysis (ssGSEA). Antitumor drugs with significant correlations between drug sensitivity and the expression of prognostic genes were identified using the CellMiner database. The distribution and expression levels of prognostic genes in immune cells were subsequently analyzed based on the TISCH database. Results This study identified eight characteristic genes that are significantly associated with LUAD prognosis and could serve as independent prognostic factors, with the low-risk group demonstrating a more favorable outcome. Additionally, a comprehensive nomogram was developed, showing a high degree of prognostic predictive value. The results from ssGSEA indicated that the low-risk group had higher immune cell infiltration. Ultimately, our findings revealed that the high-risk group exhibited heightened sensitivity to the Linsitinib, whereas the low-risk group demonstrated enhanced sensitivity to the OSI-027 drug. Conclusion The risk score exhibited robust prognostic capabilities, offering novel insights for assessing immunotherapy. This will provide a new direction to achieve personalized and precise treatment of LUAD in the future. Supplementary Information The online version contains supplementary material available at 10.1186/s13019-024-03337-y. Keywords: Telomere, Aging, Lung adenocarcinoma, Immunotherapy, Prognosis Introduction Lung cancer is one of the three most prevalent malignant neoplasms globally and is the leading cause of cancer-related mortality worldwide [[30]1]. Lung adenocarcinoma (LUAD) is the most prevalent histologic subtype of non-small cell lung cancer, representing over 50% of cases [[31]1–[32]3]. The primary etiological factor is prolonged exposure to tobacco smoke. However, it should be noted that non-smokers represent 15–20% of cases [[33]4, [34]5]. In conclusion, the development of LUAD is the result of a complex interplay between genetic and environmental factors [[35]4, [36]5]. The available treatment options for LUAD include surgery, radiation therapy, chemotherapy, and targeted therapy [[37]6, [38]7]. Despite the advent of novel diagnostic techniques and therapeutic strategies, including combination therapy and immunotherapy, which have contributed to an increase in patient survival, the five-year survival rate for LUAD patients remains low, at less than 20% [[39]8, [40]9]. Therefore, it is imperative to identify additional biomarkers for the diagnosis and prognostic assessment of LUAD. This will facilitate risk stratification and improve the survival of LUAD patients, as well as enabling patients to make informed decisions regarding the most appropriate therapeutic agent through the use of reliable prognostic assessment indicators. Telomeres are unique structures at the ends of chromosomes consisting of repetitive DNA sequences TTAGGG that play an important role in maintaining genome integrity and functional homeostasis in biological organisms [[41]10]. Changes in telomere length have been linked to a number of biological processes, including growth, proliferation, metastasis and adverse clinical outcomes in a variety of cancerous tumors, including prostate cancer [[42]11] and hepatocellular carcinoma (HCC) [[43]12]. Studies have shown that single nucleotide alterations in telomere length-related genes, including ACYP2, TERC, and TERT, are markedly correlated with the incidence of HCC (P < 0.05) [[44]12]. In addition, cell senescence is an independent risk factor for many chronic diseases and cancers, and the identification of key features of cancer cell senescence and induced senescence has been incorporated into cancer research [[45]13, [46]14]. Telomere shortening causes cellular senescence, which is a normal physiological process in biological organisms, as well as a stable cell cycle arrest mechanism in which cells participate to cope with various oncogenic stresses and inhibit the proliferation of precancerous cells [[47]15–[48]18]. Telomere - (TRGs) and aging-related genes (ARGs) play an important role in inhibiting the proliferation and metastasis of cancer cells and maintaining the normal physiological functions of biological organisms [[49]19–[50]21]. For example, telomerase reverse transcriptase protein (TERT) plays a significant role in antioxidant, anti-inflammatory, anti-aging and cell division. Recent studies have shown that TERT binds to the RNA of the oncogene RMPR, and the complex regulates the RNA level of the oncogene and regulates heterochromatin, providing a new direction in the anti-oncology aspect of combing TRGs with ARGs [[51]22]. Besides, TRGs have also been linked to prognostic assessment in several cancers, such as prostate cancer [[52]23]. In kidney renal clear cell cancer studies, risk score models associated with immune subtypes and tumor burden mutations using TRGs expression levels predict the prognosis of renal cancer patients [[53]24]. There are no reports on the combination of TRGs and ARGs in the prognostic studies of LUAD. Therefore, the objective of this study is to develop a reliable and novel prognostic model that combines the signatures of TRGs and ARGs based on a public data platform and machine learning. This model aims to provide new insights and a deeper genetic understanding of immunotherapy and drug prediction by integrating two oncogene features. Materials and methods Collection of data from multiple databases Transcriptome data of 594 samples (normal:59, tumor:535), mutations and clinical data were downloaded from The Cancer Genome Atlas (TCGA) database ([54]https://portal.gdc.cancer.gov/). The microarray data ([55]GSE72094, [56]GSE50081, [57]GSE30219) were downloaded from the Gene Expression Omnibus (GEO) database ([58]https://www.ncbi.nlm.nih.gov/geo/), and the TCGA data were merged with the [59]GSE72094 data as the training set, [60]GSE50081 and [61]GSE30219 were used as the validation set. The TRGs set was obtained from the previous study [[62]25], and the ARGs set was obtained from Aging Atlas database [[63]26]. Identification of telomeric genes associated with aging This study first used the “edger” package to analyze the difference between normal and tumor tissues of LUAD (|log Fold Change (FC)| > 1, False Discovery Rate (FDR) < 0.05) to get the differential genes, as shown in Supplementary Table [64]S1. Subsequently, the lung adenocarcinoma differential genes (LUAD-DEGs) were intersected with TRGs and ARGs to identify overlapping genes (Supplementary Table [65]S2). We then analyzed the expression levels of these overlapping genes and visualized the results with box plots. The correlations between the differential genes were calculated and displayed as a heat map using the “pheatmap” package. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was used to explore the pathways in which the overlapping genes were located, and Gene Ontology (GO) enrichment analysis was used to analyze and evaluate the overlapping genes at three levels: molecular function (MF), cellular components (CC), and biological processes (BP). Construction and analysis of a prognostic model for TRGs and ARGs signature In this study, clinical data was merged with differential gene expression to screen samples from patients with tumors surviving more than 30 days. To prevent overfitting of the model, the “glmnet” package was used to perform the Least absolute shrinkage and selection operator (LASSO) regression analysis on the candidate genes, and cross-validation was used to select the penalty parameter lambda to remove the genes with strong correlation to reduce the complexity of the model (Supplementary Table [66]S3). The “survival” package was utilized to perform multivariate regression analysis to construct the prognostic model (Supplementary Table [67]S4). The risk score formula of the model was as follows: Inline graphic where Coef(i) is the gene risk coefficient and Exp(i) is the relative gene expression. The samples were then divided into high-risk and low-risk groups based on the median risk score, and Kaplan-Meier (K-M) survival curves were used to evaluate the survival efficiency between the two risk group. Receiver operating characteristic (ROC) curves were plotted using the “timeROC” package, and the areas under the curve (AUC) values of the 1-year, 3-year, and 5-year were calculated to evaluate the performance of risk score in predicting the prognosis of patients. The risk group score distribution, survival status map and gene expression level heatmap were also plotted to evaluate the prognostic differences between the two groups. Finally, validation was performed using the validation cohort. Functional enrichment analysis of risk groups This study analyzed the pathway enrichment of the high-risk and low-risk groups using GSEA 4.3.3 software. Differential expression analysis between the high- and low-risk groups was conducted using the “edger” package (Table [68]S5, |log FC| > 1, FDR < 0.05). Additionally, the present study employed the “clusterProfiler” package to perform GO and KEGG enrichment analyses on the upregulated DEGs (logFC > 1, P < 0.05) and downregulated DEGs (logFC < -1, P < 0.05). Evaluation and construction of prognostic nomogram The present study utilized univariate and multivariate regression analyses to combine the modeled prognostic risk scores with age, gender, TMN staging, and pathological staging. Then, the nomogram and calibration curve were plotted using the “rms” and “survival” packages to predict the overall survival (OS) of patients at 1, 3, and 5 years, and the calibration curve were used to evaluate the accuracy of our nomogram. In addition, we used decision curve analysis (DCA) to assess the net benefit of nomogram and including only clinical variables. Analysis of immune landscape and immunotherapy Then, this study used “GSVA”, “ESTIMATE” package to analyze the single-sample gene set enrichment analysis (ssGSEA) for high and low risk groups. The differences between patients’ tumor microenvironments were also evaluated by plotting violin plots of immunity score, stroma score, ESTIMATE score, and tumor purity score. The level of immune cell infiltration was calculated using the CIBERSORT algorithm, and the expression of immune checkpoints was counted and box plots drawn for the high-risk and low-risk groups. Existing studies have shown that Immunophenoscore (IPS) can be used as a predictive tool for immunotherapy clinical outcomes [[69]27, [70]28], and in this study IPS data were downloaded from The Cancer Immunome Atlas (TCIA, [71]https://tcia.at/home) and we evaluated immunotherapy response in high and low risk groups. In addition, the sensitivity of immunotherapy (anti-PD-1 and anti-CTLA-4 therapies) in this study was analyzed using the tumor immune dysfunction and exclusion (TIDE, [72]http://tide.dfci.harvard.edu/) method and visualized by the TIDE score [[73]29]. P < 0.05 was considered a significant difference. Analysis of tumor mutation burden (TMB) and prediction of drug sensitivity In this study, to discriminate the mutational profiles of TMB between high-risk and low-risk groups, the samples were analyzed using the Wilcox test, the mutation data of the top 20 genes were compiled and counted, and waterfall plots were drawn using the “GenVisR” package. The CellMiner database was utilized to screen antitumor drugs significantly associated with prognostic gene expression. The “pRRophetic” package was used to predict the half maximal inhibitory concentration (IC50) of different drugs in high and low risk groups. Single-cell sequencing for analysis of prognostic genes related to TRGs and ARGs To investigate the distribution of prognostic genes in major cellular subpopulations in the tumor microenvironment (TME), this study used single-cell transcriptome datasets downloaded from the Tumor Immune Single-Cell Hub (TISCH) database ([74]http://tisch.comp-genomics.org/) for analysis. This database contains 79 high-quality single-cell transcriptome datasets from 27 tumors, mainly from the GEO and ArrayExpress databases, which can provide detailed cell type annotations at the single-cell level. And the database has the advantages of comprehensive data, easy operation, user-friendly and data visualization [[75]30]. Based on the TISCH database, the present study visualized the distribution and expression data sets of the prognostic genes related to TRGs and ARGs in [76]GSE99254 using the Uniform Manifold Approximation and Projection (UMAP) plot. In addition, this study also comparatively analyzed the expression levels of the prognostic genes in different immune cell types. Quantitative real-time transcription polymerase chain reaction (qRT-PCR) experimental validation The normal human lung epithelial cell line BEAS-2B (CL-0496) and two lung adenocarcinoma cell lines, A549 (CL-0016) and H1299 (CL-0165), were obtained from PunoSai Bio Co., Ltd. Reverse transcription was performed using the GeniuSaript III Select RT Kit for qPCR (Youji, China). Quantitative reverse transcription polymerase chain reaction (qRT-PCR) was carried out with the uGreener Fex qPCR 2X MiX kit (Youji, China). Gene expression levels were normalized to β-actin, and relative expression was quantified using the 2^−ΔΔCT method. The sequences of the target genes are provided below: Genes Forward primer (5’-3’) Reverse primer (5’-3’) PRKCQ AAAACGTGGACCTCATCTCTG TCATTAGCATTCGGCCTTGAG KLF4 CCTGGCGAGTCTGACATG AAGTCGCTTCATGTGGGAG FBP1 ATAGTGGAACCGGAGAAAAGG GACACAAGGCAATCGATGTTG FoxM1 AGAATTGTCACCTGGAGCAG TTCCTCTCAGTGCTGTTGATG TFAP2A TCTCGATCCACTCCTTACCTC GTTAATAGGGATGGCGGAGAC SNCG GAACATCGCGGTCACCTC CCTGTAGCCCTCTAGTCTCC CHEK2 CCGAACATACAGCAAGAAACAC TCCATTGCCACTGTGATCTTC GAPDH ACATCGCTCAGACACCATG TGTAGTTGAGGTCAATGAAGGG β-actin ATGTGGCCGAGGACTTTGATT AGTGGGGTGGCTTTTAGGATG [77]Open in a new tab Results Identification of analysis candidate genes related with TRGs and ARGs To identify prognostic candidate genes related to TRGs and ARGs, this study first performed differential analysis of normal and LUAD tissues (upregulated genes: 1762, downregulated genes: 3715, Fig. [78]1A). Then, this study analyzed the LUAD-DEGs (|log FC| > 1, FDR < 0.05) with TRGs and ARGs, and took the intersection to obtain 37 differential genes (DEGs), as shown in Fig. [79]1B. To further explore the DEGs, this study analyzed the expression levels of the DEGs in normal tissues and LUAD tumor tissues, and the box plots were drawn as shown in Fig. [80]1C. The results showed that BRCA1, RAD51, E2F1, FoxM1, GAPDH, LMNB1, DNB1, TRAP1, TOP2A, PCNA, PSAT1, TFAP2A, BRCA2, PARP1, CCNA2, UCHL1, RECQL4, RFC4, TERT, FEN1, CDK1, SNCG, CHEK2, FOXL2, BLM, and PRKDC were significantly upregulated in tumor tissues, and PRKCQ, ALDH2, EGR1, PCK1, JUND, PPARG, KLF4, CDKN2B, FBP1, FOS, and JUN were significantly downregulated in tumor tissues. The correlation of DEGs expression was also analyzed, and the results in Fig. [81]1D showed a strong positive correlation among the DEGs. In addition, we performed GO (Fig. [82]1E) and KEGG (Fig. [83]1F) pathway enrichment analysis of the DEGs. The results of GO analysis showed that the DEGs were enriched in the telomere organization, maintenance of telomere homeostasis, regulation of DNA metabolism, and DNA replication. The results of KEGG showed enrichment in pathways for cell cycle, cellular senescence, DNA replication, and immune cell-related functions. Fig. 1. [84]Fig. 1 [85]Open in a new tab Identification and analysis of candidate genes. (A) The volcano map shows the expression levels of the LUAD-DEGs. (B) The Venn diagram shows the intersection of the LUAD-DEGs with the TRGs and the ARGs. (C) The box plot shows the expression levels of the DEGs between the normal samples and the LUAD tumor samples. (D) The heatmap shows the correlation analysis of the DEGs. (E) GO enrichment analysis of the DEGs. (F) KEGG enrichment analysis of the DEGs Construction and validation of prognostic signature To further identify prognostic gene signatures, in this study, as shown in Fig. [86]2A-B, eleven characterized genes were obtained by LASSO analysis. The eleven characterized genes were then subjected to multivariate Cox analysis, and eight characterized genes were finally screened to construct prognostic models (Fig. [87]2C). The eight genes are GAPDH, CHEK2, TFAP2A, FBP1, SNCG, KLF4, PRKCQ, and FoxM1. Risk score based on gene expression levels and regression coefficients was calculated using the following equation: risk score = 0.264*GAPDH − 0.300*CHEK2 + 0.068*TFAP2A − 0.203*FBP1 + 0.054*SNCG + 0.097*KLF4–0.113*PRKCQ + 0.186*FoxM1. Fig. 2. [88]Fig. 2 [89]Open in a new tab Construction of a novel risk signature based on TRGs and ARGs. Cross-validation for adjustment of coefficient distributions (A) and parameter screening (B) in LASSO regression models. (C) Forest plot for multivariate Cox analysis Next, the sample (TCGA + [90]GSE72094) was subjected to risk scoring using the prognostic model. ROC analysis showed that the AUC in predicting the 1-year, 3-year, and 5-year survival of patients were 0.688, 0.708, 0.652 (Fig. [91]3A). As shown in Fig. [92]3B, the findings of the K-M curve analysis indicated that patients in the low-risk group exhibited superior OS compared to those in the high-risk group (P < 0.05). As illustrated in the risk score and survival status distribution plots, the low-risk group exhibited superior survival outcomes compared to the high-risk group. (Fig. [93]3C). Similarly, in the validation cohort, the low-risk group had a superior prognosis. We combined the data of [94]GSE50081 and [95]GSE30219 as the validation cohort, and the results of ROC analysis in Fig. [96]3D showed predicted AUC values of 0.761, 0.73, and 0.713 at 1, 3, and 5 years for patients with LUAD, respectively. Additionally, the K-M curve (Fig. [97]3E) indicated that patients with lower risk scores exhibited a higher survival rate than those with higher risk scores. According to the results of the graph of the risk score and the distribution of the survival status in Fig. [98]3F, the low-risk group had a superior survival status. Furthermore, we compared the prognostic risk model with several other models to evaluate their predictive performance and found that our model demonstrated favorable performance. This suggests that the prognostic model we have established is effective and has potential application value in predicting the prognosis of LUAD patients (Supplementary Figure [99]S1, Supplementary Table [100]S6). Fig. 3. [101]Fig. 3 [102]Open in a new tab Validation and analysis of prognostic risk signature. (A) Analysis of ROC curves in training cohort. (B) K-M curve analysis of high and low risk groups in training cohort. (C) Distribution of risk scores and survival status of high and low risk group in the training cohort. (D) Analysis of ROC curves in validation cohort. (E) K-M curve analysis in validation cohort. (F) Distribution of risk scores and survival status in the validation cohort. (G) Expression levels of prognostic signature genes between high and low risk groups in the training cohort. (H) Expression levels of prognostic signature genes in the validation cohort In this study, the gene expression levels of prognostic characterizing genes in the high-risk and low-risk groups were also examined in the training and validation cohorts, respectively. As shown in Fig. [103]3G, in the training cohort, TFAP2A, GAPDH, FoxM1, SNCG, KLF4 genes were low expressed in the low-risk group, FBP1, PRKCQ genes were low expressed in the high-risk group, and the difference of CHEK2 was not significant. In the validation cohort, TFAP2A, GAPDH, FoxM1, KLF4 genes were low expressed in the low-risk group, PRKCQ genes were low expressed in the high-risk group, and the difference of CHEK2 and SNCG was not significant (Fig. [104]3H). Meanwhile, in this study, the survival rates of eight prognostic signature genes were analyzed by K-M curves, and it was found that the survival rates of FoxM1, GAPDH, KLF4, and TFAP2A genes with different expression levels were significantly different (P < 0.05, Supplementary Figure [105]S2 A-H). Functional enrichment analysis To further explore the potential biological mechanisms between high and low risk groups, this study used GSEA to analyze the functional enrichment of genes within the two risk groups. As shown in Fig. [106]4A-C, the high-risk group was enriched in cell cycle, DNA replication, and p53 signaling pathway. The results in Fig. [107]4D-F showed that the low-risk group was enriched in alpha linolenic acid metabolism, ether lipid metabolism, and fatty acid metabolism pathways. In addition, the present study employed the GO and KEGG enrichment analysis to elucidate the biological processes and pathways that were significantly altered between the high-risk and low-risk groups. Figure [108]4G GO enrichment analysis showed that differentially upregulated genes were enriched in the intermediate filament, cytoskeleton organization, regulation of mitosis, regulation of nuclear division, and intermediate filament. Figure [109]4H shows that differentially downregulated genes were enriched in the cilium movement, humoral immune response, motile cilium and endopeptidase inhibitor activity. GO enrichment analysis revealed that differential genes in the high and low risk groups play an important role in cell cycle regulation, cell metabolism, and activation of cellular immune function. KEGG enrichment analysis showed that differentially upregulated genes were enriched in the neuroactive ligand-receptor interaction, retinol metabolism, cell cycle, and porphyrin metabolism pathways (Fig. [110]4I). KEGG enrichment analysis of differentially downregulated genes showed that most were enriched in the protein digestion and absorption, cAMP signaling pathway, and linoleic acid metabolism pathways (Fig. [111]4J). KEGG analysis results suggest that differential gene signatures may play a critical role in the genesis and development of LUAD through neuromodulatory, cell cycle and proteasome pathways. Fig. 4. [112]Fig. 4 [113]Open in a new tab The representative results of functional enrichment analysis. A representative result of enrichment in the high-risk group is (A), (B) and (C). A representative result of enrichment in the low-risk group is (D), (E) and (F). (G) Results of GO enrichment analysis of upregulated differential genes. (H) Results of GO enrichment analysis of downregulated differential genes. (I) Results of KEGG enrichment analysis of upregulated differential genes. (J) Results of KEGG enrichment analysis of downregulated differential genes Development and evaluation of the prognostic nomogram To investigate the relationship between prognostic models and clinical characteristics, this study first analyzed the differences between gender, clinical risk stage, TNM stage and risk score using Wilcoxon test. The results showed significant differences between gender (Fig. [114]5A, P < 0.05), T stage (Fig. [115]5B, P < 0.05), N stage (Fig. [116]5C, P < 0.05), stage (Fig. [117]5D, P < 0.05) and risk score. In univariate regression analysis, we observed that M stage (hazard ratio (HR) : 1.928, confidence interval (CI) : 1.085–3.428, P = 0.025), stage (HR : 2.584, CI : 1.809–3.692, P < 0.001), T stage (HR : 2. 453, CI : 1.6110–3.738, P < 0.001),N stage (HR : 2.514, CI : 1.780–3.551, P < 0.001), riskScore (HR : 1.378, CI : 1.263–1.504, P < 0.001) were significantly associated with patient prognosis (Fig. [118]5E). In multivariate regression analysis, T stage (HR : 2.081, CI : 1.273-3.400, P = 0.003), N stage (HR : 1.812, CI : 1.208–2.716, P = 0.004), riskScore (HR : 1.333, CI : 1.211–1.468, P < 0.001) were found to be independent predictors (Fig. [119]5F). Then, a nomogram was developed to provide quantitative predictions of 1-year, 3-year, and 5-year OS probabilities in patients (Fig. [120]5G). The results of the calibration curves in Fig. [121]5H show a high degree of agreement between predicted and actual OS values. The present study then used DCA curves to examine the predictive probability of the nomogram for LUAD patients at 1, 3, and 5 years. The results, Fig. [122]5I, show that such a comprehensive nomogram provides more net benefit and is more favorable for clinical management. Fig. 5. [123]Fig. 5 [124]Open in a new tab Nomogram development and evaluation. Violin plots show a significant correlation between (A) gender, (B) T stage, (C) N stage, (D) stage and risk score. (E) Univariate regression of clinicopathologic indicators and gene signatures. (F) Multivariate regression of clinicopathologic indicators and gene signatures. (G) Nomogram for predicting the survival probability of LUAD patients. (H) Calibration curves of the nomogram at 1-year, 3-year, 5-year intervals. (I) DCA curves of clinicopathologic indicators and nomogram Analysis of immune landscape and immunotherapy response Next this study examined the differences in immune-related functions and immune cell infiltration between the high-risk and low-risk groups using ssGSEA analysis. As shown in Fig. [125]6A, immune-related functions and immune cell infiltration were higher in the low-risk group than in the high-risk group. The results of ESTIMATE analysis showed that the low-risk group had higher ESTIMATEScore (Fig. [126]6B, P < 0.01), ImmuneScore (Fig. [127]6C, P < 0.001), StromalScore (Fig. [128]6D, P < 0.05) and lower TumorPurity (Fig. [129]6E, P < 0.01) score than the high-risk group. According to the results of the CIBERSORT algorithm in Fig. [130]6F, the infiltration levels of T cells CD4 memory resting, monocytes, dendritic cells resting and mast cells resting were significantly higher in the low-risk group than in the high-risk group, and the infiltration levels of T cells CD8, T cells CD4 memory activated, NK cells resting, macrophages M0, macrophages M2 and infiltration levels of neutrophils were significantly higher in the high-risk group than in the low-risk group. In addition, we evaluated the expression of immune checkpoint genes, and the results in Fig. [131]6G showed that most of the immune checkpoint genes, such as BTLA, CCL19, ITGAL, INFSF15, had higher expression levels in the low-risk group. To further evaluate the immunotherapy response, we analyzed the IPS (Fig. [132]6H-K), TIDE (Fig. [133]6L), in the high-risk and low-risk groups. IPS results show favorable therapeutic benefit of checkpoint inhibitor treatment in the low-risk group, with patients in the low-risk group having a favorable therapeutic response to anti-CTLA4 and anti-PD-1 immunotherapy. The results showed that the low-risk group had lower TIDE than the high-risk group (P < 0.05), suggesting that patients in the low-risk group had lower chances of anti-tumor immune escape, higher response rates to immune checkpoint blockade (ICB) therapy, and were more likely to benefit from ICB therapy. Fig. 6. [134]Fig. 6 [135]Open in a new tab Analysis of immune landscape and response to immunotherapy. (A) Heatmap shows ssGSEA analysis of differences in immune cell infiltration between the two risk groups. The violin diagrams show the ESTIMATE score (B), immune score (C), stromal score (D), tumor purity (E) of the high and low risk groups. (F) Comparison of the infiltration fractions of the immune cells in the two risk groups. (G) Comparison of differences in immune checkpoint expression between the two risk groups. (H-K) Comparison of IPS in the two risk groups. The violin diagrams show the comparison of TIDE (L) in the two risk group Analysis and comparison of TMB and drug sensitivity Tumor cells evade the effects of immunosurveillance by elevating TMB, and thus TMB may serve as a predictor of immune response [[136]31]. As shown in Fig. [137]7A, the TMB in the low-risk group was significantly lower than that in the high-risk group (P < 0.05). The TMB results of the top 20 mutations in the high-risk group showed that TP53 and TTN mutations were highly significant (Fig. [138]7B). In the high-risk group, 222 (95.28%) of 233 samples were mutated, but in the low-risk group, the mutation rate was 213 (91.81%) of 232 samples (Fig. [139]7C). It further indicates that the TMB in the low-risk group was lower than in the high-risk group and that there was significant variability in the immune landscape between the two risk groups. Fig. 7. [140]Fig. 7 [141]Open in a new tab Analysis of TMB and drug sensitivity between high and low risk groups. (A) The violin plot shows the variability of TMB between high-risk and low-risk groups. The TMB of TOP20 genes in the high-risk group (B) and in the low-risk group (C). (D-G) Correlation of characterized genes with drug sensitivity. (H) Correlation analysis of high and low risk groups and drug Lisitinib sensitivity. (I) Correlation analysis of high and low risk groups and drug OSI-027 sensitivity To further investigate the clinical applicability of TRGs and ARGs for precision treatment of LUAD patients, this study evaluated the correlation between characterized genes and drug sensitivity. The results showed negative correlation between PRKCQ and OSI-027 (Fig. [142]7D, -0.494, P<0.001), negative correlation between KLF4 and Arsenic trioxide (Fig. [143]7E, -0.487, P<0.001), Melphalan (Fig. [144]7G, -0.450, P<0.001), and positive correlation between FBP1 and Linsitinib (Fig. [145]7F and F, 0.470, P<0.001). In addition, this study evaluated the therapeutic effects of commonly used chemotherapeutic agents in different risk groups. The results showed that the sensitivity to Lisitinib was higher in the high-risk group than in the low-risk group (Fig. [146]7H, P < 0.05), and the sensitivity to OSI-027 was higher in the low-risk group than in the high-risk group (Fig. [147]7I, P < 0.05). Correlation of the TRGs and ARGs signature with single-cell properties In recent years, single-cell RNA sequencing (scRNA-seq) has become an important scientific tool to reveal differences in cell populations and to characterize heterogeneous cell populations. To further explore the relationship between characteristic genes and TME, we analyzed data from the TISCH database for scRNA-seq of [148]GSE99254. As shown in the results of Fig. [149]8A, this UMAP plot showed six major cell clusters, namely CD4 Tconv, CD8T, CD8Tex, Mono/Macro, Tprolif, and Treg. By further analysis, the result showed that FoxM1 (Fig. [150]8B and K) and CHEK2 (Fig. [151]8C and M) were mainly distributed in Tprolif and Mono/Macro cell clusters. FBP1 (Fig. [152]8G), TFAP2A (Fig. [153]8I), SNCG (Fig. [154]8J) and KLF4 (Fig. [155]8L) were mainly distributed in Mono/Macro cell clusters, while PRKCQ (Fig. [156]8D and F) and GAPDH (Fig. [157]8E and H) were distributed in almost all cell clusters. The findings suggest that the characteristic signature is closely associated with TME and that the signature has the potential to serve as a biomarker for predicting the efficacy of immunotherapy in LUAD patients. Fig. 8. [158]Fig. 8 [159]Open in a new tab Relevance of the characteristic genes to the properties of a single cell. (A) UMAP plot of six major cell clusters in the LUAD TME. The distribution of CHEK2 (B), FoxM1 (C), PRKCQ (D), GAPDH (E) in cell subsets. Violin plot of PRKCQ (F), FBP1 (G), GAPDH (H), TFAP2A (I), SNCG (J), FOXM1 (K), KLF4 (L), CHEK2 (M) expression at the single cell level Experimental validation of the expression of prognostic biomarkers in LUAD To validate the expression levels of the selected feature genes, this study utilized qRT-PCR analysis in LUAD cell lines. As shown in Fig. [160]9, the experimental results were consistent with the findings from bioinformatics analyses. Specifically, the expression levels of PRKCQ, KLF4, and FBP1 were significantly downregulated in A549 and H1299 cells. Conversely, the expression levels of FoxM1, GAPDH, TFAP2A, SNCG, and CHEK2 were significantly upregulated in tumor cells. Fig. 9. Fig. 9 [161]Open in a new tab The findings of qRT-PCR validation of the expression of prognostic feature genes in LUAD. *P < 0.05 Discussion In recent years, the emergence of multidisciplinary treatment modalities and new therapeutic drugs has led to significant breakthroughs in the treatment of LUAD disease [[162]32]. However, the survival of LUAD patients has not yet been significantly prolonged, and the high mortality rate of LUAD patients remains a challenging problem to solve [[163]6, [164]7]. We still need to explore reliable biomarkers to guide clinical treatment decisions. By integrating transcriptomic data and clinical information from multiple databases, this study developed and validated a risk signature associated with TRGs and ARGs for survival prediction, immune microenvironment evaluation, and antitumor drug sensitivity assessment. From the results of this study, TRGs and ARGs signatures plays an important role in survival outcome, immune microenvironment and ICB treatment efficacy in LUAD patients. In this study, we screened and analyzed eight characteristic signatures (GAPDH, CHEK2, TFAP2A, FBP1, SNCG, KLF4, PRKCQ, and FoxM1). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) is a well-characterized housekeeping protein with a recognized role in cellular metabolism [[165]33]. It is noteworthy that studies have demonstrated a correlation between elevated GAPDH expression and a poor prognosis in cancer patients [[166]34]. The K-M survival analysis in this study further demonstrated that LUAD patients exhibiting high GAPDH expression had a significantly poorer prognosis (P < 0.01). Forkhead box protein M1 (FoxM1) was identified as an oncogene, with high expression observed in numerous cancer types [[167]34]. Du et al. reported that E2F2 suppresses LUAD expression by downregulating the FoxM1/B-Myb axis [[168]35]. Elevated FoxM1 expression inhibits the downstream gene PTPN13, which in turn reduces the sensitivity of LUAD cells to the drug Gefitinib [[169]36]. The checkpoint kinase 2 (CHEK2) gene plays a pivotal role in the DNA damage pathway, regulating both the cell cycle and apoptosis [[170]37]. Transcription factor AP‑2 alpha (TFAP2A) is a key regulator of cell growth and is overexpressed in a variety of tumor tissues, including bladder cancer [[171]38] and LUAD [[172]39]. TFAP2A is aberrantly overexpressed in LUAD tissues, where it activates HMGA1 to promote glycolysis, thereby driving the progression of LUAD [[173]40]. Fructose-1,6-bisphosphatase (FBP1) is an important oncogenic factor in malignant tumors and is a rate-limiting enzyme in gluconeogenesis [[174]41]. Krüppel-like Factor 4 (KLF4), a member of the evolutionarily conserved family of zinc finger transcription factors, plays important roles in embryonic development, inflammation and malignancy, and is an important clinical biomarker and therapeutic target for several types of tumors [[175]42]. In LUAD, the abnormal expression of DNMT1 inhibits KLF4 expression through DNA methylation, thereby promoting cell proliferation and metastasis [[176]43]. PRKCQ/PKCθ is a serine/threonine kinase, a member of the novel PKC family. PRKCQ is predominantly distributed in T cells, NK cells, and it has been reported that PRKCQ regulates T cell survival by modulating the expression of pro-apoptotic and anti-apoptotic Bcl2 family members [[177]44, [178]45]. Previous studies have provided strong evidence of a close association between the signature genes and the progression of LUAD. Our findings further suggest that these eight signature genes may serve as prognostic biomarkers for LUAD, highlighting the importance of exploring their roles in LUAD. This could potentially identify novel therapeutic targets for the treatment of LUAD. We performed functional enrichment analysis to explore the potential molecular mechanisms underlying the differences between high- and low-risk groups. Our results highlighted the significant roles of the cell cycle and proteasome pathways in inhibiting LUAD progression. The cancer immunity cycle is amplified by T-cell-mediated tumor cell killing, effectively inhibiting the malignant progression of tumor cells [[179]46]. This mechanism plays a pivotal role in the efficacy of immunotherapy. Recent discoveries of proteasome inhibitors and immunomodulatory drugs have brought new hope for cancer treatment [[180]47]. Proteasome inhibition potentiates the anticancer efficacy of therapeutic agents by downregulating key apoptotic proteins, including TNF-α and NF-κB [[181]48]. Additionally, our scRNA-seq analysis revealed that FoxM1 and CHEK2 are predominantly expressed in Tprolif and Monon/Macro cells. Importantly, these cells typically exhibit high levels of immune exhaustion markers [[182]49]. qRT-PCR validation confirmed that FoxM1 and CHEK2 are significantly upregulated in LUAD. These findings suggest that the immune-suppressive microenvironment may be linked to the elevated expression of FoxM1 and CHEK2. Targeting these prognostic feature genes may offer substantial clinical value for enhancing immunotherapy in LUAD. A comprehensive understanding of the immune infiltration of TME is necessary to elucidate the molecular mechanisms and improve clinical outcomes to provide innovative immunotherapeutic approaches [[183]50]. In vivo depletion of CD4^+ T cells promotes tumor progression and can regulate the TME by mediating cytokine and co-stimulatory signals [[184]51, [185]52]. Quiescent CD4 memory T cells are differentiated and endowed with various functions, such as assisting CD8^+ T cells in performing anti-tumor functions [[186]53]. Liu et al. reported that increased numbers of M0 macrophages have been reported to correlate with poor prognosis in early LUAD [[187]54]. Our results showed higher levels of T cell CD4 memory cell infiltration in the low-risk group, but T cell CD8, T cell CD4 memory activated and M0 macrophage cells had higher immune cell infiltration in the high-risk group than in the low-risk group. These results suggest that there may be differences in the mechanisms of the immune response to the tumor between the two groups. Immune checkpoints are another key component. The immune checkpoint-related gene CCL19 plays a crucial role in suppressing tumor progression by facilitating lymphocyte recruitment to tertiary lymphoid structures (TLSs) [[188]55]. Furthermore, studies have demonstrated that CCL19 enhances T-cell infiltration and promotes CAR-T cell survival in mouse models [[189]56]. Our findings revealed that CCL19 expression is significantly higher in the low-risk group, indicating that LUAD patients with lower risk scores may experience greater benefits from immunotherapy. On the other hand, elevated TNFSF9 expression appears to have a detrimental impact on patients. Wu et al. reported that TNFSF9 is overexpressed in pancreatic cancer, where it drives M2 macrophage polarization via the Src/FAK/p-Akt/IL-1β signaling pathway, ultimately promoting cancer cell migration [[190]57]. Blocking TNFSF9 signaling has shown potential anti-tumor effects by reshaping the TME. Collectively, these findings highlight the critical roles of immune checkpoint-related genes in tumor suppression and offer valuable perspectives for advancing LUAD immunotherapy research.Furthermore, the majority of immune checkpoint-related genes, including BTLA, CD28, and TNFRSF14, are significantly upregulated in the low-risk score group. The significant expression of immune cell infiltration and immune checkpoint-related genes in LUAD reveals the potential therapeutic value of immunotherapy in LUAD. To provide more precise clinical guidance, this study examine the correlation between drug sensitivity and risk scores. Lisitinib is a potent and selective, orally active dual inhibitor of IGF-1 and the insulin receptor (IR) [[191]58]. Studies have shown that lisitinib inhibits the proliferation of a variety of tumor cell lines, including NSCLC, prostate cancer cells and colorectal cancer (CRC) cell lines [[192]59, [193]60]. This study demonstrated that the high-risk group exhibited greater sensitivity to Lisitinib, suggesting that this drug is more clinically advantageous for the treatment in the high-risk group. OSI-027 is a potent, selective, orally active and ATP-competitive inhibitor of mTOR kinase activity [[194]61, [195]62]. The results of the study showed that drug sensitivity was higher in the low-risk group than in the high-risk group, and OSI-027 may be a good choice for LUAD patients with lower risk scores. In conclusion, although the LUAD transcriptomic data utilized in this study are derived from public databases and further prospective research is necessary, GAPDH, CHEK2, TFAP2A, FBP1, SNCG, KLF4, PRKCQ, and FoxM1 have been identified as key prognostic biomarkers for LUAD. Targeting these genes holds considerable clinical potential for enhancing the treatment of LUAD. Furthermore, the present study successfully constructed a novel prognostic model combined TRGs and ARGs signatures. The construction of risk model provides a new and potentially effective methods for individualized survival prediction and clinical outcome assessment of LUAD patients. And the analysis of immunotherapy, immune cell infiltration level and drug sensitivity based on the risk model provides an effective direction to guide the precise treatment of LUAD patients. Electronic supplementary material Below is the link to the electronic supplementary material. [196]Supplementary Material 1^ (2.7MB, xls) [197]Supplementary Material 2^ (21KB, xls) [198]Supplementary Material 3^ (19.5KB, xls) [199]Supplementary Material 4^ (20.5KB, xls) [200]Supplementary Material 5^ (280.5KB, xls) [201]Supplementary Material 6^ (20.5KB, xls) [202]Supplementary Material 7^ (522.1KB, docx) Author contributions Dr Z.Y, YW.H and TT.C contributed to the study design. YY.W conducted the literature search. Z.Y, YW.H, TT.C and YY.W acquired the data. Z.Y and YY.W wrote the article. YW.H performed data analysis. TT.C drafted. Z.Y, YW.H and YY.W revised the article and gave the final approval of the version to be submitted. All authors read and approved the final manuscript. Funding Not applicable. Data availability No datasets were generated or analysed during the current study. Declarations Ethics approval and consent to participate Not applicable. Consent for publication Not applicable. Competing interests The authors declare no competing interests. Footnotes Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. References