Introduction

Pancreatic cancer poses a significant challenge within the field of oncology.¹ Due to its subtle onset, aggressive behavior, and frequent late-stage detection, it becomes a formidable adversary, causing patients and healthcare providers to grapple with its complexities.^2–5 Although pancreatic cancer is relatively uncommon, it exhibits a high mortality rate, ranking 12th among global cancer cases, with approximately 495,000 new instances and 466,000 annual fatalities.⁶ Predominantly affecting older adults, the majority of diagnoses occur in individuals aged 65 and older.⁷ Regional variations in incidence are notable, with elevated rates in North America and Europe.

The diagnostic process for pancreatic cancer typically encompasses a combination of imaging tests, including CT scans, MRI, and endoscopic ultrasound (EUS), to visualize the pancreas and evaluate tumor size and location.^8–10 Blood tests, like CA 19–9, may serve as tumor markers.¹⁰ Ultimately, confirmation of the diagnosis often necessitates a biopsy, obtainable through endoscopic or surgical procedures.¹¹ Treatment modalities for pancreatic cancer hinge on the disease’s stage.¹² Surgery stands as the primary curative option for localized tumors, yet many cases receive diagnoses at advanced stages when surgery is unviable.¹³ In such scenarios, treatment options encompass chemotherapy, radiation therapy, targeted therapy, and immunotherapy.¹⁴ Palliative care frequently enters the picture to manage symptoms and enhance the quality of life for patients with advanced-stage disease.¹⁵ Therefore, addressing the treatment quandary of pancreatic cancer underscores the need to devise novel prognostic stratification strategies and explore potential biomarkers within its tumor microenvironment.

Cell senescence denotes a condition characterized by a stable and irreversible halt in the cell cycle, resulting in a loss of the cell’s capacity for division and replication.¹⁶ Dr. Leonard Hayflick initially elucidated this phenomenon during the 1960s, and it has since become a pivotal concept in the realm of biology, particularly in the context of senescence and cancer.¹⁷ In the context of tumors, cell senescence carries a dual significance. On one hand, it operates as an initial defense mechanism, effectively arresting the proliferation of impaired cells and erecting a barrier against cancer development.¹⁸ This process, known as oncogene-induced senescence (OIS), thwarts the transformation of normal cells into malignant counterparts.¹⁹ Conversely, senescent cells present within a tumor can stimulate tumor progression by releasing a pro-inflammatory mixture termed the senescence-associated secretory phenotype (SASP). This secretion can foster an environment conducive to tumor expansion, angiogenesis, tissue remodeling, and evasion of the immune system.²⁰ Consequently, the role of cell senescence within tumors manifests as intricate and context-dependent, exerting influence on both tumor suppression and promotion.

In our research, we integrated single-cell sequencing data with bulk transcriptome sequencing data to investigate the diversity and expression patterns of genes associated with cell senescence in pancreatic cancer. Additionally, we developed a prognostic model that is pertinent to assess the outlook for individuals diagnosed with pancreatic cancer. Our study offers valuable insights that can serve as a reference point for early detection and prognosis evaluation in pancreatic cancer cases.

Materials and Methods

Single Cell Sequencing Data Download and Processing

The GEO database houses numerous single-cell sequencing datasets. Specifically, we downloaded the single cell sequencing datasets GSE155698 and GSE154778, which are related to pancreatic cancer.^21,22 GSE155698 comprises 17 tumor samples and 3 normal tissue samples, while GSE154778 encompasses 16 tumor samples (Supplemental Table 1).^21,22 For processing and analyzing these single-cell data, we utilized the “Seurat” package, version 4.3.0.1. Quality control measures were implemented as follows:Genes expressed in fewer than 3 cells were excluded; Cells with less than 200 gene expressions and those with over 9000 gene expressions were removed; Cells with mitochondrial gene expression exceeding 15% were excluded; Cells with a total gene expression surpassing 200,000 were removed; Samples with fewer than 200 cells were excluded.

To standardize the data, the “NormalizeData” function was employed. Highly variable genes, set at 3000, were identified using the “vst” method. The “SCT” method was utilized to integrate samples and eliminate mitochondrial genes and cell cycle effects. To reduce data dimensionality, PCA analysis was initially conducted. Subsequently, we employed the “tSNE” method with a dimensionality reduction setting of 20. For cell clustering, the “KNN” method was employed, setting dimensions to 20, resolution to 0.4, and random. seed to 2023. Cell annotations were based on marker genes obtained from previously published pancreatic cancer single-cell sequencing articles. Results of single cell analysis were visualized using a “tsne” chart. Additionally, the “Seurat” package’s “FindMarkers” function was used to analyze differences between the two groups. In this study, the “PercentageFeatureSet” function was utilized to calculate the cell senescence score based on the expression of cell senescence phenotype genes in cells.

Bulk Transcriptome Data Download and Processing

Three cohorts related to pancreatic cancer were acquired, encompassing not only transcriptome data but also clinical data. The TCGA transcriptome data was sourced from the UCSC Xena website (https://xenabrowser.net/).²³ Meanwhile, the PAAD-AU and PAAD-CA transcriptome data were retrieved from the ICGC website (https://dcc.icgc.org/) (Supplemental Tables 2 and 3). To ensure uniformity, all data underwent standardization. Samples that included both transcriptome data and clinical information were retained by aligning transcriptome data with clinical data.

Acquisition of Cell Senescence-Related Gene Sets

The 125 cell senescence-related genes were collected in the published article titled “a new gene set identifies senescent cells and predicts senescence-associated pathways across tissues”.²⁴

Single Sample Gene Set Enrichment Analysis (ssGSEA)

Based on the difference values obtained through the GSEA R package analysis, the order, arranged from large to small, was utilized to depict the fold change pattern of gene set expression between the two groups. Typically, when a gene set is enriched towards the top, it indicates up-regulation, while enrichment towards the bottom suggests down-regulation. For this analysis, the “clusterProfiler” package was employed, and it utilized the reference gene set “H.ARI.v7.0.entrez.GMT” to perform these assessments.

Weighted Co-Expression Network Analysis (WGCNA)

In this research, WGCNA (Weighted Gene Co-Expression Network Analysis) was employed to investigate gene sets closely associated with cell senescence phenotypes in pancreatic cancer. The soft threshold was explored across a range, specifically from 1 to 10 with a step size of 1, as well as from 12 to 20 with a step size of 2. To determine the appropriate soft threshold, the “WGCNA” package’s “pickSoftThreshold” function was utilized. Several parameters were set for the analysis: a minimum of 200 genes per module, “deepSplit” was set to 2 for enhanced clustering, and genes were organized into modules. To consolidate modules that exhibited similarity, this study applied a cut-off value of 0.3 and correlated them with cell senescence phenotypes.

Unsupervised Consistent Cluster Analysis

The “ConsensusClusterPlus” package was used to perform cluster analysis based on the expression pattern of cell senescence genes. The maximum cluster number was set to 7, clusterAlg = “hc”, innerLinkage = “average”, finalLinkage = “average”.

Construction of the Prognostic Signature

Within the TCGA cohort, the initial step involved utilizing Univariate COX regression to perform an initial screening of prognostic genes, with a significance threshold of p<0.05. Following this, LASSO regression was employed to further filter variables and construct a prognostic signature. The “glmnet” R package was utilized for this purpose, and the “maxit” parameter was set to 1000. Patients were subsequently categorized into two groups based on the median senescence signature risk score (SENRs) values, resulting in a high SENRs group and a low SENRs group. Prognostic analysis, decision curve analysis, and PCA (Principal Component Analysis) were then conducted separately for these two groups to explore their prognostic implications and assess the decision-making utility of the model.

Validation of the Prognostic Signature

PAAD-AU and PAAD-CA cohorts as validation cohort. Patients were also divided into high SENRs group and low SENRs group based on the median SENRs value. In the validation cohort, prognostic analysis, decision curve analysis and PCA analysis of high and low SENRs groups were explored.

Tumor Mutation Analysis

The “maftools” package was used to download mutation data for pancreatic cancer. By matching, samples containing both transcriptome and mutation data are retained, and a waterfall map is used to show the mutant landscape.

Immunomicroenvironment Correlation Analysis

We got through TIMER2 website (http://timer.cistrome.org/) to download a variety of methods based on the TCGA data to calculate the immune cells infiltrating data, and use the waterfall figure to show results. The expression of immune checkpoint genes, tumor necrosis genes, and leukocyte antigen genes was investigated among model groups. The “estimate” package was used to calculate the sample immunity score. In this study, the TIDE website was used to predict the immunotherapy effect of the sample as well as the TIDE score.

Pseudo-Time Series Analysis

Pseudo-time series analysis is the use of specific algorithms to analyze changes in gene expression during cell differentiation state transitions, which in turn place each cell in the appropriate position in the trajectory. In this study, the “monocle2” package was used to perform pseudo-time series analysis and explore the changes in the expression of specific genes during this differentiation process.

Clinical Correlation Analysis

The “meta” package was used for the integrated analysis. The “rpart” package was used to perform decision curve analysis based on patient clinical data and models, and was divided into different groups. In this study, the “regplot” package was used to plot nomogram based on patient data. The STRING web site (https://cn.string-db.org/) was used for protein interaction analysis.

Cell Lines, Culture Conditions and Cell Transfection

Capan-1 and PANC-1 cells were procured from the Shanghai Institute of Biological Sciences in Shanghai, China. These cells were cultured in DMEM (Dulbecco’s Modified Eagle Medium) supplemented with 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin solution. All cell cultures were maintained at a temperature of 37°C and in an environment containing 5% CO2. For transfection, Lipofectamine 3000 from Thermo Fisher Scientific in Waltham, MA, USA was employed. Cells were transfected with siRNA (RiboBio, Guangzhou, China) following the manufacturer’s instructions. Three target siRNAs were used for CAV1, which included:siCAV1-1: 5’-CCACCTTCACTGTGACGAA-3’; siCAV1-2: 5’-GCACATCTGGGCAGTTGTA-3’ siCAV1-3: 5’-GCGACCCTAAACACCCAA-3’. These siRNAs were utilized for gene knockdown experiments.

CCK-8 Assay, Colony Formation Analysis, Migration Assays

These experimental methods have been reported in our previous studies.²⁵

PCR Tests of CAV1 in Clinical Samples

For this study, 8 pancreatic cancer tissues and corresponding adjacent control tissues were obtained. Total RNA extraction from cell lines followed the manufacturer’s protocol, employing TRIzol reagents (Invitrogen, CA, USA). Subsequently, cDNA synthesis was conducted with the PrimeScript RT Reagent Kit (Takara, Nanjing, China). Real-time quantitative PCR (qRT-PCR) was carried out using AceQ Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China). The primers utilized in this investigation are detailed below:

Forward primerGGGTGTGGAAACCTCGTCTT

Reverse primerGGTGTAGAGATGTCCCTGCG

Statistical Analysis

The initial step in the prognostic gene screening involved conducting univariate COX analysis. To compare gene expression between the two groups of the model, a rank sum test was employed. Prognostic analysis between these two groups was executed using the KM (Kaplan-Meier) survival curve method. To analyze the correlation between model values, specifically SENRs, and TIDE values, the “spearman” method was applied. The analysis was conducted using R software version 4.0.5, with a defined p-value as the significance threshold. A p-value less than 0.05 was considered statistically significant.

Results

Cell Senescence Analysis at the Single-Cell Level

The research process is shown in Figure 1.

Figure 1 The flow chart.

Cell senescence phenotypes were investigated at the single-cell level in pancreatic cancer. As shown in Figure 2A–B, through quality control and sample integration, a total of 29 tumor and 3 normal samples were retained, and there was heterogeneity between tumor and normal samples. As shown in Figure 2C–E, tumor samples were clustered into 19 clusters, and a total of 11 cell types were annotated according to the expression of cell type marker genes. Subsequently, the cell senescence phenotype fraction was quantified in each cell, and this group was divided into High-Senescence group and Low-Senescence group according to the median value, as shown in Figure 2F. The High-Senescence group was mainly distributed in Monocytes, Macrophages, Granulocytes and Fibroblastes. As shown in Figure 2G, the 10 genes closely related to cell senescence are RPS8, RPL11, RPL22, UQCRH, TACSTD2, CAPZB, NFRSR1B,C1QA, MXRA8 and C1QC. As shown in Figure 2H–L, activation pathways associated with cell senescence in pancreatic cancer include Apoptosis, Epithelial Mesenchymal Transition, Hypoxia, Inflammatory Response and TGF Beta. We compared the expression of senescence related genes in pancreatic cancer with normal tissues at the single-cell level, and the results are shown in Figure S1.

Figure 2 Single cell sequencing analysis. (A and B) Quality control and sample integration. (C) All cells were clustered into 19 clusters. (D) Markers for cell clustering. (E) Annotation on cell types. (F) All cells were divided into high and low groups according to senescence score. (G) Analysis of differentially expressed genes between high - and low-senescence groups, in the form of volcano map. (H–L)GSEA between high and low groups showed that activation pathways associated with cell senescence in pancreatic cancer include Apoptosis, Epithelial Mesenchymal Transition, Hypoxia, Inflammatory Response and TGF Beta.

To Search for Genes Closely Related to Cell Senescence

In pancreatic cancer, a variety of methods were used in this study to obtain genes associated with cellular senescence. First, WGCNA was conducted in the TCGA cohort. As shown in Figure 3A, when the optimal soft threshold is 7, the data conforms to the power-law distribution, and Mean connectivity tends to be stable with the increase of the soft threshold. As shown in Figure 3B–C, the genes were clustered into 9 non-gray modules, of which the green module was most associated with the cell senescence phenotype (cor = 0.53 and p <0.05). We further explored the inter-gene correlation in the green module, as shown in Figure 3D. There was a positive correlation between Module membership in green module and Gene significance for body weight (cor = 0.64 and p <0.05). As shown in Figure 3E–F, patients are divided into two clusters according to cell senescence genes, and Cluster2 patients have poor prognosis (p<0.05), as shown in Figure 3G. Differential expressed gene analysis was carried out through package DEseq2 and screening criteria were set |LogFC| >1 and adjust p-value<0.05. A total of 6300 genes were obtained. As shown in Figure 3H, the intersection of 6295 genes related to cell senescence obtained by single-cell analysis, 1441 genes obtained by WGCNA analysis, and 6300 genes obtained by differential expressed gene analysis yielded 303 genes closely related to cell senescence.

Figure 3 Bulk RNA sequencing analysis. (A–D) Weighted coexpression network analysis(WGCNA). The green module was found to be most significantly associated with senescence phenotypes in pancreatic cancer. (E–F) Unsupervised clustering according to cell senescence genes showed that cluster2 had significantly worse prognosis. (G)Differential expressed gene analysis. A total of 6300 genes were obtained. (H) The intersection of 6295 genes related to cell senescence obtained by single-cell analysis, 1441 genes obtained by WGCNA analysis, and 6300 genes obtained by differential expressed gene analysis.

The Construction and Validation of the Prognostic Model

Starting with the 303 genes identified through the intersection of single-cell analysis, Weighted Gene Co-expression Network Analysis (WGCNA), and differential expression analysis, we aimed to narrow down this gene set to those most relevant for prognosis. To achieve this, we first conducted univariate Cox regression analysis on the TCGA cohort, setting a significance threshold of p < 0.05. This step enabled the selection of genes significantly associated with patient survival, thereby refining our set to 21 prognostic genes with a strong correlation to patient outcomes(Figure 4A).

Figure 4 Construction and validation of the prognostic model. (A) Univariate COX analysis. A total of 21 prognostic genes in the TCGA cohort and PAAD-AU cohort were obtained. (B and C)LASSO regression analysis. (D–H)Prognosis validation in the TCGA cohort. The prognosis of patients in SENRs_High group was poor (p<0.05). The prediction effect of the model was better than that of other clinical indicators, and patients could benefit more than other clinical indicators by taking appropriate clinical measures according to the model. SENRs is a risk factor for patient prognosis (p<0.05). This model can distinguish patients well. (I–M)Prognosis validation in the PAAD-AU cohort. (N–R)Prognosis validation in the PAAD-CA cohort.

To further ensure the robustness and specificity of our prognostic model, we applied Least Absolute Shrinkage and Selection Operator (LASSO) regression on these 21 genes(Figure 4B–C). LASSO regression is particularly useful for managing high-dimensional data by reducing the number of variables, effectively eliminating weaker predictors. After tuning the regularization parameter (lambda) to an optimal value of 0.05, we obtained a final prognostic model comprising six genes: CAV1, BIRC3, DCBLD2, CD109, IL1RAP, and SP100. These six genes were each assigned a coefficient, allowing us to construct a risk score that quantifies the patient’s prognosis.

This refined selection process ensured that the model retained only the most predictive genes from the initial 303-gene set, improving its predictive accuracy and clinical relevance. As shown in Figure 4D–H, patients with high-risk scores based on this model had significantly poorer prognoses (p < 0.05), indicating the model’s potential clinical utility. The prediction effect of the model was better than that of other clinical indicators, and patients could benefit more than other clinical indicators by taking appropriate clinical measures according to the model. SENRs is a risk factor for patient prognosis (p<0.05). This model can distinguish patients well. In order to further verify the accuracy of the model, we also conducted analysis in the AU and CA cohort, as shown in Figure 4I–R, SENRs_High group had poor prognosis (p<0.05), and the prediction effect of the model was better than that of other clinical indicators. Moreover, patients who took appropriate clinical measures according to the model benefited more than other clinical indicators. SENRs were all risk factors for patient prognosis (p<0.05), and the model could distinguish patients well.

Mutation and Immune Correlation Analysis

This study further explored the differences in mutation and immune infiltration between the model groups in order to better identify patients who are more suitable for immunotherapy. As shown in Figure 5A–B, in the SENRs_Low group, the mutation incidence was 74.39%, and the top 5 most mutated genes were KRAS, TP53, SMAD4, CDKN2A and MUC16. In the SENRs_High group, the mutation incidence was 90.12%, and the top 5 most mutated genes were KRAS, TP53,CDKN2A and TTN. As shown in Figure 5C, patients in SENRs_High group had a higher degree of immune infiltration, such as B cells and T cells. As shown in Figure 5D–F, most of the leukocyte antigen genes, immune checkpoint genes, and tumor necrosis genes were highly expressed in the SENRs_High group, such as HLA-E, HLA-A, CTLA4,CD86, CXCL10,EIF2A, etc. As shown in Figure 5G–J, the immune score(p<0.05), stromal score(p<0.001) and total score(p<0.01) of SENRs_High group were higher, while the tumor purity was lower (p<0.01). This may be related to the higher degree of immune cell infiltration but also the higher degree of immune escape in the high SENRs group. As shown in Figure 5K, in patients who responded to immunotherapy, the model value SENRs was lower(p<0.001). And the TIDE score was higher (p<0.001)was higher in SENR_high group, and SENRs was significantly positively correlated with TIDE (cor 0.35 and p< 0.001)(Figure 5L–M). This study further explored the activation of the pathway in the SENRs_High group and found that the TGF beta pathway and Epithelial Mesenchymal Transition pathway were significantly activated as shown in Figure 5N–O.

Figure 5 Analysis of tumor microenvironment. (A and B) Mutant landscape of SENR_low group and SENR_high group. (C)Different landscape of immune cell infiltration between the two groups. (D–F)Differences in the expression of leukocyte antigen genes, immune checkpoint genes, and tumor necrosis genes between the two groups. (G–J) the immune score, stromal score and total score of SENRs_High group were higher, while the tumor purity was lower. (K–M)Prediction of response to immunotherapy in two groups. (N and O)GSEA in SENRs_High group.(*p<0.05, **p<0.01, ***p<0.001).

Gene Expression Mapping and Pseudo-Time Series Analysis

Subsequently, the expression of the six genes in the model was analyzed in the cell type. As shown in Figure 6A–F, BIRC3 is mainly expressed in T cells, CAV1 is mainly expressed in Epithelial cells, Fibroblast cells and Pericyte cells, and CD109 and DCBLD2 are mainly expressed in Epithelial cells. IL1RAP is mainly expressed in monocyte and macrophage cells, while SP100 is highly expressed in different cells. In this study, the immune cells were screened out and simulated time sequence analysis was performed. As shown in Figure 6G–L, cells gradually differentiated from dark blue to light blue, with a total of three differentiation states. cluster1 is the main cluster in differentiation state 1, cluster11 and cluster14 are the main cluster in differentiation state 2, cluster4 and cluster6 are the main cluster in differentiation state 3. Moreover, monocytes and macrophage are early differentiated, Mast cells and plasma cells are in the middle state, and NK cells, T cells and Dendritic cells are late differentiated. The model value SENRs in each cell was calculated according to the model calculation formula. As shown in Figure 6K, SENRs_Low was mainly in the early and middle differentiation state, while SENRs_High group was mainly in the late differentiation state. As shown in Figure 6L, with the differentiation of immune cells, the expressions of BIRC3, CAV1 and DCBLD2 showed a trend of first increasing and then decreasing, SP100 mainly showed an increasing trend, while CD109 and IL1RAP mainly showed a decreasing trend.

Figure 6 Gene expression mapping and pseudo-time series analysis. (A–F) The expression and localization of 6 genes in different cells. (G–L) pseudo-time series analysis of immune cells and quasi-temporal analysis of 6 model genes.

Clinical Correlation Analysis

For clinical guidance, we conducted clinical relevance analysis. As shown in Figure 7A, we found that model values were prognostic risk factors in both the three cohorts and the combined data (p<0.05). As shown in Figure 7B–D, patients were divided into 3 groups according to the median SENRs value and N stage of the model. Group1 and Group2 corresponded to SENRs_Low group, while Group3 corresponded to SENRs_High group. Survival curve analysis suggested that Group3 had a worse prognosis than Group1. As shown in Figure 7E, by combining the model values of the patient’s TCGA-FB-A5VM and clinical data, the 1, 3 and 5-year mortality of the patient was predicted to be 0.226, 0.676 and 0.787 in the nomogram. In order to further verify the accuracy of the prediction results of Nomograph, continuous prognostic ROC curve analysis was performed in this study. As shown in Figure 7F, Nomogram predicted the prognosis of patients with AUC fluctuation of 0.75, which was higher than other clinical indicators. In order to further find the core genes among the 6 genes in the model, CAV1 interacts with CD109, DCBLD2 and BIRC3, as shown in Figure 7G. Therefore, this study further verified the gene by experiment.

Figure 7 Clinical correlation analysis. (A) The model values were prognostic risk factors in both the three cohorts and the combined data (p<0.05). (B–D)Clinical decision analysis. (E and F) Construction and validation of the nomogram. (G)Interaction networks between proteins encoded by model genes.

The Function of CAV1 in Pancreatic Cancer Cell Lines Was Verified by Cell Experiments

Because CAV1 plays a pivotal role in the protein-protein interaction network of six model genes, we conducted further cell line experiments to verify the function of CAV1. Firstly, CAV1 was knockdown with 3 siRNA in pancreatic cancer Capan-1 cell lines and PANC-1 cell lines. The results show that the knockdown efficiency of si-CAV1-2 is higher(Figure 8A). Then CCK-8 experiment was performed(Figure 8B). The results showed that the activity of Capan-1 and PANC-1 cell lines decreased significantly after CAV1 gene knockdown. And Colony formation assay showed that the proliferation activity of Capan-1 and PANC-1 cell lines decreased significantly after CAV1 gene knockdown(Figure 8C). Transwell assay showed that the migration ability of Capan-1 and PANC-1 cell lines decreased significantly after CAV1 gene knockdown(Figure 8D).

Figure 8 The function of CAV1 gene in pancreatic cancer cell line was verified by cell function assay. (A) CAV1 knockdown with 3 siRNA in pancreatic cancer Capan-1 cell lines and PANC-1 cell lines. The results show that the knockdown efficiency of si-CAV1-2 is higher. (B) CCK-8 experiment. The results showed that the activity of Capan-1 and PANC-1 cell lines decreased significantly after CAV1 gene knockdown. (C) Cloning experiments. The results showed that the proliferation activity of Capan-1 and PANC-1 cell lines decreased significantly after CAV1 gene knockdown. (D) Transwell experiment. The results showed that the migration ability of Capan-1 and PANC-1 cell lines decreased significantly after CAV1 gene knockdown.(*p<0.05, **p<0.01, ***p<0.001).

PCR Tests of CAV1 in Clinical Samples

PCR experiments showed that CAV1 was significantly up-regulated in pancreatic cancer tissues compared with adjacent normal tissues (Figure 9).

Figure 9 PCR test. PCR experiments showed that CAV1 was significantly up-regulated in pancreatic cancer tissues compared with adjacent normal tissues. (***p<0.001).

Discussion

Our study is the first to explore the significance of cellular senescence through single-cell sequencing analysis and bulk transcriptome analysis in pancreatic cancer. Through single-cell sequencing analysis, we can intuitively identify the heterogeneity of cell senescence distribution in pancreatic cancer cells. The pathways associated with cell senescence such as Apoptosis, Epithelial Mesenchymal Transition, Hypoxia, Inflammatory Response and TGF Beta were also found. At the same time, WGCNA was used to explore the gene modules most associated with cellular senescence in pancreatic cancer. The green module was found to be highly correlated with cellular senescence. This is not only a way to identify gene sets, but also a means to stratify the prognosis of pancreatic cancer patients. With unsupervised clustering, patients could be divided into two clusters based on the module genes obtained from WGCNA, with cluster2 having a significantly worse prognosis. This also provides a reference for us to evaluate the prognosis of pancreatic cancer patients.

The further construction and validation of the prognostic model is more clinical significance. The initial screening of these model genes was obtained by the intersection of genes obtained by single-cell sequencing analysis and genes obtained by WGCNA, with a total of 303 genes. These genes were double-screened, so their relevance to cellular senescence was more pronounced. Further COX and LASSO regression analyses narrowed the gene range further and assigned coefficients to each model gene so that the patient’s risk profile could be quantified by the amount of these genes expressed. Each group of patients could be divided into a high risk group and a low risk group by median risk value, and the high risk group had a significantly worse prognosis. This provides a better method for stratifying the prognosis of pancreatic cancer. The subsequent immune correlation analysis and mutation analysis showed the immune landscape and mutation landscape of pancreatic cancer, and analyzed the correlation between risk score and immunity and mutation, which could provide references for immunotherapy of pancreatic cancer. Overall, the results of our study can provide reference for the prognosis assessment and immunotherapy of pancreatic cancer.

At present, some bioinformatics studies have initially revealed the significance of senescence in the tumor microenvironment. Lv et al identified three molecular subgroups of osteosarcoma through differentially expressed age-related genes in prognosis. These subgroups have differences in immune invasion, signaling pathway activation, prognosis, and drug sensitivity, which can provide considerable references for the diagnosis and treatment of osteosarcoma.²⁶ Yue et al constructed senescence-related risk profiles in colorectal cancer, where the high risk score was associated with poorer prognosis and positively correlated with the abundance of immune cell infiltration in the tumor microenvironment.²⁷ Hua et al constructed a 5-gene senescence-related prognosis signature from high-throughput sequencing data to assess patient prognosis in pancreatic cancer.²⁸ However, compared to their results, we used single-cell sequencing data to explore the heterogeneity of senescence in pancreatic cancer cells, adding confidence to the data. Moreover, we also used weighted coexpression network analysis, which further increased the correlation of the data. Moreover, the prognostic model we constructed has a higher AUC value and is more accurate. We also used cell experiments to verify the function of the key gene CAV1 in pancreatic cancer cells.

In the realm of molecular biology and cancer research, the caveolin-1 gene, often abbreviated as CAV1, has emerged as a key player with profound implications in tumorigenesis and cancer progression.²⁹ CAV1 is a gene that encodes a critical structural protein known as caveolin-1, which is primarily associated with the formation and maintenance of caveolae, small invaginations of the cell membrane.³⁰ Beyond its structural role, CAV1 has been recognized as a multifunctional molecule involved in various cellular processes, including signal transduction, lipid homeostasis, and cellular trafficking. In recent years, an increasing body of evidence has underscored the pivotal role of the CAV1 gene in cancer biology, making it a subject of intense investigation and a promising target for cancer therapy. The significance of CAV1 in pancreatic cancer has been preliminarily elucidated by many studies. Hu et al found that cav-1 in pancreatic cancer cells stimulates IGF1R/IR and glycolysis in cancer cells and triggers a malignant state in tumor carriers.³¹ Yamao et al found that in pancreatic cancer, down-regulation of CAV1 in cancer-associated fibroblasts (CAFs) reduced the invasiveness of PC cells.³² Therefore, CAV1 of CAF may be a new target for the treatment of PC. Our study revealed the correlation between CAV1 and senescence in pancreatic cancer, and verified the function of CAV1 through cell experiments, providing a potential target for pancreatic cancer.

A limitation of our study is that our single-cell sequencing data and bulk transcriptome data were obtained from public databases, and we lacked sequencing data from the real world. Moreover, our experiments were carried out in cell lines, lacking corresponding patient tissues and animal models for verification. We will explore that in the future.