Prognostic significance and immune infiltration of microenvironment-related signatures in pancreatic cancer.

Zhong Chen and Qian Lu co-designed the study. Qian Lu, Yu Zhang collected the Data and analyzed Identification and Functional enrichment of DEGs. Qian Lu, Xiaojian Chen, Weihong Gu constructed the Protein-Protein Interaction (PPI) Network. Qian Lu, Xiaojian Chen, Xinrong Ji analyzed Overall Survival and TIMER database. All authors wrote this munscript together. Conceptualization: Qian Lu, Zhong Chen. Data curation: Qian Lu, Yu Zhang, Xiaojian Chen, Weihong Gu, Xinrong Ji. Formal analysis: Qian Lu, Yu Zhang, Xiaojian Chen, Weihong Gu, Xinrong Ji. Funding acquisition: Zhong Chen. Investigation: Qian Lu, Yu Zhang, Weihong Gu, Xinrong Ji. Methodology: Qian Lu, Yu Zhang, Weihong Gu, Xinrong Ji. Project administration: Qian Lu. Resources: Qian Lu. Software: Qian Lu. Writing – original draft: Qian Lu, Yu Zhang, Weihong Gu, Xinrong Ji. Writing – review & editing: Zhong Chen.

Among the 140 PDAC patients enrolled in the present study, 65 were female and 75 were male. Calculated by the ESTIMATE algorithm, the immune scores ranged from −1559.87 to 2648.04 while the stromal scores from −1843.32 to 2138.93, and ESTIMATE scores ranged from -3187.36 to 4161.76. The distributions of immune, stromal, and ESTIMATE scores did not associate with T stage (Fig. 1 A–1C), N stage (Fig. 1 D–1F), or M stage (Fig. 1 G–1I). Immune, stromal, and ESTIMATE scores were not correlated with T stage, N stage, or M stage. Distribution of immune scores plotted against T stage (A), N stage (D), and M stage (G). Distribution of stromal scores plotted against T stage (B), N stage (E), and M stage (H). Distribution of ESTIMATE scores plotted against T stage (C), N stage (F), and M stage (I). ESTIMATE = estimation of stromal and immune cells in malignant tumor tissues using expression data. According to immune, stromal, and ESTIMATE scores, 140 PDAC patients were divided into low and high score groups to detect potential relationship between immune/stromal/Estimate score and prognosis. For immune scores, the low score group contained 57 patients, while the high score group contained 83 patients. Kaplan- Meier survival curves showed that elevated immune scores associated with better overall survival ( P  = .033) (Fig. 2 A). By contrast, stromal scores ( P  = .656) and ESTIMATE scores ( P  = .496) were not correlated with overall survival (Fig. 2 B, 2C). Association of immune, stromal, and ESTIMATE scores with overall survival. (A) Elevated immune scores correlated with a better prognosis. (B) Stromal scores were not associated with overall survival. (C) ESTIMATE scores were not associated with overall survival. ESTIMATE = estimation of stromal and immune cells in malignant tumor tissues using expression data. We examined transcriptional microarray data of 140 PDAC cases to identify DEGs with immune and/or stromal scores. Figure 3 A and 3B show Volcano plot of gene expression profile in high vs low immune scores groups and stromal scores. Figure 3 C shows a heat map of 324 DEGs with high or low immune scores, and Figure 3 D shows a heat map of 318 DEGs with high or low stromal scores. According to immune score comparison, 320 genes were up-regulated and 4 genes down-regulated in the high score than the low score group. For high stromal score group compared with low score group, 297 up-regulated genes and 21 down-regulated genes were detected. A total of 86 DEGs were synchronously upregulated (Fig. 3 E) in the high stromal score and high immune score groups, and 4 genes were commonly downregulated (Fig. 3 F) using Venn algorithm. Comparison of gene expression profile with immune scores and stromal scores in PDAC. (A) Volcano maps show the distribution of DEGs based on immune scores. (B) Volcano maps show the distribution of DEGs based on stromal scores. Red and green dots represent high and low expressed genes, respectively. (C) Heatmaps of DEGs of immune scores. (D) Heatmaps of DEGs of stromal scores. (E) Venn diagrams showing the number of commonly upregulated DEGs in stromal and immune. (F) Venn diagrams showing the number of commonly downregulated DEGs in stromal and immune. DEGs = differential expressed genes. GO and KEGG enrichment analyses were performed to explore the biological classification of DEGs via DAVID online tool. GO analysis results indicated that biological processes of DEGs were mainly enriched in lymphocyte proliferation, mononuclear cell proliferation, regulation of lymphocyte proliferation (Table 1 ). Changes in cell component (CC) of DEGs were significantly enriched in immunological synapse, T cell receptor complex, ficolin-1-rich granule membrane (Table 1 ). Changes in molecular function of DEGs were mainly enriched in the CD4 receptor binding, MHC protein complex binding and sialic acid binding (Table 1 ). KEGG pathway analysis showed that the DEGs were ignificantly enriched in PPAR signaling pathway, Hematopoietic cell lineage and IL-17 signaling pathway (Table 1 ). GO and KEGG pathway enrichment analysis of DEGs with immune and stromal score in PDAC. GO = Gene Ontology, KEGG=Kyoto Encyclopedia of Genes and Genomes, PDAC=Pancreatic ductal adenocarcinoma We used Cytoscape to construct the PPI network of intersection DEGs (Fig. 4 A) and obtain the most significant module (Fig. 4 B-C). This network composed of 4 modules, the top 2 significant modules were selected for further analysis. These modules were named as Modules B, and C, respectively. In Module B, 13 nodes with 35 edges were formed in the network, including ITGAM, GZMA, HCST, CCL17, CD19, CD74, CCL4, TCL1A, CD83, TNFRSF1B, CD22, CD69, MS4A1 (Fig. 4 B). Module C contained 13 edges involving 7 nodes: SPN, ZAP70, IL2RB, CD1C, KLRB1, ITGAX, BTK (Fig. 4 C). PPI network and the most significant module of DEGs. (A) The PPI network of DEGs was constructed using Cytoscape. (B, C) The top two significant module was obtained from PPI network using MCODE. Upregulated genes are marked in light red; downregulated genes are marked in light blue. DEGs = differential expressed genes, MCODE = molecular complex detection, PPI = protein–protein interaction. To explore the relationship between commonly DEGs and overall survival, we performed Kaplan-Meier survival curves by using the TCGA database. The Kaplan-Meier results indicated that decreased CA9, MUC5AC, P2RY8 mRNA expression significantly associated with poor OS ( P  < .05), and elevated CD36, EBI3, PLA2G2D, SPOCK2, WDFY4, CD1D, CCL22 mRNA expression was significantly correlated with shorter OS for PDAC patients ( P  < 0.05) (Fig. 5 ). Then, we validate the prognostic value of survival related DEGs using Kaplan-Meier plotter online analysis tool. The results show CA9, EBI3, SPOCK2, WDFY4, CD1D, CCL22 mRNA expression was significantly correlated with OS ( P  < .05) (Fig. 6 ). Correlation of expression of individual DEGs with overall survival in TCGA. KM survival curves were plotted for selecting DEGs with prognostic value by comparison of groups of high (red line) and low (blue line) gene expression. DEGs = differential expressed genes. Validation of correlation of DEGs extracted from the TCGA database with overall survival by Kaplan-Meier plotter online analysis tool. KM survival curves were plotted for selecting DEGs with prognostic value by comparison of groups of high (red line) and low (black line) gene expression. DEGs = differential expressed genes. As shown in Table 2 , EBI3, P2RY8, CA9, CD36, SPOCK2, CCL22, CD1D expression (ref. low), T stage (ref. T1-T4) and N stage (ref. N0-1) were demonstrated as independent prognostic indicators for PDAC patients ( P  < .05). Multivariate Cox analysis revealed that poor OS was significantly associated with CA9 expression (ref. low; HR=2.751, p = 0.008), SPOCK2 expression (ref. low; HR = 3.310, p = 0.016), CD1D expression (ref. low; HR = 0.212, P  = .034) and N stage (ref. N0-1; HR = 2.849, P  = .029). Univariate and multivariate Cox logistic regression analysis of overall survival in TCGA cohort. CA9 = Carbonic anhydrase 9, CD1D = Cluster of Differentiation 1D, SPOCK2 = SPARC (osteonectin), cwcv and kazal-like domains proteoglycan 2. After validating prognostic value of CA9, SPOCK2, and CD1D, we performed correlation analysis between CA9, SPOCK2 and CD1D and immune infiltration level for PDAC in Figure 7 . Scatter plots were generated with partial Spearman's correlation and statistical significance. CD1D expression were significantly associated purity (correlation=-0.283). In addition, SPOCK2 and CD1D significantly associated with B cell, CD8+ T cell, CD4+ T cell, macrophage, neutrophil, and dendritic cell infiltration, CA9 significantly associated with CD8+ T cell, macrophage, neutrophil, and dendritic cell infiltration ( P  < .05). To further explore the association between genomic metrics of CA9, SPOCK2 and CD1D and the immune cell enrichment, PDAC with SCNAs were divided into four levels, including arm-level deletion, diploid/normal, arm-level gain, and high amplification. We examined infiltration of six immune cell types, including B cell, CD8+T cell, CD4+T cell, macrophage cell, neutrophil cell, and dendritic cell. The results showed that the extent of immune infiltration was largely different in PDAC with different SCNAs of CA9, SPOCK2, and CD1D (Fig. 8 ). Enrichment of B cell, CD4+T cell displayed significant downregulation in PDAC with particular SCNA of CA9. In addition, decreasing of B cell, CD8+T cell, CD4+T cell, macrophage cell and neutrophil cell enrichments was also detected in PDAC with particular SCNA of SPOCK2 and CD1D. Therefore, the results revealed that genomic alterations of CA9, SPOCK2 and CD1D were associated with the extent of immune infiltration in PDAC. Immune infiltration of CA9, SPOCK2, and CD1D. After identifying prognostic value of CA9, CD1D, and SPOCK2, we performed correlation analysis between CA9, CD1D, and SPOCK2 and immune infiltration level for PDAC. Scatter plots were generated with partial Spearman's correlation and statistical significance. (A) CA9 expression has a significant positive correlation with infiltrating levels of CD8+ T cell, macrophage, neotrophil, and dendritic cells in PDAC. (B) SPOCK2 expression has a significant positive correlation with infiltrating levels of B cell, CD8+ T cell, CD4+, T cell, macrophage, neotrophil, and dendritic cells in PDAC. (C) CD1D expression has a significant positive correlation with infiltrating levels of purity, B cell, CD8+ T cell, CD4+ T cell, macrophage, neotrophil, and dendritic cells in PDAC. The association of SCNAs of CA9, SPOCK2, and CD1D with immune infiltration in PDAC. In the genomic datasets, SCNAs of CA9, SPOCK2 and CD1D were divided into into 4 levels, including arm-level deletion, diploid/normal, arm-level gain, and high amplification. The infiltration of 6 immune cell types, including B cell, CD 8+T cell, CD4+T cell, macrophage cell, neutrophil cell and dendritic cell, was analyzed in PDAC. (A) Enrichment of B cell, CD4+T cell displayed significant downregulation in PDAC with particular SCNA of CA9. (B) Enrichment of B cell, CD 8+T cell, CD4+T cell, macrophage cell and neutrophil cell displayed significant downregulation in PDAC with particular SCNA of SPOCK2. (C) Enrichment of B cell, CD 8+T cell, CD4+T cell, macrophage cell and neutrophil cell displayed significant downregulation in PDAC with particular SCNA of CD1D.

Immune scores and stromal scores of PDAC patients were obtained from ESTIMATE ( http://bioinformatics. mdanderson.org/estimate/ ) designed by Yoshihara et al. [ 16 ] Gene expression profile (level 3 data) and clinical information for PDAC patients was collected from the TCGA data portal ( https://tcga-data.nci.nih.gov/tcga/ ). Only samples in accordance with the following Inclusion criteria in the present study: available immune and stromal scores; available gene expression information; the OS days ≥30 days. X-tile software was used to identify cutoff criteria for the ESTIMATE scores, immune score and stromal score. According to the cutoff point, we stratify PDAC patients into low-score and high-score sub-groups. Data analysis was performed by use of package limma [ 17 ] in R software (version 3.6.0). The threshold of DEGs were: |log2fold change (log2FC)| >1.5 and false discovery rate <0.05. DAVID ( http://david.ncifcrf.gov , Version 6.8) online analysis tool was undertaken to assess the role of intersection genes in PDAC. P -value <.05 was considered to be statistically significant. The Gene Ontology (GO) analysis shows the DEG function at the level of 3 main categories (molecular function, biological process and cellular component), and the KEGG analysis reveals the pathway enrichment of DEGs. The PPI network of DEGs was constructed using retrieved through the Search Tool for the Retrieval of Interacting Genes ( https://string-db.org/ ) [ 18 ] and Cytoscape software (version 3.5) [ 19 ] were used for the reconstruction of the PPI network. The Molecular Complex Detection (version1.4.2) plugin of Cytoscape was applied to identify the most significant module based upon topology to locate densely connected regions. The selection criteria were as follows: maximum depth = 100, degree cut-off = 2, node score cut-off = 0.2, and k-core = 2. The relationship between the OS of PDAC patients and the expression of DEGs was illustrated by Kaplan-Meier plots using R. Kaplan-Meier plotter ( http://kmplot. com/analysis/ ) online analysis tool was applied to verify the prognostic value of survival related DEGs. The relationship was determined by log-rank test. P-value<0.05 was considered to be statistically significant. We analyzed comprehensive correlation between prognostic genes expression and tumor-infiltrating immune cells gene markers via TIMER ( https://cistrome. shinyapps.io/timer/ ). Moreover, the association between somatic copy number alterations (SCNAs) of prognostic genes and immune cell enrichments in PDAC was also analyzed by TIMER.

PDAC is 1 of the highly fatal and most aggressive types of malignancies and approximately accounts for 90% of all cases of PC. [ 2 ] Numerous studies have reported that the tumor microenvironment (TME) was significantly correlated with the oncogenesis, progress, and prognosis of various malignancies. [ 20 – 22 ] The Interplay between tumor cells and the Immune System in TME reflect the plasticity of immune system and the tumor. [ 23 ] Tumor development is closely related to TME, and any changes of the component of TME may affect the progression of malignancies. [ 23 ] Exploring the alterations may help the advance of therapeutic strategies. Immune cells and stromal cells are the important compositions of TME, which play a key role in the progress of cancers. [ 24 ] The ESTIMATE algorithm based on gene expression profile has been applied for the exploration of breast cancer, prostate cancer, colorectal cancer and shown good clinical application value. [ 12 – 14 ] However, it has not been used for the investigation of PDAC. In the present work, we used ESTIMATE algorithm to calculate immune scores, which indicate the intensity of immune cells infiltration in tumor. [ 16 ] Our result indicated that high immune score was related with improved OS. Recent researches have similar findings. Heeren et al [ 25 ] found that the number of CD4(+) T-cell is significantly higher in lymph node-negative samples than in lymph node- positive ones. Moreover, Punt et al [ 26 ] reported that TCL1A expression of B cell associated with better survival in Cervical cancer. These results may suggest that immune cells infiltrating tumor tissue affect tumor progression. In order to explore the potential TME related gene in PDAC, we performed bioinformation analysis to identify intersection DEGs in both stromal and immune cells. Finally, 86 upregulated commonly genes and 4 downregulated commonly genes were obtained between the stromal and immune score groups. These 90 intersection genes were correlated with biological processes in the TME, including lymphocyte proliferation, mononuclear cell proliferation, regulation of lymphocyte proliferation. These processes may play a significant role in tumor development and metastasis and significantly influence OS. [ 27 ] Among molecular functions, intersection genes enriched in the CD4 receptor binding, MHC protein complex binding and sialic acid binding. CD4+ T cells constitute an important component of tumor infiltrating lymphocytes and exert a key role in tumor immune-surveillance. [ 28 ] Moreover, it has been demonstrated that interactions with sialic acid-binding receptors can significantly affect cancer progression. [ 29 ] Next, we performed Kaplan-Meier survival curves to detect prognostic immune associated genes in PDAC and validate by Kaplan-Meier plotter online analysis tool. The results indicated that CA9, EBI3, SPOCK2, WDFY4, CD1D and CCL22 were significantly correlated with OS in patients with PDAC. Of the 6 genes identified, CA9, [ 30 ] EBI3 [ 31 ] have been proved to be involved in PDAC progression or significant in predicting OS, suggesting that our bioinformation analysis based on TCGA has prognostic values. The remaining 4 genes include SPOCK2, WDFY4, CD1D, CCL22 have not or rarely been reported in PDAC prognosis previously, and may become potential biomarkers for PDAC. Subsequently, the expression of CA9, EBI3, SPOCK2, WDFY4, CD1D and CCL22 were enrolled in multivariate analysis for OS in PDAC. Importantly, CA9 and SPOCK2 and CD1D were demonstrated as independent prognostic indicators for PDAC patients. CA9, Carbonic anhydrase 9, is a cancer-related cell surface enzyme catalyzing the reversible conversion of carbon dioxide to bicarbonate ion and proton. [ 32 ] Cells will accelerate metabolism and produce excess acid in a hypoxic microenvironment, which activates the catalytic activity of CA9 to regulate the intracellular and extracellular pH perturbations. [ 33 ] Moreover, CA9 expressed in many tumor types also plays an important role in cell adhesion and spreading. [ 34 ] Therefore, in hypoxia or acidosis environments, CA9 can raise the survival advantage of tumor cells and enhance their ability of migration, invasion and metastasis. [ 35 ] A recent research showed that in response to hypoxia, PDAC cells expressing activated KRAS increase expression of CA9, via stabilization of HIF1A and HIF2A, to regulate pH and glycolysis. [ 30 ] Combined with our analysis, it highlights the possibility of CA9 as a targeted therapy for PDAC. The full name of SPOCK2 is SPARC (osteonectin), cwcv and kazal-like domains proteoglycan 2, which is a glycoprotein composed of a protein domain, 2chondroitin side chains, and heparan sulfate. [ 36 ] A recent study showed that hypermethylation of the SPOCK2 in colon cancer is significantly higher than in normal mucosal tissues adjacent to the cancer. [ 37 ] In addition, Ren F et al revealed that hypermethylation of SPOCK2 is associated with endometriosis-related ovarian endometrioid adenocarcinoma [ 38 ] and that SPOCK2 regulates biological behavior of endometrial cancer cells through regulating the protein expression of MT1-MMP and MMP2 and the activation of MMP2. [ 39 ] These evidences indicated that the abnormality of SPOCK2 may contribute to tumorigenesis and the progression of tumors. However, the role of the SPOCK2 has not been researched in PDAC. Our analysis indicates that SPOCK2 may be an important gene in the development of PDAC, which should be confirmed furtherly. CD1D, a member of the CD1 family, is a class of non-polypeptide transmembrane glycoprotein molecules, which is widely expressed in thymocyte B cells, epidermal Langham cells, dendritic cells, and intestinal epithelial cells. [ 40 ] With function of antigen presentation similar to MHC I, CD1D can specifically present antigens to NKT cells, make them activate and secrete a variety of cytokines, and directly or indirectly participate in the body's immune response. [ 41 ] Recently, a research about CD1d-binding glycolipid for iNKT-cell-based therapy showed a potent and effective treatment against human breast cancer. [ 42 ] Another clinical phase II study showed that intravenous administration of α-GalCer-pulsed antigen-presenting cells (APCs) prolonged overall survival for patients with advanced or recurrent non-small cell lung cancer. [ 43 ] (α-GalCer, an exogenous glycolipid antigen extracted from sponges, specifically binding to CD1D molecule can activate NKT cells to exert immune effects. [ 44 ] ) The effect of CD1D-restricted invariant natural killer T cells in PDAC immunotherapy are unclear. However, our analysis suggests that CD1D may have an important role in PDAC. Therefore, we considered CD1D as a promising target for PDAC immunotherapy. In summary, our findings provide a more comprehensive insight into TME and show a list of prognostic immune associated genes in PDAC. CA9, SPOCK2 and CD1D exhibited significant prognostic potential, and CD1D were significantly related with immune cell biomarkers for PDAC patients, thus indicating the relevance of monitoring and regulating the TME for PDAC prognosis and precision immunotherapies. However, further studies on these genes need to be performed to gain additional understanding of the association between TME and prognosis in PDAC.

Pancreatic cancer (PC) is 1 of the highly fatal and most aggressive types of human malignant tumors worldwide. In 2018, PC was newly diagnosed in ∼458 million people and approximately 432 million people died of PC in the whole world, ranking seventh among all cancer-associated mortalities. [ 1 ] Pancreatic ductal adenocarcinoma (PDAC) is an epithelial malignancy which originates from the ductal epithelium of pancreas and accounts for 90% of all cases of PC. [ 2 ] Due to a lack of distinct symptoms, PDAC has often already extensively metastasized at the time of diagnosis, and radical resection is almost impossible. [ 3 , 4 ] In spite of extensive efforts over the past few years, the conventional therapies for advanced PDAC, including chemotherapy, radiotherapy and targeted therapies, continues to have the poorest outcome among all the malignancies and the median overall survival (OS) remains approximately six months. [ 5 – 7 ] Thus, it becomes necessary to identify early diagnostic biomarkers and therapeutic targets of PDAC. The tumor microenvironment (TME) is a complex and synergistic mixture, comprised of fluids, immune cells, stromal cells, extracellular matrix molecules, inflammatory mediators and chemokines. The various constituents of the TME play significant roles in the tumor immune escape, tumor proliferation and metastasis, and it have been shown to be related with diagnosis and clinical outcomes. [ 8 , 9 ] Numerous clinical trials of targeting TME showed encouraging results. [ 10 , 11 ] Thus, exploring the molecular function and composition of TME is significant to effectively intervene tumor progression and immune response. Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) was created to evaluate tumor purity based on the molecular biomarker expression of stromal and immune cells. The ESTIMATE algorithm has been used to breast cancer, [ 12 ] prostate cancer, [ 13 ] colorectal cancer, [ 14 ] and other studies about tumor microenvironment. [ 15 ] However, prognostic value of ESTIMATE algorithm in PDAC is still unclear. In the present study, we applied the ESTIMATE algorithm to identify prognostic biomarkers and explore the immune infiltration of prognostic genes in PDAC.