Prognostic implications and characterization of tumor-associated tertiary lymphoid structures genes in pancreatic cancer

Abstract

Background

Pancreatic ductal adenocarcinoma (PDAC) is among the most aggressive cancers, with rising incidence and limited responsiveness to immunotherapy due to its highly suppressive tumor microenvironment (TME). Tertiary lymphoid structures (TLS), ectopic formation structures of immune cells, are linked to better prognosis and improved immunotherapy responses in PDAC. Understanding TLS’s role in PDAC could enhance immunotherapy effectiveness.

Methods

This study integrated transcriptomic and clinical data from 310 PDAC patients in GEO database. We performed consensus clustering using tumor-associated TLS (TA-TLS) genes, identifying three distinct molecular subtypes. Single-sample gene set enrichment analysis (ssGSEA) was then employed to calculate a TLS score for each patient, allowing for TLS-based evaluation. Key prognostic genes were identified using an iterative LASSO method, leading to the construction of a risk assessment model, which was validated across independent cohorts. We further analyzed the TLS score using single-cell RNA sequencing (scRNA-seq), visualized key gene expression, and validated protein expression through immunohistochemistry (IHC). Additionally, we explored the effects of DNASE1L3 on cell proliferation and migration, and its immune-related functions using Gene Set Enrichment Analysis (GSEA) and multiplex cytokine analysis.

Results

Consensus clustering revealed three PDAC molecular subtypes with significant differences in prognosis, TA-TLS gene expression, and TME features. The TLS score effectively stratified patients into high and low groups, correlating with survival outcomes and TME characteristics. Our risk model, validated across cohorts, reliably predicted patient outcomes. Validation studies showed lower expression of DNASE1L3 and IL33 in tumor tissues. scRNA-seq confirmed TLS score associations with immune cells. DNASE1L3 overexpression inhibited PDAC cell proliferation and migration, with cytokine analysis indicating increased immune activity.

Conclusions

This study elucidated the expression profile of TA-TLS genes in PDAC, constructed a TLS gene-based scoring system, and developed a related risk model. We also explored the functions and potential antitumor mechanisms of key genes, providing evidence and new insights for enhancing TLS-targeted immunotherapy strategies in PDAC.

Background

PDAC is one of the most challenging malignancies to treat, with the worst prognosis [1]. Its incidence has been increasing in recent years. Despite the diversification of treatment strategies for PDAC, overall therapeutic outcomes remain extremely limited. On one hand, PDAC progresses rapidly, with many patients experiencing distant metastases shortly after surgery [2]. On the other hand, the unique biological behavior of PDAC, characterized by its complex and disordered TME, contributes to its high invasiveness and metastatic potential [3, 4]. These challenges make it imperative to develop new therapeutic approaches and to precisely stratify patients by risk to maximize treatment benefits.

Immunotherapy has shown exciting efficacy in the treatment of various solid tumors, especially with the advent of immune checkpoint inhibitors (ICIs) [5]. However, due to the immunosuppressive state of the PDAC tumor microenvironment, the response to immunotherapy is generally low, and it remains difficult to improve patient prognosis [6, 7]. Key research directions in immunotherapy include identifying PDAC patients who are sensitive to immunotherapy and overcoming the TME’s immunosuppressive state to enhance PDAC’s immunogenicity and improve immune response.

In recent years, research on TLS within the tumor microenvironment has garnered significant attention. TLS, similar to secondary lymphoid organs, are pathological lymphoid structures that emerge in tissues with chronic inflammation and in tumors [8, 9]. These structures are primarily composed of CD20 + B cells at their center, surrounded by CD3 + T cells, and include various other immune cells and structural components such as dendritic cells (DCs), follicular dendritic cells (FDCs), fibroblastic reticular cells (FRCs), stromal cells, and high endothelial venules (HEVs). Immune cells and cytokines within TLS drive anti-tumor immune responses [10]. Both antibody responses and T-cell responses can be initiated within TLS, where B cells and T cells contribute to humoral and cellular immunity, respectively, exerting localized anti-tumor effects that inhibit tumor growth [11]. Clinically, the presence of TLS in most solid tumors is associated with better prognosis and a stronger response to immunotherapy [12, 13]. TLS are emerging as potential biomarkers for predicting clinical outcomes and the effectiveness of immunotherapy. Understanding the mechanisms by which TLS participate in anti-tumor immune responses and identifying key regulatory factors involved in TLS formation are crucial for developing related immunotherapies.

In PDAC, the presence of TLS is believed to correlate with better prognosis, and genes associated with TLS are becoming a focal point for inducing TLS formation to enhance the effectiveness of immunotherapy in PDAC [14, 15]. Beyond the histological identification of TLS components and characteristics, comparing and analyzing the genetic features of TLS offer an alternative perspective to elucidate the clinical value and underlying mechanisms of TLS [16,17,18]. Previous studies have summarized TLS-related genes in solid tumors, gathering TA-TLS genes and revealing strong correlations with the immune microenvironment [19]. Further analysis of these genes is instrumental in explaining the role of TLS and their potential research value from a genomic perspective. However, these genetic characteristics have not yet been further analyzed or validated in PDAC. Therefore, identifying key TLS-related genes in PDAC and conducting a genomic quantitative analysis of TLS-associated features will help to further clarify the mechanisms by which TLS influence PDAC prognosis.

Our study provides a comprehensive analysis of the expression profile of TA-TLS genes in PDAC (pancreatic ductal adenocarcinoma) and integrates TLS scoring and risk prognostic models for multidimensional bioinformatics analysis, filling the gap in the understanding of TLS gene characteristics and the tumor microenvironment in PDAC. By constructing a scoring system based on TLS-related genes, the study successfully stratified PDAC patient cohorts, revealing a close relationship between TLS genes and tumor microenvironment features, and identifying key genes with significant prognostic value. Furthermore, our research employs multi-omics data integration to elucidate the potential mechanisms of TLS genes in PDAC. In order to clearly present the ideas of our work, the flowchart of the research is presented here (Fig. 1). Through these analyses, we provide new insights and directions for future immune therapeutic strategies, clinical diagnostics, and prognostic evaluation in PDAC, offering a novel perspective for advancing tumor immunity research related to TLS.

Materials and methods

Bulk RNA sequencing data of PDAC and processing

We extracted bulk RNA sequencing data and corresponding survival prognosis information, primarily overall survival (OS), from four pancreatic cancer cohorts (GSE21501 [20], GSE57495 [21], GSE62452 [22], GSE71729 [23]) in the Gene Expression Omnibus (GEO) database. We then merged the expression matrices and removed batch effects between datasets using the ‘Combat’ method from the ‘SVA’ R package. Additionally, we collected PDAC transcriptomic data and clinical information from the TCGA database and other GEO datasets, including GSE85916. To ensure the consistency and correspondence of patient samples, we filtered out samples lacking complete survival data.

Single-cell RNA sequencing data analysis

scRNA-seq analysis was performed on previously published data (GSE212966 [24]) using the ‘Seurat’ package (version 4.3.1). To ensure relative accuracy, mitochondrial genes were excluded, and the top 2000 variable features of the dataset were identified. Batch effects were eliminated using the ‘FindIntegrationAnchors’ and ‘IntegrateData’ methods. The integrated dataset underwent further dimensionality reduction by principal component analysis (PCA) and uniform manifold approximation and projection (UMAP) for visualization. Cell annotation was based on the original descriptions in the literature.

Consensus clustering and gene expression analysis

Based on the characteristic TA-TLS genes summarized by Fridman et al. [19], we performed consensus unsupervised clustering analysis on the PDAC samples using the ‘ConsensusClusterPlus’ R package to classify them into different molecular patterns [25]. Differential gene expression analysis was conducted using the ‘limma’ package [26].

TLS score estimation using ssGSEA

We calculated the TLS score for each patient using single-sample gene set enrichment analysis (ssGSEA [27]) based on the characteristic TA-TLS gene identified through univariate Cox regression analysis. According to the scoring results, the PDAC cohort was divided into two groups: TLS_H and TLS_L.

Construction of risk score using iterative LASSO

We developed a gene expression signature model for predicting prognosis from a pancreatic cancer cohort of 310 samples. The conditions set included a variance in expression level greater than 0.5 and a univariate Cox proportional hazards analysis P-value (Wald test for predictive potential) less than 0.5. The gene set was subjected to penalized multivariate Cox proportional hazards survival modeling using a variable selection algorithm based on L1 regularization (LASSO) [28].

Support vector machine (SVM)-based survival analysis

A prognostic model using SVM was constructed to validate the significance of DNASE1L3, RAMP2, and IL33. RNA-seq data and clinical survival information were used as inputs, with z-score normalization ensuring comparability. The Survival-SVM algorithm with a linear kernel was applied, optimizing hyperparameters via cross-validation [29]. Predictive accuracy was evaluated using the concordance index (C-index) and time-dependent ROC analysis. Patients were stratified into high- and low-risk groups based on median risk scores, with Kaplan-Meier survival analysis and log-rank tests assessing group differences. Analyses were implemented using scikit-survival and lifelines in Python.

Drug sensitivity prediction

We predicted the chemotherapy response for each sample using the Genomics of Drug Sensitivity in Cancer 2 (GDSC2), the largest publicly available pharmacogenomics database (https://www.cancerrxgene.org/). The prediction was carried out with the ‘oncoPredict’ R package, employing the calcPhenotype function to estimate tumor sample drug sensitivity based on cell line expression profiles [30].

Immunohistochemistry (IHC)

Using three key genes identified through the iterative LASSO method, we performed IHC identification and analysis on 20 pairs of PDAC tumor and matched para-cancerous paraffin-embedded samples. All samples were provided by the Hepatobiliary and Pancreatic Surgery Department of Peking University First Hospital following guidelines from the Ethics Committee of Peking University First Hospital. Each patient provided written informed consent prior to participation in the study. The paraffin samples for immunohistochemistry underwent deparaffinization and antigen retrieval, were washed, and then incubated with the primary antibody overnight at 4 °C. The next day, they were incubated with HRP-conjugated secondary antibodies at room temperature for 30 min. Finally, the slides were sealed with neutral resin, covered with coverslips, and the expression of target proteins was analyzed by two independent pathologists. Staining intensity was determined by the proportion of stained cells (0 = 0%, 1 ≤ 25%, 2 = 26 to 50%, 3 = 51 to 75%, 4 = 76 to 100% positive). The antibodies used were DNASE1L3 (1:400 dilution), RAMP2 (1:200 dilution), and IL-33 (1:800 dilution), all purchased from Proteintech (67041-1-Ig, 13223-2-AP, 12372-1-AP).

Cell culture and transfection

Human pancreatic cancer cell lines, including PANC-1, MiaPaCa-2, AsPC-1, T3M4, Patu8988, and normal pancreatic duct cells (hTERT-HPNE), were purchased from the American Type Culture Collection (ATCC). All cells were cultured in DMEM (Gibco, USA) with 10% fetal bovine serum (Gibco, USA) and 1% penicillin-streptomycin solution (Gibco, USA). The cells were maintained under normoxic conditions in a cell culture incubator (ThermoFisher, USA) with 5% CO2 and 20% O2. To construct DNASE1L3 overexpression stable cell lines, overexpression lentivirus and a negative control were acquired from Genechem (Shanghai, China). The DNASE1L3 gene (NM_004944.4) was sourced from Genechem’s cDNA library. The lentiviral vector plasmid GV492 (Ubi-MCS-3FLAG-CBh-gcGFP-IRES-puromycin), along with the DNASE1L3 gene sequence, was digested using AgeI and BamHI restriction enzymes, followed by cloning through the In-fusion recombination method. The recombinant vector was confirmed by DNA sequencing.

Lentivirus production involved transfecting the viral vectors and two helper plasmids, psPAX2 and pMD2.G, into 293T cells using Lipofectamine 2000 (Invitrogen; Thermo Fisher Scientific, Inc.). Seventy-two hours post-transfection, infectious lentiviruses were harvested, centrifuged to remove cell debris, and filtered through 0.45 μm cellulose acetate filters. The virus titer, determined by fluorescence-activated cell sorting of GFP-positive 293T cells, was approximately 2.5 × 10^8 transducing units (TU)/mL, and stored at -80 °C for further use.

Western blot

For protein extraction, cells were washed with PBS and lysed using RIPA buffer (Solarbio, Beijing, China) containing PMSF. The lysate was centrifuged at 12,000 g for 5 min at 4 °C to remove cell debris. Equal amounts of protein were then loaded onto SDS-PAGE gels and transferred to 0.22 μm PVDF membranes (Millipore, Massachusetts, USA). The membranes were blocked with skim milk, incubated overnight at 4 °C with primary antibodies, and washed three times with TBS-T. Afterward, they were incubated at room temperature for 2 h with HRP-conjugated secondary antibodies. Protein bands were visualized using the enhanced chemiluminescence (ECL) method. The primary antibodies used were DNASE1L3 (Proteintech, #67041-1-Ig, 1:2000, USA) and β-tubulin (Proteintech, # 10094-1-AP, 1:5000, USA).

Cell proliferation and migration experiment

Cell viability was assessed using colony formation and cell proliferation assays. In the cell proliferation assay, control and transfected cells were cultured in 96-well plates (3000 cells/well), with five replicates per group. Viability was measured every 24 h. After adding 10 µl of CCK-8 solution to each well, the plates were incubated at 37 °C for 1.5 h, and absorbance was measured at 450 nm. For the colony formation assay, control and transfected cells were seeded into six-well plates and cultured for 2 weeks in medium containing 10% FBS. Colonies were then fixed with paraformaldehyde and stained with crystal violet (Solarbio, #G1062, China).

In the Transwell migration assay, 1 × 10^5 cells were seeded into the upper chamber with 200 µL of serum-free medium (8.0 μm pore size, Corning, USA). The lower chamber contained 700 µL of medium with 10% FBS as a chemoattractant. After 24 h of incubation, non-invading cells were removed from the upper chamber, which was then fixed with paraformaldehyde, stained with crystal violet, and photographed.

Luminex assay

We employed a previously developed and validated method to detect cytokine levels in cell culture supernatants using Luminex liquid chip technology (Luminex 200 system). This assay detects 48 cytokines, including β-NGF, CCL27, CCL11, FGF-basic, G-CSF, GM-CSF, GRO-α (Gro-α/KC/CXCL1), HGF, IFN-α2, IFN-γ, IL-1α, IL-1β, IL-1Rα, IL-2, IL-3, IL-4, IL-5, IL-6, IL-7, IL-8/CXCL8, IL-9, IL-10, IL-12(p40), IL-12(p70), IL-13, IL-15, IL-16, IL-17 A, IL-18, IP-10/CXCL10, LIF, MCP-1/CCL2, MCP-3/CCL7, MIP-1α/CCL3, MIP-1β, MIF, PDGF-BB, RANTES, SCF, SCGF-β, SDF-1α, TRAIL, TNF-α, TNF-β, and VEGF-A. The supernatant was diluted twofold, and 50µL was added to each well for analysis. To ensure reproducibility, masked samples were included on each plate and tested in triplicate. The quality criteria for the experiment were: (1) a standard curve with an r² > 0.99; (2) recovery within the linear range of the standard curve (recovery = measured concentration/set concentration * 100%) falling within 100 ± 30%; and (3) a coefficient of variation (CV = SD/mean) within the linear range of the standard curve not exceeding 25%. The fluorescence values obtained were then applied to the standard curve equation to calculate the sample concentrations.

Fluorescent multiplex immunohistochemistry

Tissue slides were processed as follows: First, they were dewaxed with xylene and then rehydrated through a graded ethanol series. Endogenous peroxidase activity was blocked using 3% hydrogen peroxide (Sinopharm Chemical Reagent Co., China, #73113760) for 10 min. The slides were then washed with PBS, subjected to antigen retrieval by heating, and subsequently washed again with PBS. To prevent nonspecific binding, 5% BSA (Sigma, Shanghai, China, #B2064) was applied to each slide, and incubation was performed at room temperature for 20 min. Primary antibodies were then added to the slides in a volume of 100 µL each: DNASE1L3 (Proteintech, #67041-1-Ig, 1:400, USA), CD20 (Proteintech, #60271-1-Ig, 1:200, USA), and CD3 (Proteintech, #81324-1-RR, 1:200, USA). The slides were incubated overnight at 4 °C.The following day, the slides were washed with PBS and incubated with a labeled secondary antibody (Abcam, Goat anti-rabbit IgG H&L (HRP), 1:2000, ab205718) at 37 °C for 30 min. After another PBS wash, each section was treated with 100 µL of Try-488 Tyramine Conversion Reagent (Runnerbio, Bry-try488) and incubated for 10–30 min at room temperature. Finally, the sections were mounted with an anti-fluorescence quenching sealant containing DAPI (Beyotime Biotechnology, China, P0131).

Co-culture experiment

Jurkat cells were seeded in the upper chamber of a 0.4 μm pore size transwell, while DNASE1L3-overexpressing PANC-1 cells and NC PANC-1 cells were seeded in the lower chamber. The cells were co-cultured for 48 h in DMEM medium containing 10% fetal bovine serum and 1% penicillin-streptomycin, maintained at 37 °C in a 5% CO2 environment.

Statistical analysis

We conducted R software (version 4.3.3) and GraphPad Prism (version 10) for data analysis and visualization. Statistical significance for normally distributed quantitative data was assessed using student’s t-test, while the Wilcoxon test were applied for non-parametric comparisons. For all survival analyses, we utilized the Kaplan-Meier method from the R package “Survminer” to evaluate patient prognosis. P-value < 0.05 was considered statistically significant in this study.

Results

Molecular subtyping and characterization of tertiary lymphoid structures genes in pancreatic cancer at the transcriptomic level

We conducted a transcriptomic analysis to subtype and TA-TLS genes in pancreatic cancer. Utilizing the 39 TA-TLS genes summarized by Fridman [19] (Additional file 1: Table S1), we performed unsupervised consensus clustering on 310 PDAC clinical samples from the GEO-meta datasets (GSE21501, GSE57495, GSE62452, GSE85916). Our objective was to explore the expression profiles and prognostic value of TLS-related genes in PDAC. The consensus clustering analysis stratified the 310 samples into three distinct molecular subtypes (Cluster number is set to 3): Cluster A (83 samples), Cluster B (160 samples), and Cluster C (67 samples) (Fig. 2A, Additional file 1: Table S2). Principal component analysis (PCA) was employed to visualize the characteristics of each molecular subtype, revealing distinct boundaries between the three subtypes and indicating marked differences in the molecular features among the TLS clusters (Fig. 2B). Clinical prognosis analysis across the three groups demonstrated that TLSclusterC was associated with the poorest prognosis, in contrast to TLSclusterA and TLSclusterB (Fig. 2C). Correspondingly, a heatmap depicting TA-TLS gene expression in the three subtypes showed that TLSclusterC, which had the poorest prognosis, exhibited the lowest expression of TA-TLS genes. In contrast, TLSclusterA and TLSclusterB, associated with better prognoses, exhibited progressively higher expression levels (Fig. 2D). Furthermore, we quantified the immune cell infiltration characteristics of each patient using the ssGSEA algorithm to analyze the relationship between TLS gene characteristic subgroups and the tumor immune microenvironment comprehensively (Fig. 2E, Additional file 1: Table S2). TLSclusterA and TLSclusterB, associated with better prognoses, exhibited higher immune cell scores across most immune cell types, including various T cells, B cells, monocytes, and macrophages, suggesting a more robust anti-tumor immune function in these clusters. Notably, dendritic cell scores were significantly elevated in TLSclusterC, indicating a potentially detrimental role of DCs in this TLS molecular pattern. Subsequently, we conducted a differential analysis of biological processes based on previously summarized TME characteristics [31] for the TLS molecular patterns (Fig. 2F). The results indicated that positive tumor immune features, such as antigen processing machinery, CD8 + T effector cells, immune checkpoints, TMEscoreA, and TMEscoreB, were significantly elevated in TLSclusterA and TLSclusterB. In contrast, epithelial-mesenchymal transition (EMT) characteristics were markedly increased in TLSclusterC. These findings suggest that the tumor microenvironment corresponding to the three distinct TLS molecular patterns exhibits significant heterogeneity.

Establishment of a scoring model and analysis of tumor microenvironment characteristics

To further quantify and analyze the heterogeneity across different TLS molecular patterns, we identified differentially expressed genes among the three patterns and visualized them using a Venn diagram, revealing a total of 161 intersecting genes (Fig. 3A, Additional file 1: Table S4). Next, we conducted univariate Cox regression analysis on these intersecting genes based on prognostic outcomes, identifying 20 signature genes (Fig. 3B, Additional file 1: Table S5). Gene Ontology (GO) enrichment analysis of these 20 prognostic differential genes indicated significant enrichment in processes such as glycoprotein regulation, antigen processing and presentation, and dendritic cell immune regulation (Fig. 3C). Utilizing these 20 signature genes, we applied the ssGSEA method to score the 310 PDAC samples from the GEO-meta cohort. Each sample was assigned a TLSscore, allowing for the stratification of patients into high-score (TLS_H) and low-score (TLS_L) groups (Fig. 3D, Additional file 1: Table S6). Survival analysis demonstrated that the TLS_H group exhibited significantly better prognostic outcomes compared to the TLS_L group (Fig. 3E). For the TLSscore, we further validated the expression of 39 TA-TLS genes between the high-score and low-score patient groups. The results demonstrated that most TA-TLS genes exhibited significantly higher expression levels in the TLS_H group compared to the TLS_L group (Figure S1A). This finding indicates a strong concordance between the scoring system based on 20 signature genes and the expression of TA-TLS genes.

The relationship between TLSscore and the different TLS clusters has been further analyzed. TLSclusterC, associated with the poorest prognosis, exhibited the lowest TLSscore, while TLSclusterA, linked to the best prognosis, had the highest score, followed by TLSclusterB (Fig. 4A). A Sankey diagram was used to illustrate the relationships between TLS clusters, TLSscores, and clinical outcomes. This visualization revealed that TLSclusterA predominantly corresponded to the TLS_H group, whereas TLSclusterC was primarily associated with the TLS_L group. Notably, most patients in the TLS_L group had deceased, while the TLS_H group had a slightly higher proportion of surviving patients (Fig. 4B). Using the ‘ESTIMATE’ package, we performed single-sample gene set enrichment analysis to determine tumor purity in PDAC patient expression data. This analysis generated stromal scores (indicating the presence of stromal components in the tumor tissue), immune scores (reflecting the extent of immune cell infiltration), and ESTIMATE scores (inferring tumor purity) (Additional file 1: Table S7). A heatmap of these results revealed that various immune-related signatures were significantly higher in the TLS_H subgroup compared to the TLS_L subgroup (Fig. 4C). To further illustrate the scoring results, violin plots were used to compare the Stromal Score, Immune Score, and ESTIMATE Score, demonstrating that relevant tumor microenvironment scores were significantly elevated in the TLS_H subgroup relative to the TLS_L subgroup (Fig. 4D).

Development of prognostic signature genes of TLS

We used 310 PDAC samples as the training set for Lasso-penalized multivariate Cox proportional hazards modeling. After 500 iterations, genes with the highest frequency of survival prediction (greater than 50) were selected. The gene set corresponding to the peak area under the curve (AUC) was identified as the final signature, resulting in the selection of three key prognostic genes (Fig. 5A, B). The formula for the prediction model is as follows:

Signature = DNASE1L3× (− 0.21624559) + RAMP2× (− 0.19559646) + IL33× (− 0.03507206).

We generated Kaplan-Meier survival curves for the high- and low-risk groups within the training set, which demonstrated that the 310 PDAC patients could be effectively stratified into distinct risk groups with significant differences in prognosis (Fig. 5C, Additional file 1: Table S8). The predictive power of the model was further validated in additional datasets, including the TCGA-PDAC cohort and the GSE85916 cohort, where it consistently stratified patients into high- and low-risk groups with significant prognostic differences (Fig. 5D, E). Furthermore, we constructed a prognostic model using the aforementioned key prognostic genes through an SVM-based approach to validate their significance (Additional file 1: Table S9). The SVM prognostic model effectively stratified patients; however, its performance in ROC analysis was inferior to that of the Lasso-penalized multivariate Cox proportional hazards model (Figure S1B, C). We also analyzed the expression levels of immune checkpoint-related genes across the two risk groups, based on our prior investigation into tumor immunity. The analysis revealed that most immune checkpoint genes were expressed at lower levels in the high-risk group and higher levels in the low-risk group, suggesting a potential correlation between risk stratification and the corresponding tumor immune response (Fig. 5F). Given the critical role of chemotherapy in treating pancreatic cancer, we also sought to predict chemotherapy response within the risk groups identified through our stratification model. Using the ‘oncopredict’ package, we trained models to assess drug sensitivity. Two-sample testing identified several drugs with significant differences in sensitivity between the groups (Additional file 1: Table S10). Based on logFC values for drug sensitivity and their statistical significance, we highlighted six chemotherapy drugs with the most relevance: Vinblastine_1004, Nutlin-3a_1047, Mitoxantrone_1810, Vinblastine_1818, Docetaxel_1819, and Vinblastine_2048 (Fig. 5G).

Single-cell analysis reveals TA-TLS gene expression characteristics in the PDAC tumor microenvironment

To further explore the molecular characteristics and heterogeneity associated with the TLSscore and the identified prognostic signature genes, we conducted an analysis of scRNA-seq data. This approach allowed us to examine the expression patterns of the TLSscore and the prognostic signature genes (DNASE1L3, RAMP2, and IL33) across various cell populations within the tumor microenvironment at a single-cell resolution. We integrated scRNA-seq data from six PDAC samples and performed stringent cell filtering and data quality control on the integrated dataset (Fig. 6A), and resulting in a comprehensive cell-gene matrix encompassing 31,503 cells. After rigorous quality control using Seurat and harmony tools, we performed clustering analysis and visualized the data using UMAP. The cells were annotated as T cells (CD3D), fibroblasts (LUM), ductal cells (KRT19), macrophages (CD68), B cells (MS4A1), stellate cells (RGS5), endothelial cells (CDH5), mast cells (TPSAB1), neutrophils (S100A8), plasma cells (MZB1), Schwann cells (S100B), and endocrine cells (CHGA) (Fig. 6B). The cell-type composition within each sample was also displayed, revealing significant variability among samples but consistently showing a predominance of T cells, fibroblasts, ductal cells, and macrophages (Fig. 6C). To assess TLS activity at the single-cell level, we calculated the TLSscore for each cell by integrating the TLSscore with the single-cell data. Using AUCell, we evaluated the gene sets comprising the TLSscore, finding the highest scores in endocrine cells, followed by T cells, B cells, plasma cells, and macrophages, with the lowest scores observed in fibroblasts and ductal cells (Fig. 6D). A similar evaluation using ssGSEA showed concordant results, with endocrine cells achieving the highest scores, followed by plasma cells, Schwann cells, macrophages, T cells, and B cells, while fibroblasts and ductal cells again exhibited the lowest scores (Fig. 6E). We also analyzed the expression levels of the three prognostic signature genes—DNASE1L3, RAMP2, and IL33—at the single-cell level (Fig. 6F-H). We statistically analyzed the expression levels of DNASE1L3, RAMP2, and IL33 in each lineage cell cluster and found that DNASE1L3, RAMP2, and IL33 were mainly expressed in the endothelial cell cluster (Figure S1A-C). This expression pattern underscores the potential roles of these genes in the endothelial and immune compartments of the PDAC TME.

Analysis and identification of TLSscore-related PDAC prognostic genes

We conducted a comprehensive analysis to identify and characterize the PDAC prognostic genes DNASE1L3, RAMP2, and IL33, focusing on their correlation with immune cells, prognostic significance, and expression profiles in the context of TLSscore. Using ssGSEA, we evaluated the relationship between these genes and immune cell infiltration in PDAC transcriptome data. All three genes demonstrated positive correlations with the majority of immune cells (Fig. 7A-C). Specifically, DNASE1L3 showed the strongest correlation with activated B cells, RAMP2 with mast cells, and IL33 with eosinophils. Further detailed analysis revealed significant positive correlations between DNASE1L3 and activated B cells, as well as between RAMP2 and IL33 with mast cells (Fig. 7D-F). To explore the prognostic relevance of these correlations, we conducted survival analysis in the GEO-meta cohort (Fig. 7G-I). The analysis indicated that higher expression levels of DNASE1L3, RAMP2, and IL33 were associated with significantly better prognoses compared to lower expression levels, suggesting these genes may play tumor-suppressive roles in PDAC. Next, we examined the expression levels of these genes in the GSE71729 pancreatic cancer cohort. The results showed that DNASE1L3 and IL33 were significantly more expressed in normal tissues than in tumor tissues, while RAMP2 exhibited no significant difference in expression between normal and tumor tissues (Fig. 8A-C). We utilized the transcriptomic data of pancreatic cancer from TCGA to classify DNASE1L3 expression into High and Low groups. Statistical tests were conducted on associated pathological information, including age, gender, and staging. The results revealed that the age in the DNASE1L3-High group was significantly lower than that in the DNASE1L3-Low group (Figure S2A). Regarding gender, there was no significant difference between the DNASE1L3-High and Low groups (Figure S2B). In the pathological staging analysis, although the overall difference between the DNASE1L3-High and Low groups was not statistically significant, it was observed that early-stage patients (IA + IB + IIA) were relatively fewer in the DNASE1L3-High group, and all stage III and IV patients were present in the DNASE1L3-High group (Figure S2C). To further validate these findings, we collected 20 pairs of pancreatic cancer tissues and matched adjacent normal tissues from the Hepatobiliary and Pancreatic Surgery Department at Peking University First Hospital (Beijing, China). Immunohistochemical staining was performed for DNASE1L3, RAMP2, and IL33, followed by quantification of staining intensity (Fig. 8D, E). The results confirmed that DNASE1L3 and IL33 were expressed at lower levels in tumor tissues compared to normal tissues, whereas RAMP2 showed no significant difference (Figure S2D-F). These findings reinforce the potential significance of DNASE1L3 and IL33 as prognostic markers in PDAC.

DNASE1L3 suppresses PDAC proliferation and migration and is closely linked to tumor immunity

Our integrative analysis of survival data, bulk RNA expression, and IHC results from clinical samples identified DNASE1L3, a gene strongly associated with TA-TLS, as a potential tumor suppressor in PDAC. Given its lower expression in tumor tissues compared to normal tissues, we proceeded to validate DNASE1L3 expression levels in various cell lines. Western blot analysis confirmed that DNASE1L3 is significantly more highly expressed in the human normal hTERT-HPNE compared to several PDAC cell lines (Fig. 9A, Figure S2G). To elucidate the functional role of DNASE1L3 in PDAC, we established DNASE1L3 overexpression models in the PANC-1 and Mia PaCa-2 cell lines using LV-DNASE1L3 lentivirus. Validation at the protein level confirmed successful overexpression of DNASE1L3 (Fig. 9B). We then evaluated the impact of DNASE1L3 overexpression on PDAC cell proliferation through cell proliferation assays and colony formation assays. The results demonstrated that increased DNASE1L3 expression led to a significant reduction in the growth rate of PDAC cells (Fig. 9C) and a marked decrease in the number of colonies formed (Fig. 9D), indicating a substantial decrease in the cells’ proliferative capacity. Moreover, Transwell migration assays revealed that DNASE1L3 overexpression significantly impaired the migration ability of PDAC cells (Fig. 9E). To further understand the mechanisms underlying the anti-tumor effects of DNASE1L3, we performed GSEA. The analysis identified several immune-related signaling pathways that are closely associated with DNASE1L3 (Additional file 1: Table S11). The top three pathways, ranked by p-value, include the Chemokine signaling pathway, the cGMP-PKG signaling pathway, and the Cell adhesion molecules pathway (Fig. 9F). These findings suggest that DNASE1L3 exerts its tumor-suppressive effects not only by inhibiting cell proliferation and migration but also by engaging key immune-related pathways in the PDAC tumor microenvironment. To further investigate the role of this key gene in tumor immunity, we employed Luminex liquid chip technology for high-throughput detection of secreted cytokines in the culture supernatants of two DNASE1L3-overexpressing cell lines and their corresponding negative controls, and developed a standard curve for each cytokine. (Fig. 9G, Figure S3). The results revealed that, compared to the negative control group, the supernatants of the DNASE1L3-overexpressing cell lines exhibited elevated levels of numerous cytokines (Additional file 1: Table S12). Notably, cytokines that were consistently elevated in both DNASE1L3-overexpressing cell lines included HGF, IL-2Rα, IL-1α, GM-CSF, CCL27, IL-12, CCL11, MIG, IFN-α2, TRAIL, IL-4, SCF, IL-1Rα, FGF-basic, CCL7, G-CSF, CXCL10, TNF-α, and IL-1β. Other cytokines showed an increase in only one of the overexpressing cell lines, while a few exhibited minimal changes. Additionally, we examined the expression of DNASE1L3 in pancreatic cancer tissues from patients with and without TLS using multicolor immunofluorescence staining. The analysis revealed significantly higher expression of DNASE1L3 in pancreatic cancer tissues containing TLS compared to those without TLS (Figure S4A). To further investigate its functional impact, we co-cultured PANC-1 pancreatic cancer cells overexpressing DNASE1L3 and their corresponding negative control cells with Jurkat cells in a transwell chamber for 48 h. This co-culture experiment demonstrated notable alterations in immune exhaustion markers, with a marked reduction in the expression of LAG3 and TIM3 proteins in Jurkat cells co-cultured with DNASE1L3-overexpressing PANC-1 cells (Figure S4B).

Discussion

The prognosis and clinical translational research of pancreatic cancer remain focal points and challenges in current oncology studies. While immune therapies, particularly those targeting PD-1 and PD-L1, have expanded treatment options and significantly improved outcomes for various solid tumors, the benefits for pancreatic cancer patients have been limited [32,33,34]. Only a small subset of patients with microsatellite instability-high (MSI-H) pancreatic cancer have shown notable responses to immunotherapy, leaving the majority of patients with suboptimal outcomes [35]. Recently, the discovery of TLS in tumors, and their association with favorable clinical outcomes, has spurred a wave of clinical and experimental research. Emerging evidence suggests that TLS presence in pancreatic cancer is linked to better prognoses [36]. Moreover, the spatial distribution of TLS whether intratumoral or peritumoral might also influence clinical outcomes [14]. Notably, mature TLS are increasingly recognized as potential hubs for intratumoral immune responses in PDAC patients, highlighting the urgent need to elucidate the role of TLS in PDAC [37]. This includes detailed investigations at the genetic level and predictive studies. Despite some advances in understanding the pathological features of TLS, research in this area faces significant limitations, such as the lack of comprehensive expression profiling of TLS signature genes in PDAC and a limited understanding of their prognostic implications. This underscores the necessity of further clarifying the role of TLS in PDAC, particularly in determining the prognostic value of their signature genes and exploring deeper therapeutic strategies.

In this study, we initially compiled a set of 39 TA-TLS genes based on Fridman’s summary and used these to classify pancreatic cancer patients in the GEO-meta cohort into three clusters. These clusters exhibited significant differences in survival outcomes, gene expression profiles, and TME characteristics. Specifically, the TLSclusterC, which had the poorest prognosis, corresponded to the lowest expression levels of TA-TLS genes. This finding aligns with the established understanding that the presence of TLS in tumors is associated with improved clinical outcomes, a conclusion supported by the gene expression profiles in our integrated pancreatic cancer cohort. To further investigate the underlying mechanisms by which TA-TLS genes impact PDAC, we conducted a detailed analysis of immune cell characteristics and TME-related features across these molecular clusters. We observed that TLSclusterA and TLSclusterB, which are associated with better prognoses, generally had higher immune cell scores. Mature TLSs, as potential sites of tumor immune responses, are closely linked to favorable prognoses, with B cells and T cells—the main components of TLS—showing higher scores in these clusters, while TLSclusterC had the lowest scores, particularly for immature B cells, activated B cells, CD4 + T cells, and CD8 + T cells. Correspondingly, in the analysis of TME signature scores, the CD8 + T effector score was significantly higher in TLSclusterA and TLSclusterB, with the associated immune checkpoint score exhibiting a similar pattern. Conversely, TLSclusterC, with the poorest prognosis, showed greater activity in the EMT score. These results suggest that the expression of TLS-related genes can influence the production and activation of immune cells in PDAC. Moreover, our molecular classification based on TA-TLS genes effectively correlates with relevant characteristics of the tumor immune microenvironment. To further explore the immunological roles and prognostic potential of TA-TLS genes, we performed differential gene expression analysis across the three molecular clusters, identifying a total of 161 DEGs. We then selected 20 genes that were strongly associated with prognosis and analyzed their functions, finding them to be closely related to immune regulation, immune activation, and antigen presentation. Examples include pathways such as antigen processing and presentation, regulation of antigen processing and presentation, and the BCL-2 family protein complex.

Based on these findings, we proposed using these 20 prognostically significant genes as a gene set to evaluate the TLS score in PDAC patients. Using this gene set as the basis for ssGSEA calculation, we assessed the existing GEO-meta cohort and classified patients into two groups: TLS_H and TLS_L. Then we compared the survival outcomes, correlations with TLS molecular clusters, and TME characteristics between the two groups. The results revealed that individuals with higher TLS scores exhibited better prognoses, with their molecular profiles predominantly aligning with TLSclusterA and TLSclusterB. These individuals also displayed stronger immune cell infiltration characteristics and higher TME scores. This indicates that our TLSscore effectively correlates with survival outcomes and TME features, providing a valuable tool for evaluating the heterogeneity in the TME landscape driven by differences in TLS gene expression. Subsequently, we performed iterative LASSO regression on the gene set within the GEO-meta dataset, ultimately identifying three core genes including DNASE1L3, RAMP2, and IL33 around which we constructed a clinical prognostic model. Previous studies have reported associations between these genes and cancer. DNASE1L3 has been shown to enhance tumor immunity and inhibit tumor progression in various cancers. It exerts its antitumor effects through multiple mechanisms, such as repairing DNA damage and activating immune cells, as well as regulating downstream oncogenes via epigenetic modifications [38,39,40]. RAMP2 is likely associated with endothelial proliferation, vascular remodeling, and favorable cancer prognosis [41,42,43]. IL33 acts as a mediator in several tumors, playing a role in various complex aspects of tumor immunity [44,45,46]. However, research on these genes remains insufficient in pancreatic cancer, necessitating further investigation. Our prognostic model demonstrated effective stratification of survival outcomes within the GEO-meta cohort, successfully dividing the population into high- and low-risk groups. This model was also validated in two independent PDAC datasets, where it continued to show significant results. To validate the significance of TLS-related genes, we performed a side-by-side comparison using the SVM modeling approach. The results demonstrated that the SVM model, based on three key prognostic genes, effectively stratified patient risk. However, its ROC curve performance was slightly inferior to the iterative LASSO method. Nonetheless, this analysis provides strong support for the importance of the identified TLS-related genes. Additionally, we performed drug sensitivity screening based on the risk stratification provided by this model, identifying 179 chemotherapeutic agents that could be leveraged for personalized treatment following risk assessment.

Building on the risk prediction model that effectively evaluates the survival prognosis of PDAC patients, we conducted a detailed analysis and validation of the core genes DNASE1L3, RAMP2, and IL33, which are integral to the model. First, we performed single-cell level scoring analysis on the 20-gene panel used to construct the TLSscore. Both AUC and ssGSEA scores demonstrated the mapping of this gene panel at the single-cell level, with high-scoring cell populations predominantly concentrated in endothelial cells, plasma cells, T cells, and B cells. Given that TLS-associated cells and structures include T cells, B cells, and HEVs, these findings suggest that the 20 TLS-related prognostic genes we identified may play a significant role at the single-cell level. Additionally, we assessed the expression density of the three key genes, revealing that all three are highly expressed in endothelial cells, with DNASE1L3 also showing expression in macrophages and B cells. This indicates that these genes may exert certain immunological effects within endothelial cells in PDAC. Next, we examined the correlation between these three genes and immune cells at the transcriptome level. All three genes were positively correlated with a wide range of immune cells, particularly B cells and T cells. In survival analyses, DNASE1L3, RAMP2, and IL33 each exhibited a strong association with favorable clinical outcomes in PDAC. Expression analysis at the RNA level revealed that DNASE1L3 and IL33 were significantly downregulated in PDAC tissues compared to normal tissues, while RAMP2 showed no significant difference. To further validate these findings, we conducted protein-level assessments in our clinical tissue sample cohort. The results were consistent with the RNA findings: DNASE1L3 and IL33 were more highly expressed in normal pancreatic tissues, while their expression was reduced in tumor tissues. These observations suggest that the core TLSscore genes DNASE1L3 and IL33 may function as potential tumor suppressors in PDAC. To explore the functional roles of these core genes, we focused on DNASE1L3, which had a high weight in the prognostic model. Our experimental results demonstrated that DNASE1L3 overexpression inhibited the proliferation and migration of PDAC cells, further supporting its role as a tumor suppressor in PDAC. We further investigated the potential signaling pathways associated with this gene and found a strong connection with the Chemokine signaling pathway. Subsequent analysis of the supernatants from NC and DNASE1L3-overexpressing cell lines revealed a significant increase in the levels of numerous cytokines following DNASE1L3 overexpression. These cytokines play crucial roles in immune response, cell proliferation, activation, inflammation, and tissue repair. For instance, IL-2Rα, IL-1α, CCL27, IL-12, MIG, IL-4, CXCL10, and IL-1β are key regulators of T cell and B cell immune responses; GM-CSF and G-CSF promote myeloid cell production; and TRAIL and TNF-α are tumor necrosis factor-related cytokines. These findings suggest that DNASE1L3 may exert its anti-tumor effects by modulating cytokine production, thereby influencing the recruitment, regulation, and activation of immune cells. Additionally, we validated the correlation between TLS and DNASE1L3 using multicolor immunohistochemistry. We found that DNASE1L3 expression was relatively higher in patients with TLS compared to those without TLS, suggesting that DNASE1L3 may be associated with the presence of TLS. We also confirmed that overexpression of DNASE1L3 in pancreatic cancer cell lines can partially reverse the immune exhaustion of Jurkat cells, providing further insight into the immunological function of DNASE1L3. Previous studies in liver and colorectal cancers have also investigated and validated the anti-tumor effects of DNASE1L3, suggesting that DNASE1L3 may exert its effects by enhancing tumor immune responses [38, 39, 47]. In this study, we performed analogous validation in pancreatic cancer and confirmed the immune-related functions of DNASE1L3.

Based on these promising results, we conducted a comprehensive analysis of TA-TLS gene expression in PDAC, developing the TLSscore and associated prognostic model. This not only provided broad insights into PDAC prognosis but also indirectly confirmed the potential tumor-suppressive effects of TLS-related genes in PDAC. The formation of TLS and the induction of more TLS are emerging as promising strategies in immunotherapy. Identifying the key genes involved in TLS formation and understanding the underlying mechanisms are critical to advancing these therapeutic approaches [48, 49]. Our findings offer valuable experience for conducting multi-omics studies on TLS in PDAC. However, we must acknowledge the limitations of our approach. First, the accuracy of the scoring model based on TA-TLS genes requires further refinement and should be complemented by clinical guidelines for improvement. Additionally, the mechanisms by which key genes derived from TA-TLS regulate immune responses need to be validated through extensive experimental verification, particularly in clinical samples. It is crucial to identify methods that might induce TLS formation within PDAC tumor tissues, aiming to reshape the tumor microenvironment and enhance the response to immunotherapy. We are actively integrating relevant clinical samples to further elucidate the clinical value of TLS-related genes in predicting patient prognosis and their underlying mechanisms of action in tumors.

Conclusions

In conclusion, our study extends valuable insights into the expression patterns and roles of TLS related genes in pancreatic cancer. TA-TLS signature genes can serve as reference points for prognostic prediction in pancreatic cancer and offer potential entry points for personalized clinical decision-making, particularly in the context of future immunotherapy strategies for pancreatic cancer.

Fig. 1

A flowchart of this study, illustrating the research design and process

Fig. 2

(A) Consensus clustering of PDAC in the GEO-meta cohort based on TA-TLS genes, showing results for k (cluster number) set to 2, 3, and 4. (B) Principal component analysis (PCA) demonstrating distinct differences in transcriptomes among the three clusters. (C) Kaplan–Meier curves illustrating overall survival (OS) differences among the three clusters (P = 0.045). (D) Heatmap showing the expression profiles of TA-TLS genes across the three clusters. (E) Boxplot depicting ssGSEA scores of immune cells among the three clusters. (F) Boxplot displaying TME signatures across the three clusters in the GEO-meta cohort. (* p < 0.05; ** p < 0.01; *** p < 0.001; **** p < 0.0001)

Fig. 3

(A) Venn diagram showing the differentially expressed genes (DEGs) across the three clusters. (B) Forest plot illustrating 20 key prognostic genes among the DEGs. (C) GO pathway enrichment results for the prognostic genes. (D) Clustering based on TLSscore obtained from ssGSEA scoring, resulting in TLS_H and TLS_L groups. (E) Kaplan–Meier curves analyzing the prognosis between the TLS_H and TLS_L groups

Fig. 4

(A) Correlation analysis between TLSscore and TLScluster using the Kruskal-Wallis test. (B) Sankey diagram illustrating the relationship among TLScluster, TLSscore, and survival status. (C) Heatmap of tumor microenvironment characteristics and immune cell expression between TLS_H and TLS_L groups. (D) Violin plot showing differences in TME scores between TLS_H and TLS_L groups. (* p < 0.05; ** p < 0.01; *** p < 0.001)

Fig. 5

(A, B) AUC curves for the risk assessment model constructed using key prognostic genes. (C) Survival analyses for high and low Riskscore groups in the GEO-meta cohort. (D) Validation of survival analyses for the risk assessment model in the TCGA-PDAC cohort. (E) Validation of survival analyses for the risk assessment model in the GSE85916 cohort. (F) Boxplot of differential expression of immune checkpoint genes between high and low Riskscore groups. (G) Boxplot of drug sensitivity between high-risk and low-risk groups. (* p < 0.05; ** p < 0.01; *** p < 0.001)

Fig. 6

(A) Boxplot showing filtering criteria during scRNA-seq processing. (B) UMAP plot depicting the major cell types in PDAC, with each dot representing a single cell and cell types color-coded. (C) Proportions of different cell types across PDAC patient samples. (D) UMAP plot of AUC scores based on the gene set used to construct the TLSscore. (E) Boxplot of ssGSEA scores based on the gene set used to construct the TLSscore. (F-H) Density plots showing single-cell expression levels of the three key genes from the risk scoring model: DNASE1L3, RAMP2, and IL33

Fig. 7

(A-C) Correlation analysis of DNASE1L3, RAMP2, and IL33 with immune infiltration. (D-F) Scatter plots showing the correlation between DNASE1L3, RAMP2, and IL33 with their most associated immune cells. (G-I) Kaplan-Meier survival curves analyzing the expression levels of DNASE1L3, RAMP2, and IL33

Fig. 8

(A-C) Violin plots depicting the RNA expression levels of DNASE1L3, RAMP2, and IL33 in normal samples and PDAC tissues. (D) Immunohistochemistry (IHC) images showing the protein expression levels of DNASE1L3, RAMP2, and IL33 in normal pancreatic tissues and PDAC tumors

Fig. 9

(A) Western blot analysis showing DNASE1L3 expression levels across various cell lines. (B) Confirmation of DNASE1L3 overexpression in PANC-1 and Mia PaCa-2 cells via Western blot. (C) CCK-8 assay measuring the proliferation of PDAC cells with DNASE1L3 overexpression compared to NC-transfected PDAC cells. (D) Representative images and quantification of the colony formation assay for PANC-1 and Mia PaCa-2 cells with DNASE1L3 overexpression and siNC. Magnification: 200×. (E) Transwell migration assay showing the effect of DNASE1L3 on the migration ability of PDAC cells. Magnification: 200×. (F) Top three signaling pathways identified by GSEA analysis of DNASE1L3. (G) Heatmap of cytokine levels in the culture supernatants of PDAC cells with DNASE1L3 overexpression and NC. (* p < 0.05; ** p < 0.01; *** p < 0.001; **** p < 0.0001)