- 1Department of Pathology, Shandong Cancer Hospital and Institute, Shandong First Medical University and Shandong Academy of Medical Sciences, Jinan, China
- 2Geneplus-Beijing Institute, Beijing, China
- 3Geneplus-Shenzhen Clinical Laboratory, Shenzhen, China
- 4Xiangya School of Pharmaceutical Sciences, Central South University, Changsha, China
- 5Key Laboratory Experimental Teratology of the Ministry of Education, Department of Physiology, School of Basic Medical Sciences, Cheeloo College of Medicine, Shandong University, Jinan, China
- 6School of Tropical Medicine and Laboratory Medicine, Hainan Medical University, Haikou, China
- 7Geneplus-Shenzhen Institute, Shenzhen, China
- 8Department of Thoracic Surgery, The Second Affiliated Hospital, Hengyang Medical School, University of South China, Hengyang, China
Introduction: The connection between aging and cancer is complex. Previous research has highlighted the association between the aging process of lung adenocarcinoma (LUAD) cells and the immune response, yet there remains a gap in confirming this through single-cell data validation. Here, we aim to develop a novel aging-related prognostic model for LUAD, and verify the alterations in the genome and immune microenvironment linked to cellular senescence.
Methods: We integrated a comprehensive collection of senescence genes from the GenAge and CellAge databases and employed the least absolute shrinkage and selection operator (LASSO) Cox analysis to construct and validate a novel prognostic model for LUAD. This model was then utilized to examine the relationship between aging, tumor somatic mutations, and immune cell infiltration. Additionally, we explored the heterogeneity of senescence and intercellular communication within the LUAD tumor microenvironment (TME) through single-cell transcriptomic data analysis.
Results: By exploring the expression profiles of 586 cellular senescence-related genes in 428 LUAD patients, we constructed an aging-related genes (ARGs) risk model included 10 ARGs and validated it as an independent prognostic predictor for LUAD patients. Notably, patients with low aging scores (LAS group) exhibited better survival, lower tumor mutation burden (TMB), lower somatic mutation frequency, lower tumor proliferation rate, and an immune activated phenotype compared to patients with high aging scores (HAS group). While the HAS group was enriched in tumor cells and showed a lower infiltration of CD8-CCR7, CD8- CXCL13, CD8-GNLY, FCGR3A NK cells, XCL1 NK cells, plasma cell (PC) and other immune subsets. Furthermore, the SPP1 and TENASCIN pathways, associated with tumor immune escape and tumor progression, were also enriched in the HAS group. Additionally, our study also indicated that senescence levels were heterogeneous in the LUAD tumor microenvironment (TME), especially with tumor cells in the LAS group showing higher age scores compared to those in the HAS group.
Conclusions: Collectively, our findings underscore that ARRS through ARGs serves as a robust biomarker for the prognosis in LUAD.
1 Introduction
Cancer with complex molecular characteristics (1), remains a significant global health challenge, accounting for a substantial number of deaths and impacting life expectancy worldwide. Amid the array of cancer types, lung cancer emerges as the second most prevalent contributor to cancer-related mortality, marked by a discouraging 5-year relative survival rate of just 23% (2). Lung adenocarcinoma (LUAD), the predominant histological subtype within non-small cell lung cancer, constitutes over 40% of all lung cancer cases (3). Notably, LUAD continues to rise in incidence among current smokers, former smokers, and even non-smokers, and its five-year survival rate remains dishearteningly low at approximately 15%, as a significant majority of patients are typically diagnosed at advanced stages of the disease (4). Hence, there is still a compelling need to formulate a novel prognostic model for predicting the outcomes of LUAD to advance more potent strategies for diagnosis and treatment.
Aging is a ubiquitous biological process that results in a progressive and irreversible decline in physical function across all organ systems, which presents with genomic instability, telomere attrition, epigenetic alterations, loss of proteostasis, disabled macroautophagy, deregulated nutrient-sensing, mitochondrial dysfunction, stem cell exhaustion, chronic inflammation, altered intercellular communication, cellular senescence, and dysbiosis (5–8). Cellular senescence refers to the essentially irreversible arrest of cell proliferation (growth) that occurs when cells experience potentially oncogenic stress (damage to DNA, strong mitogenic signals, damage or disruptions to the epigenome, and ectopic expression of certain tumor suppressors) (9, 10). Several evidences have shown that cellular senescence plays a double-edged role in initiation, growth, and progression of tumor (11, 12). Senescent tumor cells wield influence over the tumor microenvironment (TME) via the senescence-associated secretory phenotype (SASP). On one hand, by emitting pro-inflammatory cytokines, chemokines, growth factors, and proteases like IL-6, IL-8, and TGF-β, senescent cells can trigger paracrine senescence, transforming neighboring non-senescent cells into senescent counterparts. This process recruits and activates immune cells within the TME, leading to outcomes that can either hinder or foster tumor growth. M1 macrophages and natural killer cells, for instance, can eliminate tumor cells and foster their senescence through the secretion of IFN-γ and TNF-α, thereby restraining tumor expansion. On the other hand, senescent tumor cells may activate myeloid-derived suppressor cells and M2 macrophages via SASP, affecting the clearance of senescent tumor cells, in turn, driving tumor progression and vascularization (9, 13, 14). Given the role of cellular senescence in constraining tumor development, it emerges as a potential target for tumor therapy. Hence, unraveling the impact of senescence in tumorigenesis is paramount importance.
In recent years, several studies have focused on the role of senescence in LUAD (15–20). For example, Lin et al. constructed a cellular senescence-related signature (SRS) by leveraging senescence-related genes. They found that SRS involved in the regulation of the tumor immune microenvironment through SASP was a robust biomarker for the immunotherapeutic response and prognosis in LUAD (15). In a similar vein, another research by Lin et al. explored cellular senescence patterns within LUAD by analyzing mRNA expression profiles of 278 cellular senescence-related genes, demonstrating the association between cellular senescence patterns and tumor immune infiltration in LUAD (16). Besides, Liu et al. developed a 12-gene signature for LUAD using 91 cancer-related senescence genes to assess survival outcome (19). Nonetheless, prior investigations were marked by limitations. Firstly, all focused on only a subset of senescence genes. Secondly, the assessment of the TME was largely confined to the bulk transcriptomic level. As a result, the role of senescence in LUAD has yet to undergo systematic evaluation, and the intricate interplay between senescence and LUAD prognosis has remained obscure.
This current study seeks to overcome these limitations by integrating a comprehensive collection of 586 senescence genes sourced from the GenAge and CellAge databases. Employing the least absolute shrinkage and selection operator (LASSO) Cox analysis, a novel prognostic model for LUAD was constructed and validated. This model was further investigated the relationship between aging and tumor somatic mutation or immune cell infiltration. Furthermore, this study delved into the senescent heterogeneity and intercellular communication of various cells within the LUAD TME through the analysis of single-cell transcriptomic data. In summary, this study enriches our understanding of the profound impact of cell senescence on the survival outcomes of patients with LUAD, which unravels the complex associations between senescence, the immune landscape, and the intricate genetic makeup of the tumor, ultimately illuminating novel avenues for therapeutic interventions and prognostic assessments.
2 Materials and methods
2.1 Data source and processing
In the training cohort, bulk RNA sequencing (RNA-seq) data, somatic mutation data and clinical information for LUAD were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) (21). After excluding non-primary cases and patients with incomplete follow-up information, we analyzed 428 patients from the TCGA dataset as the training set. For the validation cohort (GSE31210, GSE50081, and GSE30219) (22–25), transcriptome data were obtained from data series in the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) (26). Single-cell RNA-seq (scRNA-seq) data (GSE189357) comprising nine patients with LUAD was also download from the GEO database (27). Fragments per kilobase million (FPKMs) values or raw gene expression counts were normalized to transcripts per kilobase million (TPMs) in both the training and validation cohorts. Genes that were not expressed in more than half of the samples were excluded from the expression profiles. The clinical features of 428 patients are listed in Table 1.
2.2 Aging gene set and screening
The set of aging-related marker genes (ARGs) was obtained from two databases, GenAge and CellAge. Initially, 279 ARGs were selected from CellAge (https://genomics.senescence.info/cells/) (28), and an additional 307 ARGs were obtained from GenAge (https://genomics.senescence.info/genes/index.html) (29) (Supplementary Table S1). Univariate Cox analysis was conducted by survival (version 3.3-1) packages to preliminarily identify ARGs associated with the overall survival (OS) of LUAD patients in the TCGA cohort (30), resulting in a final gene set comprising 102 ARGs (Supplementary Table S2).
2.3 Construction and validation of an ARGs risk model
We utilized the “glmnet” (version 4.1-8) package in R software (version 4.1.2) to perform the LASSO Cox regression analyses (family=“cox”) to screen out the prominent genes (31, 32). The “lambda.1se” value, determined through tenfold cross-validation, was employed as the lambda for model fitting. Ten genes were ultimately selected to construct the risk model. The prognostic capability of the ten genes was assessed using Kaplan-Meier survival curves generated with the survminer (version 0.4.9) and survival (version 3.3-1) R packages (30). Subsequently, we calculated a risk score for each sample, as a linear combination of gene expression levels within the signature set, weighted by their respective LASSO Cox regression coefficients, using the following formula:
Here, “Coef (genei)”, signifies the LASSO Cox regression coefficient, “Expre (genei)”, represents the expression level of each gene, and “n” denotes the number of genes included in the model. In addition, the R package “survival” (version 3.3-1) was used to construct multiple multivariate Cox analysis to determine the independent prognostic factor in LUAD patients (30).
In the TCGA training cohort, LUAD patients were classified into high aging score group (HAS group) and low aging score group (LAS group) based on the median value of ARRSs. The prognostic capability of the risk model in terms of OS and progression-free survival (PFS) was assessed using Kaplan-Meier survival curves generated with the survminer (version 0.4.9) and survival (version 3.3-1) R packages (30). Additionally, we also compared the clinicopathological characteristics of TCGA-LUAD patients between the HAS group and the LAS group using Fisher’s Exact Test.
To validate the ARGs Risk Model, we calculated the risk score for patients in the validation cohort (GSE31210, GSE50081, and GSE30219) using the same formula as applied to the TCGA-LUAD cohort. Patients in the validation cohort were also categorized into high and low-risk groups based on the median value of ARRS. Kaplan-Meier curves were generated to assess the relationship between ARRS and OS in the validation cohort.
2.4 Functional enrichment analysis of differentially expressed genes based on HAS and LAS groups
We used the “DESeq2” (version 1.36.0) R package to calculate fold-changes and identify differentially expressed genes (DEGs) based on the two risk groups (false discovery rate (FDR) < 0.05 and |Log2FC| > 1) (33). Subsequently, we conducted Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses on these DEGs using the “clusterProfiler” (version 4.7.1.002) R package (34). Pathways with adjusted p-values less than 0.05 were considered significant.
2.5 Immune infiltration between the HAS-group and LAS-group from TCGA-LUAD cohort
The “estimate” R package, a powerful tool for quantifying the immune stromal, and ESTIMATE scores, which was based on the expression of related molecular biomarkers in immune and stromal cells, to predict the TME (35). The “xCell” is a robust algorithm that analyzes the infiltration levels of 64 immune and stroma cell types, including extracellular matrix cells, epithelial cells, hematopoietic progenitors, innate and adaptive immune cells (36). Herein, we utilized the R package estimate (version 1.0.13) and xCell (version 1.1.0) to evaluate the immune infiltration score and immune cell infiltration in each patient between HAS and LAS subgroups. Additionally, the T cell-inflamed gene expression profile (GEP) was calculated as a weighted sum of standardized expression values of 18 genes (CCL5, CD27, CD274, CD276, CD8A, CMKLR1, CXCL9, CXCR6, HLA-DQA1, HLA-DRB1, HLA-E, IDO1, LAG3, NKG7, PDCD1LG2, PSMB10, STAT1, TIGIT) as described in previous literature (37–39). The single sample gene set enrichment analysis (ssGSEA) algorithm in “gsva” (version 1.42.0) R package was performed to compare differences in 13 gene sets associated immune function and 4 gene sets related to angiogenesis, matrix, matrix remodeling, and tumor proliferation rate from previous studies (40–42). Box plots were developed using ggplot2 software (version 3.4.3) in R to display the differences between the two groups (43).
2.6 The genetic landscapes of HAS-group and LAS-group
Genetic landscapes were analyzed and visualized using the “maftools” (version 2.12.0) R package (44). Tumor Mutation Burden (TMB) was defined as the number of somatic, non-silent, protein-coding mutations in the coding regions per megabase (mut/Mb) and counted using ‘maftools’ (version 2.12.0). The mutated samples of tumor-related and DNA damage repair (DDR) pathways in HAS and LAS groups were compared using Fisher’s exact test (with p <0.05 indicates a significant difference) and visualized using “ggradar” (version 2.12.0) and ggplot2 (version 3.4.3) R packages (43).
2.7 Single-cell RNA-seq analysis
Raw matrix data were obtained from the GEO database for subsequent analysis (27). Initially, cells with low quality were filtered out based on the following criteria: 1) fewer than 200 expressed genes, 2) total molecule count per cell less than 800, and 3) greater than 10% of reads mapped to the mitochondrial genome. Additionally, the “DoubletFinder” R package (45) was utilized to identify and remove doublet cells using default parameters.
The “Seurat” package (version 4.3.0) (46) was employed to normalize the single-cell gene expression data using the “NormalizeData” and “ScaleData” functions, respectively. Subsequently, the top 2,000 highly variable genes for each sample were selected using the “FindVariableFeatures” function. Principal component analysis (PCA) was performed using the “RunPCA” function, and the first 20 principal components were used for Uniform Manifold Approximation and Projection (UMAP) analysis with the “RunUMAP” function. Following UMAP analysis, cells were clustered using an unsupervised method with a resolution parameter of 1.5 employing the “FindNeighbors” function. Differential expression analysis was conducted on the original log-normalized data by comparing cells within each cluster to all other clusters using the “FindAllMarkers” function. Clusters were annotated based on the expression of well-known markers and differentially highly expressed genes.
Subgroup analysis of each cell group, including T/NK cells, B cells, and myeloid cells, was performed using the standard Seurat pipeline. Specific markers were used for grouping and are listed in Supplementary Table S3. Bar plots were generated to illustrate the percentage of cells between the two groups. Additionally, cell occupancy differences were assessed using Fisher’s exact test. The cytotoxic and exhausted scores for T cell subgroups, as well as the hallmark pathways compared between HAS and LAS groups, were calculated using the ssGSEA algorithm in the “gsva” package (version 1.42.0) based on different sets of genes (42).
2.8 Identification of cancer cells
To identify cancer cells, we utilized the inferCNV (version 1.13.0) tool (https://github.com/broadinstitute/inferCNV), as previously described in studies by Liu, He, et al. and Chen et al. (47, 48). The inferCNV package compares gene expression profiles of each cell to reference gene expression profiles from other cells. Initially, raw count data and cell type annotations for all cells were extracted from the Seurat object. Immune cells and stromal cells were chosen as reference cells. A gene ordering file was generated from the human GRCh38 assembly, containing chromosomal start and end positions for each gene. These files were used to create an inferCNV object using the “CreateInfercnvObject” function, followed by running inferCNV with default parameters. The calculated copy number variation (CNV) signal was defined as the mean square of CNV estimates across all genomic locations. CNV R-scores were calculated as the Pearson correlation coefficient between each cell’s CNV pattern and the average CNV pattern of the top 5% of cells from the same tumor based on CNV signal. Cells with CNV R-scores ≥0.3 were classified as tumor cells.
2.9 Aging-related risk scores based on pseudo-bulks
The Seurat object was transformed into a “SingleCellExperiment” object, followed by the computation of pseudo-bulks. Pseudo-bulks, which represent the sum of counts, were calculated using aggregation-based methods in the muscat (version 1.10.1) R package (https://github.com/HelenaLC/muscat). The ARRSs were then derived using the previously described formula based on the pseudo-bulks. Patients were stratified into two groups, HAS and LAS, based on the median value of ARRSs. Additionally, age scores for each cell were calculated based on ten ARKGs at the single-cell level using the ssGSEA algorithm.
2.10 Cell-cell interactions
CellChat (version 1.5.0) is an open-source R package (https://github.com/sqjin/CellChat) utilized for the analysis, comparison, and visualization of single-cell RNA sequencing data intercellular communication (49). In this study, CellChat was employed to infer cell-cell interactions across 24 immune subgroups, fibroblasts, normal epithelial cells, tumor cells, and endothelial cells for both the HAS and LAS groups. Subsequently, major signaling changes between the HAS and LAS groups were computed.
2.11 Statistical analysis
The Wilcoxon test was conducted to examine differences in variables between two groups, while the Kruskal-Wallis test was used to assess differences among groups greater than two. Gene mutation differences between the HAS and LAS groups were determined using Fisher’s exact test.
3 Results
3.1 Construction and validation of aging-related risk score
The workflow of the whole study was graphically presented in Figure 1A. We compiled a comprehensive list of 586 aging-associated genes sourced from the CellAge and GenAge databases. Among these genes, 102 were significantly associated with clinical survival (p < 0.05) based on univariate Cox analysis (detailed results shown in Supplementary Table S2), conducted on the expression matrix and clinical survival information of 428 LUAD samples obtained from the TCGA dataset. Subsequently, to construct the ARGs risk model, we performed LASSO Cox regression analysis on the aforementioned 102 genes and the gene expression profiles of the training cohort (Figures 1B, C). Through this analysis, we successfully identified 10 aging-related key genes (ARKGs), including. BRCA2, CSNK1E, EEF1E1, GAPDH, IGFBP3, IL1A, PSEN1, XRCC5, XRCC6, and YWHAZ. And low RNA expression for the 10 ARKGS was correlated with longer survival time in LUAD (Supplementary Figure S1). Utilizing these ten ARKGs and their corresponding risk coefficients, we established an aging risk signature. The risk score of every patient was calculated using this formula. Patients in the training cohort were stratified into two groups: the high aging score group (HAS group) and the low aging score group (LAS group) based on median values of ARRSs. Upon investigating the expression levels of the ten ARKGs, we found that they were significantly higher in HAS group than LAS group (Supplementary Figure S2, Supplementary Table S4).
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1347770/fimmu-15-1347770-HTML/image_m/fimmu-15-1347770-g001.jpg)
Figure 1. Identification of ARKGs related to prognosis in the TCGA training cohort. (A) The workflow of the present study. (B) Selection of optimal candidate genes in the LASSO model. (C) LASSO coefficients of prognosis-associated ARKGs. (D, E) Kaplan-Meier curves for overall survival (D) and progression free survival (E) of the TCGA-LUAD cohort in the HAS and LAS groups. (F) Forest plots showing results of multivariate Cox regression analysis between Risk score, clinical information and overall survival. (G-I) Kaplan-Meier curves for overall survival of validation cohorts in the high and low groups: GSE50081 (G), GSE30219 (H), GSE31210 (I).
We compared the clinicopathological characteristics, including age, gender, race, OS, TNM tumor grade, and AJCC tumor grade, of TCGA-LUAD patients between the HAS group and the LAS group (Table 1). The results showed significant differences in OS status (P = 0.002), T grade (P = 0.002), N grade (P = 0.001), and AJCC tumor grade (P = 0.004) between the groups. Survival analysis demonstrated that the LAS-group exhibited significantly improved overall survival (OS) (hazard ratio (HR) = 0.45, 95% confidence interval (CI) = 0.31 – 0.66, P = 0.000044) and progression-free survival (PFS) (HR =55, 95% CI = 0.41 – 0.74, P = 0.000052) compared to the HAS-group (Figures 1D, E). Upon integrating age, gender, TNM tumor grade, and AJCC tumor grade into the multivariate Cox regression analysis, ARRSs emerged as the sole significant survival-related risk factor (HR = 9.32, 95% CI = 4.50 – 19.29, P = 0.0000000018) (Figure 1F), suggesting that ARRSs was an independent prognostic factor for LUAD.
To validate the prognostic roles of the above risk model, we applied the same stratification method to three independent datasets form the GEO database. Consistent with the findings from the training cohort, patients with in the high ARRSs group displayed significantly worse survival outcomes compared with the low ARRSs group in all three cohorts, namely GSE50081 (HR = 0.32, 95% CI = 0.18 – 0.56, P = 0.000081), GSE30219 (HR = 0.52, 95% CI = 0.29 – 0.95, P = 0.038), and GSE31210 (HR = 0.36, 95% CI = 0.19 – 0.71, P = 0.005) (Figures 1G-I).
3.2 The genetic characteristics of HAS-group and LAS-group
To explore the genetic features in LUAD with different ARRSs, we further investigated the genomic differences between the HAS group and the LAS group based on somatic mutation data in the TCGA-LUAD cohort (Figure 2A; Supplementary Figures S3A-C). We observed that HAS group had a higher mutation frequency than the LAS group, particularly in the top 20 genes such as, TP53, TTN, CSMD3, ZFHX4, RYR3, CSMD2, SI, LRRC7, and PAPPA2 (detailed P values shown in Supplementary Table S5) between HAS and LAS groups (Figure 2B). Additionally, the HAS group displayed a higher tumor mutation burden (TMB) but a lower occurrence of co-occurring mutations between genes, indicating distinct genomic alteration patterns (Figure 2C; Supplementary Figure S3D). Further analysis of ten tumor-related pathways revealed significantly higher mutation frequencies in the Hippo (P = 0.011), NOTCH (P = 0.013), and TP53 (P = 0.011) pathways in the HAS-group compared to the LAS-group (Figures 2D, F; Supplementary Figure S4B). Similarly, higher mutation rates were observed in the HAS group among the eight DDR pathways, with five of them being statistically significant (Figures 2E, G; Supplementary Figure S4A).
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1347770/fimmu-15-1347770-HTML/image_m/fimmu-15-1347770-g002.jpg)
Figure 2. Genomic alterations differences between the HAS and LAS group from the TCGA-LUAD cohort. (A) Genomic alterations landscape between the HAS (left) and LAS (right) group. (B) Mutation frequency differences of the top 20 mutation genes in the HAS group compared to the LAS group. The asterisk to the right of the gene indicates that the mutations in the gene were significantly different in the two groups, as determined by Fisher’s exact test. (C) The TMB between HAS and LAS groups. The HAS group had a higher TMB (2.66 (IQR: 0.04, 11.985) compared to the LAS group (1.58 (IQR: 0.02, 6.9)) with P value = 0.000065 compared by the Wilcoxon test. The frequency of mutated genes in each tumor-related pathway (D, F) and DDR pathway (E, G) difference between two groups. The asterisks in (D, E) denote significant differences of mutated genes in different pathways identified by Fisher’s exact test which showed in (F, G).
3.3 ARRSs is associated with cell proliferation and immune function
Differential expression analysis of gene expression data based on the HAS group and LAS group identified a total of 1664 differentially expressed genes (DEGs) under a threshold of adjusted p < 0.05, comprising 707 up-regulated and 957 down-regulated genes (Figure 3A). GO enrichment analysis for DEGs revealed that in the HAS-group, biological processes were predominantly enriched in cell cycle, cell division, and cell development, indicating a potential involvement in regulating normal cell function and organismal development (Figure 3B). Furthermore, based on gene sets from Bagaev, et al. (40), we found that the tumor proliferation rate, and matrix remodeling of the HAS group were significantly higher than those of the LAS group (Figure 3F, detailed P values were shown in the Supplementary Table S4).
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1347770/fimmu-15-1347770-HTML/image_m/fimmu-15-1347770-g003.jpg)
Figure 3. Transcriptomic differences between the HAS and LAS group from the TCGA-LUAD cohort. (A) Volcano Plot of DEGs between the HAS and LAS group. (B, C) Top 20 biological processes of GO enrichment results between the HAS (B) and LAS (C) group. (D) Stromal score, immune score and ESTIMATE score between the two groups. (E) GEP score between the two groups. (F, G) Boxplots of gene sets related to tumor proliferation (F) and immune-related functions (G). (H, I) Box plot of T cells (H), B cells (H), and myeloid cells (I) infiltration in “Xcell” between the two groups. (J) Box plot of cell death between the two groups. "ns" indicates p > 0.05, * indicates p ≤ 0.05, ** indicates p ≤ 0.01, *** indicates p ≤ 0.001, and **** indicates p ≤ 0.0001. The actual P determined by the Wilcoxon test, and the medians (IQR) in Figures 2D-F were all displayed in Supplementary Table S4. All abbreviations presented in Figure 3 showed as following: GEP, T cell-inflamed gene expression profile; CCR, cytokine and cytokine receptor; HLA, human leukocyte antigen; MHC, major histocompatibility complex.
In contrast, the LAS-group exhibited enrichment in immune response mechanisms, encompassing cell activation, signal transduction, and production of immune mediators (Figure 3C). Based on another gene set related to immune functions (41), we observed that nine of the 13 immune function gene sets had significantly higher ssGSEA scores in all LAS groups than the HAS group (Figure 3G, detailed P values were shown in the Supplementary Table S4), especially type II IFN response, T cell co-stimulation, and HLA. Immune estimations for LUAD patients within the training set (TCGA-LUAD) showed notably increased StromalScore, ImmuneScore, ESTIMATEScore, and GEP score in the LAS group when compared to the HAS group (Figures 3D, E). Xcell analysis revealed the immune infiltration of TME (36). The results indicated that LAS group had an activated TME, with significantly increased numbers of T cells, such as CD8+ T cells, CD8+ Tcm, CD4+ Tem, and CD4+ Tcm, and significantly decreased numbers of Th1 and Th2 (Figure 3H, detailed p values were shown in the Supplementary Table S4). Additionally, B cells such as plasma cells (Figure 3H), and myeloid cells such as Mast cells, and various DCs (Figure 3I, detailed p values were shown in the Supplementary Table S4), were also significantly increased in the LAS group. Furthermore, we explored the relationship between ARRSs and various cell death pathways. The findings revealed that significantly elevated scores for Alkaliptosis, Cuproptosis, and Oxeiptosis in the HAS-group, whereas Autophagy, Lysosome-dependent cell death, Necroptosis, and Parthanatos scores were markedly higher in the LAS-group (Figure 3J; Supplementary Table S4).
3.4 The single cell alta of HAS-group and LAS-group
To further investigate whether the ARRSs is heterogeneous in the TME, we utilized a single-cell dataset (GSE189357) containing over 10,000 cells from 9 patients. Initially, the single-cell dataset was converted to pseudo-bulks and then ARRSs were calculated. Subsequently, the 9 patients were divided into HAS (n = 5) and LAS (n = 4) groups based on the median value of ARRSs. Notably, two of the three invasive adenocarcinoma (IAC) samples were categorized into the HAS group, exhibiting significantly higher aging scores compared to the LAS group (Figures 4A, B). Employing the standard pipeline in Seurat (46), we identified six major cell types, including T/NK cells, B cells, myeloid cells, fibroblasts, endothelial cells, and epithelial cells (Figures 4C, D). Subsequently, the epithelial cells were further subdivided into tumor cells and normal epithelial cells (Figure 4E). Interestingly, we observed an enrichment of tumor cells and endothelial cells in the HAS (P = 6.44E-66, odds ratio (95% CI) = 1.57 (1.49, 1.66), Supplementary Table S6) and LAS (p = 0, odds ratio (95% CI) = 4.56 (4.2, 4.95)) groups, respectively (Figure 4F).
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1347770/fimmu-15-1347770-HTML/image_m/fimmu-15-1347770-g004.jpg)
Figure 4. Single cell atlas and cell-cell interactions between the HAS (n = 5) and LAS (n = 4) group. (A) ARRSs based on single cell pseudo-bulks differed between the HAS and LAS group. The HAS group had a higher ARRS 0.4496 (IQR: 0.4304, 0.463) compared to the LAS group 0.4147 (IQR: 0.3969, 0.4228) with P value = 0.0159 compared by the Wilcoxon test. (B) Alluvial diagram showed the grouping of the nine samples. (C) UMAP plot for cells displaying the six major cell types from patients. (D) Dot plot depicting mean expression levels and percentage of cells expressing signature genes across the six major cell types. (E) Distribution of normal and tumor cells in epithelial cells from LUAD. (F) The composition of the cell compartment, displaying the average frequencies of cell subsets for HAS and LSA groups. (G-I) The UMAP plot and the average frequencies of different T cell, B cell and myeloid cell subgroups. (J, K) Comparison of the significant ligand-receptor pairs of SPP1 signaling (J) and TENASCIN signaling (K) for the HAS group. Dot color reflects communication probabilities and dot size represents computed p-values. Empty space means the communication probability is zero. p-values are computed from one-sided permutation test. (L, M) Circle plot showed cell–cell communication mediated by CD70-CD27 (L) and BAG6-NCR3 (M) in the LAS group. All abbreviations presented in Figure 4 showed as following: ARRS, aging related risk score; IQR, interquartile range; AIS, lung adenocarcinoma in situ; MIA, minimally invasive adenocarcinoma; IAC, invasive adenocarcinoma; tumor, tumor cells; normal, normal epithelial cells; Fib, fibroblasts; End, endothelial cells.
Furthermore, we conducted subtype annotation specifically for immune cells including T/NK cells, B cells, and myeloid cells (Figures 4G-I; Supplementary Table S2). T/NK cells were subdivided into eight T cell subpopulations and two NK cell subpopulations (Figure 4G). Functional scoring of T-cell subsets revealed that FCGR3A NK cells (T09) and CD8-GNLY (T08) had the highest cytotoxic scores, while CD8-CXCL13 (T06) had the highest exhausted score (Supplementary Figure S5). We compared the cellular infiltration in the HAS and LAS groups and found that the T and NK cell subpopulations were significantly differed between the HAS and LAS groups (Supplementary Table S6). Specifically, CD4-CCR7 (T01, P = 1.53E-134, odds ratio (95% CI) = 1.81 (1.72, 1.9)), and CD4-FOXP3 (T03, P = 4.55E-21, odds ratio (95% CI) = 1.41 (1.31, 1.52)) were enriched in the HAS group, whereas CD8-CCR7 (T05, P = 6.42E-14, odds ratio (95% CI) = 1.34 (1.24, 1.45)), CD8-CXCL13 (T06, P = 2.35E-64, odds ratio (95% CI) = 5.66 (4.51, 7.16)), CD8-GNLY (T08, P = 2.94E-23, odds ratio (95% CI) = 1.39 (1.3, 1.48)), FCGR3A NK cells (T09, P = 4.29E-98, odds ratio (95% CI) = 1.95 (1.83, 2.08)), and XCL1 NK cells (T10, P = 2.65E-34, odds ratio (95% CI) = 2.02 (1.8, 2.27)) were enriched in the LAS group. For B cell subsets, naive and memory B cells were more prevalent in the HAS group, whereas plasma cell (PC) subsets (B03 P = 1.33E-24, odds ratio (95% CI) = 2.3 (1.96, 2.69); B04, P = 3.51E-48, odds ratio (95% CI) = 2.96 (2.55, 3.43)) and stressed PC (B05, P = 7.24E-13, odds ratio (95% CI) = 3.02 (2.21, 4.12)) were more prevalent in the LAS group. The Mast cells (M01, P = 7.44E-292, odds ratio (95% CI) = 2.77 (2.62, 2.92)) showed a tendency to increase in the LAS group compared to the HAS group, while neutrophils (M02, P = 4.84E-164, odds ratio (95% CI) = 5.25 (4.55, 6.08)), S100B DC (M06, P = 1.08E-66, odds ratio (95% CI) = 2.05 (1.88, 2.23)), TXN DC (M07, P = 9.61E-15, odds ratio (95% CI) = 2.47 (1.93, 3.19)), and proliferation myeloid cells (M09, P = 1.51E-16, odds ratio (95% CI) = 1.87 (1.6, 2.2)) were significantly more prevalent in the HAS group. These results provide further evidence of heterogeneity in immune cell infiltration between groups with differing ARRs at the single-cell level, especially the LAS enriched more cytotoxic T/NK cells and antibody-secreting B cells.
3.5 Inference of cell-cell interactions
Given that senescence alters intercellular communication, we conducted a comparative analysis of intercellular communication between the HAS and LAS groups for each cell subset based on single-cell data. Significant differences were observed in several signaling networks between the HAS and LAS group (Supplementary Figure S6). Notably, SPP1 was exclusively present in in the HAS group (Supplementary Figure S6; Figure 4J). Especially, the interaction of SPP1-CD44 has been reported to inhibit T-cell activation and promote tumor immune evasion (50, 51). Additionally, TENASCIN was frequently observed in the HAS group, with tumor cells in this group interacting with other cells, including tumor cells themselves, via TNC - SDC1/SDC4 or TNC - ITGA8_ITGB1/ITGAV_ITGB6 (Supplementary Figure S6; Figure 4K). TNC is an extracellular matrix glycoprotein known to contribute to tumor progression, and increased TNC expression in LUAD tissues correlates with an unfavorable clinical outcome for patients (52). Conversely, certain pathways were exclusively or more frequently observed in the LAS group (Supplementary Figure S6). For example, the secreted signaling BAG, and CD70 pathways were uniquely found in the LAS group (Supplementary Figure S6). The BAG6-NCR3 interaction targeting T09 might trigger NK cell cytotoxicity (Figure 4L). Furthermore CD70-CD27 interaction was observed between B02 and PC or between B02 and T cells. CD27 receptor activation provides a costimulatory signal promoting T cell and B cell activity to enhance anti-tumor and anti-infection immunity (Figure 4M) (53).
3.6 Cellular senescence heterogeneity in the tumor microenvironment
Using single-cell data, we evaluated the senescence levels of individual cells and compared the senescence levels among different cell subpopulations (Figure 5A). We observed lower age scores in B01, B03, B05, M02, and endothelial cell subpopulations, while M03 and M05 exhibited higher age scores (Figure 5A). Subsequently, we compared the senescence levels of cell subpopulations between the HAS and LAS groups (Figure 5B). Most T cell subsets (e.g., T05, T07) displayed higher age scores in the HAS group than in the LAS group (Figures 5B, C). Moreover, endothelial and fibroblast cells exhibited higher age scores in the HAS group, whereas normal epithelial cells and tumor cells showed higher age scores in the LAS group (Figures 5B, C). Age scores for different subpopulations of myeloid and B cells varied between the HAS and LAS groups (Figures 5B, C). For instance, the age scores of B01, B02, M05 and M06 were significantly lower in the HAS group than in the LAS group, while B04 and M02 showed higher scores in the HAS group (Figures 5B, C). As cellular damage caused by reactive oxygen species (ROS) is a major trigger for senescence (54), we assessed and compared the “reactive oxygen species pathway”. Our results revealed higher scores for this pathway in the HAS group for B01, B02, M05, M06, normal epithelial cells, and tumor cells, whereas the HAS group for B04, M02, T05, T07, endothelial, and fibroblast cells exhibited lower scores (Figure 5D), consistent with the trend observed in age scores (Figure 5C).
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1347770/fimmu-15-1347770-HTML/image_m/fimmu-15-1347770-g005.jpg)
Figure 5. Single cell age score. (A) The age score for each cell. (B) The medians of age scores for each cell type in the HAS and LAS groups. (C) The box plots demonstrating between-group differences in the HAS and LAS groups for age scores for specific cell types. We used the following notation for statistical significance: “ns” indicates p > 0.05, ** indicates p ≤ 0.01, *** indicates p ≤ 0.001, and **** indicates p ≤ 0.0001. The actual P determined by the Wilcoxon test, and the medians (IQR) in (C) were all displayed in Supplementary Table S4. (D) Heatmap showed the activity of hallmarks between the HAS and the LAS groups for different cell types. All abbreviations presented in Figure 5 showed as following: tumor, tumor cells; normal, normal epithelial cells; Fib, fibroblasts; End, endothelial cells.
4 Discussion
Cellular senescence involves the cessation of cell-cycle and the release of inflammatory cytokines with autocrine, paracrine and endocrine activities (55). The SASP represents a significant feature of senescent cells, encompassing the release of various cytokines, chemokines, growth factors and proteases (56). The impact of cellular senescence on cancer is intricate, displaying both advantageous and detrimental effects. Nevertheless, the extent to which the senescent heterogeneity of immune infiltration cells within tumors, as well as the interplay between tumor senescence and immune infiltration in LUAD, remains incompletely documented. In the current investigation, we constructed an independent prognostic model based on cellular senescence-related genes, and comprehensively analyzed the role of aging in genomic alterations and immune landscape in LUAD, which might hold the potential to facilitate the development of personalized immunotherapy.
This study successfully identified a novel and independent prognostic risk model incorporating ten significantly upregulated genes in LUAD. Ten genes were selected from a comprehensive list of 586 aging-associated genes obtained from the CellAge and GenAge databases. These genes also have been previously reported as positive regulators of tumor development. For example, CSNK1E, a member of the serine/threonine protein kinase family, controls circadian rhythms, which is closely related to the animals longevity (57). Inhibition of CSNK1E has been show to selectively inhibit tumor cell development (58), and elevated CSNK1E expression is associated with poor prognosis in patients with ovarian cancer and malignant melanoma (59, 60). EEF1E1, a tumor suppressor, plays a role in ATM/ATR-mediated p53 activation (61), and serves as a poor prognosis predictor in lung cancer (62). Overexpression of EEF1E1 in transgenic mice resulted in a significantly shorter mean lifespan (63). GAPDH directly participates in tumor progression, invasiveness, and metastasis (64), and conditions such as oxidative stress impair GAPDH catalytic activity, leading to cellular aging and apoptosis (65). Increased expression of PSEN1 in colorectal cancer is associated with enhanced tumor development through heightened EGFR signaling via NOTCH1 processing and activation of the COX-2-PGE2 pathway (66). PSEN1-null mice die shortly after birth (67), although PSEN1’s role in human aging remains largely unknown. YWHAZ is an adapter protein implicated in several signal transduction pathways (68) and interacts with numerous proteins associated with aging, such as the INS/IGF1 pathway (69, 70). YWHAZ has also been shown to mediate lung cancer malignancy and β-catenin protein through its complex with β-catenin (71). IL1A, a pivotal inflammatory cytokine, is thought to be one of the critical upstream regulators of other SASP-related genes (72, 73) and drives tumor growth and metastasis (74). IGFBP3, a member of the insulin-like growth factor-binding protein (IGFBP) family, regulates IGF1 and IGF2 by altering the interaction of IGFs with their cell surface receptors. Interestingly, the cell growth regulator IGFBP3 exhibits a unique pattern, as elevated levels are associated with a good prognosis in patients with advanced NSCLC (75). BRCA2, XRCC5, and XRCC6 are all DDR related genes, involved in DNA damage and repair. Mice deficient for BRCA2 and XRCC5 have a reduced lifespan (76, 77). XRCC5/6 are associated with poor prognosis and can be used as diagnostic and prognostic biomarkers for LUAD (78). BRCA2’s role in cancer well-established, as elevated BRCA2 expression is associated with a significantly reduced number of stromal cells and high infiltration of both beneficial and detrimental immune cells in breast cancer (79). BRCA2 has also been demonstrated to exhibit increased mRNA levels and poor prognosis in lung cancer (80). These findings collectively provide compelling evidence that this newly proposed prognostic risk model has the potential to reflect LUAD prognosis by considering genomic alterations and the immune landscape.
Genetic instability is a common characteristic of both aging and cancer (81), encompassing changes in DNA damage, DNA damage response and repair, mutations, replication stress, transposition, chromosome aberrations, telomere shortening, micronuclei, and DNA fragments (82). In our study, we found that the HAS group exhibited more frequent gene mutations and higher TMB, indication the presence of an unstable genome and immunogenic potential in patients with HAS. Furthermore, the mutation frequency of the Hippo, NOTCH, TP53, and DDR pathways in the HAS group were also significantly increased. Hippo is an important pathway regulating differentiation, stem cell renewal, and oncogenic transformation (83). In cancer research, the activated Hippo pathway is considered as a tumor suppressor pathway due to its ability to impede cell proliferation and facilitate apoptosis (84). Similarly, NOTCH (85) and TP53 (86) pathway mutations have also been reported to associate with unfavorable prognosis in lung cancer. DNA damage response plays a significant role in maintaining genomic integrity and closely associated with lung cancer progression and treatment (87, 88). These researches provide additional insights into our observed outcomes that patients with HAS experience poorer survival when compared to those with LAS patients.
Cellular senescence functions as a stress response characterized by a halt in proliferation and heightened secretion of pro-inflammatory cytokines (89). Senescent cells recruit immune cells, facilitating their own immune clearance, thereby restoring tissue homeostasis. In the context of cancer, various stressors such as oncogenic signaling, replication stress, hypoxia, reactive oxygen species, nutrient deprivation, and exposure to cytokines within the tumor microenvironment can trigger senescence. This underscores the significant link between tumor cell senescence and immune cell infiltration. Through a bulk-transcriptome analysis, we observed that senescence-associated genes exert a strong influence on the immune microenvironment in LUAD. Specifically, the LAS group showed an activated TME, this manifested as a noteworthy increase in the quantities of CD8+ T cells, CD8+ Tcm, CD4+ Tem, CD4+ Tcm, plasma cells, mast cells and DC, alongside heightened ImmuneScore, GEP score and type II IFN response, T cell co-stimulation, and HLA scores, in addition to enriched immune response pathways. These findings were further corroborated though single-cell analysis, which revealed that CD8-CCR7 (T05), CD8-CXCL13 (T06), CD8-GNLY (T08), FCGR3A NK cells (T09), XCL1 NK cells (T10), plasma cell sets (B03, B04, B05), and mast cells (M01) were more enriched in the LAS group (Figure 4). In contrast, the HAS group displayed an immunosuppressive microenvironment with lower immune function scores and a higher tumor proliferation rate (Figure 4). Additionally, based on the cellular communication results, we identified some signaling pathways specific to the HAS group, such as SPP1 and TENASCIN (Figure 4), which contribute to tumor immune escape and tumor progression (50–52). These results suggest that the HAS group might promote tumor cell invasion by evading immune surveillance, enhancing proliferation and immune escape, leading to poor prognosis in LUAD.
In addition to bulk-level senescence assessment, we also compared senescence at the single-cell level and found significant heterogeneity in cellular senescence. Interestingly, we found that the age scores for tumor cells in the HAS group were significantly lower than that in the LAS group (Figure 5C), suggesting that senescence at the bulk-level is not the same as senescence at the cellular level. Senescent tumor cells might augment the immune response against tumors (90), which is entirely consistent with the highly senescent tumor cells and activated immune microenvironment in the LAS group. However, it’s worth noting that these senescent cells could also reinforce the tumor’s resistance to immunotherapy through potent immunosuppressive mechanisms (91, 92). Therefore, more in-depth studies at the cellular level remain essential.
Herein, we also explore the relationship between senescence and other modes of cell death. Patients in the HAS-group demonstrated a propensity for Alkaliptosis and ROS cell death mechanisms such as Oxeiptosis (93) and Cuproptosis. These endogenous damages, coupled with certain exogenous factors, induced a wide array of genetic injuries, including point mutations and deletions (94), ultimately leading to significantly higher TMB in the HAS-group compared to the LAS-group. To counteract DNA damage, the HAS-group employed a series of intricate DNA repair and maintenance mechanisms associated with cell proliferation and differentiation, ensuring the preservation of proper chromosomal structure and stability (8, 94). Conversely, the interactions among lysosome-dependent cell death, autophagy, and apoptosis played a more significant role in the LAS-group. Meanwhile, the LAS-group exhibited immunological functions in response to cellular senescence, engaging in tissue repair through immune cell recruitment and immune clearance of senescent cells.
More novel analyses were added to our study, although studies related to senescence in LUAD have been reported (15–18, 20). Firstly, although previous studies have also compared differences between aging subgroups in terms of mutations, or TMB (15–18). Patients with higher risk scores had noticeably increased TMB and mutated more frequently for TP53 (15, 16, 18), which is consistent with the results we found. Furthermore, our study was the first to compare at the pathway level which showed significant differences in patients with different ARRs. Second, existing researches related to senescence in LUAD have found that the lower risk scores group embodies an immune-activated microenvironment. Lin, et al., 2023 showed that the ASRS was positively correlated with most immunomodulator-related mRNAs, including chemokines, and immune inhibitors, and receptors (18). This study collected a previously reported set of 13 immune-related gene sets (41) and comprehensively compared the immunity of different subgroups. We found that nine of the 13 immune function gene sets were positively correlated with ARRS score, including APC to stimulation, cytokine and cytokine receptor (CCR), Check-point, cytolytic activity, inflammation-promoting, HLA, T cell co-stimulation, T cell co-stimulation, and type II IFN response (Figure 3). Thirdly, previous studies based on different datasets and different methods have been performed to show the association between immune infiltration and senescence. However, sometimes inconsistent results were obtained by different software. Our study evaluates the association between immune infiltration and senescence for the first time at the single cell level, and using scRNA-seq, this study compared cellular communication between different senescence groups, revealing possible alterations in cellular communication caused by senescence (Figure 4). Finally, we assessed senescence at the cellular level for the first time and found significant inter-cellular heterogeneity for senescence. In particular, we found an opposite trend between the overall senescence score and the tumor cell senescence score. This study still had some limitations, the limited availability of single-cell samples and immune cohort samples may introduce some bias in our model validation. Although we validated the aging score model using several external independent public datasets, prospective clinical trials verification of our model is still necessary. Nevertheless, we hope that this model can contribute to the comprehension of the molecular mechanisms of cellular senescence and TME in LUAD.
In conclusion, our study identified and validated a senescence-related signature based on 10 senescence-related genes as an independent prognostic significance for patients with LUAD, indicating that the senescence levels are heterogeneous in LUAD immune microenvironment, and the HAS group might promote tumor cell invasion by evading immune surveillance, enhancing proliferation and immune escape, leading to poor prognosis in LUAD.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article.
Author contributions
KR: Writing – original draft, Conceptualization. LC: Writing – original draft, Conceptualization, Data curation. CW: Investigation, Methodology, Software, Writing – original draft. XT: Data curation, Investigation, Methodology, Software, Writing – original draft. WA: Writing – review & editing, Methodology. QW: Writing – review & editing. YM: Investigation, Methodology, Software, Writing – original draft. YH: Project administration, Supervision, Writing – original draft. XX: Writing – original draft. JB: Project administration, Supervision, Writing – original draft. XL: Writing – original draft. XFX: Supervision, Writing – original draft. MZ: Supervision, Writing – original draft.
Funding
The author(s) declare financial support was received for the research, authorship, and/or publication of this article. National Natural Science Foundation of China (No. 82071035,82371165). Natural Science Foundation of Shandong Province (No. ZR2022LZL001). Hunan Provincial Department of Science and Technology Clinical Medical Technology Innovation Guidance Project Foundation (No. 2021SK51711).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2024.1347770/full#supplementary-material
Abbreviations
ARGs, aging-related marker genes; ARKGs, aging-related key genes; ARRSs, aging-related risk scores; CI, confidence interval; CNVs, copy number variations; DDR, DNA damage repair; DEGs, differentially expressed genes; HR, hazard ratio; IQR, interquartile range; LUAD, lung adenocarcinoma; SASP, senescence-associated secretory phenotype; SRS, senescence-related signature; TMB, tumor mutation burden; TME, tumor microenvironment.
References
1. Hanahan D. Hallmarks of cancer: new dimensions. Cancer Discovery. (2022) 12:31–46. doi: 10.1158/2159-8290.CD-21-1059
2. Siegel RL, Miller KD, Wagle NS, Jemal A. Cancer statistics, 2023. CA: Cancer J Clin. (2023) 73:17–48. doi: 10.3322/caac.21763
3. Duma N, Santana-Davila R, Molina JR. Non-small cell lung cancer: epidemiology, screening, diagnosis, and treatment. Mayo Clinic Proc. (2019) 94:1623–40. doi: 10.1016/j.mayocp.2019.01.013
4. Spella M, Stathopoulos GT. Immune resistance in lung adenocarcinoma. Cancers. (2021) 13:384. doi: 10.3390/cancers13030384
5. Guo J, Huang X, Dou L, Yan M, Shen T, Tang W, et al. Aging and aging-related diseases: from molecular mechanisms to interventions and treatments. Signal Transduct Target Ther. (2022) 7:391. doi: 10.1038/s41392-022-01251-0
6. Mittelbrunn M, Kroemer G. Hallmarks of T cell aging. Nat Immunol. (2021) 22:687–98. doi: 10.1038/s41590-021-00927-z
7. Cho SJ, Stout-Delgado HW. Aging and lung disease. Annu Rev Physiol. (2020) 82:433–59. doi: 10.1146/annurev-physiol-021119-034610
8. López-Otín C, Blasco MA, Partridge L, Serrano M, Kroemer G. Hallmarks of aging: An expanding universe. Cell. (2023) 186:243–78. doi: 10.1016/j.cell.2022.11.001
9. Campisi J. Aging, cellular senescence, and cancer. Annu Rev Physiol. (2013) 75:685–705. doi: 10.1146/annurev-physiol-030212-183653
10. Campisi J, d’Adda di Fagagna F. Cellular senescence: when bad things happen to good cells. Nat Rev Mol Cell Biol. (2007) 8:729–40. doi: 10.1038/nrm2233
11. Schmitt CA, Wang B, Demaria M. Senescence and cancer - role and therapeutic opportunities. Nat Rev Clin Oncol. (2022) 19:619–36. doi: 10.1038/s41571-022-00668-4
12. Anczuków O, Airhart S, Chuang JH, Coussens LM, Kuchel GA, Korstanje R, et al. Challenges and opportunities for modeling aging and cancer. Cancer Cell. (2023) 41:641–5. doi: 10.1016/j.ccell.2023.03.006
13. Hernandez-Segura A, Nehme J, Demaria M. Hallmarks of cellular senescence. Trends Cell Biol. (2018) 28:436–53. doi: 10.1016/j.tcb.2018.02.001
14. Calcinotto A, Kohli J, Zagato E, Pellegrini L, Demaria M, Alimonti A. Cellular senescence: aging, cancer, and injury. Physiol Rev. (2019) 99:1047–78. doi: 10.1152/physrev.00020.2018
15. Lin W, Wang X, Wang Z, Shao F, Yang Y, Cao Z, et al. Comprehensive analysis uncovers prognostic and immunogenic characteristics of cellular senescence for lung adenocarcinoma. Front Cell Dev Biol. (2021) 9:780461. doi: 10.3389/fcell.2021.780461
16. Lin W, Wang X, Xu Z, Wang Z, Liu T, Cao Z, et al. Identification and validation of cellular senescence patterns to predict clinical outcomes and immunotherapeutic responses in lung adenocarcinoma. Cancer Cell Int. (2021) 21:652. doi: 10.1186/s12935-021-02358-0
17. Zhang W, Li Y, Lyu J, Shi F, Kong Y, Sheng C, et al. An aging-related signature predicts favorable outcome and immunogenicity in lung adenocarcinoma. Cancer Sci. (2022) 113:891–903. doi: 10.1111/cas.15254
18. Lin T, Wang H, Liu Y, Zhao F, He X. Gene set variation analysis-based aging and senescence score as a prognostic indicator and therapeutic guide in lung adenocarcinoma. Front Genet. (2023) 14:1176292. doi: 10.3389/fgene.2023.1176292
19. Liu X, Lin L, Cai Q, Sheng H, Zeng R, Zhao Y, et al. Construction and validation of a prognostic model based on novel senescence-related genes in non-small cell lung cancer patients with drug sensitivity and tumor microenvironment. Adv Biol (Weinh). (2023) 7:e2300190. doi: 10.1002/adbi.202300190
20. Xu Q, Chen Y. An aging-related gene signature-based model for risk stratification and prognosis prediction in lung adenocarcinoma. Front Cell Dev Biol. (2021) 9:685379. doi: 10.3389/fcell.2021.685379
21. Collisson EA, Campbell JD, Brooks AN, Berger AH, Lee W, Chmielecki J, et al. Comprehensive molecular profiling of lung adenocarcinoma. Nature. (2014) 511:543–50. doi: 10.1038/nature13385
22. Okayama H, Kohno T, Ishii Y, Shimada Y, Shiraishi K, Iwakawa R, et al. Identification of genes upregulated in ALK-positive and EGFR/KRAS/ALK-negative lung adenocarcinomas. Cancer Res. (2012) 72:100–11. doi: 10.1158/0008-5472.CAN-11-1403
23. Yamauchi M, Yamaguchi R, Nakata A, Kohno T, Nagasaki M, Shimamura T, et al. Epidermal growth factor receptor tyrosine kinase defines critical prognostic genes of stage I lung adenocarcinoma. PloS One. (2012) 7:e43923. doi: 10.1371/journal.pone.0043923
24. Der SD, Sykes J, Pintilie M, Zhu CQ, Strumpf D, Liu N, et al. Validation of a histology-independent prognostic gene signature for early-stage, non-small-cell lung cancer including stage IA patients. J Thorac Oncol: Off Publ Int Assoc Study Lung Cancer. (2014) 9:59–64. doi: 10.1097/JTO.0000000000000042
25. Rousseaux S, Debernardi A, Jacquiau B, Vitte AL, Vesin A, Nagy-Mignotte H, et al. Ectopic activation of germline and placental genes identifies aggressive metastasis-prone lung cancers. Sci Trans Med. (2013) 5:186ra66. doi: 10.1126/scitranslmed.3005723
26. Barrett T, Wilhite SE, Ledoux P, Evangelista C, Kim IF, Tomashevsky M, et al. NCBI GEO: archive for functional genomics data sets--update. Nucleic Acids Res. (2013) 41:D991–5. doi: 10.1093/nar/gks1193
27. Zhu J, Fan Y, Xiong Y, Wang W, Chen J, Xia Y, et al. Delineating the dynamic evolution from preneoplasia to invasive lung adenocarcinoma by integrating single-cell RNA sequencing and spatial transcriptomics. Exp Mol Med. (2022) 54:2060–76. doi: 10.1038/s12276-022-00896-9
28. Avelar RA, Ortega JG, Tacutu R, Tyler EJ, Bennett D, Binetti P, et al. A multidimensional systems biology analysis of cellular senescence in aging and disease. Genome Biol. (2020) 21:91. doi: 10.1186/s13059-020-01990-9
29. Tacutu R, Thornton D, Johnson E, Budovsky A, Barardo D, Craig T, et al. Human Ageing Genomic Resources: new and updated databases. Nucleic Acids Res. (2018) 46:D1083–d1090. doi: 10.1093/nar/gkx1042
30. Therneau TM, Grambsch PM. Modeling Survival Data: Extending the Cox Model. New York: Springer. (2000).
31. Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Softw. (2010) 33:1–22. doi: 10.18637/jss.v033.i01
32. Simon N, Friedman J, Hastie T, Tibshirani R. Regularization paths for cox’s proportional hazards model via coordinate descent. J Stat Softw. (2011) 39:1–13. doi: 10.18637/jss.v039.i05
33. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. (2014) 15:550. doi: 10.1186/s13059-014-0550-8
34. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. (2012) 16:284–7. doi: 10.1089/omi.2011.0118
35. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. (2013) 4:2612. doi: 10.1038/ncomms3612
36. Aran D, Hu Z, Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. (2017) 18:220. doi: 10.1186/s13059-017-1349-1
37. Ayers M, Lunceford J, Nebozhyn M, Murphy E, Loboda A, Kaufman DR, et al. IFN-γ-related mRNA profile predicts clinical response to PD-1 blockade. J Clin Invest. (2017) 127:2930–40. doi: 10.1172/JCI91190
38. Cristescu R, Mogg R, Ayers M, Albright A, Murphy E, Yearley J, et al. Pan-tumor genomic biomarkers for PD-1 checkpoint blockade-based immunotherapy. Sci. (2018) 362:eaar3593. doi: 10.1126/science.aar3593
39. Wang Y, Weng W, Liang R, Zhou Q, Hu H, Li M, et al. Predicting T cell-inflamed gene expression profile in hepatocellular carcinoma based on dynamic contrast-enhanced ultrasound radiomics. J Hepatocell Carcinoma. (2023) 10:2291–303. doi: 10.2147/JHC.S437415
40. Bagaev A, Kotlov N, Nomie K, Svekolkin V, Gafurov A, Isaeva O, et al. Conserved pan-cancer microenvironment subtypes predict response to immunotherapy. Cancer Cell. (2021) 39:845–865.e7. doi: 10.1016/j.ccell.2021.04.014
41. He Y, Jiang Z, Chen C, Wang X. Classification of triple-negative breast cancers based on Immunogenomic profiling. J Exp Clin Cancer Res: CR. (2018) 37:327. doi: 10.1186/s13046-018-1002-1
42. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf. (2013) 14:7. doi: 10.1186/1471-2105-14-7
44. Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. (2018) 28:1747–56. doi: 10.1101/gr.239244.118
45. McGinnis CS, Murrow LM, Gartner ZJ. DoubletFinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Syst. (2019) 8:329–337.e4. doi: 10.1016/j.cels.2019.03.003
46. Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM 3rd, et al. Comprehensive integration of single-cell data. Cell. (2019) 177:1888–1902 e21. doi: 10.1016/j.cell.2019.05.031
47. Liu Y, He S, Wang XL, Peng W, Chen QY, Chi DM, et al. Tumour heterogeneity and intercellular networks of nasopharyngeal carcinoma at single cell resolution. Nat Commun. (2021) 12:741. doi: 10.1038/s41467-021-21043-4
48. Chen K, Wang Y, Hou Y, Wang Q, Long D, Liu X, et al. Single cell RNA-seq reveals the CCL5/SDC1 receptor-ligand interaction between T cells and tumor cells in pancreatic cancer. Cancer Lett. (2022) 545:215834. doi: 10.1016/j.canlet.2022.215834
49. Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun. (2021) 12:1088. doi: 10.1038/s41467-021-21246-9
50. Klement JD, Paschall AV, Redd PS, Ibrahim ML, Lu C, Yang D, et al. An osteopontin/CD44 immune checkpoint controls CD8+ T cell activation and tumor immune evasion. J Clin Invest. (2018) 128:5549–60. doi: 10.1172/JCI123360
51. Nallasamy P, Nimmakayala RK, Karmakar S, Leon F, Seshacharyulu P, Lakshmanan I, et al. Pancreatic tumor microenvironment factor promotes cancer stemness via SPP1-CD44 axis. Gastroenterology. (2021) 161:1998–2013.e7. doi: 10.1053/j.gastro.2021.08.023
52. Schlensog M, Ruehlmann AC, Haeberle L, Opitz F, Becher AK, Goering W, et al. Tenascin-C affects invasiveness of EGFR-mutated lung adenocarcinoma through a putative paracrine loop. Biochim Biophys Acta Mol Basis Dis. (2023) 1869:166684. doi: 10.1016/j.bbadis.2023.166684
53. Liu W, Maben Z, Wang C, Lindquist KC, Li M, Rayannavar V, et al. Structural delineation and phase-dependent activation of the costimulatory CD27:CD70 complex. J Biol Chem. (2021) 297:101102. doi: 10.1016/j.jbc.2021.101102
54. Liu H, Lv R, Song F, Yang Y, Zhang F, Xin L, et al. A near-IR ratiometric fluorescent probe for the precise tracking of senescence: a multidimensional sensing assay of biomarkers in cell senescence pathways. Chem Sci. (2024) 15:5681–93. doi: 10.1039/D4SC00595C
55. Huang W, Hickson LJ, Eirin A, Kirkland JL, Lerman LO. Cellular senescence: the good, the bad and the unknown. Nat Rev Nephrol. (2022) 18:611–27. doi: 10.1038/s41581-022-00601-z
56. Basisty N, Kale A, Jeon OH, Kuehnemann C, Payne T, Rao C, et al. A proteomic atlas of senescence-associated secretomes for aging biomarker development. PloS Biol. (2020) 18:e3000599. doi: 10.1371/journal.pbio.3000599
57. Hurd MW, Ralph MR. The significance of circadian organization for longevity in the golden hamster. J Biol Rhythms. (1998) 13:430–6. doi: 10.1177/074873098129000255
58. Yang WS, Stockwell BR. Inhibition of casein kinase 1-epsilon induces cancer-cell-selective, PERIOD2-dependent growth arrest. Genome Biol. (2008) 9:R92. doi: 10.1186/gb-2008-9-6-r92
59. Rodriguez N, Yang J, Hasselblatt K, Liu S, Zhou Y, Rauh-Hain JA, et al. Casein kinase I epsilon interacts with mitochondrial proteins for the growth and survival of human ovarian cancer cells. EMBO Mol Med. (2012) 4:952–63. doi: 10.1002/emmm.201101094
60. Yang J, Jiang Q, Liu L, Peng H, Wang Y, Li S, et al. Identification of prognostic aging-related genes associated with immunosuppression and inflammation in head and neck squamous cell carcinoma. Aging. (2020) 12:25778–804. doi: 10.18632/aging.v12i24
61. Park BJ, Kang JW, Lee SW, Choi SJ, Shin YK, Ahn YH, et al. The haploinsufficient tumor suppressor p18 upregulates p53 via interactions with ATM/ATR. Cell. (2005) 120:209–21. doi: 10.1016/j.cell.2004.11.054
62. Hassan MK, Kumar D, Naik M, Dixit M. The expression profile and prognostic significance of eukaryotic translation elongation factors in different cancers. PloS One. (2018) 13:e0191377. doi: 10.1371/journal.pone.0191377
63. Oh YS, Kim DG, Kim G, Choi EC, Kennedy BK, Suh Y, et al. Downregulation of lamin A by tumor suppressor AIMP3/p18 leads to a progeroid phenotype in mice. Aging Cell. (2010) 9:810–22. doi: 10.1111/j.1474-9726.2010.00614.x
64. Sirover MA. Pleiotropic effects of moonlighting glyceraldehyde-3-phosphate dehydrogenase (GAPDH) in cancer progression, invasiveness, and metastases. Cancer Metastasis Rev. (2018) 37:665–76. doi: 10.1007/s10555-018-9764-7
65. Nicholls C, Li H, Liu JP. GAPDH: a common enzyme with uncommon functions. Clin Exp Pharmacol Physiol. (2012) 39:674–9. doi: 10.1111/j.1440-1681.2011.05599.x
66. Gamez-Belmonte R, Mahapatro M, Erkert L, Gonzalez-Acera M, Naschberger E, Yu Y, et al. Epithelial presenilin-1 drives colorectal tumour growth by controlling EGFR-COX2 signalling. Gut. (2023) 72:1155–66. doi: 10.1136/gutjnl-2022-327323
67. Shen J, Bronson RT, Chen DF, Xia W, Selkoe DJ, Tonegawa S. Skeletal and CNS defects in Presenilin-1-deficient mice. Cell. (1997) 89:629–39. doi: 10.1016/S0092-8674(00)80244-5
68. Yaffe MB, Rittinger K, Volinia S, Caron PR, Aitken A, Leffers H, et al. The structural basis for 14-3-3:phosphopeptide binding specificity. Cell. (1997) 91:961–71. doi: 10.1016/S0092-8674(00)80487-0
69. Ogihara T, Isobe T, Ichimura T, Taoka M, Funaki M, Sakoda H, et al. 14-3-3 protein binds to insulin receptor substrate-1, one of the binding sites of which is in the phosphotyrosine binding domain. J Biol Chem. (1997) 272:25267–74. doi: 10.1074/jbc.272.40.25267
70. Tatar M, Bartke A, Antebi A. The endocrine regulation of aging by insulin-like signals. Sci (New York NY). (2003) 299:1346–51. doi: 10.1126/science.1081447
71. Chen CH, Chuang SM, Yang MF, Liao JW, Yu SL, Chen JJ. A novel function of YWHAZ/β-catenin axis in promoting epithelial-mesenchymal transition and lung cancer metastasis. Mol Cancer Res. (2012) 10:1319–31. doi: 10.1158/1541-7786.MCR-12-0189
72. Orjalo AV, Bhaumik D, Gengler BK, Scott GK, Campisi J. Cell surface-bound IL-1alpha is an upstream regulator of the senescence-associated IL-6/IL-8 cytokine network. Proc Natl Acad Sci United States America. (2009) 106:17031–6. doi: 10.1073/pnas.0905299106
73. Leon KE, Buj R, Lesko E, Dahl ES, Chen CW, Tangudu NK, et al. DOT1L modulates the senescence-associated secretory phenotype through epigenetic regulation of IL1A. J Cell Biol. (2021) 220:e202008101. doi: 10.1083/jcb.202008101
74. Garlanda C, Dinarello CA, Mantovani A. The interleukin-1 family: back to the future. Immunity. (2013) 39:1003–18. doi: 10.1016/j.immuni.2013.11.010
75. Wang YA, Sun Y, Palmer J, Solomides C, Huang LC, Shyr Y, et al. IGFBP3 modulates lung tumorigenesis and cell growth through IGF1 signaling. Mol Cancer Res. (2017) 15:896–904. doi: 10.1158/1541-7786.MCR-16-0390
76. Donoho G, Brenneman MA, Cui TX, Donoviel D, Vogel H, Goodwin EH, et al. Deletion of Brca2 exon 27 causes hypersensitivity to DNA crosslinks, chromosomal instability, and reduced life span in mice. Genes Chromosomes Cancer. (2003) 36:317–31. doi: 10.1002/gcc.10148
77. Vogel H, Lim DS, Karsenty G, Finegold M, Hasty P. Deletion of Ku86 causes early onset of senescence in mice. Proc Natl Acad Sci United States America. (1999) 96:10770–5. doi: 10.1073/pnas.96.19.10770
78. Fan Y, Gao Z, Li X, Wei S, Yuan K. Gene expression and prognosis of x-ray repair cross-complementing family members in non-small cell lung cancer. Bioengineered. (2021) 12:6210–28. doi: 10.1080/21655979.2021.1964193
79. Satyananda V, Oshi M, Endo I, Takabe K. High BRCA2 gene expression is associated with aggressive and highly proliferative breast cancer. Ann Surg Oncol. (2021) 28:7356–65. doi: 10.1245/s10434-021-10063-5
80. Yan B, Xie B, Huang M, Guo J, Sun J, Chen J, et al. Mutations and expressions of breast cancer 1/2 in lung cancer. Thorac Cancer. (2023) 14:1753–63. doi: 10.1111/1759-7714.14920
81. López-Otín C, Pietrocola F, Roiz-Valle D, Galluzzi L, Kroemer G. Meta-hallmarks of aging and cancer. Cell Metab. (2023) 35:12–35. doi: 10.1016/j.cmet.2022.11.001
82. Bao H, Cao J, Chen M, Chen M, Chen W, Chen X, et al. Biomarkers of aging. Sci China Life Sci. (2023) 66:893–1066. doi: 10.1007/s11427-023-2305-0
83. Lian I, Kim J, Okazawa H, Zhao J, Zhao B, Yu J, et al. The role of YAP transcription coactivator in regulating stem cell self-renewal and differentiation. Genes Dev. (2010) 24:1106–18. doi: 10.1101/gad.1903310
84. Pan D. The hippo signaling pathway in development and cancer. Dev Cell. (2010) 19:491–505. doi: 10.1016/j.devcel.2010.09.011
85. Wang CX, Yan J, Lin S, Ding Y, Qin YR. Mutant-allele dispersion correlates with prognosis risk in patients with advanced non-small cell lung cancer. J Cancer Res Clin Oncol. (2023) 149:8545–55. doi: 10.1007/s00432-023-04801-3
86. Hernández Borrero LJ, El-Deiry WS. Tumor suppressor p53: Biology, signaling pathways, and therapeutic targeting. Biochim Biophys Acta Rev Cancer. (2021) 1876:188556. doi: 10.1016/j.bbcan.2021.188556
87. Qin C, Fan X, Sai X, Yin B, Zhou S, Addeo A, et al. Development and validation of a DNA damage repair-related gene-based prediction model for the prognosis of lung adenocarcinoma. J Thorac Dis. (2023) 15:6928–45. doi: 10.21037/jtd
88. Zhao Y, Qing B, Xu C, Zhao J, Liao Y, Cui P, et al. DNA damage response gene-based subtypes associated with clinical outcomes in early-stage lung adenocarcinoma. Front Mol Biosci. (2022) 9:901829. doi: 10.3389/fmolb.2022.901829
89. Chen HA, Ho YJ, Mezzadra R, Adrover JM, Smolkin R, Zhu C, et al. Senescence rewires microenvironment sensing to facilitate antitumor immunity. Cancer Discovery. (2023) 13:432–53. doi: 10.1158/2159-8290.CD-22-0528
90. Marin I, Boix O, Garcia-Garijo A, Sirois I, Caballe A, Zarzuela E, et al. Cellular senescence is immunogenic and promotes antitumor immunity. Cancer Discovery. (2023) 13:410–31. doi: 10.1158/2159-8290.CD-22-0523
91. Shahbandi A, Chiu FY, Ungerleider NA, Kvadas R, Mheidly Z, Sun MJS, et al. Breast cancer cells survive chemotherapy by activating targetable immune-modulatory programs characterized by PD-L1 or CD80. Nat Cancer. (2022) 3:1513–33. doi: 10.1038/s43018-022-00466-y
92. Prieto LI, Sturmlechner I, Goronzy JJ, Baker DJ. Senescent cells as thermostats of antitumor immunity. Sci Trans Med. (2023) 15:eadg7291. doi: 10.1126/scitranslmed.adg7291
93. Holze C, Michaudel C, Mackowiak C, Haas DA, Benda C, Hubel P, et al. Oxeiptosis, a ROS-induced caspase-independent apoptosis-like cell-death pathway. Nat Immunol. (2018) 19:130–40. doi: 10.1038/s41590-017-0013-y
Keywords: cellular senescence, lung adenocarcinoma, tumor microenvironment, heterogeneity, machine learning
Citation: Ru K, Cui L, Wu C, Tan XX, An WT, Wu Q, Ma YT, Hao Y, Xiao X, Bai J, Liu X, Xia XF and Zhao MQ (2024) Exploring the molecular and immune landscape of cellular senescence in lung adenocarcinoma. Front. Immunol. 15:1347770. doi: 10.3389/fimmu.2024.1347770
Received: 01 December 2023; Accepted: 08 August 2024;
Published: 29 August 2024.
Edited by:
Kevin M. Mcbride, University of Texas MD Anderson Cancer Center, United StatesReviewed by:
Petar Ozretić, Rudjer Boskovic Institute, CroatiaIvana Samarzija, Rudjer Boskovic Institute, Croatia
Copyright © 2024 Ru, Cui, Wu, Tan, An, Wu, Ma, Hao, Xiao, Bai, Liu, Xia and Zhao. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Miao Q. Zhao, zhaomqsd@163.com
†These authors have contributed equally to this work