Exploring prognostic and immunological characteristics of pancreatic ductal adenocarcinoma through comprehensive genomic analysis of tertiary lymphoid structures and CD8 + T-cells

Purpose Tertiary lymphoid structures (TLSs) and CD8 + T-cells are potential prognostic indicators for pancreatic ductal adenocarcinoma (PDAC). We established a novel scoring system for evaluating the risk for PDAC based on TLS- and CD8 + T-cell-related genes. Methods We analyzed single-cell sequence data from PDAC patients in the Genome Sequence Archive. Bioinformatics and machine algorithms established and validated a scoring method (T-C score) based on PDAC survival-related genes highly expressed in TLSs and CD8 + T-cells. Patients were stratified into the low- and high-T-C score groups. Differences in survival, pathway enrichment, mutation status, immune cell infiltration, expression of immune checkpoint-associated genes, tumor stemness, and response to antitumor therapy were compared through computer simulation methods. Results Overall survival differed significantly between the training and validation cohorts’ low- and high-T-C score groups. The low-T-C score group correlated with lower tumor mutation burden and lower levels of tumor stemness compared with the high-T-C score group. Patients with lower T-C scores exhibited advantages in immunotherapeutic responses and might be more sensitive to the chemotherapeutic regimen and multi-kinase inhibitors. Conclusion The T-C score could serve as an effective model for predicting the survival and therapeutic responses of patients with PDAC. Supplementary Information The online version contains supplementary material available at 10.1007/s00432-024-05824-0.


Introduction
Neoadjuvant and adjuvant treatments, including chemo-, radio-, and immunotherapy, achieve limited success due to the heterogeneity of the pancreatic ductal adenocarcinoma (PDAC) tumor microenvironment (TME).CD8 + T-cells are the main component of tumor-infiltrating lymphocytes in the TME of PDAC; they harbor immune-monitoring functions and can detect antigens secreted by malignant cells B-cell centers and dendritic cells surrounded by a rim of T-cells (Workel et al. 2019).TLSs usually occur at sites of chronic inflammation as well as various cancer types, including pancreatic cancer (Gunderson et al. 2021).However, unlike primary and secondary lymphoid structures, antigenic stimulation is crucial during the development of TLSs, and the formation of TLSs represents an ongoing adaptive immune response (Lutz et al. 2014).Tumor-associated TLSs, which are usually correlated with a good prognosis, are also thought to be a site for immune cell recruitment and activation (Helmink et al. 2020).The presence of TLSs is closely related to the ICB response in melanomas, sarcomas, and urothelial cancers, and TLS infiltration can be increased by ICB treatment (Cabrita et al. 2020;Helmink et al. 2020;Petitprez et al. 2020).In addition, the formation of TLSs is relevant to the CD8 + T-cell infiltration; TLSs with significant plasma cell infiltration are correlated with improved CD8 + T-cell infiltration and overall survival (OS) in ovarian cancer (Kroeger et al. 2016).Recently, a positive correlation between tumor-infiltrating CD8 + T-cells and TLSs in PDAC was shown by Takeshi Tanaka et al. (Tanaka et al. 2023).However, the relationship between the degree of TLSs and CD8 + T-cell infiltration and the ICB therapeutic responses in PDAC remains unclear.In addition, whether high TLSs and CD8 + T-cell infiltration contribute to the OS and therapeutic response in patients with PDAC and the biological mechanism behind them are still poorly illustrated.
In this study, we aimed to establish a novel scoring system for evaluating the risk of PDAC based on TLSs-and CD8 + T-cell-related genes.To do this, we characterized TLSs and CD8 + T-cell-related gene clusters to identify the potential biological functions of TLSs and CD8 + T-cells in PDAC.Firstly, weighted gene co-expression network analysis (WGCNA) was performed to identify CD8 + T-cellassociated genes in PDAC.Secondly, unsupervised cluster analysis was carried out to identify two distinct TLSs and CD8 + T-cells related subtypes.Furthermore, the TLSs-CD8 + T-cell (T-C) score was developed and then tested for its ability to predict survival and responses of patients to chemotherapeutic agents and ICBs.This scoring system can potentially serve as an effective model for predicting the survival and therapeutic responses of patients with PDAC.

Materials and methods
A preprint has previously been published (Hu et al. 2023).The procedures were described in the Supplementary File in detail.

Identification of survival-related TLSAGs and CTCRGs (T-C genes)
CTCRGs in PDAC were identified (please see the Supplementary File, Figures S1-S3).Tumor-infiltrating CD8 + T-cells and TLSs are positively correlated in PDAC [15].So, it is likely that the degree of TLSs, together with CD8 + T-cell infiltration, relate to the overall ICB therapeutic responses in PDAC.However, this has not yet been established.Therefore, for this study, we identified TLSAGs and CTCRGs related to survival.As we described before, 37 TLSAGs were obtained from previously published articles (Table S1), and 31 of these TLSAGs were commonly expressed across databases (the median expression levels of these genes are also displayed in Table S1).These 31 TLSAGs were combined with the previously obtained CTCRGs for univariate Cox regression analysis.As a result, eight TLSAGs and seven CTCRGs with p < 0.01 were significantly related to survival, so we named these genes T-C genes.The univariate Cox analysis of the T-C genes showed that four genes were protective factors with HR < 1, and 11 genes were risk factors with HR > 1 for PDAC prognosis (Fig. 1a).In addition, we detected widespread CNVs in the T-C genes.Copy number gains were most frequent, and MYC, MCL1, AIM2, CD40, LAMP3, BCL6, USHBP1, RGS10, MKI67, HMGB2, and BHLHE40 showed extensive CNVs in amplifications, whereas FAM107A, APO-BEC3G, and GFRA2 showed copy number losses (Fig. 1b).The chromosomal locations of the T-C genes with CNVs are shown in Fig. 1c.Through the GSCA Lite database, we found the most frequent CNVs associated with the T-C genes were heterozygous.MYC, AIM2, MCL1, and CD40 were characterized by high amplification, while GFRA2, CLEC3B, FAM107A, APOBEC3G harbored more deletions (Fig. 1d-e).

Identification of TLS and CD8 + T cell related subtypes (T-C cluster) in PDAC
To explore the relationship between TLSs, CD8 + T Cell, and tumorigenesis, we conducted unsupervised clustering upon the 174 PDAC patients of TCGA-PAAD samples based on the mRNA expression level of the 15 T-C genes.We tried k from 2 to 5, and k = 2 produced the best results in terms of clustering and OS analysis (Fig. 2a, Figure S4a).The two subgroups were designated as T-C cluster-1 and T-C cluster-2, and the patients in T-C cluster-1 had significantly better survival (p = 0.002, Fig. 2b).Furthermore, different TME and mutational features were observed in the two clusters (Supplementary File, Figure S5).S3).Then, the PDAC samples were stratified into high or low T-C score groups based on the median T-C score of 1.121.Figure 4a predicting 1-year,  2-year, 3-year, 4-year, and 5-year survival were 0.774,  0.734, 0.783, 0.717, and 0.661, respectively (Fig. 4d).

Prognostic value of the T-C score for PDAC patients
Next, we validated the prognostic value of the T-C score.The relationship between clinical features and T-C score in the TCGA-PAAD cohort is depicted in Fig. 5a-f.No significant differences were observed in age, gender, and N stage (Fig. 5a-c).At the same time, PDAC, which was better differentiated, smaller in tumor size, or early in the TNM stage, had a significantly smaller T-C score (Fig. 5d-f).Then, the univariate and multivariate Cox regression analyses were

Construction of the T-C score
As we had established two patient clusters with differences in survival, we could then use this information to fully investigate DEGs between the clusters and to identify genes critical for the prognosis of patients with PDAC.Based on the 'limma' package, differential expression analysis between the two clusters was performed, and 1049 DEGs were identified (Fig. 3a).According to these DEGs, the unsupervised clustering method was utilized to divide the TCGA-PAAD cohort into three gene clusters with distinct survival (Fig. 3b).The Kaplan-Meier analysis revealed that compared with other gene clusters, the patients in gene cluster-2 had significantly longer OS (log-rank p = 0.011).In contrast, the patients in gene clusters-1 and − 3 had relatively poor prognoses (Fig. 3c).Based on the T-C gene signatures A and B (Table S2), the dimension reduction was performed by the Boruta algorithm.The transcriptomic profile of the 78 most abundant DEGs identified among the gene clusters was displayed as a heatmap (Fig. 3d).Based on the 78 most abundant DEGs, we calculated the T-C score for each PDAC sample in the TCGA-PAAD cohort using

Functional enrichment analysis between the T-C score groups
Additionally, GSVA analysis was performed to elucidate the difference in bioinformatic functions between the T-C score groups.As depicted by the bar plot, the high T-C score group showed significant enrichment in tumorigenesis-related pathways, such as KEGG_P53_ SIGNALING_PATHWAY, KEGG_WNT_SIGNALING_ PATHWAY, KEGG_NOTCH_SIGNALING_PATHWAY, KEGG_CELL_CYCLE, KEGG_HOMOLOGOUS_ RECOMBINATION, KEGG_TIGHT_JUNCTION, KEGG_THYROID_CANCER, KEGG_SMALL_CELL_ LUNG_CANCER, KEGG_BLADDER_CANCER, KEGG_PANCREATIC_CANCER.In contrast, the low T-C score group was mainly enriched in pathways related to benign diseases (Figure S6a).

The correlation between the T-C score and TIME
To further illuminate the intrinsic immune status that contributed to the distinct survival pattern, the correlation performed in the TCGA-PAAD cohort (Fig. 6a S6d).Similar results could be found in immune analyses using xCELL, EPIC, and MCP-Counter algorithms (Figure S7b-d), which indicated that the low T-C score group could be an immune-active phenotype.The higher levels of B-cells and monocytes may also indicate a stronger signal for TLS infiltration.Moreover, the expression levels of the ICPs-related genes also showed significant differences between the two groups.Among the 28 differentially expressed ICPs-related genes, 21 genes were significantly upregulated in the low T-C score group (Figure S6e), which suggested that the low T-C score group might benefit more from ICB.  analysis (Figure S8e-g).Moreover, patients in the high T-C score group exhibited significantly higher gene alterations involved in amplifications, deletions, and TMB (Figure S9a-c).The results demonstrated that patients with high T-C scores had worse survival, which coincided with the notion that the higher mutation accumulation in cancer negatively correlated with the OS of patients (Pleasance et al. 2020).Furthermore, the T-C score is related to PDAC tumor stemness (Supplementary File, Figure S10), and relationships were observed between the T-C score and the sensitivity to treatments (Supplementary File, Figure S11).

Discussion
In this study, we demonstrated the presence of TLSs was coupled with higher infiltration of tumor-associated CD8 + T cells from a single-cell perspective and then established a scoring method (T-C score) based on the consensus clustering of the survival-related TLSs-and patients in the high T-C score group exhibited a higher frequency of genomic alterations in general (98.81% in the high T-C score group vs. 63.51% in the low T-C score group).As the most mutated genes in PDAC, the mutation frequencies of KRAS and TP53 in the high T-C score group were 88% and 75%, respectively, which were remarkably higher than those in the low T-C score group (31% and 41%, respectively) (Figure S8c).The alteration frequencies of RNF213 and ADGRL2 also significantly differed between the high and low T-C score groups (Figure S8c).Furthermore, extensive co-occurrence was found among the top 25 mutated genes, including between TP53 and KRAS, CDKN2A and KRAS, CDKN2A and TP53, CACNA1B and SMAD4, RREB1 and TTN, RREB1 and MUC16, ADGRL2 and RNF213, PCDH9 and ADAMTS16 (Figure S8d).(Connor et al. 2019;Notta et al. 2016).KRAS, TP53, CDKN2A, and SMAD4 are among the earliest genetic variants in PDAC (Bazzichetto et al. 2020).In the current research, though not all significant, higher mutation frequencies in KRAS (88% vs. 31%), TP53 (75% vs. 41%), CDKN2A (24% vs. 12%), and SMAD4 (26% vs. 19%) were found in the high T-C score group, the mutation analysis revealed a positive correlation between the TMB and T-C score, as well as between CNV and T-C score.Notably, most patients in the T-C cluster-2 were in the high T-C score group.These findings suggested that gene mutations in PDAC might drive the profile of TLSs and CD8 + T-cell infiltration toward a pattern of high T-C score in PDAC.Stem-like cell phenotypes are closely related to many key steps in PDAC progression and the cancer metastatic cascade, as cancer cells expressing those features exhibit increased cell motility and invasiveness (Yang et al. 2022;ZhangLiu et al. 2022;Zhou et al. 2022).Therefore, we estimated PDAC stemness using the OCLR algorithm.Higher levels of stemness score were found in the high T-C score group, and a positive correlation was found between T-C score and stemness score, indicating more mutations in genes that encode oncogenes and epigenetic modifiers, frequent perturbations in specific mRNA/ miRNA transcriptional networks, and more deregulation of signaling pathways for PDAC patients in the high T-C score group (Malta et al. 2018).
CD8 + T-cell-associated gene expression profiles in PDAC patients.
Previous studies mainly focused on the role of T cells, highlighting the association of tumor-infiltrating T lymphocytes with patients' prognosis in numerous tumors, including PDAC (Balachandran et al. 2017;Carstens et al. 2017;Fridman et al. 2017;Germain et al. 2014;Rouanne et al. 2021;Sautes-Fridman et al. 2016;Wirsing et al. 2014).However, recently, B-cells and TLSs were also discovered to be critical components in the TIME in terms of patient prognosis and response to immunotherapy (Cabrita et al. 2020;Germain et al. 2014;Helmink et al. 2020;Petitprez et al. 2020;Thommen et al. 2018).Inducing the production of TLSs has been suggested as a new direction for individualized immunotherapy for patients with pancreatic cancer (Zhou et al. 2021).TLSs represent a special TIME status and serve as a beacon for antigen presentation and T-cell activation.Previous publications also demonstrated that different types of TLSs may exist in individual tumors independent of the spatial location of the TLS (Cabrita et al. 2020;Horeweg et al. 2022).In this study, the investigation of 31 TLSAGs identified eight significantly related to survival and included in the T-C score, supporting the previous research to some degree.
This study identified 638 CTCTGs closely correlated with CD8 + T cell abundance in PDAC samples.A close relationship between CD8 + T-cells and TLSs is found in many cancers.Cabrita et al. (Cabrita et al. 2020) found that the presence of TLSs was coupled with increased tumorassociated CD8 + T-cell infiltration in melanoma, and the combination of infiltration of both TLSs and CD8 + T-cells was associated with better survival outcomes compared to CD8 + T-cells infiltration alone.A positive correlation between CD8 + T-cell infiltration and TLS presence was also observed in PDAC by Gunderson et al. (Gunderson et al. 2021).
Although the roles of TLSs and CD8 + T-cells in tumor biology and antitumor immunity have been widely studied, the potential clinical translation is still hampered by the heterogeneity of TLSs and their complex relationship with CD8 + T-cells.By focusing on developing a scoring method that utilized the gene signatures of both TLSs and CD8 + T-cells in PDAC, we aimed to better stratify pancreatic cancer.In this study, the survival-associated TLSAGs and CTCRGs were selected by univariate Cox analysis, and based on that, PDAC samples were initially stratified into two T-C clusters based on consensus clustering.The patients in T-C cluster-1 had relatively better survival and lower TMB compared to T-C cluster-2 patients, which was consistent with the notion that mutation accumulation in cancer has a negative correlation with the OS of patients (Pleasance et al. 2020).The T-C cluster-1 also had higher Nevertheless, this study proposes a novel score for the prognosis of PDAC.Part of the score can be obtained from the pathological examination of the surgical specimen.The other part of the score relies upon advanced gene amplification and quantification single-cell techniques.We agree that such techniques might not be accessible to hospitals yet, but such technologies are becoming more reliable and available and less expensive.Hence, this score could be determined by any laboratory possessing adequate equipment and could contribute to improving the prognosis of PDAC.

Conclusion
In summary, we integrated transcriptome data and the expression of TLSs-and CD8 + T-cell-related genes to establish a new independent prognostic marker for PDAC, the T-C score.This score could be used for the prediction of OS, immune infiltration status, tumor stemness, tumor mutation, and sensitivity to immune therapy and chemotherapeutics for PDAC.Therefore, the T-C score may assist with diagnosing and treating patients with PDAC.
The TIME is closely related to the clinical prognosis of patients with tumors (Chen et al. 2020).Different types of immune cell populations, such as lymphocytes and myeloid-originated cells, comprise a major part of the TIME (Zhang et al. 2020).Through ssGSEA and CIBER-SORTx, we recognized that the low T-C score group was a more immune-activated phenotype with significantly higher infiltrations of naïve B cells, CD8 + T cells, activated NK cells, and monocytes.In contrast, the high T-C score group was heavily infiltrated by M0 macrophages, which was found to be an independent risk factor for various cancers (Farha et al. 2020;Huang et al. 2020;Jairath et al. 2020;ZhangZou et al. 2022).Recent advances have indicated that the active pre-treatment immune status could positively impact immunotherapy (Chabanon et al. 2021;D'Alise et al. 2021).More tumor-infiltrating CD8 + T cells in the low T-C score group would be expected to have a marked influence on immunotherapy, as the classical mechanism behind the anti-PD-1/PD-L1 treatment is that the mutual suppression effect conducted by tumor cells and tumor-infiltrating CD8 + cytotoxic T cells (CTLs) could be abolished by ICBs.So, CD8 + CTLs could regain their cytotoxicity against tumor cells (Ishii et al. 2020).In light of the previously reported capability of TIME status in predicting sensitivity to ICB treatment in patients with cancer (Ouyang et al. 2021), we analyzed the ICP gene expression, and a higher level was found in the low T-C score group.In accordance with this result, a significantly higher response rate to ICBs was observed in the same group in TIDE analysis.Though classically, a high TMB has been regarded as a predictive factor for the response to ICB treatment (Chen et al. 2019;Cristescu et al. 2018), a different result was recently published regarding all solid cancer types (McGrail et al. 2021).Thus, the T-C score may be an alternative marker for predicting therapeutic response to ICBs.Emerging evidence has exhibited that multi-kinase inhibitors alone or synergistically used with HDAC inhibitors could suppress PDAC tumor growth (Dent et al. 2021;Yu et al. 2022).In the current study, the IC50 and CMap analysis demonstrated a possible higher sensitivity for using multi-kinase and HDAC inhibitors in the low-risk group, which provided another therapeutic possibility for these PDAC patients.
This study had limitations.Bioinformatics-based investigations are subjected to the limitations of the algorithms used.A batch effect may occur even after correction due to different RNA sequencing methods and platforms.The ICB therapeutic response was predicted, and the clinical data in the ICGC-CA, E-MTAB-6134, and TCGA-PAAD cohorts did not exactly match, which may limit the validation process.Large-scale prospective studies are required to verify the value of the T-C score in predicting clinical benefits for PDAC patients.

Fig. 2
Fig. 2 Stratifying the PDAC patients in the TCGA-PAAD cohort according to the expression profile of T-C genes and the mutational landscape between T-C clusters.(a) Clustering heatmap of the T-C genes in the TCGA-PAAD cohort (n = 174).(b) Survival analysis of the two T-C clusters.(*p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001; Ns, not significant).TMB, tumor mutation burden Fig. 1 Identification of T-C genes and genetic mutational landscape of the T-C genes in PDAC.(a) Eight TLSAGs and seven CTCRGs were recognized by the univariate Cox analysis as T-C genes with p < 0.01, forest plot showing the results of the univariate Cox regression analysis.(b) Frequencies of CNV gain, loss, and non-CNV among T-C depicts the distribution of samples among T-C clusters, gene clusters, and T-C score groups.Compared with the patients in the low T-C score group, those with high T-C scores indicated a poorer prognosis (log-rank p < 0.001, Fig. 4b).The scatter plot depicted the survival status of individual patients and the T-C score distribution (Fig. 4c).The AUC values for

Fig. 3
Fig. 3 The construction of the T-C score.(a) The volcano map shows the differentially expressed genes (DEGs) between the two T-C clusters.(b) Consensus matrix heatmap defining three gene clusters according to the DEGs.(c) Survival analysis of the PAAD patients in three gene clusters (log-rank p = 0.011).(d) Unsupervised clustering of the DEGs further stratified PAAD patients into three gene clusters.Heat map showing the clinical features of the three gene cluster -b), demonstrating that the T-C score was an independent risk factor for PDAC.According to time-dependent AUC values, the T-C score displayed good predictive ability for OS in comparison with age, gender, differentiation grade, or stage in the TCGA-PAAD cohort (Fig. 6c).Furthermore, we calculated T-C scores for our validation cohort, including ICGC-CA, E-MTAB-6134, and the joint cohort of GSE71729 and GSE85916, and in each cohort, the median T-C score was used to stratify patients into the high or low T-C score groups.Then Kaplan-Meier analyses (Figure S4c-e) and univariate and multivariate Cox regression analyses were carried out (Fig. 4g-h), and the same trend as in the TCGA cohort was found throughout all validation cohorts.The survival status of individual patients and the distribution of the T-C score for each validation cohort are displayed in Figure S4f-h.ROC curves for each validation set are depicted in Figure S4i-k, and Fig. 6d shows the comparison of AUC values for the TCGA-PAAD, ICGC-CA, and E-MTAB-6134 cohorts.

Fig. 4
Fig. 4 The construction of the T-C score.(a) A Sankey plot was used to reveal the correlation between T-C clusters, gene clusters, T-C score groups, and survival of the PDAC patients in the TCGA-PAAD cohort.(b) Survival analysis of the PDAC patients of TCGA-PAAD in Low and High T-C score groups (logrank p < 0.001).(c) Distribution of the T-C score and survival status of PDAC patients.(d) ROC curves for the 1-, 2-, 3-, 4-, and 5-year survival times based on the T-C score Given the tight correlation between genomic alterations and tumor immunity, we implemented somatic mutation and CNV analyses to investigate the genomic alterations between T-C score groups.The top 25 genes with the most frequent alteration between T-C score groups in the TCGA cohort were visualized in oncoplots (FigureS8-a, b).The between the T-C score and the PDAC TIME was investigated.The ESTIMATE algorithm was applied to the TCGA-PAAD cohort, and the correlation analyses indicated that the T-C score was significantly and negatively correlated with the ESTIMATE score (R=-0.23,p = 0.002; tau=-0.170,p = 8.95e-04), immune score (R=-0.24,p = 0.002; tau=-0.171,p = 8.31e-04), and stromal scores (R=-0.19,p = 0.01; tau=-0.144,p = 4.86e-03) (Figure S6b).Meanwhile, the low T-C score group demonstrated significantly higher levels of ESTIMATE score, immune score, and stromal score (Figure S6b).Furthermore, the relationship between tumor-infiltrated immune cells and the T-C score was evaluated.The ssGSEA heatmap (Figure S6c) showing the infiltration of 28 immune cell types indicated the low T-C score group was more immune infiltrated compared to the high T-C score group.In the boxplot (Figure S7a), the low T-C score group had significantly higher levels of infiltration in 14 out of 16 immune cell types that were differently infiltrated between the groups.The CIBERSORTx web tool was also used to predict immune status between the T-C score groups, which exhibited significantly escalated infiltration of naive B-cells, CD8 + T-cells, and monocytes.It activated NK cells in the

Fig. 5
Fig. 5 The prognostic value of the T-C score in PDAC.Comparison of the T-C score among clinical features including age (a), gender (b), N stage (c), histological grade (d), T stage (e), and TNM stage (f).OS, overall survival; AUC, area under curve

Fig. 6
Fig. 6 The prognostic value of the T-C score in PDAC.Both univariate (a) and multivariate (b) Cox regression analyses suggested that the T-C score was a significant prognostic factor for OS in the training (TCGA-PAAD) and validation (ICGC-CA and E-MTAB-6134) cohorts.(c) The time-dependent AUC values of gender, age, histological grade,