Consensus Clustering for Immunosenescence Subtype Identification and DEGs Pathway Analysis

The “limma” program was utilized to identify 1740 DEGs between PAAD and normal tissues. Subsequently, 37 IMSRGs were pinpointed (as depicted in Fig. 1A and B). The expression levels of these 37 IMSRGs exhibited significant variation between the normal tissues and PAAD samples, as illustrated in Fig. 1C.

Subsequently, employing consensus clustering, two clusters of PAAD associated with immunosenescence were delineated. Utilizing the TCGA dataset, these two clusters demonstrated distinct immunosenescence gene expression profiles following k-means clustering analysis, as depicted in Fig. 2A-B. Furthermore, survival analysis suggested that these immunosenescence-derived subgroups exhibited markedly distinct clinical outcomes. Notably, the C2 subgroup exhibited a better survival rate compared to the C1 subgroup, as illustrated in Fig. 2C.

In this research, pivotal DEGs and critical signaling pathways within the two immunosenescence subgroups were uncovered to elucidate the molecular mechanisms influencing prognosis. A total of 1640 dysregulated genes were identified, as depicted in Fig. 2D. These genes were found to be enriched in functions of the immune system, such as cytokine-cytokine receptor interactions, chemokine signaling pathways, cell adhesion molecules (CAMs), hematopoietic cell lineage, phagosome, immune system mechanisms, immune reactions, modulation of immune system processes, defense mechanisms, and cell activation, as represented in Fig. 2E-F.

Fig. 1

Screening and expression of differential genes. (A) The diverse expression of DEGs between PAAD and normal tissue. (B) Venn diagram showing common genes of IMSRGs and DEGs. (C) Heatmap illustrating the expression of 37 immunosenescence genes in the PAAD and normal tissue specimens utilizing the TCGA database

Fig. 2

Identification of Immunosenescence-associated subgroups through consensus clustering. (A)The delta area curve for consensus clustering illustrates the comparative differences in the area beneath the cumulative distribution function (CDF) curve for k values ranging from 2 to 10; (B) The consensus matrix of the clustering analysis via k-means clustering (k = 2); (C) KM curves depicting patient OS in the C1 and C2 subcategories; (D) Distribution of the DEGs between the C1 and C2 subgroups in the TCGA dataset; and (E-F) the KEGG and GO enrichment analyses of signaling pathways

Establishment and Validation of the Immunosenescence Risk Signature

In this study, we have developed an innovative prognostic model predicated on the IMSRGs. Our analysis revealed that three IMSRGs exhibited a significant correlation with OS as determined by the Cox univariate regression analysis, which is depicted in Fig. 3A. Subsequently, the Lasso regression analysis was employed to assess these three IMSRGs, and they were consequently selected to form the predictive model, as illustrated in Fig. 3B-C. By setting the lambda value at 0.009976, we were able to pinpoint three key genes. The predictive model is encapsulated in the following formula: Riskscore = (-0.4629 * NCAM1) - (0.1211 * SESN1) + (0.1518 * MAFF).

Utilizing the median risk score (RS), we categorized our study population into a high-risk cohort (H) and a low-risk cohort (L). We utilized the TCGA dataset as our training cohort and the GSE85916 dataset for validation purposes. In the training cohort, the receiver operating characteristic (ROC) curves were developed with area under the curve (AUC) values of 0.65 for the 1-year survival, 0.73 for the 3-year survival, and 0.69 for the 5-year survival. Similarly, in the validation cohort, the ROC curves were generated with AUC values of 0.69 for the 1-year survival, 0.81 for the 3-year survival, and 0.77 for the 5-year survival (Fig. 3D). Furthermore, the connection between the RSs and OS status was evaluated. The findings revealed that a significantly lower proportion of patients in the H cohort survived, both in the training and validation cohorts, with a more pronounced effect in the H cohort (Fig. 3E). The Kaplan-Meier (KM) survival assessment corroborated that individuals in the H cohort experienced shorter survival times compared to those in the L cohort, as evidenced in both the TCGA and GSE85916 cohorts (Fig. 3F).

Fig. 3

Constructing and confirming the Immunosenescence risk profile. (A) Univariate Cox regression examination of IMSRGs; (B-C) Lasso Cox regression evaluation; (D) ROC curves with the AUC values; (E) RS allocation, OS condition of each subject, and the heat diagrams illustrating the prognostic utilization of the 3-gene indicator; and (F) KM plots of the subject’s OS in the elevated-risk cohort and reduced-risk cohort classifications

KEGG, GO, and GSEA Analysis Related to Immunosenescence Risk Signature

In our research, we identified a total of 1,533 DEGs when comparing the H and L cohorts, as illustrated in Fig. 4A. These genes were notably involved in various biological processes, including pancreatic secretion, cell adhesion molecules (CAMs), insulin secretion, protein digestion and absorption, hematopoietic cell development, and extracellular matrix components such as the extracellular region, vesicles, and extracellular space. Additionally, they were associated with responses to external stimuli (Fig. 4B-C).

GSEA was employed to further illuminate the operational distinctions between the H and L cohorts. This analysis highlighted significant differences in pathways such as neuroactive ligand-receptor interactions, primary bile acid synthesis, beta-alanine metabolism, vascular smooth muscle function, neurotransmitter transport mechanisms, intercellular signaling in multicellular organisms, regulation of transporter activity, modulation of neurotransmitter levels, and synaptic transmission specifically involving glutamate (Fig. 4D-E). These findings underscore the complex biological underpinnings that distinguish the H and L cohorts at the molecular level.

Fig. 4

DEGs, KEGG analysis, GO analysis, and GSEA analysis. (A) The allocation of the DEGs between the H and L cohorts; (B-C) KEGG and GO enrichment analyses of signaling pathways; (D-E) The GSEA pinpoints the core signaling pathways

Somatic Mutations and the Tumor Microenvironment in Diverse Immunosenescence Risk Cohorts

In this analysis, we explored the disparities in the TME across two distinct patient cohorts. Figure 5A presents a comprehensive overview of the immune cell infiltration patterns among 176 PAAD patients, as documented in the TCGA dataset. Employing the Cibersort algorithm and the lm22 gene signature matrix, we assessed the presence of 22 distinct immune cell types in each cohort.

Our findings revealed that individuals in the H cohort displayed a notable reduction in the prevalence of native B cells, plasma cells, CD8 + T cells, stimulated CD4 + T cell memory, and monocytes, as opposed to those in the low-risk (L) cohort. Conversely, the H cohort demonstrated increased levels of memory B cells, M0 macrophages, and eosinophils (Fig. 5B).

Furthermore, the L cohort demonstrated superior immune activity, as evidenced by higher ImmuneScore, StromalScore, and EstimateScore, indicating a more robust immune response within the TME (Fig. 5C). Additionally, the L cohort exhibited enhanced expression of various immune checkpoint molecules and human leukocyte antigen (HLA) genes, suggesting a more active immune surveillance mechanism in comparison to the H cohort (Fig. 5D-E). These observations highlight the complex interactions among immune cell dynamics and the risk stratification of PAAD patients.

Fig. 5

Immunological characteristics of distinct Immunosenescence Risk Cohorts. (A) Comparative distribution of immune cell infiltration; (B) Notable disparities in immune cell infiltration between cohorts; (C) Stromal, Immune, and Estimate scores for L and H cohorts; Differential expression patterns of various immune checkpoint genes (D); and HLA genes (E) across the two categories. (* P<0.05;** P<0.01;*** P<0.001;****P<0.0001)

Connection Between the Immunosenescence Risk Signatures and the TME

Our investigation revealed a notable correlation between the RS and the immune cell composition. In particular, there was a positive correlation observed among individuals with a higher RS, notably with memory B cells, M0 macrophages, and regulatory T cells (Tregs). This suggests a potential immunosuppressive role of these cell types in the context of increased risk. Conversely, a negative correlation was identified with M2 macrophages, stimulated mast cells, monocytes, stimulated natural killer (NK) cells, and CD8 T cells, indicating that these cell types might be protective or their reduced presence could be indicative of a more aggressive disease phenotype. These relationships were systematically depicted across multiple figures (Fig. 6A-H), providing a nuanced view of the immunological landscape in relation to risk stratification.

Fig. 6

The association between the TME and Immunosenescence risk signature. (A–H) The correlation between RS and the immune-infiltrating cell populations

Correlation Analysis Between IMSRG Expression Level and the Infiltrating Immune Cells

Utilizing Lasso-cox regression analysis, we identified three IMSRGs. We further explored the correlation of the expression of these genes with the infiltration of six distinct immune cell populations—neutrophils, macrophages, B cells, dendritic cells, CD4 + T cells, and CD8 + T cells—leveraging the TIMER database. The outcomes are illustrated in Fig. 7.

Our analysis suggested that MAFF expression demonstrated a notable positive association with the presence of CD8 + T cells, neutrophils, and dendritic cells. Additionally, NCAM1 levels exhibited a positively associated with an increased infiltration of CD8 + T cells, CD4 + T cells, macrophages, neutrophils, and dendritic cells. Intriguingly, SESN1 expression was observed to have a positive association with the infiltration degrees of all six immune cell types examined. These findings highlight the possible influence of these genes in modulating the immune response within the tumor microenvironment.

Fig. 7

Correlation between IMSRGs abundance and immune cell infiltration as analyzed by TIMER. (A-C) Expression patterns of MAFF, NCAM1, and SESN1 in relation to immune cell infiltration levels

Drug Sensitivity of the IMSRGs in PAAD

Employing GSCA, we observed that the IMSRGs were associated with resistance to a range of therapeutic agents, as depicted in Figs. 8A-B. Through comprehensive drug response profiling (CTRP) sensitivity assays, it was revealed that NCAM1 and SESN1 exhibited resistance to several medications, while MAFF displayed a heightened sensitivity to these treatments, notably to drugs like apicidin, axitinib, and doxorubicin (Fig. 8A).

Further corroboration was provided by the Genomics of Drug Sensitivity in Cancer (GDSC) drug sensitivity tests, which indicated a reciprocal relationship in drug sensitivity between MAFF and the other two genes, SESN1 and NCAM1. Specifically, MAFF was found to be particularly responsive to drugs such as docetaxel, trametinib, tubastatin A, and navitoclax, in contrast to the resistance patterns observed for SESN1 and NCAM1 (Fig. 8B). These findings highlight the differential drug response profiles linked to the expression of these key genes, potentially offering insights into personalized therapeutic strategies.

Fig. 8

Drug responsiveness of IMSRGs in pancreatic ductal adenocarcinoma. (A) Association between gene expression and susceptibility to CTRP compounds (Top 30). (B) Interplay of gene expression levels and reactivity to GDSC agents (Top 30)

Single Cell Analysis of MAFF

In our current research, MAFF was identified as a key factor significantly linked to the prognosis in PAAD. To delve deeper into the interplay between MAFF expression and the tumor immune landscape, we utilized the TISCH. Our analysis within the PAAD_CRA001160 dataset revealed that MAFF expression was notably elevated in various cell types, encompassing ductal cells, endothelial cells, endocrine cells, and stellate cells, as visualized in Fig. 9A-B. Similarly, in the PAAD_GSE111672 dataset, MAFF was found to be prominently expressed in malignant cells and ductal cells, as illustrated in Fig. 9C-D. These outcomes indicate a profound role for MAFF in the cellular composition and immune dynamics of PAAD.

Fig. 9

Single cell analysis of MAFF in PAAD. (A-B) Cellular-level mapping unveiled the allocation of MAFF across diverse immune cell populations in PAAD_CRA001160 and (C-D) PAAD_GSE111672

Multilevel Expression Verification and in Vitro Functional Exploration of MAFF

We conducted a search using GEPIA2.0, and the findings indicated that MAFF exhibited elevated expression in PAAD, with high MAFF levels correlating to an unfavorable prognosis (Fig. 10A-B). Then we found that the HPA database showed that MAFF was highly expressed in PAAD tissues compared with normal pancreatic tissues (Fig. 10C). To examine the increased MAFF expression in PAAD, we employed W-B analyses. The results revealed that MAFF displayed notable upregulation in pancreatic cancer cell lines, specifically PANC-1, and SW1990, in comparison to the normal pancreatic cell line H6C7 (Fig. 10C-D).

Fig. 10

The prognostic value and expression of MAFF. (A) Differential analysis of mRNA expression of MAFF in PAAD in GEPIA2.0 database. (B) Kaplan–Meier curve for OS between the high- and low-MAFF expression cohort in PAAD in GEPIA2.0 database. (C) IHC images of MAFF in PAAD tissues and normal pancreatic tissues from HPA database. (D-E) The expression level of MAFF was examined by W-B. (** P<0.01;*** P<0.001)