Immunological Classification of Pancreatic Carcinomas to Identify Immune Index and Provide a Strategy for Patient Stratification

Background Cancer immunotherapy has produced significant positive clinical effects in a variety of tumor types. However, pancreatic ductal adenocarcinoma (PDAC) is widely considered to be a “cold” cancer with poor immunogenicity. Our aim is to determine the detailed immune features of PDAC to seek new treatment strategies. Methods The immune cell abundance of PDAC patients was evaluated with the single-sample gene set enrichment analysis (ssGSEA) using 119 immune gene signatures. Based on these data, patients were classified into different immune subtypes (ISs) according to immune gene signatures. We analyzed their response patterns to immunotherapy in the datasets, then established an immune index to reflect the different degrees of immune infiltration through linear discriminant analysis (LDA). Finally, potential prognostic markers associated with the immune index were identified based on weighted correlation network analysis (WGCNA) that was functionally validated in vitro. Results Three ISs were identified in PDAC, of which IS3 had the best prognosis across all three cohorts. The different expressions of immune profiles among the three ISs indicated a distinct responsiveness to immunotherapies in PDAC subtypes. By calculating the immune index, we found that the IS3 represented higher immune infiltration, while IS1 represented lower immune infiltration. Among the investigated signatures, we identified ZNF185, FANCG, and CSTF2 as risk factors associated with immune index that could potentially facilitate diagnosis and could be therapeutic target markers in PDAC patients. Conclusions Our findings identified immunologic subtypes of PDAC with distinct prognostic implications, which allowed us to establish an immune index to represent the immune infiltration in each subtype. These results show the importance of continuing investigation of immunotherapy and will allow clinical workers to personalized treatment more effectively in PDAC patients.


INTRODUCTION
Pancreatic ductal adenocarcinoma (PDAC) is one of the most invasive and lethal malignancies. For patients diagnosed with PDAC, the prognosis remains extremely poor, with less than a 10% survival rate (1,2). Most patients with PDAC are diagnosed as advanced or metastatic, and the median survival time is less than 1 year. Conventional approaches using chemotherapy and radiotherapy have only moderately increased the overall survival time of patients (3). In recent years, immune checkpoint blockade (ICB) treatment based on programmed cell death protein 1 (PD-1) and its programmed death ligand 1(PD-L1) have produced promising clinical results in many types of cancers (4). However, PDAC is widely considered to be a poorly immunogenic type of cancer, and ICB treatment has thus far shown a low response rate in the treatment of PDAC. Indeed, the objective response rate in unselected PDAC patients is only 3.1% in clinical trials (5). Therefore, acquiring a complete understanding of the heterogeneity of the immune responses of patients and the mechanisms that underline the effectiveness of immunotherapy treatment for PDAC is of great importance.
In the past decade, substantial progress in research on PDAC molecular subtypes has increased our understanding of molecular pathogenesis (6,7). Given the advancement, many attempts have been made to evaluate the expression of biomarkers to identify patients most likely to benefit from ICB. However, low patient numbers have limited observation opportunities for studying the associations between treatment response and PD-L1 expression, or microsatellite instability status. Ines de Santiago et al. performed an integrative metaanalysis of 353 patients from four different studies to derive the PDAC classification based on immunological parameters (8). The results of this study revealed that the expression characteristics related to tumor-infiltrating immune cells in different PDAC subtypes may assist in guiding immunotherapy treatments. The success of cancer immunotherapy treatment depends on multiple key steps that are involved in immune activation during the cancer immune cycle (9). The low response rates to this kind of treatment and the limited number PDAC patients who have benefitted from ICB have been attributed to low immunogenicity and multiple immunosuppressive mechanisms (10,11). Balachandran et al. highlighted the key barriers that limit immunotherapy efficacy, including cytokines, i m m u n e c e l l t y p e s , a n d c e l l u l a r c o m p o n e n t s o f immunosuppression (12). Thus, more detailed immune characterization of PDAC is needed in the investigation of novel therapeutic strategies.
In this study, we implemented the ssGSEA approach to evaluate immune characteristics based on the marker genes of 119 immune cells, then classified PDAC into different immune subtypes (ISs) based on the immune scores. Subsequently, the response patterns of different ISs to immunotherapy were analyzed to verify their reproducibility in different datasets. We then established the immune index to reflect the different degrees of immune infiltration in the patients. Finally, based on the coexpression network analysis, we identified potential prognostic markers related to the immune index and functionally validated them in vitro.

Data Source
Gene expression profiles and clinical follow-up information data on pancreatic cancer were downloaded from the International Cancer Genome Consortium (ICGC) data portal, which contains a total of 237 samples after excluding probes with empty values and samples without clinical data. The TCGA-PAAD gene expression data (n = 177) and the corresponding clinical data were obtained from the UCSC Xena (http://xena.ucsc.edu/) website. An additional 250 PDAC gene expression microarray data and clinical data were derived from the NCBI Gene Expression Omnibus (GEO) database as validation cohorts, including GSE28735 (n = 42) (13), GSE57495 (n = 63) (14), GSE62452 (n = 66) (15), and GSE85916 (n = 79). Using the "Combat" algorithm in R package "sva," the batch effects were removed, and the expression values were quantile-normalized across the different samples. The scaled estimate values derived from RNA-seq by Expectation Maximization (RSEM) were converted to transcripts per million (TPM) values by multiplying them by one million. A total of 100 immune gene signatures were obtained from the R package "IOBR" (Table S1) (16), and after which we calculated the immune scores using the single-sample gene set enrichment analysis (ssGSEA) algorithm in the "GSVA" package (17).

Clustering Analysis
Based on the normalized enrichment score of the immune characteristics, the R package "ConsensusClusterPlus" (18) was used to construct a consistency matrix to classify the samples via clustering to generate the immune subtypes in PDAC. The "PAM" algorithm and "1-Pearson Correlation" were used to measure distance, and bootstraps were performed 500 times. Each bootstrapping process included 80% of the patients in the training set. The number of clusters was set from 2 to 10, and the optimal classification was determined through calculating the consistency matrix and the consistency cumulative distribution function.

Immunoprofiling
To investigate the relative abundance of tumor-infiltrating immune cells, we performed CIBERSORT (19) to quantify the distribution of the 22 immune cell types in individual specimens of the ICGC cohort. Immune-related signatures of PDAC, including the Th1/IFNg gene, cytolytic immune activity (CYT), and 47 immune checkpoints, were obtained from previous publications (20)(21)(22). Additionally, we utilized the Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) tool to determine the score that represents the proportion of immune and stromal cells (23). We used the "TIDE" software (24) with default parameters to analyze the differences in immune efficacy. The anti-CTLA4, anti-PD1, and anti-PD-L1 treatment profile data were obtained from previous studies (25,26).

LDA and Construction of the Immune Index
As different subtypes have different gene signatures, we conducted a linear discriminant analysis (LDA) for dimensionality reduction and established a subtype classification feature index to better quantify the immune characteristics of patients in different sample cohorts (27,28). In order to do this, z-transformation was performed on each feature across 14 prognostic-related immune features, and the centroid of each group was dispersed as much as possible based on Fisher's LDA optimization standards. The goal was to find one linear combination, "A," that could maximize variance across all classes and ensure that the first two features of the model could clearly distinguish samples of the different subtypes. The LDA score was computed for the discriminative functional markers by adding the first two linear discriminants, LD1 and LD2, from the LDA of the combined normalized data to produce an immune index. The classification performance of the immune index in different subtypes was determined by the area under the curve (AUC) of multiclass receiver operating characteristic (ROC) curves.

Weighted Correlation Network Analysis
The R software package "WGCNA" (29) was used to identify the co-expression modules in immune genes. Genes with the top 5,000 standard variations were retained and subjected to the clustering analysis in the TCGA cohort. The co-expression network conformed to the scale-free network with beta values ranging between 1 and 20. Meanwhile, the linear model was established via logarithms of the adjacency degree of a node (log k) and the appearance probability of the node [log(p(k))] with a correlation greater than 0.85. To ensure a scale-free network, the nearest soft threshold was selected and used to filter the coexpression module. In the next step, the expression matrix was converted into an adjacency matrix and then converted into a topological matrix. Based on the topological overlap matrix (TOM), we used the average-linkage hierarchical clustering method to cluster genes in order to maintain the minimum number of genes in each module of base 30 according to the standard of the hybrid dynamic cut tree. After this, we calculated the eigengenes of each module to perform a cluster analysis on the modules, then merged the modules that were closer to each other into a new module with the following characteristics: height = 0.25, DeepSplit = 2, and minModuleSize = 30. The gene networks were built and visualized by the software Cytoscape version 3.8.0 (30).

Prognostic and Functional Enrichment Analysis
The "coxph" function of the "survival" package in R was used to perform univariate Cox analysis on significant gene modules. The prognostic performance was determined by analyzing the HR score of each gene module based on its degree of significance with p <0.05. Moreover, the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and gene ontology (GO) analyses were carried out using the "clusterProfiler" package in R (31).

Cell Culture
The human pancreatic ductal cell line PANC-1 and pancreatic carcinoma cell HPDE6-C7 were purchased from the American Type Culture Collection (ATCC, Manassas, VA, USA). The PANC-1 and HPDE6-C7 cells were maintained in DMEM medium (Gibco, Grand Island, USA) with 100 U/ml of penicillin, 100 mg/ml of streptomycin, and 10% fetal bovine serum (Gibco, Grand Island, USA) in an incubator containing 5% CO 2 at 37°C. The cell lines were subcultured every 2 to 3 days following digestion with 1 ml 0.25% trypsin and 0.02% EDTA (Sigma-Aldrich, MO, USA). The viability was measured by the number of surviving cells as a percentage of the total number of cells. The average viability of over 95% was determined by Trypan Blue staining at 37°C in an incubator containing 5% CO 2 .

Small Interfering RNA and Transfection
The small interfering RNAs (siRNAs) were purchased from GEMA Gene Company (Pudong, Shanghai, China), and their sequences are listed in Table S2. We designed triplicate sets for each gene, and the cells were seeded into six−well plates that measured 3 × 10 5 cells per well.

Cell Invasion Assay
The cell invasion assay was performed using 8-mm Matrigelcoated Transwell inserts (Corning Costar, NY, USA) to evaluate invasion capacities in vitro. Dulbecco's modified Eagle's medium (DMEM) containing 10% FBS was added to the lower chamber.

Statistical Analysis
The unpaired Student's t-test was used for comparison between two continuous variables and a normally distributed variable. Non-normally distributed variables were analyzed using the Wilcoxon rank-sum test. To compare three or more groups, analysis of variance (ANOVA) and the Kruskal-Wallis test were performed on both the parametric method and the nonparametric method.

Identification of Immune Subtypes in PDAC Based on Immune Gene Signatures
First, we implemented the ssGSEA approach to calculate the absolute enrichment scores of 110 immune features in the "IOBR" package from the ICGC and TCGA cohorts. Next, we used a total of 250 PDAC samples with complete follow-up and overlapping immune characteristics from GSE28735, GSE57495, GSE62452, and GSE85916 as validating cohorts and pretreated them by removing batch effects, thus leading to an independent cohort. Data before and after normalization were inspected using principal component analysis (PCA), which revealed that the batch effect was successfully removed using the "ComBat" algorithm ( Figures 1A, B). We then performed a univariate Cox analysis, which showed that 30 immune features in the ICGC cohort, 40 immune features in the TCGA cohort, and 10 immune features in the GEO cohort were significantly associated with PDAC prognosis ( Table S3). As shown in Figure 1C, there were few intersections across these three cohorts, indicating that the immune features were inconsistent among the datasets of the different platforms. Therefore, we selected 14 immune prognostic-related risk features over at least two cohorts for subsequent analysis (p < 0.05). The distribution and abundance of these 14 immune characteristics in the three cohorts are shown in Figure 1D. We then classified 237 PDAC samples, which we used as a training cohort, from the ICGC cohort according to the 14 immune features.
Using consensus clustering cumulative distribution functions (CDF) and the CDF Delta area curve, we determined that the optimal cluster number is three, at which point relatively stable clustering results can be applied and result in three PDAC ISs ( Figures 1E-G). Upon further analysis of the prognostic characteristics of the three ISs, we observed that IS1 has the poorest prognosis, while IS3 has favorable prognosis in all three cohorts ( Figures 1H-J), suggesting that the three immune subtypes showed consistency in different PDAC cohorts. Consequently, we compared the levels of the TNM staging system, clinical stage, and grade among the three ISs in the TCGA cohort, which was revealed to be in line with the survival data. IS3 patients account for the lowest proportion of M1 (metastasis) and the highest proportion of low-grade patients (G1, G2) ( Figure S2).

The Differences in the Innate and Adaptive Immune Signatures Among the Three ISs
It is well-known that the cyclic GMP-AMP synthase-stimulator of interferon genes (cGAS-STING) signaling pathway is an innate immune pathway that can induce the release of type I IFNs and other inflammatory factors by recognizing foreign cytosolic DNAs (cDNAs) to promote innate immunity. In addition, this pathway also functions as a detector of self-DNA released from tumor cells and dying cells (32,33). Activation of this pathway is significantly associated with tumor progression, and its role in cancer immunotherapy has been well identified in recent years (34)(35)(36)(37)(38). This includes pancreatic cancer, as it exhibits a strong connection to the type I interferon pathway (39)(40)(41)(42)(43). We compared the expression levels of four key genes in the cGAS-STING signaling pathway across the three ISs, including cGAS (encoding cGAS protein), TMEM173, TBK1, and IRF3 (encoding STING protein). As shown in Figure S1, there were no significant differences in the gene expressions of the four genes among the three ISs. We also found that IFNg, a type II interferon that is predominantly produced by T helper (TH) CD4 and CD8 cytotoxic T lymphocyte (CTL) effector T cells during antigen-specific immunity, was significantly upregulated in IS3 (Kruskal-Wallis test, p = 3.4e−06) ( Figure 2A). In line with these results, the cytolytic immune activity (CYT) score, which was based on the cytotoxic T-cell (CTL) markers proposed by Rooney et al. (21), consistently showed the identical trend that IS3 held the highest cytotoxic immunity (Kruskal-Wallis test, p = 2.6e−13) ( Figure 2B) and was associated with positive prognosis. All together, these results suggest that our classification system accurately reflects the differences in adaptive immunity in pancreatic cancer patients rather than innate immunity.

Difference in Immune Profiles Among the Three ISs
In the ICGC cohort, we performed CIBERSORT analysis to evaluate the proportions of the 22 immune cells in each sample ( Table S4). The distribution of the immune cell proportions in all samples is shown in Figures S3A, B, while the proportions across the three ISs are shown in Figures S3C, D. These results exhibit the significant differences in the immune characteristics of different subgroups. Notably, IS3 has the highest proportion of CD8 + T cells, resting memory CD4 T cells, and B cells, as well as the lowest proportion of M0 and M2 macrophages. Immune infiltration analysis showed that IS3 has the highest immune infiltration, followed by IS2, while IS1 had the lowest stromal and immune score (Figures 3A-C; Table S5). All together, we observed that IS3 patients presented a preferable cell-mediated immune response and humoral immune response. Subsequently, we implemented the TIDE algorithm to evaluate the potential clinical effects of immunotherapy in the three ISs we defined. As shown in Figure S4A, there was apparently no difference in the TIDE scores of the three ISs. At the same time, IS3 had the highest T-cell dysfunction score compared with IS1 (p < 0.0001) and IS2 (p < 0.0001), but the lowest score in T-cell exclusion (p < 0.0001) (Figures S4B, C; Table S6). These results may partially explain why T-cell infiltration was the highest in the IS3 patients, but the objective response rates of ICB in clinical trials were extremely low. Therefore, there are likely other underlying mechanisms that lead to T-cell dysfunction due to sustained antigen exposure. We used submap analysis to compare the similarity between the gene signatures of the three ISs and melanoma patients treated with anti-PD1 and anti-CTLA4 (44), which revealed that IS3 patients were more sensitive to PD-1 inhibitors than the other two subtypes ( Figure 3D). We also analyzed the responsiveness of traditional chemotherapy drugs, including gemcitabine, cisplatin, erlotinib, and fluorouracil (5-FU), in all three ISs. IS1 was more sensitive to gemcitabine (p < 0.01 and p < 0.0001, respectively), cisplatin (p < 0.0001), and erlotinib (p < 0.0001) compared with IS2 and IS3, while IS3 was more sensitive to 5-FU (p < 0.0001) ( Figures 3E-H). Thus, our overall computational calculations predicated that IS3 has the highest immune cell infiltration. Therefore, stratification might be helpful in patient selection for ICB trials or further studies on drug resistance, regardless of the wide ineffectiveness with PDAC treatment so far.

Construction of the Immune Index and Its Related Immune Characteristics
To further assess the differences in the gene signatures across the three ISs, as well as quantify the immune characteristics of individual patients, we performed LDA analysis. As shown in Figure 4A, the first linear discriminants (LD) of the model clearly distinguish between the samples of different subtypes. Following this, we calculated the LDA score of each patient in the three cohorts, which make up the immune index. The differences in the scores were significant across the three ISs (Table S7; Figures 4B-D). The classification performance of the immune index in the different subtypes is depicted in Figure 4E with a multiclass comprehensive prediction AUC of 0.89. After applying the immune index to the TCGA and GEO cohorts, we observed a similar performance to the ICGC dataset (multiclass AUC = 0.84 and 0.78, respectively) ( Figures 4F,  G). Together, a high index indicated a higher immune infiltration as seen in IS3 and vice versa as seen in IS1. Additionally, to observe the relationship between the immune index and immune cell characteristics, the immune scores in 28 immune cells were calculated using ssGSEA analysis (Table S8) (45). Most of the immune cells in IS3 had significantly higher enrichment scores than IS1, including activated B cells, activated CD8 T cells, activated CD4 T cells, T follicular helper cells, myeloid-derived suppressor cells (MDSCs), and natural killer cells, all of which were significantly positively correlated with the highest immune index (Figures S5A-C).
Recent studies have revealed that the recruitment of tumorassociated immune cells that occurs during inflammation that accompanies tumor progression has been shown to promote tumor growth and contribute to angiogenesis, invasion, and metastasis (46)(47)(48). Therefore, we also analyzed the distributions of inflammation-related genes of the ISs in the ICGC cohort, the results of which revealed that the monocyte/myeloid lineage (HCK), T cells (LCK), major histocompatibility class II complex (MHC-II) molecules, and gene clusters are all highly expressed in IS3 compared with IS1 and IS2 ( Figure 5A; Table S9). All seven inflammation-related metagenes were then quantified with the enrichment score using the ssGSEA approach as shown in Figure  S6A, which revealed that the enrichment score of the five metagenes in IS3 are significantly higher than in IS1 and IS2 (p < 0.0001), with the exception of MHC-I and interferon. The relationship between the LDA score and the enrichment score was identified using Pearson's correlation coefficient. Based on this, we found that antitumor immune metagenes, including HCK, LCK, MHC-II, and STAT1 were significantly positively correlated with the immune index (p < 0.001) ( Figure S5B). Figure 5B. Notably, the LDA score has a significant positive correlation with many critical inhibitory immune checkpoints [i.e., CTLA4, BTLA, IDO1, PDCD1 (PD-1), IDO2, LAG3, etc.]. Of those, CTLA4 (r = 0.72, p < 0.0001) and PDCD1 (r = 0.54, p < 0.0001) show high correlation, whereas no significant correlation between the immune index and CD274 was discovered (r = 0.04, p = 0.530) (Figures 5B-D). We then analyzed the difference between immune index and melanoma patients treated with anti-PD1 and anti-CTLA4 as well as urothelial carcinoma patients treated with PD-L1 under different responses. The LDA score was significantly upregulated in patients with an immunotherapy response compared with SD or PD patients in anti-CTLA4 and anti-PD1 treatments ( Figures 5F, G). However, there was no big difference in anti-PD-L1 patients across the three ISs ( Figure 5H). These results reinforce the results of the submap analysis in Figure 3D that indicate that immune index can be a good representation of the gene profile. It also further suggests that patients with high immune infiltration (IS3) might benefit from anti-PD-1 therapy.

Identification of the Co-Expression Gene Module of the Immune Index
According to the weighted correlation network analysis (WGCNA) method described previously, samples were clustered ( Figure S6A) and immune gene co-expression modules were identified with a soft threshold of 5 for the scale-free network (Figures 6A, B). As a result, a total of 18 gene modules were generated based on a dissimilarity measure (1-TOM), where the gray module was a collection of genes that could not be gathered into other modules (Figures 6C, S6B;  Table S10). Furthermore, we analyzed the correlation between each module and the age, gender, clinical stage, grade, IS1, IS2, IS3, and LDA score of the patient. As shown in Figure 6D, the tan module showed a significant positive correlation with IS1 and a negative correlation with IS3 and its LDA score. On the other hand, the light cyan module showed the opposite result. The genes in these two modules are significantly associated with gene significance and module membership (IS1 is tan and IS3 is light cyan) (Figures 6E, F). Next, we investigated the correlation between each of these gene modules and the immune index. In this analysis, each of these 13 modules were significantly associated with the LDA score ( Figure 6G) where tan is shown as a risk factor (p < 0.0001, HR = 1.4) and light cyan is shown as a favorable factor (p = 0.002, HR = 0.76) ( Figure 7A; Table S11), as expected. We subsequently screened the co-expressed genes of prognosticrelated modules with co-expression weights greater than 0.35 and identified 11 hub genes ( Figure 7B; Table S12). These hub genes reflected the immune index and showed up as potential candidates for biomarkers. Functional enrichment analysis revealed that genes within tan modules were enriched for biological terms (GO) or KEGG pathways with strong significance in homologous recombination, the p53 signaling pathway, DNA replication, and the cell cycle ( Figure S7). Meanwhile, the light cyan module was significantly associated with immune-related processes including positive regulation of lymphocyte activation, positive regulation of cell activation, regulation of T-cell activation, etc. (Figures 7C, D). Among these hub genes, high expressions of CSTF2, TPX2, FANCG, and ZNF185 were significantly associated with poor outcomes (log rank test, p = 0.009, p < 0.0001, p = 0.0357, and p = 0.006, respectively) ( Figures 7E and S8).
These data suggest that the expression levels of CSTF2, TPX2, FANGG, and ZNF185 may be negatively correlated with the immune index and warrant validation of their biological functions in vitro to find potential biomarkers indicative of the immune index.
Given that in biological and therapeutic systems upregulated genes are more likely to be manipulated than downregulated genes, we decided to focus our validation study on investigating the functions of CSTF2, TPX2, FANCG, and ZNF185 in PDAC. Firstly, we examined the expression levels of these genes in the human pancreatic ductal cell line PANC-1 as well as the pancreatic carcinoma cell HPDE6-C7. As presented in  Figure 8A, increased expression of ZNF185, CSTF2, and FANCG (p < 0.01), but not TPX2, was observed in the HPDE6-C7 cell line (p > 0.05). These results indicate that ZNF185, CSTF2, and FANCG were overexpressed in PDAC and might be considered as biomarkers of poor prognosis. Next, for each gene (ZNF185, CSTF2, and FANCG), we designed triplicate sets of siRNAs. We downregulated ZNF185, CSTF2, and FANCG in HPDE6-C7 through transfection with three siRNAs using the empty vector as a negative control. As a result, gene expressions significantly dropped after transfection ( Figure S10), where si-ZNF185-3, siCSTF2-2, and si-FANCG-3 exhibited the best transfection efficiency among each group (p < 0.01, p < 0.001, and p < 0.001, respectively). These siRNAs were subsequently selected for functional cell experiments. Weakened cell viabilities were observed in ZNF185, CSTF2, and FANCG downregulated HPDE6-C7 cells compared with NC-transfected  Figure 8B). Additionally, the migration assay showed that downregulation of ZNF185, CSTF2, and FANCG in HPDE6-C7 cells can significantly inhibit cell invasion capability (p < 0.01; Figure 8C). Collectively, these findings indicate that ZNF185, CSTF2, and FANCG are involved in the proliferation and metastasis of PDAC.

DISCUSSION
During the last decade, the emergence of immunotherapy has presented a significant advantage for certain malignancies, including melanoma, lung cancer, and other tumor patients (49)(50)(51)(52). This has resulted in the improvement of the prognosis of treated patients, suggesting that there could be substantial benefits associated with immunotherapy. Even though immunotherapy has achieved impressive clinical benefits, its applications have shown limited efficacy due to multiple immune escape mechanisms. There is an urgent need to better understand the tumor immune microenvironment to design more effective strategies and improve the response and outcomes of immunotherapy. In this study, we reevaluated available OMICS data from an immunological perspective to understand the immune response of PDAC and provide clinical implications for personalized immunotherapy. Pancreatic tumors avoid immune responses through different m e c h a n i s m s . F i r s t , p a n c r e a t i c T M E i s r i c h i n immunosuppressive cells such as MDSCs, regulatory T cells (Tregs), tumor-associated macrophages (TAMs), and immunosuppressive antigen-presenting cells (APCs), which lead to refractory ICB (12,53,54). In our study, unsupervised consensus clustering showed that there is substantial variation in the immune characteristics of PDAC. Among all subtypes in our study, IS3 is associated with the best prognosis, whereas IS1 indicates the poorest survival probability. Cyclic GMP-AMP synthase (cGAS) is a nucleotidyl transferase that is critical for the recognition of double-stranded DNA in the cytosol and generates second messenger cyclic GMP-AMP (cGAMP). It promotes the activation of STING to elicit kinase TBK1 and its downstream substrate, the transcription factor IRF3. The phosphorylated IRF3 then translocates into the nucleus and stimulates the production of type I interferons (IFNs) and other cytokines, leading to immunity regulation (55,56). Type I IFNs are antitumor cytokines and modulators of innate immune activation that have recently been shown to cause radiosensitization in pancreatic cancer, and they all bind to the cell surface receptor complex (IFNAR) (39,43,57). The deactivation of the IFN1-IFNAR1 pathway by cancerassociated fibroblasts (CAFs) leads to tumorigenic effects and tumor growth in pancreatic cancer (41,43). This pathway affects the pancreas by regulating the production of T cells so that when STING is activated in human cells, the T-cell infiltration can be decreased (41). However, our findings have not shown a significant difference in gene expressions of these four key genes within cGAS-STING pathways among the three ISs, suggesting that our classifications have little connection with innate immunity in pancreatic cancer. Furthermore, it is well known that IFN-g has been shown to exhibit antitumor activity both in vitro and in vivo, and deficiencies in IFN-g (IFN-g −/− ) or the IFN-g receptor (IFN-gR −/− ) promote tumor development in mice (58,59). The CYT score has been reported to be associated with improved prognosis in glioma, gastric cancer, and hepatocellular carcinoma, among other cancers (21,(60)(61)(62). In our study, IS3 patients represent the highest score in both IFN-g and CTY. Their maintenance of strong immune activity may suggest a good response to immunotherapy. Since most of the expression levels of inhibitory immune checkpoints were upregulated more in IS3 than in other ISs [i.e., CTLA4, BTLA, IDO1, PDCD1 (PD-1), IDO2, LAG3, etc.], we can conclude that different subgroups have varied therapeutic responses to ICB, which indicates that IS3 may be sensitive to ICB and vice versa.
We also demonstrated stromal and immune phenotypes in IS3, including a high proportion of T cells CD8, T cells CD4 memory resting, and B cells as well as a low proportion of macrophages M0 and macrophages M2 based on CIBERSORT and ESTIMATE algorithms. The presence of M2-type macrophages is associated with poor clinical outcomes in various types of cancers due to certain cellular interactions within metastatic sites. Nevertheless, the low CD8 T cells and the high M2-type macrophages were infiltrated in IS1, with unsurprisingly minimal immune infiltration. Our study revealed that IS3 is more significantly correlated with melanoma patients who respond well to anti-PD1 treatment than IS1 and IS2. However, the high levels of T-cell dysfunction in IS3 might explain why some PDAC patients have high T-cell infiltration but are generally tolerant to ICB. Additionally, IS1 patients were more sensitive to traditional chemotherapeutic agents, including gemcitabine, cisplatin, and erlotinib, but more chemoresistant to 5-FU. These findings might provide a strong rationale in patient selection for systemic therapy optioncombination regimens.
Unlike PCA, which summarizes the total variation in a dataset, LDA derives synthetic variables from a linear combination of features looking for maximum separation of two or more classes of objects. This allows us to derive an immune index to characterize immunity based on LDA scores. Consistent with the highest immune infiltrations, IS3 represented the highest immune index, which showed as significantly positively associated with the immune cells engaged in the antitumor response. High LCK has been identified with a good prognosis in breast cancer, endometrial carcinoma, and melanoma (63)(64)(65). It was in this way that we first identified the potential roles of LCK in PDAC. Strikingly, Tiziana et al. reported that MHC-II was highly correlated with the LCK and HCK metagene (66), which reflects the infiltration of T cells and monocytes (67). Indeed, IS3 shows high immune infiltration and plays critical roles in antitumor immune responses through activating MHC-II, LCK, and HCK as well as augmentation of its expression levels. Pancreatic tumors remain refractory to ICB as it is characterized by its T-cell exclusion and an immunosuppressive tumor microenvironment (68)(69)(70). However, while the LDA score is significantly associated with the expression levels of anti-CTLA4 and anti-PD1, it is not associated with PD-L1. The high immune index represents a good clinically predicted response to both anti-CTLA4 and anti-PD1 therapy. Actually, our results strengthen the findings that patients with a high immune index (IS3) and low Tcell exclusion might benefit from PD-1 therapy; however, the complications resulting from severe T-cell dysfunction needs to be further investigated in the future.
We found several interesting components within the hub gene modules related to IS3 as well overall survival: ZNF185, FANCG, and CSTF2. These candidates have been functionally validated, and their potential roles have also been identified by previous studies. ZNF185, an actin cytoskeleton-associated protein from the LIM family of Zn-finger proteins, has been demonstrated to be a bona fide p53 target gene following DNA damage, which is consistent with the enrichment results of the GO biological process in the tan module-p53 single signaling pathway (71). In addition, the overexpression of ZNF185 has been shown to promote chemoresistance, tumor proliferation, and inhibition of apoptosis by downregulating SMAD4 in PDAC (72). Furthermore, inherited mutations of the DNA repair genes FANCG have been thoroughly demonstrated to be responsible for accelerating genomic instability due to the loss of the G1/S checkpoint, leading to increased risk in PDAC (73,74). Even though the functional role of CSTF2 has not been conclusively proven, the overexpression results in the shortening of 3′ UTRs and promotes pathogenesis and poor prognosis in several types of cancer, including endometrial carcinoma, lung cancer, and bladder cancer (75)(76)(77). The results indicate that ZNF185, FANCG, and CSTF2 cannot be effective biomarkers for risk assessments and reflect the immune index for personalized treatment decision-making in PDAC. Nevertheless, the immune-related roles of these genes will have to be further verified in future studies.
Although computational approaches for analyzing bulk RNA sequencing data from public databases have been well established, this study shares some limitations with previous bioinformatics studies. Given the recalcitrant nature of pancreatic cancer, there are rarely research cohorts regarding immunotherapy for it. The TCGA, ICGC, and GEO cohorts were collected prior to the immunotherapy. The projection of pancreatic cancer by assessing the similarities between molecular features in gene expression data with melanoma immunotherapy cohorts may cause some bias. Furthermore, given that we used these three gene signatures to represent the immune infiltration, we might be unable to comprehensively profile the immune landscape of the PDAC sample. The true predictive effect of these three gene signatures (ZNF185, FANCG, CSTF2) on immune infiltration needs to be evaluated in future studies. Although we predict that the IS3 group may benefit from anti-PD1 therapy and show low levels of T-cell exclusion, whether the true response rate is higher than in other groups requires further validation in other cohorts receiving anti-PD1 treatment.
In summary, we identified three immunological subtypes of PDAC, each with distinct immune infiltration and prognostic characteristics. We further established an immune index to quantify the abundance of immune infiltration in different ISs. Finally, based on the co-expression network analysis, potential biomarkers (ZNF185, FANCG, CSTF2) related to the immune index were identified and functionally validated in vitro. Overall, our results provide new insights into the stratification and selection of patients for personalized immunotherapy assessment in PDAC.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.