Mature tertiary lymphoid structures tumor microenvironment-based risk model to assess patients with pancreatic ductal adenocarcinoma

Introduction

Pancreatic ductal adenocarcinoma (PDAC), which comprises over 90% of all instances of pancreatic cancer, is not only the most common but also one of the deadliest types of the illness. It is characterized by an unfavorable prognosis and a significant death rate (1).

Tertiary lymphoid structures (TLS) are atypical aggregations of immune cells that develop in extra-lymphatic tissues. Typically, these formations develop in non-lymphoid tissues as a result of chronic inflammation, as in the case of malignancy. As TLS promotes the infiltration of immune cells into the tumor microenvironment (TME) of solid tumors, its presence may be correlated with patient survival. Prior studies have established that anomalous lymph nodes, which are TLS, contain tumor-infiltrating B cells. The function of these B cells is to facilitate tumor immune surveillance by delivering tumor antigens to T cells and producing anti-tumor antibodies (2). B cells significantly influence the action of TLS in PDAC. Long-term patient survival is associated with successful T cell infiltration and activation in TME (3). Hence, additional research on techniques to manipulate the infiltration and activation of immune cells in the TME of PDAC is imperative and could have a significant impact on the prognosis of PDAC patients. Therefore, exploring the role of immune cells in the TME plays an important role in immunotherapy. A more thorough characterization of TLS can investigate the condition of the cells and help to improve their value as prognostic and predictive indicators, and research into the functional state of TLS is critical for the development of therapeutic targets (4,5).

TLS play a central role in the induction and maintenance of tumor-specific immunity. To investigate the manifestation of the immune role of TLS in PDAC, the immune role was elucidated by screening for differential genes and mapping the immune landscape in the TME. Further, screening was conducted for drug targets of action and revealing drugs that may have potential therapeutic effects of drugs.

Among the 39 genes identified in TLS, CCL2, 3, 4, 5, 8, 18, 19, and 21, and CXCL9, 10, 11, and 13 were identified as chemokine marker genes in a recent study (6). Therefore, we established a scoring method for 11 chemokine markers and investigated the differential genes and molecular mechanisms in different TLS contents by establishing this scoring classification method at the transcriptomic level, studied the immune infiltration of PDAC in the TME, and investigated the potential therapeutic drug profile by taking intersections with drug targets in the DrugBank database and Therapy Target Database, which may be informative for current clinical trials and drug studies. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1601/rc).

Methods

Retrieval and processing of PDAC datasets

Publicly available datasets of PDAC include The Cancer Genome Atlas Pancreatic Adenocarcinoma (TCGA-PAAD), and Gene Expression Omnibus (GEO; GSE226840, GSE130563), which included RNA-Seq, clinical information, and treatment. TCGA sequencing data [fragments per kilobase million (FPKM) values] were downloaded from TCGA database (https://portal.gdc.cancer.gov), and GEO was downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/gds). The datasets were analyzed by the R software (R Foundation for Statistical Computing, Vienna, Austria) using the ID transformation package. The IMvigor210 immunotherapy cohort was acquired from the R software using the “IMvigor210CoreBiologies” R utility. GSE226840 was classified into four groups according to the expression of TLS biomarker CD23 in the GEO database: normal, CD23high, CD23low, and CD23neg, and the expression was analyzed using the GEO 2R tool. For prognostic analysis, 177 patients with PDAC were enrolled in the TCGA-PAAD; patients lacking clinical data were eliminated randomly. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Establishment of the PDAC-related TLS signature

From the TCGA database, standard data normalized by log2 scale and robust multi-array averaging (RMA) were extracted, the microenvironment cell populations (MCP) counting method was used to predict the number of lymphocytes, and we designated a set of 11 chemokines—CCL2, CCL3, CCL4, CCL5, CCL8, CCL18, CCL19, CCL21, CXCL10, CXCL11, and CXCL13—from the TCGA-PAAD data focused on transcriptome level-specific assessment of TLS in PDAC by classifying TLS expression into TLS-high and TLS-low. Such a stratification approach enables determination of the different gene expression profiles of TLS and the predictive role of TLS features in PDAC prognosis. This approach was established based on the following rationale: the aforementioned chemokines play a central role in mediating the homing and recruitment of lymphocytes (such as T cells, B cells, and dendritic cells) as well as the initiation and formation of TLS. In particular, CXCL13, CCL19, and CCL21 are widely recognized as key drivers of TLS formation and function. Therefore, the overall expression level of this chemokine cluster serves as an effective molecular surrogate for TLS presence and activity within the TME. This gene expression-based quantitative approach provides an objective, reproducible metric compared to histopathological assessment, which is particularly suited to analyzing large-scale transcriptomic datasets (e.g., TCGA). It has been successfully applied in multiple cancer studies to infer TLS status (7,8). After calculating the lower quartile, median, and upper quartile, samples exceeding the median were designated as TLS-high, whereas those falling below the median were designated as TLS-low.

The GSE226840 dataset in the GEO database was classified into four groups on account of CD23 expression in TLS as normal, TLS-CD23high, TLS-CD23low, and TLS-CD23neg. Mature tertiary lymphoid structures (mTLS) are well-structured lymphoid formations that consist of B lymphocytes and CD23+ follicular dendritic cells (FDCs) (9). Immunohistochemical methods were used to classify PDAC-associated TLS into three stages of maturation: (I) the initial stage, where no CD23 staining is seen, embodied as TLS-CD23neg; (II) the intermediate stage, where dispersed CD23+ cells assume a dendritic morphology, labeled as TLS-CD23low; and (III) the final stage of the CD23+ mfdc network, labeled as TLS-CD23high. CD23high. normal, uninflamed neighboring tissue was used as a control and labeled as normal. Based on the differences at the molecular level, it is possible to study the trends in gene expression in TLS during maturation.

Exploration of TLS differentially expressed genes (DEGs)

Using the edgeR package and the PDAC in TCGA data, differential genes, as well as up-regulated and down-regulated genes in the TLS-high and TLS-low cohorts, were identified, and the GSE226840 dataset from the GEO database was categorized into four groups based on the expression settings of the TLS marker CD23: normal, CD23high, CD23low, and CD23neg, respectively. Differential genes were analyzed using the GEO 2R tool and NetworkAnalyst website (https://www.networkanalyst.ca/), and the limma (volcano plot) algorithm produced a volcano plot (Padj less than 0.05) to illustrate the relationship between mathematical significance (−log10 P value) and change magnitude (log2-fold change). The genes that were down-regulated were indicated by the color blue (P<0.05, fold change ≤−2), and gray indicated non-significant genes. Venn diagrams were generated using the limma (vennDiagram) algorithm to investigate the overlapping relationships of significant genes between multiple groups, and analyzed for gene co-expression for the four groups of samples (inclusion criterion was Padj less than 0.05). Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment was performed to analyze differential gene enrichment in signaling pathways, and networks were mapped to show how these pathways are linked and to demonstrate pathway-related influences and disease conditions. Based on the KEGG dataset, it is possible to visualize the gene set network in terms of the connections between signaling pathways, performing over representation analysis (ORA). The heat maps and volcano maps demonstrated the expression of differential genes. The red dots represent genes that were considerably elevated. These genes exhibited a positive logFC value and possessed a false discovery rate (FDR) below a specific threshold at 0.05. Green dots represent genes that had been down-regulated to a significant extent. The expression levels of these genes exhibited negative log-fold change values. Gene expression distribution was also depicted for each of the four samples. Gene Ontology (GO) and KEGG enrichment analyses were separately conducted for the DEGs using the R programming language. The graph uses color to represent the q-value, which denotes the adjusted significance level. The color red signifies a lower q-value, indicating higher relevance. Conversely, the color blue signifies a higher q-value, indicating lower significance.

GO, KEGG, gene set enrichment analysis, weighted gene co-expression network analysis (WGCNA), and unsupervised cluster analysis

The clusterProfiler package was utilized to apply GO enrichment (P<0.05) and KEGG analysis (P<0.05, FDR <0.05) to the PDAC sample information in the TCGA database to respectively identify DEGs. The outcomes were displayed with the ggplot2 and circlize packages. Gene set enrichment analysis (GSEA) was also performed using the clusterProfiler package and was performed separately for up- and down-regulated genes. The WGCNA package was employed to create gene co-expression networks for the cohort with the aim of detecting co-expressed gene modules, exploring the association between gene networks and phenotypes, and identifying the pivotal genes in the networks by constructing weighted gene co-expression networks. Furthermore, we condensed the co-expression patterns of the characterized genes and assessed the association between each characterized gene and immune cells. The WGCNA-R package was utilized to generate co-expression networks for every gene in the dataset. The soft threshold was established at 6, the mean value of expression matrix filtering was maintained at 0.5, the P value for immune-infiltration data was deemed less than 0.05, and the weighted neighbor-joining matrix was transformed to assess network connectivity. A clustering tree was constructed using the hierarchical clustering method. The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database was utilized to analyze the DEGs module (www.string-db.org/). A protein-protein interaction (PPI) network was developed using Cytoscape version 3.10.1 software (https://cytoscape.org/), which enabled the visualization of complex networks and the integration of attribute data of any nature. In accordance with their expression patterns, genes were categorized based on their weighted correlation coefficients; genes exhibiting comparable patterns were consolidated into a single module. The ConsensusClusterPlus package was employed for unsupervised clustering within the TLS-high and TLS-low cohorts, enabling further categorization of these genes.

Screening for risk genes and constructing TME risk models for scoring

Risk genes were screened and scored in the TCGA expression matrix using the limma R package to draw heat maps. We also compared the relationship between risk scores of TLS-high and TLS-low groups and visualized this relationship by plotting a histogram. At the same time, using the maftools package, various data (histological mutations, copy number variations, clinical data, expression data, etc.) were summarized and stored as MAF objects, and based on the MAF files, risk genes were plotted in a maftools waterfall plot (Oncoplot), which demonstrated the mutation information in the MAF files. Oncoplot waterfall plots are used in TME studies to visualize the distribution of gene mutations. Oncoplot waterfall plots are created on the basis of mutation data provided by the TCGA database in the MAF format employing the maftools package for the examination of the mutation status of the genes associated with risk. The waterfall plot allows us to quickly identify the genes with the highest number of mutations in all samples (right bar) and quickly determine which mutations account for the most significant percentage of a sample (top bar).

Cox-PH and prognostic analysis

Cox survival regression models were employed to evaluate the significance of each parameter as a prognostic factor for PDAC through univariate and multivariate analyses. Separate analyses were conducted to assess the prognostic impact of clinical characteristics, risk genes, and immune-infiltrating cells. The effects of these factors as covariates on prognosis were individually examined. DEGs were subjected to univariate Cox analysis using the survival package to obtain the hazard ratio (HR) and z-value (Wald statistic). The z-value was determined by dividing the regression coefficient (coef) by its standard error (se(coef)), expressed as z = coef/se(coef), and the P value was calculated. Calculations were made for HR regarding overall survival (OS), and independent prognostic indicators were subsequently integrated into a multivariate Cox proportional hazard model to predict PDAC patient survival. The influence of prior malignancies on OS was assessed through HR and their respective 95% confidence intervals (CIs).

The significance of covariate inclusion was determined by P values. Kaplan-Meier curves were generated to illustrate the impacts of immune-infiltrating cells, risk genes, and clinical characteristics on patients’ prognosis. Additionally, the TLS-high group and TLS-low group were compared in high- and low-risk cohorts, evaluating their associated prognostic outcomes. Furthermore, the roles of TLS-high and TLS-low groups in relation to prognosis were compared.

Investigation of immune infiltration and TME characteristics

We analyzed immune infiltration in PDAC samples retrieved from the TCGA database. Using transcriptome data, to estimate the composition and amount of immune cells in mixed cell populations, we employed a deconvolution technique available on the Cell-type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) website (http://cibersortx.stanford.edu/). CIBERSORT is an online program used to separate expression matrices of different types of human leukocytes, providing a set of gene expression signatures for 22 leukocyte subtypes. Immune infiltration histograms, heatmaps, and bar charts were plotted separately to show the characteristics of immune infiltration. Scatter plots were employed to demonstrate the relationship between immune cell distribution and TME risk scores, box plots to demonstrate the distributions of immune cells in the risk variance groups in the TME risk model situation, and a lollipop chart showed the linear relationship between immune cells and risk model scores. The Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) algorithm was used to compute the ImmuneScore, StromalScore, and ESTIMATE score, normalized using the limma package, and computes the TumorPurity, comparing the difference between the different TLS expression groups.

Anti-programmed cell death 1 (PD-1) inhibitor immunotherapy response prediction in TLS

We further explored an independent immunotherapy cohort (IMvigor210) based on 348 anti-PD-1 inhibitor pretreatment samples to predict the response in variant expressed TLS level. Clinical information collected included best confirmed overall response, binary response, and immune phenotype. Immune phenotype included excluded, inflamed, and desert. We plotted the survival curve to predict the survival after immune therapy in the variant-expressed TLS group as validation.

Statistical analysis

The figures in this study were plotted and data processing was conducted using R Software version 4.3.1 (https://www.rstudio.com) and GraphPad Prism 9.5.0 (GraphPad Software, San Diego, CA, USA). The statistical analysis included the Mantel-Cox test, Gehan-Breslow-Wilcoxon test, and unpaired two-sample t-test. Statistical significance was denoted by both an adjusted q value and a P value below 0.05. For immune score, stromal score, and ESTIMATE score, we used two-sample t-test to compare the means of the two sets of data for any significant differences. Analysis of variance (ANOVA) summary and Brown-Forsythe test were used for unsupervised clustering of C1, C2, and C3 respectively. The violin map and the sandbox map of immune cells were subjected to Mann-Whitney U test to demonstrate the differences between cells. The survival curves used log-rank (Mantel-Cox test) test and the Gehan-Breslow-Wilcoxon test to compare the differences in survival.

Results

Establishment of TLS signature and molecular expression differences in PDAC

We constructed a new gene signature consisting of 11 chemokine genes (Figure 1A). These chemokines were associated with different proportions of prognostic risk (Figure 1B). The expression of differential genes was analyzed in each of the four groups, and the top 10 groups with the highest differential gene significance were PLA2G1B, IL7R, CCL19, CXCR4, CD24, AMBP, LTBR, SPP1, MS4A1, and APP. The PDAC samples were categorized into two categories based on their TLS levels: TLS-high and TLS-low, and a heat map (Figure 1C) and volcano maps (Figure 1D) were plotted. The volcano map demonstrated that there were more up-regulated than there were down-regulated genes. We explored the differences in PLA2G1B on immunohistochemical images in normal tissue, TLS-high, and TLS-low groups. Figure 1E refers to PLA2G1B in normal tissue; Figure 1F refers to PLA2G1B in tumor tissue. Immunohistochemical images showed that PLA2G1B expression was higher in normal tissues than in tumor tissues. DEGs were more significant in the normal (normal) and TLS-CD23high groups (Figure 2A) relative to the other expression groups (CD23neg vs. normal, CD23high vs. CD23low, CD23low vs. CD23neg) (Figure 2B-2D), suggesting that the final stage of the CD23+ mfdc network is more significant and high CD23 expression in TLS was associated with enhanced immune responses and may play an important role in the TME which supporting the importance of TLS as a regulatory center for anti-tumor immune responses. Upregulated and downregulated genes were more significant, among which the up-regulated genes included AMBP, APP, CD24, CD63, CDH1, CXCL12, MUC1, and IL22RA1, and the down-regulated genes included APOE, BTLA, CCL19, CD22, CD48, CD53, CD63, CD79A, CXCR4, and IL7R. The Wayne diagram (Figure 2E) showed that CD23neg versus normal and normal versus TLS-CD23high genes intersected the most, with 39 genes. Pathways showing high significance in KEGG enrichment results included cytokine-cytokine receptor interaction, NF-kappa B signaling pathway, and intestinal immune network for immunoglobulin A (IgA) production (FDR <0.01) (Figure 2F). Bar and bubble plots were obtained by GO enrichment (Figure 2G,2H), and the first three items of the enrichment results were positive regulation of lymphocyte activation, positive regulation of leukocyte activation leukocyte cell-cell adhesion, and positive regulation of cell activation leukocyte cell-cell adhesion. The GO enrichment indicated that these genes likely had significant functions in the TME and immune response, particularly in areas related to immunoglobulins and immune cell receptors. For the KEGG signaling pathways ranking based on differential expression (DE) method used (P value), the network centered on the intestinal immune network for IgA production was more significant, as can be seen from the visualization plot (Figure 2I). The pathways connected to it included cytokine-cytokine receptor interaction, Th1 and Th2 cell differentiation, and inflammatory bowel disease (IBD). The metabolic pathways included purine metabolism and pyrimidine metabolism. The p53 signaling pathway existed independently of other signaling pathways.

Figure 1 TLS signature establishment and function enrichment analysis. (A) The heatmap of the establishment of TLS signature. (B) Forest plot for different signature chemokines. (C) The heatmap of gene expression in TCGA. (D) Volcano plot in TCGA cohort. In the figure, the color scale indicates the direction of change: red represents up-regulation (or increased expression) and green represents down-regulation (or decreased expression). This is consistent with common conventions in the field for visualizing such differential analysis. (E) PLA2G1B protein expression in normal group, acquired from the Human Protein Atlas database (https://www.proteinatlas.org/ENSG00000170890-PLA2G1B/cancer/pancreatic+cancer#img). The right panel of this image is presented at a 4× magnification, showing the spatial expression and localization characteristics of PLA2G1B in normal pancreatic parenchymal and stromal tissues. (F) PLA2G1B protein expression in pancreatic cancer group, acquired from the Human Protein Atlas database (https://www.proteinatlas.org/ENSG00000170890-PLA2G1B/tissue/pancreas). The right panel of this image is presented at a 4× magnification, displaying the expression pattern and distribution of PLA2G1B in pancreatic cancer tissue with malignant morphological characteristics. The tissue microarray sections were stained by immunohistochemistry using the antibodies HPA047822 and HPA060803. The images were acquired at 4× magnification. FC, fold change; FDR, false discovery rate; HR, hazard ratio; TCGA, The Cancer Genome Atlas; TLS, tertiary lymphoid structures.

Figure 2 Exploration of TLS DEGs and functional enrichment in the GEO database. (A) Volcano plot between normal and CD23-high. (B) Volcano plot between CD23-neg and normal. (C) Volcano plot between CD23-high and CD23-low. (D) Volcano plot between CD23-neg and CD23-low. Red represents up-regulation and green represents down-regulation. This is consistent with common conventions in the field for visualizing such differential analysis. (E) Venn diagram in GSE226840 cohort. Venn diagram showing the overlap of DEGs among the indicated pairwise comparisons (normal, CD23-high, CD23-low, and CD23-neg groups) in dataset GSE226840. Numbers in the diagram indicate the number of DEGs in each region, identified using limma with Benjamini-Hochberg adjusted P<0.05 (both up- and down-regulated genes included). (F) KEGG function analysis. (G) GO enrichment barplot. (H) GO enrichment bubble diagram. (I) Enrichment pathway network diagram. (J) GSEA enrichment of up-regulated genes. (K) GSEA enrichment of down-regulated genes. (L) GSEA enrichment in TCGA. BP, biological process; CC, cellular component; DEGs, differentially expressed genes; GEO, Gene Expression Omnibus; GO, Gene Ontology; GSEA, gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; TCGA, The Cancer Genome Atlas; TLS, tertiary lymphoid structures.

The GSEA plots (Figure 2J) showed differential gene enrichment in the motor proteins (enrichmentScore189), cell cycle (enrichmentScore157), and carbon metabolism (enrichmentScore115). The upward pathways included citrate cycle [tricarboxylic acid (TCA) cycle], fructose and mannose metabolism, mucin type O-glycan biosynthesis, nitrogen metabolism, and pentose and glucuronate interconversions (Figure 2K). The downward pathways included allograft rejection, asthma, autoimmune thyroid disease, graft-versus-host disease, and primary immunodeficiency (Figure 2L). The GSEA enrichment results indicated that in particular, pathways related to immune response and metabolism were significantly enriched in the sequenced dataset. These significantly enriched sets of genes suggested that related biological processes and metabolic pathways may play an important role under specific experimental conditions or disease states. WGCNA employs hierarchical clustering analysis (HCA) to assign genes to various modules and represents each module with distinct colors. Figure 3A showed the results of soft thresholding power selection in the construction of gene co-expression network to ensure the network has scale-free topology, which represented the scale independence and mean connectivity under different soft thresholds. The results showed that between soft thresholds of 6 and 10, the R2 value exceeded 0.8, indicating that the network conformed better to the scale-free properties at these thresholds. The correlation between the brown plate co-expression and seven lymphocyte species was calculated in 40 modules (Figure 3B). The module-feature relationship plots showed that the brown plate co-expression was positively associated with T cells CD8 (R=0.49, P=1.9e−82), T cells CD4 naive (R=0.47, P=4e−75), T cells CD4 memory activated (R=0.31, P=1.9e−31), T cells follicular helper (R=0.25, P=1.1e−20), and T cells CD4 memory resting (R=0.36, P=1.4e−42). Meanwhile, there was a negative correlation with T cells regulatory (Tregs) (R=−0.22, P=2.9e−16) and T cells gamma delta (R=−0.26, P=2.7e−22). The module–feature relationship plot and the differential genes in the TLS-high group and TLS-low group showed that the co-expressed brown module was a key module of PDAC related to the TLS and was used for further analysis. The heatmap (Figure 3C) illustrates the relative abundance or expression levels of various T cell types in numerous samples. A significant presence of CD8+ T cells may suggest an anti-immune response, whereas a large proportion of Tregs may indicate immunosuppression in the sample. The presence of different colors in specific regions indicated the presence of gene groups that share comparable patterns of gene expression (Figure 3D-3F).

Figure 3 WGCNA related immune filtration and gene expression. (A) Scale independence. (B) Module-trait relationships. (C) Sample dendrogram and trait heatmap. (D) Dynamic tree cut. (E) Merged dynamic. (F) Clustering of module eigengenes. (G) PPI network constructed by TLS differentially expressed genes in the STRING database. (H) Cluster selected by the PPI network. (I) PPI network constructed by Cytoscape. PPI, protein-protein interaction; TLS, tertiary lymphoid structures; WGCNA, weighted gene co-expression network analysis.

A network was constructed for the 1065 genes in the co-expressed brown panel using the STRING database and the Molecular Complex Detection (MCODE) plugin in Cytoscape (Figure 3G,3H). The minimum required interaction score was set in the STRING database: highest confidence 0.900, average node degree: 21.8, PPI enrichment P value: <1.0e−16; the results showed that CD8A was the gene with the highest degree in the PPI network, with a degree of 96 (Figure 3I). GO enrichment was performed in the enrichment function that comes with the STRING database and KEGG enrichment in the STRING database, the local network cluster showed that the most significant clustering was MHC II, with a strength of 1.27 and an FDR of 1.36e−07. The consensus matrix (Figure 4A,4B) plotted the clustering results at k=3, showing the consistency of the clustering of the samples across subsamples, which could be divided the PDAC cohort into three clusters (C1, C2, and C3), and the three clusters were categorized into a high immune checkpoint inhibitor (ICI) score group, C3, a medium ICI score group, C1, and a low ICI score group, C2, by means of the immunity score (Immunity Score). In unsupervised clustering, a consensus cumulative distribution function plot (Figure 4C) showed that the optimal k value was 3, delta area plot (Figure 4D) and the tracking plot (Figure 4E). Distinct immunization scores indicated various stages of TLS. In addition, a clinicopathological information heat map was created based on the TCGA sequence, visually linking tumor staging stage, T-classification, N-classification, M-classification, gender, TLS, and risk (Figure 4F).

Figure 4 Cluster analysis of differential genes, impact of immune cell and PDAC characteristics on the riskiness of patient clinical outcomes. (A) Cluster heatmap showing the consensus matrix k=3. (B) Consensus matrix legend. (C) Consensus CDF. (D) Relative change in area under CDF curve. (E) Samples tracking plot. (F) Heatmap of clinical characteristics related with patients outcome. *, P<0.05; **, P<0.01. CDF, cumulative distribution function; M, metastasis; N, node; PDAC, pancreatic ductal adenocarcinoma; T, tumor.

Modeling TLS risk scores in relation to TME

The results obtained from univariate Cox analysis were further subjected to multivariate Cox analysis of the effects of prognosis-related genes on clinical outcomes of PDAC patients, and finally seven genes, including CXCL11, CASC8, REEP2, TNNT1, SLC16A11, DUSP26, and CHGA, were identified, which were further applied to construct the prognostic risk model. CXCL11, CASC8, TNNT1, and DUSP26 were identified as high-risk, and REEP2, SLC16A11, and CHGA were identified as low-risk. The formula for calculating the risk score was as follows:

riskScore = 0.35418 × CXCL11 + 0.48872 × CASC8 + 0.18693 × TNNT1 + 0.72090 × DUSP26 − 0.48344 × REEP2 − 0.67860 × SLC16A11 − 0.25736 × CHGA.

The risk scores obtained from the Cox regression analysis of all patients were arranged in order from smallest to largest to make a substate plot (Figure 5A) to assess the prognosis of the patients. On the plot, the X-axis denotes patient serial numbers, whereas the Y-axis represents survival time. Green markers indicate survival and red markers indicate death. Examination of the plot revealed that the high-risk group demonstrated a reduced duration of survival (as evidenced by a descending pattern of scattered dots) and a greater incidence of mortality (with a larger number of red dots compared to green dots). The gene expression heatmap (Figure 5B) showed the relationship between risk score and risk genes. The bar graph (Figure 5C) showed the distribution of risk scores, and it can be seen that the TLS-low group had a higher risk score than the TLS-high group. The principal component analysis (PCA) visualization chart (Figure 5D) depicted the level of risk associated with specific genes. The forest plot related to risk genes is shown in Figure 5E. In the oncoplot (Figure 5F-5H), the horizontal coordinate represented the sample, the vertical coordinate the gene, and there was a color to show that the gene had mutated in the sample, the closer to black the higher the mutation frequency. The top bar graph represents the overall mutation status of these genes, and the different colors reflected the specific mutation status of the sample. The waterfall plots demonstrated mutations in 431 samples. The majority of mutations were missense mutations (green), followed by shift deletions (blue), and nonsense mutations (red). WNK4 had the highest mutation rate of 5%, followed by PTPRN2, KCNJ3, and ABCC8, with a mutation rate of 4%. Figure 5I shows the pathological image between CXCL11, REEP2, DUSP26, and CHGA.

Figure 5 Establishment of risk model constructed by DEGs, mutation state of risk genes, and expression of risk genes. (A) The curve diagram and scatter plot of patient risk scores and survival time. Each point represents one patient. Red indicates high‑risk patients, and green indicates low‑risk patients. (B) The heatmap of risk score correspond with samples in TCGA cohort. (C) The risk score between TLS-high and TLS-low group. (D) PCA analysis of patient risk. (E) The forest plot of risk genes related with clinical outcomes. (F) The mutation landscape of alteration in 108 of 431 samples (25.06%). (G) The mutation landscape of alteration in 101 of 431 samples (23.43%). (H) The mutation landscape of alteration in 99 of 431 samples (22.97%). (I) The expression of top risk genes in normal and pancreatic cancer, all retrieved from the Human Protein Atlas (HPA) database (CXCL11: https://www.proteinatlas.org/ENSG00000169248-CXCL11/tissue/pancreas; REEP2: https://www.proteinatlas.org/ENSG00000132563-REEP2/tissue/pancreas; DUSP26: https://www.proteinatlas.org/ENSG00000133878-DUSP26/tissue/pancreas; CHGA: https://www.proteinatlas.org/ENSG00000100604-CHGA/tissue/pancreas). All images in this panel are captured at a 4× magnification, showing the distinct spatial expression and localization patterns of the four molecules in normal pancreatic parenchymal/stromal tissues and malignant pancreatic tumor tissues for direct comparative analysis. CXCL11 was detected using antibody CAB009921. REEP2 was detected using antibody HPA031813. DUSP26 was detected using antibody HPA018221. CHGA was detected using antibody HPA017369. DEGs, differentially expressed genes; HR, hazard ratio; NA, not applicable; PCA, principal component analysis; TCGA, The Cancer Genome Atlas; TLS, tertiary lymphoid structures.

Cox-PH and prognostic analysis

Univariate and multivariate Cox regression analyses of clinically relevant factors were performed, and the results showed that tumor stage (HR =1.102, 95% CI: 0.748–1.624) and T-classification (HR =1.359, 95% CI: 0.874–2.111) had non-significant effects on risk (P>0.05) (Figure 6A,6B) and N-classification (HR =1.801, 95% CI: 1.048–3.097) was significant risk factor (P<0.05). In multivariate Cox regression, TLS (HR =0.633, 95% CI: 0.412–0.947), with HR <1, indicated that TLS was a protective factor favoring the prognosis of PDAC patients, whereas the effect of gender on risk was not significant (HR =0.755, 95% CI: 0.494–1.154) (P>0.05). Univariate and multivariate Cox regression analyses were performed on the immune cells analyzed in CIBERSORT to investigate the effect of immune cells on prognostic outcomes. In univariate Cox analysis (Figure 6C), B cells naive, B cells memory, plasma cells, T cells CD8, monocytes, and macrophages M0 were significantly associated with prognosis (P<0.05), with B cells naive (HR =0.6, 95% CI: 0.36–0.99), plasma cells (HR =0.59, 95% CI: 0.36–0.99), T cells CD8 (HR =0.55, 95% CI: 0.32–0.94), monocytes (HR =0.49, 95% CI: 0.30–0.81) shown to be protective factors, with HR less than 1, favoring patients’ prognostic recovery. These types of immune cells were positively correlated with the function of TLS; macrophages M0 (HR =1.93, 95% CI: 1.18–3.17), B cells memory (HR =2.24, 95% CI: 1.24–4.04) with HR >1 shown to be risk factors, increasing the risk of occurrence and decreasing the survival rate. In multivariate Cox analysis (Figure 6D), monocytes, macrophages M2, and dendritic cells resting were significantly associated with prognosis (P<0.05); monocytes (HR =0.49, 95% CI: 0.26–0.90) were shown to be a protective factor, and dendritic cells resting (HR =2.02, 95% CI: 1.06–3.86) were a risk factor and unfavorable to patient prognosis.

Figure 6 Cox regression analyses of clinical characteristics and immune cells in PDAC prognosis. (A) Uni-Cox analysis of PDAC clinical characteristics in forest plot. (B) Multi-Cox analysis of PDAC clinical characteristics. (C) Uni-Cox analysis of immune cells related with patients clinical outcomes in forest plot. (D) Multi-Cox analysis of immune cells. HR, hazard ratio; N, node; NK, natural killer; PDAC, pancreatic ductal adenocarcinoma; T, tumor; TLS, tertiary lymphoid structures.

The Kaplan-Meier survival curves demonstrated that the low-risk group exhibited a more favorable prognosis compared to the high-risk group in the risk model (Figure 7A). The prognostic model’s sensitivity and specificity were assessed by computing the area under the receiver operating characteristic (ROC) curve (AUC), utilizing ROC curves across time. The sensitivity and specificity of the prognostic model were evaluated by calculating the AUC, which was 0.757 for the 1-year curve, 0.792 for the 3-year curve, and 0.839 for the 5-year curve, which was seen to be the largest AUC at 5 years, and the predictive accuracy of the model was better (Figure 7B,7C). Kaplan-Meier survival curves of risk genes were plotted to elucidate the effects of risk genes on the survival status of high-risk and low-risk groups. The high expression of BEX1 (Figure 7D), MIR7-3HG (Figure 7E), PEMT (Figure 7F), and SMIM6 (Figure 7G) showed better survival than low expression. CXCL10 had a better prognosis for the low-risk group, and the high-risk group had the lowest survival rate at 3 years (Figure 7H). To study the TLS expression and patient prognosis, the Kaplan-Meier survival curve showed better survival in TLS-high group compared to TLS-low group (Figure 7I), log-rank test showed a P value of 0.01, Chi-squared of 6.599, and a significant difference in survival curve. In the male group (Figure 7J) log-rank test, Chi-squared was 5.454, and the P value was 0.02, which was significant difference, and the Kaplan-Meier survival curve demonstrated superior survival rates in the TLS-high group as compared to the TLS-low group. In the female group, using Wilcoxon test, Chi-squared was 3.958, and the P value was 0.047, which was significant, and the TLS-high group had a better survival rate than the TLS-low group (Figure 7K). In the Stage I/II group (Figure 7L), the Kaplan-Meier survival curves showed that the TLS-high group had a better survival rate than the TLS-low group; using the Wilcoxon test, the Chi-squared was 4.816, P value was 0.03, and the survival curves were significantly different.

Figure 7 Kaplan-Meier survival curve and prognostic models predict patient survival. (A) Kaplan-Meier survival curve between high- and low-risk groups. (B) ROC curve calculating the true positive rate and optimal threshold. (C) ROC curve calculating the 1-, 3-, and 5-year AUC. (D) Kaplan-Meier survival curve of BEX1 related with patient prognosis. (E) Kaplan-Meier survival curve of miR7-3HG related with patient prognosis. (F) Kaplan-Meier survival curve of PEMT. (G) Kaplan-Meier survival curve of SMIM6. (H) Kaplan-Meier survival curve of CXCL10. (I) Kaplan-Meier survival curve between the TLS-high and TLS-low groups. (J) Kaplan-Meier survival curve between the male TLS-high and TLS-low groups. (K) Kaplan-Meier survival curve between the female TLS-high and TLS-low groups. (L) Kaplan-Meier survival curve between the TLS-high and TLS-low groups in stage I/II. AUC, area under the curve; ROC, receiver operating characteristic; TLS, tertiary lymphoid structures.

The landscape of immune infiltration created by TLS

We performed immune infiltration analysis on PDAC samples. The heatmap showed the distribution of TLS score and the immune cells type (Figure 8A). Originally, box plots were designed to illustrate the distribution of various immune cells in each sample. Notably, the TLS-high group generally exhibited a higher percentage of immune cells compared to the TLS-low group. This trend was observed in the expression levels of CD8 T cells, activated memory CD4 T cells, activated NK cells, and naive B cells within both groups. In the sandbox figure (Figure 8B), B cells naive, B cell memory, T cells CD8, T cells CD4 memory resting, T cells CD4 follicular helper, macrophages M0, macrophages M1, macrophages M2, and neutrophils, these immune cells were significantly different between the TLS-high and TLS-low groups. At the same time, the expression of B cells naive, B cell memory, CD8 T cells, macrophages M1, dendritic cells resting have higher expression in TLS-high than TLS-low group. The macrophages M0, macrophages M2, and T cells CD4 memory resting were more highly expressed in the TLS-low group than they were in the TLS-high group. In the violin figure (Figure 8C), there were significant differences in B cells naive, B cell memory, T cells CD8, T cells CD4 memory activated, T cells CD4 follicular helper, macrophages M1, dendritic cells resting, and neutrophils. These cells were significantly higher in the TLS-high group. However, T cells CD4 memory resting, natural killer (NK) cells resting, macrophages M0, and mast cells activated were higher in the TLS-low group. The correlation of immune cells was shown by the correlation heatmap (Figure 8D), which showed the strongest effect of naive B cells and plasma cells, with the antagonistic effect of CD8 T cells and M0 macrophages. In the immune cell correlation analysis, the association with the risk score can be analyzed. In correlation plots between immune cell expression levels and risk score (Figure 8E-8M), the expression of macrophages M1, T cells CD4 memory resting, and T cells CD4 follicular helper had a positive correlation with risk score. Lollipop plots (Figure 8N) can show correlation coefficients for different immune cells. Significantly related cell types were macrophage M1 (cor =0.31, P=0.002) and B cells naive (cor =0.24, P=0.01).

Figure 8 Immune infiltration analysis of 22 immune cells. (A) Heatmap of immune cells in TCGA cohort. (B) The sandbox figure of immune infiltration between the TLS-high and TLS-low groups. (C) The violin figure of immune infiltration between the TLS-high and TLS-low groups. (D) The correlation heatmap of immune cells. (E) The correlation between B cell naive and risk score. (F) The correlation between monocytes and risk score. (G) The correlation between macrophages M1 and risk score. (H) The correlation between NK cells activated and risk score. (I) The correlation between T cells CD4 memory resting and risk score. (J) The correlation between T cells follicular helper. (K) The correlation between T cells CD8. (L) The correlation between B cell naive and risk score expression. (M) The correlation between macrophages M1 and risk score. (N) The lollipop figure shows the correlation coefficient between immune cells. (O) The bar graph shows the correlation of the risk score and immune cells. ns, not significant; *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001. NK, natural killer; TCGA, The Cancer Genome Atlas; TLS, tertiary lymphoid structures.

The bar graph (Figure 8O) showed that naive B cells, T cells CD8, T cells regulatory (Tregs), monocytes, and NK cells activated had a negative correlation with risk score. The CIBERSORT algorithm assessed immune cell presence in the tumor, revealing distinctions between risk variance groups. Additionally, the ESTIMATE tool determined tumor purity, ESTIMATE score, ImmuneScore, and StromalScore for each PDAC sample (Figure 9A). A t-test compared TLS-high and TLS-low groups for these scores, with P values over 0.05 indicating no significant difference. The violin plot (Figure 9B-9D) showed that the TLS-high group had higher ESTIMATE, immune, and stromal scores compared to the TLS-low group. However, the TLS-low group exhibited reduced tumor purity (Figure 9E). At the same time, we compared the immune score in C1, C2, and C3 clusters (Figure 9F). The immune score was the highest between C1 and C2; C1 was higher than C2.

Figure 9 Tumor microenvironment score, immune prediction, and validation of TLS function in IMvigor210 cohort. (A) Stromal score, immune score, and ESTIMATE score between high risk and low risk group. (B) Immune score between the TLS-high and TLS-low groups. (C) ESTIMATE score between the TLS-high and TLS-low groups. (D) Stromal score between the TLS-high and TLS-low groups. (E) Tumor purity between the TLS-high and TLS-low groups. (F) Immune score between the C1, C2, and C3 groups after consensus cluster. (G) Kaplan-Meier survival curve between the TLS-high and TLS-low groups in IMvigor210 cohort. (H) Best confirmed overall response between the TLS-high and TLS-low groups. (I) Binary response after anti PD-1 treatment in the TLS-high and TLS-low groups. (J) Venn diagram of drug and disease target. (K) The cluster selected by PPI network from Cytoscape. (L) GO enrichment of possible treatment target. Data are presented as mean ± standard deviation (or SEM). ****, P<0.0001. CR, complete response; GO, Gene Ontology; PD, progressive disease; PD-1, programmed cell death 1; PPI, protein-protein interaction; PR, partial response; SD, stable disease; SEM, standard error of the mean; TLS, tertiary lymphoid structures; TME, tumor microenvironment.

Anti-PD-1 inhibitor immunotherapy response prediction

We utilized the IMvigor210 cohort to predict immunotherapy response in the variant TLS expression group. The Kaplan-Meier analysis revealed a higher survival rate among patients with high TLS levels following immunotherapy (Figure 9G). Individuals with an elevated TLS score exhibited a heightened likelihood of responding favorably to immunotherapy. The comprehensive results suggest that TLS could serve as a valuable biomarker for identifying PDAC patients suitable for immunotherapy. Figure 9H,9I show the proportion of binary response in the TLS-high and TLS-low groups.

Network pharmacology analysis of PDAC samples as potential therapeutic targets

By taking the intersection of the TLS-high and TLS-low group differential genes of PDAC in GSE226840 and TCGA, 85 DEGs were obtained; 63 drug targets were screened in the DrugBank database and Therapeutic Target Database, and 10 targets were seen by taking the intersection using the Wayne diagram (Figure 9J), which were CD96, CD37, CD2, CD52, PTPRC, ATM, CD27, CTLA-4, ITGAL, and IL2RA. Using the STRING database to obtain protein interactions, the MCODE plugin in Cytoscape analyzed disease and drug targets to establish a disease-drug target PPI network (Figure 9K), which could be observed as three clusters (red for disease targets, green for drug targets). For Cluster 1, with a score of 20.5, the highest degree target was SELL 114, followed by CD2, degree 110. For Cluster 2, with a score of 13.565, the highest degree target was CD4, degree 180, followed by PTPRC; in Cluster 3, CD4 ranked highest with a degree of 180, followed by PTPRC. A bar chart (Figure 9L) was drawn to show this difference by sorting the degree of the targets, and the results showed that CD4 was the highest ranked target. The specific drug and target information is provided in Tables S1,S2.

Discussion

PDAC is characterized by a poor prognosis and high mortality rate, with diagnosed patients having an aggregate 5-year survival rate of less than 8% (10). Staging by the level of TLS expression can be linked to patient prognosis, changes in the TME can be interpreted in terms of molecular mechanisms, and risk models can be developed to understand the impact of molecular mechanisms on patient clinical outcomes. A better understanding of tumor immune mechanisms is pivotal to improve clinical outcomes.

In this study, we analyzed the gene expression profile of the different developmental stages of TLS by examining transcriptomic data from GEO data. This study aimed to extract refined prognostic features from the complex biology of TLS through an integrated analytical process. We first used 11 widely validated chemokine signatures as molecular proxies to classify patients into TLS-high/low phenotypes, thus defining the core biological background of the study (4). Subsequently, in order to go beyond simple phenotype associations and explore the core regulatory network of the immune microenvironment that drives this phenotype, we adopted WGCNA and used immune cell infiltration levels rather than TLS grouping itself as trait inputs. This critical design focuses our analysis on the internal immune ecosystem of tumors, and the identified brown module, due to its high correlation with key immune cells such as cytotoxic T cells and B cells, as well as significant overlap with the TLS DEG set, has been identified as the TLS-related immune core co expression network. From this immune network rich in TLS biological information, we further extracted seven gene prognostic features through Cox regression. This step aims to maximize the predictive performance and clinical practicality of the model, so the final selected gene combinations represent the “signal integrator” that can best capture prognostic information, but do not represent the core components of the TLS structure. This cleverly explains why the final model gene has a low degree of direct overlap with the initial chemokine marker. Our model focuses on predicting the output, which is a broad effector immune state that is coordinated by multiple immune cells.

Although the seven gene signature itself does not directly point to the B cell pathway, there is strong evidence that its prognostic efficacy is rooted in the core biological functions of TLS, particularly B cell-mediated immunity. Firstly, as the source of this signature, the brown module’s gene set showed a highly significant enrichment in core processes such as the B cell receptor (BCR) signaling pathway (P=1.09e−24) and B cell activation (P=4.07e−17) in detailed pathway analysis. Secondly, more directly, the survival analysis forest plot based on immune cell infiltration clearly shows that naive B cells are a significant protective factor (P=0.01). In summary, our seven gene risk score, as a highly integrated indicator, relies on the immune environment shaped by a functionally robust TLS, and the presence and function of B cells as the core pillar of this environment are the underlying biological basis for the model to predict good prognosis.

In terms of the signaling pathway, apelin treatment decreases NF-κB activation in the pancreas induced by acute pancreatitis (AP) and chronic pancreatitis (CP), thereby reducing inflammatory and fibrosis responses via inhibition of neutrophil recruitment and activity of pancreatic stellate cells. It was observed that PDAC samples characterized by elevated NF-κb scores displayed significantly more pronounced T cell exhaustion signatures (11). In this way, it functions to initially to decrease the expression of the NF-κB signaling pathway. Between the GO enrichment and network diagram, we found that the intestinal immune network for IgA production is crucial to the TME of TLS in PDAC. Somatic mutation and affinity maturation are processes that take place within the germinal center (GC) and rely on the assistance of T cells. Anatomically restricted dark zone and light zone compartments comprise the GC (12). In GCs, B cells periodically shift from dark to light zones (13). In the bright zone, they interact with helper T cells and FDCs. The interaction of CXCL13 with CD39 enhances the interaction of CD4 T cells with B2 cells (14). Chronic TCR activation and TGF-β signaling affect CXCL13 expression and enriched in TLS⁺ tumors; TGF-β producing fibroblasts serve a crucial role in promoting CXCL13 production by PDAC T cells (9). After class switch recombination and somatic hypermutation (SHM), B1 cells differentiate into lgA⁺ B cells. FasL membrane-induced Fas-mediated mortality is averted by IgA BCR signaling (15). According to experimental findings, IgA is indispensable for optimal GC engagement and GC-derived memory formation, and it regulates the production of intestinal homing plasma cells (15) (Figure 10). In the GC reaction, the complexity of the gut microbiome determines affinity. As a result, there is a correlation between the gut microbiome and the PDAC TME.

Figure 10 The mechanism of intestinal immune network for IgA production in correlation with tertiary lymphoid structures’ germinal center. BCR, B cell receptor; FDC, follicular dendritic cell; GC, germinal center; IgA, immunoglobulin A; IL, interleukin; MHC, major histocompatibility complex; SHM, somatic hypermutation; TGF-β, transforming growth factor beta.

Enrichment analyses by WGCNA indicated that the brown module was correlated with TLS DEGs expression by comparing differential genes in modules. The PPI network constructed by the STRING database represents a classical human leukocyte antigen (HLA) cluster and includes HLA-drb1, HLA-dmb, and HLA-doa, among others. The majority of these genes are attributed to HLA-II. One study found that intratumoral HLA-II presentation was mostly controlled by professional antigen-presenting cells (APCs) instead of cancer cells (16). By designating which genes are granted exclusive access to the HLA-II presentation pathway, we will facilitate enhanced HLA-II-directed cancer therapies that will aim to improve the prognosis of patients with PDAC. The most significant target is CD8a. One study indicated that CD8α is essential for keeping CD8⁺ T lymphocytes in a condition of physiological quiescence in peripheral lymphoid organs. Both naive and memory CD8⁺ T cells acquired activation phenotypes spontaneously following inducible deletion of CD8α; these cells subsequently perished in the absence of specific antigen exposure (17).

Phenotype variants such as stage, gender, T staging, and N staging are factors that are intensively involved in the prognosis and diagnosis of PDAC patients. As confirmed by the Cox regression analyses, TLS constituted a distinct risk factor.

The Kaplan-Meier survival analyses indicated an adverse prognosis associated with elevated CXCL10 expression and BEX1, PEMT, SMIM6, and MIR7-3HG were related to a favorable prognosis. A previous study suggested that the signaling pathway involving CXCL10/CXCR3 in macrophages may enhance their ability to be attracted to the pancreas (18).

By comparing the immune cells’ correlation coefficient, it occurred to us that macrophages M1 increase risk (P=0.002) and B cells naive decrease the risk (P=0.01). In this way, naive B cells play a crucial role in the TME of TLS. Within the GC, these naive B cells transform into GC B cells, exhibiting SHM and immunoglobulin-like switch recombination (CSR) of the BCR in the dark zone (17). Naive B cells are believed to have a significant role in GCs of TLS, where the BCR interacts with FDCs and helper T cells. The current study suggests that microRNAs (miRNAs) play a role in the regulation of B-cell tolerance (19). A study of PDAC biomarkers suggested a role for miR-181b-5p in the regulation of B cells associated with the BCR pathway (20). Potential therapeutic targets can be uncovered by analyzing these miRNA-regulated protein-coding genes.

We conducted a comparative analysis of the distribution of several immune cell types between the TLS-high and TLS-low groups. The findings revealed a statistically significant increase in the proportion of B cell naive, CD8 T cell, and CD4 T cell memory activated in the TLS-high group. The study revealed that B cell naive, CD8 T cell, and CD4 T cell memory activation are crucial in the functioning of the TLS.

The validation cohort IMvigor210 indicates that those in the TLS-high group may exhibit heightened responsiveness to immunotherapy, leading to improved outcomes and prognosis. The cancer immune phenotypes consist of three categories: “inflamed”, “immune excluded”, and “desert”. “Inflamed” refers to tumor cells that express programmed cell death ligand 1 (PD-L1), and the TME has more infiltrating immune cells, especially lymphocytes (21). In the context of “immune excluded”, there are suppressor cells present in the immunological microenvironment that exhibit increased TGF-β signaling, myeloid inflammation, and tumor angiogenesis (22). The characteristic of “desert” manifests a lack of immune infiltration, a lack of antigen presentation (low MHCI), and high tumor growth (22,23). In unsupervised cluster analysis, C3 reflects higher immunity scores than C1 and C2. As a result, C3 represents “inflamed”, C1 represents “immune excluded”, and C2 represents “desert”. Our immunophenotyping of TLS can better guide immunotherapy. Although our model was developed in PDAC, we sought to validate its core biological implication—that a TLS-high phenotype signifies an immunotherapy-responsive TME—in the IMvigor210 cohort of metastatic urothelial cancer. The successful stratification of both prognosis and immunotherapy response in this independent cohort suggests that the immune landscape captured by our signature represents a fundamental, cross-tumor mode of anti-tumor immunity that is susceptible to checkpoint blockade, transcending specific cancer types.

The current state of immunotherapies for PDAC is limited in their ability to improve survival rates (24,25). The resistance of PDAC to immunotherapies can be attributed to two factors: a low mutational burden and a TME that is densely packed, inaccessible, and characterized by immunosuppressive, hypoxic, and fibrotic features. These characteristics collectively render the tumor immunologically inactive (26). By studying the difference in immune escape profile in different expressions of TLS, the sensitivity of specific tissues to immunotherapy can be increased to reduce the immune escape phenomenon. Network pharmacology reveals that the highest degree target is CD4, along with CD8a tertiary. It suggests that chimeric antigen receptor T-cell (CAR-T) therapy against CD4, CD8a, and protein tyrosine phosphatase receptor type C (PTPRC), also known as CD45 targets, may have a potential role in the future in developing immunotherapies, such as anti-PD-1/anti-PD-L1, anti-CTLA-4, and anti-TIGIT. The ongoing inquiry is focused on the utilization of CAR-T treatment for solid malignancies, particularly PDAC. The application of CAR-T therapy for PDAC treatment is mainly impeded by the interstitial barrier, an immunosuppressive milieu, insufficient chemotaxis, and the “ontarget, offtarget” phenomenon (27).

The constraints of our work are as outlined below: Initially, our material was acquired from publicly accessible databases, and the verification process was conducted using data from the same public databases. The validation we conducted is insufficient to encompass all possible variants of PDAC patients from all pertinent regions. Furthermore, because of the small sample size of TCGA and GEO cases, the statistical analysis may not have yielded statistically significant results. Furthermore, the hub targets we acquired were computed using an algorithm. To potentially solve this issue, it may be worth considering external experimental validation as a viable therapy target. We believe it is crucial to subject our research findings in the field of immunotherapy development to rigorous validation at both the transcriptome level and through external experimental verification, given the gravity of scientific research and the significance of this area of study.

Conclusions

This study successfully developed and validated a novel TLS-derived gene signature that serves as a robust and independent prognostic biomarker for PDAC. Beyond its predictive utility, the model provides mechanistic insights into the immune landscape of PDAC, substantiating the positive role of TLS in patient survival and immunotherapy response. Our analysis, guided by the model, underscores the critical function of IgA in the GC response within the TLS, and highlights naive B cells as a pivotal cellular component, potentially through their BCR-mediated interactions with FDCs and helper T cells. Furthermore, we identified potential therapeutic targets from the TLS-related gene network. Ultimately, this prognostic model not only offers a practical clinical tool but also establishes a foundational framework for future research aimed at elucidating the underlying mechanisms and developing targeted strategies to modulate TLS for therapeutic benefit.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1601/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1601/prf

Funding: This study was supported by the General Program of the National Natural Science Foundation of China (grant No. 72274127), the Training Program of the Major Research Plan of the National Natural Science Foundation of China (grant No. 92046015), and the R&D Program of Beijing Municipal Education Commission (grant No. KZ202010025047).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1601/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Dreyer SB, Beer P, Hingorani SR, et al. Improving outcomes of patients with pancreatic cancer. Nat Rev Clin Oncol 2025;22:439-56. [Crossref] [PubMed]
2. Senturk ZN, Akdag I, Deniz B, et al. Pancreatic cancer: Emerging field of regulatory B-cell-targeted immunotherapies. Front Immunol 2023;14:1152551. [Crossref] [PubMed]
3. Zou X, Lin X, Cheng H, et al. Characterization of intratumoral tertiary lymphoid structures in pancreatic ductal adenocarcinoma: cellular properties and prognostic significance. J Immunother Cancer 2023;11:e006698. [Crossref] [PubMed]
4. Schumacher TN, Thommen DS. Tertiary lymphoid structures in cancer. Science 2022;375:eabf9419. [Crossref] [PubMed]
5. Fridman WH, Meylan M, Petitprez F, et al. B cells and tertiary lymphoid structures as determinants of tumour immune contexture and clinical outcome. Nat Rev Clin Oncol 2022;19:441-57. [Crossref] [PubMed]
6. Sautès-Fridman C, Petitprez F, Calderaro J, et al. Tertiary lymphoid structures in the era of cancer immunotherapy. Nat Rev Cancer 2019;19:307-25. [Crossref] [PubMed]
7. Tokunaga R, Nakagawa S, Sakamoto Y, et al. 12-Chemokine signature, a predictor of tumor recurrence in colorectal cancer. Int J Cancer 2020;147:532-41. [Crossref] [PubMed]
8. Zhang K, Xie X, Zou LH, et al. Tertiary Lymphoid Structures Are Associated with a Favorable Prognosis in High-Grade Serous Ovarian Cancer Patients. Reprod Sci 2023;30:2468-80. [Crossref] [PubMed]
9. Vanhersecke L, Bougouin A, Crombé A, et al. Standardized Pathology Screening of Mature Tertiary Lymphoid Structures in Cancers. Lab Invest 2023;103:100063. [Crossref] [PubMed]
10. Chaves-Almagro C, Auriau J, Dortignac A, et al. Upregulated Apelin Signaling in Pancreatic Cancer Activates Oncogenic Signaling Pathways to Promote Tumor Development. Int J Mol Sci 2022;23:10600. [Crossref] [PubMed]
11. Somani VK, Zhang D, Dodhiawala PB, et al. IRAK4 Signaling Drives Resistance to Checkpoint Immunotherapy in Pancreatic Ductal Adenocarcinoma. Gastroenterology 2022;162:2047-62. [Crossref] [PubMed]
12. Kennedy DE, Okoreeh MK, Maienschein-Cline M, et al. Novel specialized cell state and spatial compartments within the germinal center. Nat Immunol 2020;21:660-70. [Crossref] [PubMed]
13. Inoue T, Baba Y, Kurosaki T. BCR signaling in germinal center B cell selection. Trends Immunol 2024;45:693-704. [Crossref] [PubMed]
14. Ukita M, Hamanishi J, Yoshitomi H, et al. CXCL13-producing CD4+ T cells accumulate in the early phase of tertiary lymphoid structures in ovarian cancer. JCI Insight 2022;7:e157215. [Crossref] [PubMed]
15. Raso F, Liu S, Simpson MJ, et al. Antigen receptor signaling and cell death resistance controls intestinal humoral response zonation. Immunity 2023;56:2373-2387.e8. [Crossref] [PubMed]
16. Ran Z, Mu M, Lin S, et al. Deciphering the MHC immunopeptidome of human cancers with Ligand.MHC atlas. Brief Bioinform 2025;26:bbaf314. [Crossref] [PubMed]
17. Zheng L, Han X, Yao S, et al. The CD8α-PILRα interaction maintains CD8(+) T cell quiescence. Science 2022;376:996-1001. [Crossref] [PubMed]
18. Yin H, Chen Q, Gao S, et al. The Crosstalk with CXCL10-Rich Tumor-Associated Mast Cells Fuels Pancreatic Cancer Progression and Immune Escape. Adv Sci (Weinh) 2025;12:e2417724. [Crossref] [PubMed]
19. Roy U, Raghavan SC. Regulation of B-cell development and differentiation by microRNAs during immune response and their implications in immunological disorders. J Immunol 2025;214:3199-207. [Crossref] [PubMed]
20. Şen S, Özgel MÇ, Tunçer ŞB, et al. A Pilot Study of Exploring miRNA-Protein Interaction Networks in Pancreatic Ductal Adenocarcinoma Patients: Implications for Diagnosis and Prognosis. Diagnostics (Basel) 2025;15:2479. [Crossref] [PubMed]
21. Zhu C, Ma J, Zhu K, et al. Spatial immunophenotypes predict clinical outcome in intrahepatic cholangiocarcinoma. JHEP Rep 2023;5:100762. [Crossref] [PubMed]
22. Wu B, Zhang B, Li B, et al. Cold and hot tumors: from molecular mechanisms to targeted therapy. Signal Transduct Target Ther 2024;9:274. [Crossref] [PubMed]
23. Wang T, Song W, Tang Y, et al. Breaking the immune desert: Strategies for overcoming the immunological challenges of pancreatic cancer. Biochim Biophys Acta Rev Cancer 2025;1880:189353. [Crossref] [PubMed]
24. Hu ZI, O'Reilly EM. Therapeutic developments in pancreatic cancer. Nat Rev Gastroenterol Hepatol 2024;21:7-24. [Crossref] [PubMed]
25. Sang W, Zhou Y, Chen H, et al. Receptor-interacting Protein Kinase 2 Is an Immunotherapy Target in Pancreatic Cancer. Cancer Discov 2024;14:326-47. [Crossref] [PubMed]
26. Farhangnia P, Khorramdelazad H, Nickho H, et al. Current and future immunotherapeutic approaches in pancreatic cancer treatment. J Hematol Oncol 2024;17:40. [Crossref] [PubMed]
27. Wehrli M, Guinn S, Birocchi F, et al. Mesothelin CAR T Cells Secreting Anti-FAP/Anti-CD3 Molecules Efficiently Target Pancreatic Adenocarcinoma and its Stroma. Clin Cancer Res 2024;30:1859-77. [Crossref] [PubMed]

(English Language Editor: J. Jones)