- 1Department of General Surgery, Clinical Laboratory of Integrative Medicine, The First Affiliated Hospital of Dalian Medical University, Dalian, Liaoning, China
- 2Department of Visceral, Vascular and Endocrine Surgery, Martin-Luther-University Halle-Wittenberg, University Medical Center Halle, Halle, Germany
The worldwide prevalence of pancreatic cancer has been rising in recent decades, and its prognosis has not improved much. The imbalance of substance and energy metabolism in tumour cells is among the primary causes of tumour formation and occurrence, which is often controlled by the neuroendocrine system. We applied Cox and LASSO regression analysis to develop a neuroendocrine regulation- and metabolism-related prognostic risk score model with three genes (GSK3B, IL18 and VEGFA) for pancreatic cancer. TCGA dataset served as the training and internal validation sets, and GSE28735, GSE62452 and GSE57495 were designated as external validation sets. Patients classified as the low-risk population (category, group) exhibited considerably improved survival duration in contrast with those classified as the high-risk population, as determined by the Kaplan-Meier curve. Then, we combined all the samples, and divided them into three clusters using unsupervised clustering analysis. Unsupervised clustering, t-distributed stochastic neighbor embedding (t-SNE), and principal component analysis (PCA) were further utilized to demonstrate the reliability of the prognostic model. Moreover, the risk score was shown to independently function as a predictor of pancreatic cancer in both univariate and multivariate Cox regression analyses. The results of gene set enrichment analysis (GSEA) illustrated that the low-risk population was predominantly enriched in immune-associated pathways. “ESTIMATE” algorithm, single-sample GSEA (ssGSEA) and the Tumor Immune Estimation Resource (TIMER) database showed immune infiltration ratings were enhanced in the low-risk category in contrast with the high-risk group. Tumour immune dysfunction and exclusion (TIDE) database predicted that immunotherapy for pancreatic cancer may be more successful in the high-risk than in the low-risk population. Mutation analysis illustrated a positive link between the tumour mutation burden and risk score. Drug sensitivity analysis identified 44 sensitive drugs in the high- and low-risk population. GSK3B expression was negatively correlated with Oxaliplatin, and IL18 expression was negatively correlated with Paclitaxel. Lastly, we analyzed and verified gene expression at RNA and protein levels based on GENPIA platform, HPA database and quantitative real-time PCR. In short, we developed a neuroendocrine regulation- and metabolism-associated prognostic model for pancreatic cancer that takes into account the immunological microenvironment and drug sensitivity.
Introduction
Pancreatic cancer is becoming an increasingly major health problem throughout the world, with the age-standardized incidence rates rising from 5.0/100,000 persons in 1990 to 5.7/100,000 in 2017 (1). Pancreatic cancer is projected to overtake lung cancer as the major cancer killer in the United States by 2030, moving up from its current position as the third leading contributor to cancer-associated deaths (2–4). Surgery is the only current hope for curing pancreatic cancer, which is exceedingly aggressive and has a very dismal prognosis. However, about 80-85% of patients with pancreatic cancer are diagnosed as unresectable or metastatic due to the occult onset and lack of early screening methods. Even when diagnosed to be resectable, pancreatic cancer prognosis is dismal, with just a 20% 5-year survival probability after surgery (5). Surgery and chemotherapy are common treatment options for pancreatic cancer, however, they have limited effectiveness. Researchers discovered that the median survival periods for individuals with metastatic pancreatic cancer treated with FOLFIRINOX (leucovorin, irinotecan, oxaliplatin, and fluorouracil) or gemcitabine were 11.1 months and 6.8 months, respectively, while those with resected pancreatic cancer treated with modified FOLFIRINOX or gemcitabine had median survivals of 54.4 months and 35.0 months, respectively (6, 7). Although the survival benefit of FOLFIRINOX is superior to gemcitabine, there are also more toxic side effects. Targeted therapy and immunotherapy have made enormous strides in recent years, bringing revolutionary advances to the treatment of cancer, but their effectiveness is not optimal for pancreatic cancer. A study by Hong et al. (8) illustrated that the overall response rate (complete or partial response) of ibrutinib plus durvalumab (a PD-L1-targeting antibody) in the treatment of pancreatic cancer was only 2%. Therefore, it is still extremely vital and urgent to investigate the mechanism of occurrence and development and to identify an effective therapy for pancreatic cancer.
Metabolic alteration is an important topic in cancer biology research, and metabolic reprogramming is considered to be one of the hallmarks of cancer, participating in the process of cancer occurrence, development and metastasis (9). To adapt to the microenvironment of hypoxia and nutrient deficiency, establish survival advantages and achieve rapid growth, tumour cells change their material and energy metabolism patterns, which is called metabolic reprogramming (10). The research found that tumour cells significantly increase the demand and uptake of glucose, rapidly produce ATP through the glycolysis pathway, and aerobic glycolysis is performed even in the presence of oxygen, also referred to as the “Warburg effect”. Lactic acid produced by glycolysis will accumulate in the tumour microenvironment (TME), boost tumour cells invasiveness and reduce the anti-tumour immunity (11). Glutamine is an important source of nitrogen and carbon in biosynthetic reactions. Normal cells can synthesize glutamine by themselves, while tumour cells can obtain glutamine from the microenvironment by solute carrier group in addition to their own synthesis to meet the proliferation needs (12). Lipids, mainly including fatty acids, cholesterol, phospholipids and acrylamide, are not only the basic structure of cell membranes, but also a source of signalling molecules and energy, and the lipid metabolism reprogramming of tumour cells can promote their proliferation, invasion and metastasis (13, 14). Ringel et al. (15) showed that obesity can cause the metabolic changes for fatty acid, impair the function and infiltration for CD8+T cells, and thus inhibit anti-tumour immunity. Hypoxia and nutrient deprivation in TME will lead to metabolic competition between immune and tumour cells, and the high metabolism and strong adaptability of tumour cells will further change the metabolic characteristics of the TME, causing metabolic pressure on immune cells, while continuously accumulating toxic metabolites, negatively affecting the immunity, thereby promoting immune suppression and escape (16, 17). Besides, metabolic alteration in tumour cells is also intimately linked to the sensitivity of chemotherapy, targeted therapy, immunotherapy, and radiotherapy (18–22). Obesity has been shown to increase tumour cells’ resistance to chemotherapy, radiation, and biological and endocrine-targeted treatments by altering the way fatty acids are metabolized in TME (23).
The imbalance of material and energy metabolism of tumour cells is often regulated by the neuroendocrine system. Therefore, tumorigenesis and progression are aided by neuroendocrine control. Cancers of the liver, pancreas, colorectal, breast, and uterus are all thought to be linked to neuroendocrine regulation disorders such as diabetes, obesity, and depression (24–28). An analysis of data from a large-scale cohort study in the US involving 112,818 women and 46,207 men found that those with new-onset diabetes exhibited a 2.97-fold (95% CI, 2.31-3.82) greater risk of developing pancreatic cancer than those without diabetes, while those with long-term diabetes recorded a 2.16-fold (95% CI, 1.78-2.60) higher risk (29). As a digestive organ, the pancreas has both endocrine and exocrine functions, and its lesions are often accompanied by abnormal regulation of blood glucose. Therefore, diabetes is not only a risk factor for pancreatic cancer but also one of its secondary diseases. The mechanism for diabetes causing pancreatic cancer is complex, including hyperglycemia, hyperinsulinemia, insulin resistance, chronic inflammation, and so on (30). Persons with a body mass index (BMI) of 25 to 29.9 were found to have an odds ratio (OR) for pancreatic cancer that was 1.19 (95% CI, 1.02-1.40), whereas those with a BMI of 30 to 34.9 were found to have an OR of 1.62 (95% CI, 1.19-2.21) compared with normal-weight individuals (31). A study by Rebours et al. (32) found that pancreatic fatty infiltration was related to the formation of pancreatic intraepithelial neoplasia.
Based on the above evidence, our study was the first to develop a neuroendocrine regulation- and metabolism-related prognostic model of pancreatic cancer utilizing The Cancer Genome Atlas (TCGA) database, then verified the model through Gene Expression Omnibus (GEO) database. Unsupervised clustering, t-distributed stochastic neighbor embedding (t-SNE), and principal component analysis (PCA) were further conducted to demonstrate the consistency and reliability for the prognostic model. In addition, the association between prognostic models and clinicopathological characteristics, the tumour immune microenvironment, the tumour mutation load, and treatment sensitivity was investigated. Finally, we explored and verified gene expression at RNA and protein levels based on GENPIA platform, the Human Protein Atlas (HPA) database and quantitative real-time PCR, and explored gene expression distribution in different subcellular structures and cell types via HPA database.
Materials and methods
Data acquisition and processing
The transcriptome profiles (including 178 pancreatic tumour tissues and 4 normal pancreatic tissues), somatic mutation status, copy number variation (CNV), and matching clinical data (including 185 pancreatic cancer samples) were extracted in the TCGA database (https://portal.gdc.cancer.gov/). To get the validation set, GSE28735 (including 45 pancreatic cancer samples), GSE62452 (including 69 pancreatic cancer samples) and GSE57495 (including 63 pancreatic cancer samples) were derived in the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The “SVA” R package was utilized to eliminate any batch effects that existed across various data sets (33). To elucidate the differences between pancreatic tumors and normal samples more precisely, the TCGA dataset (containing 178 pancreatic tumour tissues and 4 normal tissues) and the Genotype-Tissue Expression Project (GTEx) project dataset (containing 167 normal pancreatic tissues) were extracted in the UCSC Xena database (https://xenabrowser.ucsc.edu/datapages/) and the gene expression data were normalized using the log2(x+1) transformation.
Identification of differentially-expressed neuroendocrine regulation- and metabolism-related genes
Based on the “limma” R packages, with the filter condition: |log2FC| > 1, and adjusted p-values < 0.05, we detected differentially expressed genes (DEGs) between tumour and normal samples. The screened differential genes were visualized by the “ggplot2” R package. Neuroendocrine regulation-related and metabolism-related genes were retrieved in the GeneCards database (https://www.genecards.org/) (33), setting filter parameters: relevance score > 5. We identified differentially expressed neuroendocrine regulation- and metabolism-related genes (NMRGs) by merging DEGs and NMRGs.
Development and analysis of protein-protein interaction network
We explored the possible interactions for differentially expressed NMRGs using the STRING database (https://cn.string-db.org/) at a minimum interaction score of 0.4 (medium confidence) (34). The PPI network was created and shown utilizing Cytoscape (version 3.9.1). The PPI network’s critical modules and hub genes were then isolated utilizing the MCODE plugin, with screening conditions established as below: Max depth = 100, k-score = 2, node score cutoff = 0.2, and degree cutoff = 2.
Genetic alteration and enrichment analysis
Mutation landscape of hub genes in pancreatic cancer was analyzed and visualized utilizing the “maftools” R package. The “RCircos” R programme was applied to evaluate and illustrate the CNV of hub genes. “org.Hs.eg.db”, “clusterProfiler”, “enrichplot” and “ggplot2” R packages were applied to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses for the above differentially expressed NMRGs, based on the filter condition: q value < 0.05. GO enrichment analysis involves molecular function (MF), cellular component (CC), and biological process (BP).
Construction and verification of the prognostic model
Patients from TCGA were randomized at a ratio of 7:3 into the training and internal validation sets by utilizing the “caret” R package. The external validation set included patients from GSE28735, GSE62452, and GSE57495. Firstly, we conducted univariate Cox regression analysis to preliminarily filter the genes linked to prognosis. Next, the least absolute shrinkage and selection operator (LASSO) regression analysis was utilized to overcome overfitting with the “glmnet” R package. Lastly, a predictive model was developed utilizing multivariate Cox regression analysis. Below is the risk score equation for the prognostic model:
Patients were assigned a risk score using the approach and then categorized into high- and low-risk groups as per how their score compares to the training set’s median risk score. To contrast the prognosis between the high- and low-risk groups, survival curves were generated utilizing Kaplan-Meier method. The prognostic model was evaluated utilizing a time-dependent Receiver Operating Characteristic (ROC) curve and area under the curve (AUC).
Cluster analysis
For identifying the prognostic model’s consistency and dependability, we pooled patients from TCGA, GSE28735, GSE62452, and GSE57495 and performed an unsupervised cluster analysis using the “ConsensusClusterPlus” R package. The cumulative distribution function (CDF) curve and consensus matrix were applied to get the best possible clustering parameter values. The Kaplan-Meier technique was utilized to compare survival times across the categories. Additional verification of the clusters and prognostic model was achieved via the use of PCA and t-SNE.
Clinicopathologic correlation, independent prognostic analysis, and nomogram model construction
We integrated patients’ clinical data and risk scores categorized them as per the clinicopathological criteria, followed by performing the Wilcoxon and Kruskal-Wallis rank sum test to assess risk scores across groups. For identifying whether the risk score independently acted as a predictor of pancreatic cancer, univariate and multivariate Cox regression analyses were conducted. Clinicopathological parameters (age, stage, pathological grade, and gender) and risk scores were considered in the development of the nomogram model for predicting the prognosis of pancreatic cancer utilizing the “rms” R package. The predictive performance of the nomogram model was tested using the calibration curve (35).
Gene set enrichment analyses
To compare the variations in biological processes and metabolic pathways between low- and high-risk groups, gene set enrichment analyses (GSEA) based on gene sets “c5.go.v7.5.1.symbols.gmt” and “c2.cp.kegg.v7.5.1.symbols.gmt” was executed utilizing R package “limma”, “org.hs.eg.db”, “ClusterProfiler” and “enrichplot”, the filter criteria was set as follows: |normalized enrichment score (NES)| > 1 and q value < 0.05.
Immune analysis
To evaluate the variations in TME between low- and high-risk groups, the “ESTIMATE” algorithm was utilized to compute each patient’s stromal, immune, and ESTIMATE scores (36). Single-sample GSEA (ssGSEA) was carried out to determine the infiltration scores of 16 distinct immune cells and the activity score of 13 immune-associated pathways for each patient utilizing the “GSVA” and “GSEABase” R packages. Immunocyte infiltration scores of all TCGA tumors were retrieved from the Tumour Immune Estimation Resource database (TIMER, http://timer.cistrome.org/) to further evaluate their association with risk scores. These scores were calculated using a variety of algorithms, notably, MCPCOUNTER, EPIC, QUANTISEQ, XCELL, CIBERSORTABS, TIMER, and CIBERSORT (37). Tumour immune dysfunction and exclusion (TIDE) can measure the sensitivity of immune checkpoint blockade by simulating tumour immune evasion mechanism (38). The online TIDE platform (http://tide.dfci.harvard.edu/) was utilized to determine the TIDE score, as well as T-cell dysfunction and exclusion scores for each patient.
Mutation analysis
To compare and contrast the mutation profiles of individuals at high and low risk, we utilized the “maftools” in R package. Somatic, insertion, base substitution, coding, and deletion mutations were all included in the definition of tumour mutation burden (TMB) (39). TMB estimates are calculated as overall mutation frequency divided by 38 Mb since this value is commonly retrieved based on the length of human exons (40). The association between risk scores and TMB was analyzed with the use of the Spearman correlation test. The ideal cut-off value of TMB was also determined, and the affect of TMB on the prognosis of pancreatic cancer was assessed, with the use of the “survminer” and “survival” R packages.
Drug sensitivity analysis
To analyze the disparity in medication responsiveness between high- and low-risk groups, we used the “pRRophetic” package in R (41). Using the CellMiner database (https://discover.nci.nih.gov/cellminer/), we retrieved pertinent gene expression data and Food and Drug Administration (FDA) authorized drug sensitivity data to further assess the link between drug sensitivity and the genes in the predictive model (42). The connection between gene expression and drug responsiveness was analyzed using the Pearson test.
Gene expression and distribution exploration
Differences in the expression of GSK3B, IL18 and VEGFA at the RNA level between normal and tumor pancreatic tissues were analyzed via the GEPIA platform (http://gepia2.cancer-pku.cn/#analysis). Protein expression patterns in normal cells, tissues, as well as cancer tissues may be generated using the public HPA platform (https://www.proteinatlas.org/) (43). For protein-level confirmation of gene expression, we used immunohistochemistry pictures of normal and pancreatic cancer tissue obtained from the HPA database. We also obtained relevant data and photos to examine the distribution of gene expression in diverse subcellular structures and cell types. In contrast to traditional bulk RNA-seq that generates mixed gene expression data from tissues, single cell RNA-seq can provide a transcriptional data of individual cells (44). Tumor Immune Single-cell Hub (TISCH, http://tisch.comp-genomics.org) is a single cell RNA-seq database focusing on TME (45). The distribution of model genes in different cells in the TME of pancreatic cancer was further explored through the TISCH database.
Cell lines and culture
Cell lines of HPDE6-C7, CF-PAC1, Panc-1 and BxPC-3 were from our laboratory, which were consistent with our previous research. And the HPDE6-C7 was a human pancreatic ductal epithelium, the others were pancreatic cancer cell lines. We used the with DMEM mixed with 10% FBS (Gibco, USA) to culture HPDE6-C7, BxPC-3, and Panc-1 cell lines. And CF-PAC1 were cultured with IMDM mixed with 10% (FBS) (Procell, China). All the cell lines were incubated at Cell Incubator.
The validation of the hub RNA expression with quantitative real-time PCR
The total RNAs was isolated from HPDE6-C7 cell lines and CF-PAC1, Panc-1and BxPC-3 cell lines by extraction tool named TRIzol (Accurate Biotechnology). We utilized the Reverse Transcription Reagent to prepare the cDNAs. RT-PCR was carried out by using qPCR Kit (Accurate Biotechnology). The experiment reagents were from our laboratory. And the GAPDH was served as the control standard. The analysis and quantification of RNA expression level adopted the ΔΔCt method. All the primer sequences obtained from GenePharma (Suzhou, China) were for human, which were as follows: IL18, 5’- TCTTCATTGACCAAGGAAATCGG-3’ (Forward), 5’- TCCGGGGTGCATTATCTCTAC-3’ (Reverse); GSK3B, 5’- GCCCAGAACCACCTCCTTTGC-3’ (Forward), 5’- CACCTTGCTGCCGTCCTTGTC-3’ (Reverse); VEGFA, 5’- GCCTTGCCTTGCTGCTCTACC-3’ (Forward), 5’- CTTCGTGATGATTCTGCCCTCCTC-3’ (Reverse).
Statistical analysis
For this project, we utilized R (version 4.1.2) and GraphPad Prism 9 to conduct statistical analysis and visual representation of data. The Wilcoxon rank sum test was conducted to examine the disparities between the two groups. More than two groups were compared using the Kruskal-Wallis rank sum test. Parametric and nonparametric variables were compared using Pearson or Spearman correlations, respectively. Kaplan-Meier with log-rank test was utilized to analyze survival data. p-value < 0.05 denotes a remarkably outcome.
Results
Detection and analysis of differentially-expressed NMRGs
The whole study process is depicted in Figure 1. We identified 5552 DEGs between pancreatic tumors and normal samples by combining TCGA and GTEx databases (Figure 2A), 1173 neuroendocrine regulation-related genes and 1131 metabolism-related genes from GeneCards database (Table S1). In total, 85 differentially-expressed NMRGs were detected by taking intersection (Figures 2B, C). We constructed PPI using 85 differentially expressed NMRGs, and further identified three core modules with 45 hub genes using the MCODE plug-in (Figures 2D, E). 45 hub genes were preserved for subsequent analysis.
Figure 2 Identification and analysis of differentially-expressed neuroendocrine regulation- and metabolism-related genes (NMRGs). (A) A volcano map of differentially expressed genes. (B) A venn diagram of intersection of differentially expressed, neuroendocrine regulation-related and metabolism-related genes. (C) A heat map of differentially-expressed NMRGs. (D) PPI of differentially-expressed NMRGs. (E) Three core modules in PPI. (F) Genetic mutation of hub genes. (G) Frequencies of CNV gain and loss. (H) Location of the CNV on the chromosomes. (I) GO enrichment analysis. (J) KEGG enrichment analysis.
Genetic landscape and enrichment analysis
We analyzed 45 hub genes for mutation landscape and CNV in pancreatic cancer. The results showed that 15 (9.49%) of 158 samples had gene mutations, and the three most common mutated genes were ATM (4%), CTNNB1 (3%) and SKT11 (2%), missense mutation was the most common type, C>T accounted for the highest proportion in single nucleotide variants (SNV) (Figure 2F). CNV was present in all 45 hub genes, and the three most frequent genes are AKT2, TNFRSF11B and VEGFA (Figure 2G). Figure 2H depicted the chromosomal position of the CNV for each of the 45 hub genes.
We conducted GO and KEGG enrichment analyses on these 45 hub genes to delve into their biological roles and mechanisms. GO enrichment analysis illustrated that “regulation of small molecule metabolic process”, “epithelial cell proliferation”, “Positive regulation of kinase activity” and “regulation of polysaccharide process” were the most remarkably enriched pathways in BP, “secretory granule lumen”, “cytoplasmic vesicle lumen”, “vesicle lumen” and “platelet alpha granule lumen” were the most significantly enriched pathway in CC, “signalling receptor activator activity”, “receptor ligand activity”, “cytokine activity” and “insulin receptor binding” were the most remarkably enriched pathway in MF (Figure 2I). KEGG enrichment analysis illustrated that “PI3K−Akt signalling pathway”, “Alzheimer’s disease”, “HIF−1 signalling pathway”, “proteoglycans in cancer”, “non−alcoholic fatty liver disease”, and “thyroid hormone signalling pathway” were the most significantly enriched pathway (Figure 2J). We could find that the enriched pathways are mainly related to metabolism, neuroendocrine system diseases and tumors.
Establishment and validation of the NMRGs-related prognostic model
Patients from the TCGA cohort were randomised at a 7:3 ratio into the training and internal validation sets. Cox and LASSO regression analyses were used to establish a prognostic model with three NMRGs-related genes (Figures 3A-C). Below is the equation of the prognostic model: risk score = (1.69076150270978 * GSK3B expression) + (0.755709134258276 * IL18 expression) + (0.453448880701734 * VEGFA expression). The patient’s survival time in the low-risk group was considerably elevated as opposed to that of patients in the high-risk group, as shown by a survival analysis conducted on the training set, the internal validation set, and the whole TCGA dataset (Figures 3D-F). In comparison to the low-risk category, patients classified as the high-risk category fared worse in terms of overall survival (OS) (Figures 3G-I). The AUC value for 1, 3 and 5 years were 0.726, 0.669 and 0.787 in the training cohort (Figure 3J), 0.652, 0.735 and 0.869 in the internal validation cohort (Figure 3K), and 0.698, 0.700 and 0.840 in the whole TCGA cohort (Figure 3L), respectively. All of these pointed to the high predictive power of our prognostic model. We used GSE62452, GSE57495 and GSE28735 datasets as external validation sets to additionally illustrate the predictive significance of our prognostic model. Longer survival rates were recorded for patients classified in the low-risk category (Figures 3M-O).
Figure 3 Construction and verification of the NMRGs-related prognostic model. (A) LASSO regression analysis with coefficient path diagram. (B) LASSO regression analysis with cross-validation curve. (C) Coefficient of three genes in the prognostic model. Kaplan-Meier curve of the training set (D), internal validation set (E) and the whole TCGA dataset (F). Risk score and survival status distribution of training set (G), internal validation set (H) and the whole TCGA dataset (I). ROC curve of the training set (J), internal validation set (K) and the whole TCGA dataset (L). Kaplan-Meier survival curve of GSE62452 (M), GSE57495 (N) and GSE28735 dataset (O).
NMRG-based consensus clustering
To additionally illustrate the reliability of our prognostic model, we combined all the samples from TCGA, GSE62452, GSE57495 and GSE28735, and then divided them into three clusters using unsupervised clustering analysis (Figures 4A-C). Survival analysis showed that the survival time of Cluster 3 was remarkably higher than Cluster 1 and Cluster 2 (Figure 4D). In all the pooled datasets, the patients having a low risk had a remarkably more favorable prognosis in contrast with those identified as having a high risk (Figure 4E). The alluvial diagram displayed the patients’ distribution in the three NMRG-related clusters and two NMRG-related risk score groups, and all cluster 3 patients were mapped to the low-risk subgroup, and all high-risk patients were mapped to cluster 1 and cluster 3, which indicated that our clusters and groups were reasonable and reliable (Figure 4F). PCA and t-SNE results showed that our clusters and groups could clearly distinguish different patients, and this further demonstrates the good consistency and reliability of our prognostic model (Figures 4G-J).
Figure 4 Consensus clustering. (A) Cumulative distribution function (CDF) curve. (B) Relative change of the area under the CDF curve when cluster number k = 2–9, and the optimal k = 3. (C) Heat map of consensus matrix when k = 3. (D) Kaplan-Meier curve of three clusters. (E) Kaplan-Meier curve of the high- and low-risk score groups. (F) Alluvial diagram of changes in three clusters and two risk score groups. (G) PCA analysis of three clusters. (H) PCA analysis of the high- and low-risk score groups. (I) t-SNE analysis of three clusters. (J) t-SNE analysis of the high- and low-risk score groups.
Clinicopathologic correlation, independent prognostic analysis, and nomogram model construction
The results showed that there was no difference in risk scores between different age and gender groups (Figures 5A, B). Risk scores were higher in higher pathological grade and the difference was statistically significant (Figure 5C). Risk scores increased gradually in the higher TNM stage, and the differences were close to statistical significance (Figure 5D). The expression of IL18 was associated with higher pathological grade and TNM stage (Figure S1). Univariate and multivariate Cox regression analyses confirmed that age and risk score independently acted as prognostic markers for pancreatic cancer (Figures 5E, F). Subsequently, utilizing clinicopathological parameters and risk score, we designed a nomogram model to predict the survival of patients with pancreatic cancer (Figure 5G). The calibration curve illustrated that the 1-, 3-, and 5- years survival rate predicted by the nomogram was close to the real survival, which signified that our nomogram model has outstanding predictive significance (Figure 5H).
Figure 5 Clinicopathologic correlation, independent prognostic analysis, and nomogram model construction. Differences in risk score in age groups (A), gender groups (B), grade groups (C) and stage groups (D). (E) The forest map of univariable Cox regression. (F) The forest map of multivariable Cox regression. (G) The nomogram prediction model based on risk score and clinicopathological characteristics. (H) Calibration curve of nomogram model of predicting 1, 3, 5 years survival rate.
Gene set enrichment analyses
The GSEA was executed to investigate the variations between the high- and low-risk groups in terms of biological processes and metabolic pathways. In total, 274 pathways were considerably enriched in the gene set “c5.go.v7.5.1.symbols.gmt” (Table S2). In the high-risk group, the top 5 pathways with considerable enrichment were “epidermis development”, “keratinization”, “keratinocyte differentiation”, “skin development” and “cadherin binding” (Figure 6A). Additionally, the top 5 enriched pathways in the low-risk group were “B cell receptor signalling pathway”, “regulation of ion transport”, “signal release”, “presynapse” and “T cell receptor complex” (Figure 6B). In total, 6 pathways were considerably enriched in the gene set “c2.cp.kegg.v7.5.1.symbols.gmt”. The pathways with remarkable enrichment in the high-risk group were “cell cycle”, “pathways in cancer”, “small cell lung cancer” and “steroid hormone biosynthesis” (Figure 6C). Furthermore, the substantially enriched pathways in the low-risk were “neuroactive ligand-receptor interaction” and “primary immunodeficiency” (Figure 6D). According to the GSEA results, a remarkable enrichment in the high-risk group was mainly enriched in some pathways related to cancer, whereas the low-risk group was predominantly enriched in pathways linked to immune response.
Figure 6 Gene set enrichment analyses. The top 5 significantly enriched pathways in gene set “c5.go.v7.5.1.symbols.gmt” of the high- (A) and low- (B) risk score groups. The significantly enriched pathways in gene set “c2.cp.kegg.v7.5.1.symbols.gmt” of the high- (C) and low- (D) risk score groups.
Immune analysis
To analyze the link between immune infiltration and risk scores, we used several algorithms to contrast high- and low-risk groups. The “ESTIMATE” algorithm proved that stromal score had no significantly difference between high- and lowrisk groups (Figures 7A). The low-risk category has elevated immune and ESTIMATE scores in contrast with the high-risk population (Figures 7B, C). ssGSEA findings demonstrated that the low-risk category exhibited superior performance in immune cell infiltration and immune-related pathway in contrast with the high-risk subgroup, including NK cells, mast cells, CD8+ T cells, pDCs (Plasmacytoid dendritic cells), TIL (tumour Infiltrating lymphocyte), cytolytic activity, Type II IFN response, and T cell co-stimulation (Figures 7D, E). Immune cell infiltration analysis indicated that the risk score was significantly and inversely linked to naive CD4+ T cells, CD8+ T cells, macrophage M2, NK T cells, B cells, and T cell regulatory (Tregs) infiltration, significantly positively correlated with neutrophil and endothelial cell infiltration (Figure 7F). The expression of GSK3B was significantly inversely linked to the naive CD4+ T cells, memory B cells, NK T cells, and Tregs infiltration, and significantly positively linked to macrophage, neutrophil, cancer-associated fibroblast (CAFs), B cell plasma and activated mast cell infiltration (Figure 7G). The expression of IL18 was strongly and inversely linked to macrophage M2, mast cell and endothelial cell infiltration, and significantly positively correlated with neutrophil, macrophage M1, B cells, and CD8+ T cell infiltration (Figure 7H). The expression of VEGFA was significantly inversely linked to CD8+ T cell, B cell naïve, macrophage M2, endothelial cell infiltration, and significantly positively correlated with eosinophil and macrophage M0 infiltration (Figure 7I). Based on the calculated TIDE score, as well as T cell dysfunction and exclusion scores for each sample, we discovered that the high-risk population exhibited remarkably elevated T cell exclusion score (Figure 7J), whereas the low-risk population had significantly elevated T cell dysfunction score and TIDE score (Figures 7K, L). In addition, further analysis showed that the expression of GSK3B, IL18 and VEGFA were significantly associated with lower TIDE score (Figure S2). This illustrated that pancreatic cancer in the high-risk population had a greater likelihood of responding to immunotherapy as opposed to the low-risk population.
Figure 7 Immune analysis. Stromal (A), immune (B) and ESTIMATE (C) score in the high- and low-risk score groups. 16 immune cell infiltration scores (D) and the activity score of 13 immune-related pathways (E) in the high- and low-risk score groups. Relationship between immune cell infiltration and risk score (F), GSK3B (G), IL18 (H) and VEGFA (I) expression level. T cell exclusion (J), T cell dysfunction (K) and TIDE score (L) in the high- and low-risk score groups. *p<0.05;**p<0.01.
Tumor mutation burden
The “maftools” R package was utilized to analyze the mutation landscape in the high- and low-risk groups. A greater mutation frequency was seen in the high-risk group for the five most frequently mutated genes (TTN, CDKN2A, SMAD4, TP53, and KRAS) (Figures 8A, B). Previous studies indicated that activation mutations of the proto-oncogene KRAS and inactivation mutations of the tumour suppressor gene TP53, SMAD4 and CDKN2A were intimately linked to the occurrence, progression and dismal prognosis of pancreatic cancer (46, 47). TMB was depicted to be substantially enhanced in the high-risk group as per the Wilcoxon test (Figure 8C). Spearman test indicated a positive link between TMB and risk score (Figure 8D). We further analyzed the relationship between model genes and TMB, and the results showed that the expressions of GSK3B and VEGFA were significantly positively correlated with TMB (Figure S3). Survival analysis confirmed that TMB was linked to a worse outcome for patients with pancreatic cancer (Figure 8E). Accordingly, low risk and low TME were correlated with the best prognoses, whereas high risk and high TMB were linked to the worst prognoses (Figure 8F).
Figure 8 Mutation analysis. Genetic mutation landscape in the high- (A) and low- (B) risk score groups. (C) Tumour mutation burden (TMB) of the high- and low-risk score groups. (D) Correlation between TMB and risk score. (E) Kaplan-Meier curve of the high- and low-TMB categories. (F) Kaplan-Meier curve of the different TMB and risk score categories.
Drug sensitivity
Drug therapy is an important treatment for pancreatic cancer, especially for advanced pancreatic cancer. However, different patients have different sensitivity to different drugs. Therefore, it may be more effective and scientific to make individualized treatment plans for different patients. By applying the “pRRophetic” package in R software to predict drug sensitivity, we discovered that in high-risk group patients, 18 drugs (including BIBW2992, Bicalutamide, Gefitinib, Lapatinib, etc.) had significantly lower IC50 values, and in the low-risk group, 26 drugs (including Axitinib, Metformin, Roscovitine, Sunitinib, Vinblastine, etc.) had significantly lower IC50 values (Table 1). We selected 14 drugs shown in Figure 9A. Based on the relevant data from CellMiner database, we found that three genes in the model were associated with the sensitivity of 78 drugs (Supplementary Table 3), and the top 25 drugs with the most significant sensitivity were shown in Figure 9B. VEGFA expression was positively linked to the sensitivity of Abiraterone and Zoledronate, and inversely linked to Fludarabine, Cytarabine and Cladribine. GSK3B expression was inversely linked to Oxaliplatin and brigatinib. IL18 expression was negatively linked to Paclitaxel, VINORELBINE, Vinblastine and Sulfatinib.
Figure 9 Drug sensitivity analysis. (A) Drug sensitivity in the high- and low-risk score groups. (B) Correlation between drug sensitivity and GSK3B, IL18 and VEGFA expression level.
Gene expression verification and distribution analysis
We used the GEPIA platform and discovered that GSK3B, IL18 and VEGFA RNA expression levels were elevated in tumors than in normal tissues (Figures 10A-C). qRT-PCR suggested that GSK3B, IL18 and VEGFA RNA expression level in tumor cells was remarkably higher in contrast with that in normal cells, in line with findings based on GEPIA platform (Figures 10M-O). Immunohistochemistry images derived from the HPA database showed that GSK3B and IL18 expression at protein level in tumour tissues was elevated in contrast with that in normal tissues, in line with findings of RNA expression levels (Figures 10D, E). However, there was no remarkably variations in the protein expression level of VEGFA in tumour tissues in contrast with normal tissues (Figure 10F). Subsequently, we additionally examined the distribution of the three genes’ expression in various subcellular structures and cell types by HPA database. GSK3B was detected in the nucleoplasm and mainly expressed in pancreatic endocrine, ductal and exocrine glandular cells (Figures 10G, J). IL18 was detected in the nucleoplasm, Golgi apparatus, and cytosol, and was also predicted to be secreted extracellular and mainly expressed in mixed cell types (Figures 10H, K). VEGFA was predicted to be secreted extracellular and predominantly expressed in pancreatic endocrine and ductal cells (Figures 10I, L). In addition, the single-cell dataset CAR001160 from the TISCH platform was utilized for further exploring the distribution of model genes in different cells in the TME of pancreatic cancer. Results showed that in the tumor microenvironment of pancreatic cancer, GSK3B was mainly distributed in endothelial cells, malignant cells and B cells, IL18 was mainly distributed in dendritic cells, monocytes/macrophages and malignant cells, and VEGFA was mainly distributed in malignant cells, monocytes/macrophages and ductal cells (Figure S4).
Figure 10 Gene expression verification and distribution analysis. Differences in the expression of GSK3B (A), IL18 (B) and VEGFA (C) at RNA level between pancreatic normal tissues and tumor tissues based on GEPIA platform. Immunohistochemical images of GSK3B (D), IL18 (E) and VEGFA (F) in pancreatic normal and tumor tissues. Distribution of GSK3B (G), IL18 (H) and VEGFA (I) expression in different subcellular structures. Distribution of GSK3B (J), IL18 (K) and VEGFA (L) expression in different cell types. Differences in the RNA expression of GSK3B (M), IL18 (N) and VEGFA (O) between pancreatic normal cells and tumor cells based on RT-PCR. ns, no statistical significance; *p<0.05;**p<0.01;***p<0.001.
Discussion
Pancreatic cancer is a very challenging malignant tumour with insidious onset, rapid progression and poor prognosis (48). The only current hope for curing pancreatic cancer is via major surgery. Unfortunately, by the time most patients are diagnosed, their chances of undergoing radical surgery have already been missed, and the effect of adjuvant therapy such as chemotherapy and radiotherapy on pancreatic cancer is not obvious (49). Finding effective new treatments for pancreatic cancer is crucial. In addition to antiangiogenic therapies and immunotherapies already in clinical practice, metabolic regulation is considered another promising approach for cancer treatment (50). Cancer is characterised in part by the metabolic reprogramming of its tissues. Compared with normal cells, pancreatic cancer cells undergo a series of metabolic alterations: (1) reprogramming the metabolism of intracellular nutrients, such as lipids, amino acids, and glucose; (2) enhancing nutrient supply through scavenging and recycling; (3) microenvironmental interactions involving metabolic processes and other components (51). These metabolic changes are conducive to the survival and division of pancreatic cancer cells in an environment of hypoxia and nutrient deprivation. The metabolism of tumors is often regulated by the neuroendocrine system, and abnormal neuroendocrine regulation may cause metabolic disorders. Studies showed that obesity could lead to abnormal adipose metabolism, chronic inflammation, insulin resistance, and hyperglycemia, and further affect the secretion of different hormones, growth factors, inflammatory cytokines, adipokines, and free fatty acids, which were considered to be the risk biomarkers for cancer morbidity and mortality (25, 52). Besides obesity, neuroendocrine diseases such as diabetes, depression, and anxiety can also enhance the risk of many malignancies and are linked to dismal prognoses (53, 54).
We developed an NMRGs-related prognostic model for pancreatic cancer using data in the TCGA database and by performing Cox and LASSO regression analyses, and the prognostic model was verified by GSE62452, GSE57495 and GSE28735 datasets. Unsupervised clustering, PCA and t-SNE analysis additionally proved the prognostic model’s reliability and consistency. The prognostic model included three genes: GSK3B, IL18 and VEGFA. GSK3B is a multifunctional serine/threonine kinase, which is implicated in various biological activities such as metabolism, cell cycle, DNA damage repair, cell proliferation, and apoptosis, and is associated with diabetes, tumors, psychiatric and neurodegenerative diseases (55–59). Darrington et al. (60) discovered that GSK3B expression level was elevated in prostate cancer (PCa) tissue in contrast with that in normal prostate tissue, and GSK3B inhibitors could reduce the growth of PCa cells. Mamaghani et al. (61) and Ougolkov et al. (62) illustrated that GSK3B expression was elevated in pancreatic cancer samples in contrast to normal pancreatic samples, which was congruent with our findings. Furthermore, the GSK3B inhibitor could suppress the survival and proliferation of pancreatic cancer cells by attenuating the activity of nuclear factor-kappaB (NF-κB). Studies indicated that GSK3B may play two roles in tumors: (1) promoting cancer through induced activation of NF-κB; (2) anti-cancer effect by preventing epithelial to mesenchymal transition (EMT) and metastasis (59, 63). Histone deacetylases (HDACs) could down-regulate the expression of E-cadherin in pancreatic cancer to promote EMT and metastasis, and HDACs inhibitors could suppress the proliferative and migratory capacities of pancreatic cancer cells (64). Edderkaoui et al. (63) found that metavert, a molecule that inhibits both GSK3B and HDACs activity, could significantly reduce tumour size, prevent metastasis, increase the killing of paclitaxel- and gemcitabine-resistant pancreatic cancer cells. IL-18 is a pro-inflammatory and immunomodulatory cytokine of the IL-1 family that is converted from an inactive precursor protein (pro-IL18) by caspase-1-induced cleavage of an N-terminal fragment and may have anti-cancer and oncogenic effects depending on the tissue and cellular environment (65, 66). Liu et al. (67) found that tongue squamous cell carcinoma may be prevented from advancing if IL18 is overexpressed since it may cause apoptosis and decrease the activity of the cells. However, Li et al. (68) illustrated that individuals with colorectal cancer who had an elevated blood IL-18 level had a worse prognosis. Kim et al. (69) found that IL18 could directly enhance the migratory ability of gastric cancer cells by filamentous-actin polymerization and tensin down-regulation. Guo et al. (70) illustrated that the expression level of IL18 was remarkably elevated in pancreatic cancer patient plasma in contrast with pancreatic benign tumors, pancreatitis, and healthy human plasma, elevated in pancreatic cancer tissues in contrast with normal tissues and was linked to a dismal prognosis of pancreatic cancer. This was consistent with our study results. VEGFA is a member of the vascular endothelial growth factor family, which participates in tumour angiogenesis and is intimately linked to tumour development and metastasis, and may be employed as a possible target for tumour therapy (71–73). Our study discovered that high expression of VEGFA was linked to the poor prognosis for pancreatic cancer.
TME is an intricate and comprehensive system in which tumour cells originate and live, which consists of tumour cells, stromal cells, immune cells, and extracellular matrix. TME is intimately linked to tumorigenesis, progression, and patient prognosis (74). GSEA results illustrated that pancreatic cancer in the low-risk group was predominantly enriched in immune-associated pathways. Furthermore, the “ESTIMATE” algorithm confirmed that patients in the low-risk subgroup had an elevated immune score. ssGSEA analysis further confirmed that immunocyte infiltration scores and immune-associated functional pathway scores were elevated in the low-risk group. Immune cell infiltration analysis confirmed that the infiltration degree of CD8+ T cells and NK cells was elevated in the low-risk group, whereas the infiltration degree of CAFs and neutrophil cells was elevated in the high-risk. CD8+ cytotoxic T cells perform an instrumental function in anti-tumour immunity by killing tumour cells, and can also inhibit angiogenesis by secreting interferon-gamma (IFN-γ), which is widely believed to be linked to improved prognosis of tumour patients (75, 76). Similar to CD8+ cytotoxic T cells, NK cells also perform an integral function in anti-tumour immunity, which can produce cytotoxic effects through effector cytokines, cytotoxic molecules and Fas pathway, and then kill tumour cells (74, 77). CAFs constitute the majority of stromal cells in the TME, including antigen-presenting CAFs (apCAFs), inflammatory CAFs (iCAFs), and Myofibrotic CAFs (myCAFs) and can reshape the extracellular matrix to enhance interstitial sclerosis, promote tumour invasion, induce chemotherapy resistance, inhibit antitumor T-cell response, and promote tumour growth (78). Tumour-associated neutrophils have been recognized as key players in malignant transformation, tumour progression, anti-tumour immunity and angiogenesis, and were associated with poor prognosis for advanced cancers and poor outcomes of immune checkpoint inhibitors therapy (79–81). CXCR2, the CXC receptor expressed by neutrophils, can bind with its ligand chemokine family (CXCL1, CXCL2, CXCL3, CXCL5, CXCL7, and CXCL8) to recruit neutrophils to the TME and participate in the mobilization of tumour-associated neutrophils (81, 82). Steele et al. (46) showed that inhibition of CXCR2 could slow tumour formation, prevent metastasis, and enhance the response to chemotherapy and immunotherapy in pancreatic cancer.
Immunotherapy is a promising treatment that has revealed considerable efficacy in numerous tumors, including melanoma, non-small cell lung cancer, renal cell carcinoma, hepatocellular carcinoma, and Hodgkin’s lymphoma (83–85). However, pancreatic cancer does not appear to be sensitive to immunotherapy, with a low overall response rate. The TIDE score was utilized to predict the link between immunotherapy and risk score, and the findings revealed that patients in the high-risk group exhibited a high likelihood of responding to immunotherapy. A greater TMB is often related to a greater rate of immunotherapy response, as evidenced by a series of studies (86, 87). TMB was considerably enhanced in the high-risk, as shown by our research, suggesting a positive link between risk score and TMB. Therefore, patients with pancreatic cancer in the high-risk group may have a greater sensitivity to immunotherapy, which is consistent with the result of TIDE score prediction. Besides surgery, chemotherapy is the main treatment for pancreatic cancer, especially for advanced pancreatic cancer. Our research revealed that the high- and low-risk groups differ remarkably in their susceptibility to certain small molecular drugs and chemotherapeutic medications. The level of IL18 expression is inversely linked to Paclitaxel sensitivity. GSK3B expression was inversely linked to the sensitivity of Oxaliplatin. Paclitaxel is often used together with gemcitabine to enhance the prognosis of patients with pancreatic cancer (88). Oxaliplatin is one of the chemotherapeutic drugs in FOLFIRINOX regimen (first-line treatment for pancreatic cancer) (89). Therefore, our risk score model is helpful to develop individualized treatment plans for patients with pancreatic cancer.
To our knowledge, this is the first study to use bioinformatics to comprehensively analyze the prognostic role of NMRGs in pancreatic cancer. Nonetheless, this investigation is not without its drawbacks. First, our data are from online databases TCGA and GEO, and the real prospective clinical cohorts are needed for further validation. Secondly, basic investigations still need to be conducted to better comprehend the function of NMRGs in the etiology and progression of pancreatic cancer.
Conclusion
In summary, we established an NMRGs-related prognostic risk score model through the TCGA database, and the model was validated using GSE62452, GSE57495 and GSE28735 datasets. Unsupervised clustering analysis, PCA and t-SNE analysis further illustrated that the prognostic model has very good reliability. The prognostic risk score model contained three genes: GSK3B, IL18 and VEGFA, all of which were highly expressed in pancreatic cancer tissue and were associated with poor prognosis. In addition, our prognostic risk score model and model genes were closely linked to the immune infiltration microenvironment, TMB, and drug sensitivity, and can provide evidence for the treatment strategy of pancreatic cancer patients.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Author contributions
BiZ, QY, and BoZ conceptualized and designed the study. BiZ and BoZ collected the data from the database. BiZ and QY analyzed the data. QY and XC conducted the experiments. BiZ, SL, ZW, HL and FM wrote the article. XC and DS revised the manuscript. All authors contributed to the article and approved the submitted version.
Funding
This research was supported by the Key Projects of Overseas Training Foundation of the Higher Education Institutions of Liaoning Province, China (No. 2020GJWZD004); National Natural Science Foundation of China (No. 82000075); National Natural Science Foundation of Liaoning (No. 2020-BS-195).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fendo.2022.1078424/full#supplementary-material
Supplementary Figure 1 | Clinicopathologic correlation. Differences of GSK3B expression in grade groups (A) and stage groups (B). Differences of IL18 expression in grade groups (C) and stage groups (D). Differences of VEGFA expression in grade groups (E) and stage groups (F).
Supplementary Figure 2 | The relationship between TIDE score and model genes. TIDE score in different GSK3B (A), IL18 (B) and VEGFA (C) expression groups.
Supplementary Figure 3 | The relationship between TMB and model genes. TMB in different GSK3B (A), IL18 (B) and VEGFA (C) expression groups. Correlation between TMB and GSK3B (D), IL18 (E) and VEGFA (F) expression.
Supplementary Figure 4 | Single-cell analysis of model genes expression. (A) CAR001160 annotation of all cell types and percentage of each type. (B, C) Percentages and expressions of GSK3B, IL18 and VEGFA in different cells of pancreatic cancer tissues.
Supplementary Table 1 | 5552 differentially expressed genes between pancreatic tumour and normal tissues were identified by TCGA and GTEx databases, and 1173 neuroendocrine regulation-related genes and 1131 metabolism-related genes were obtained in the GeneCards platform.
Supplementary Table 2 | Gene set enrichment analyses (GSEA) identified 274 remarkably enriched pathways between the high- and low-risk score groups.
Supplementary Table 3 | Three genes in the model were associated with the sensitivity of 78 drugs based on the relevant data from CellMiner database.
References
1. GBD 2017 Pancreatic Cancer Collaborators. The global, regional, and national burden of pancreatic cancer and its attributable risk factors in 195 countries and territories, 1990-2017: A systematic analysis for the global burden of disease study 2017. Lancet Gastroenterol Hepatol (2019) 4(12):934–47. doi: 10.1016/S2468-1253(19)30347-4
2. Li B, Zhang C, Wang J, Zhang M, Liu C, Chen Z. Impact of genetic variants of ABCB1, APOB, CAV1, and NAMPT on susceptibility to pancreatic ductal adenocarcinoma in Chinese patients. Mol Genet Genomic Med (2020) 8(6):e1226. doi: 10.1002/mgg3.1226
3. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2019. CA Cancer J Clin (2019) 69(1):7–34. doi: 10.3322/caac.21551
4. Rahib L, Smith BD, Aizenberg R, Rosenzweig AB, Fleshman JM, Matrisian LM. Projecting cancer incidence and deaths to 2030: The unexpected burden of thyroid, liver, and pancreas cancers in the united states. Cancer Res (2014) 74(11):2913–21. doi: 10.1158/0008-5472.CAN-14-0155
5. Mizrahi JD, Surana R, Valle JW, Shroff RT. Pancreatic cancer. Lancet (2020) 395(10242):2008–20. doi: 10.1016/S0140-6736(20)30974-0
6. Conroy T, Desseigne F, Ychou M, Bouché O, Guimbaud R, Bécouarn Y, et al. FOLFIRINOX versus gemcitabine for metastatic pancreatic cancer. N Engl J Med (2011) 364(19):1817–25. doi: 10.1056/NEJMoa1011923
7. Conroy T, Hammel P, Hebbar M, Ben Abdelghani M, Wei AC, Raoul JL, et al. FOLFIRINOX or gemcitabine as adjuvant therapy for pancreatic cancer. N Engl J Med (2018) 379(25):2395–406. doi: 10.1056/NEJMoa1809775
8. Hong D, Rasco D, Veeder M, Luke JJ, Chandler J, Balmanoukian A, et al. A phase 1b/2 study of the bruton tyrosine kinase inhibitor ibrutinib and the PD-L1 inhibitor durvalumab in patients with pretreated solid tumors. Oncology (2019) 97(2):102–11. doi: 10.1159/000500571
9. DeBerardinis RJ, Chandel NS. Fundamentals of cancer metabolism. Sci Adv (2016) 2(5):e1600200. doi: 10.1126/sciadv.1600200
10. Yang LN, Ning ZY, Wang L, Yan X, Meng ZQ. HSF2 regulates aerobic glycolysis by suppression of FBP1 in hepatocellular carcinoma. Am J Cancer Res (2019) 9(8):1607–21.
11. Liberti MV, Locasale JW. The warburg effect: How does it benefit cancer cells? Trends Biochem Sci (2016) 41(3):211–8. doi: 10.1016/j.tibs.2015.12.001
12. Cluntun AA, Lukey MJ, Cerione RA, Locasale JW. Glutamine metabolism in cancer: Understanding the heterogeneity. Trends Cancer (2017) 3(3):169–80. doi: 10.1016/j.trecan.2017.01.005
13. Snaebjornsson MT, Janaki-Raman S, Schulze A. Greasing the wheels of the cancer machine: The role of lipid metabolism in cancer. Cell Metab (2020) 31(1):62–76. doi: 10.1016/j.cmet.2019.11.010
14. Fu Y, Zou T, Shen X, Nelson PJ, Li J, Wu C, et al. Lipid metabolism in cancer progression and therapeutic strategies. MedComm (2020) 2(1):27–59. doi: 10.1002/mco2.27
15. Ringel AE, Drijvers JM, Baker GJ, Catozzi A, García-Cañaveras JC, Gassaway BM, et al. Obesity shapes metabolism in the tumor microenvironment to suppress anti-tumor immunity. Cell (2020) 183(7):1848–66.e26. doi: 10.1016/j.cell.2020.11.009
16. Le Bourgeois T, Strauss L, Aksoylar HI, Daneshmandi S, Seth P, Patsoukis N, et al. Targeting T cell metabolism for improvement of cancer immunotherapy. Front Oncol (2018) 8:237. doi: 10.3389/fonc.2018.00237
17. Lin D, Fan W, Zhang R, Zhao E, Li P, Zhou W, et al. Molecular subtype identification and prognosis stratification by a metabolism-related gene expression signature in colorectal cancer. J Transl Med (2021) 19(1):279. doi: 10.1186/s12967-021-02952-w
18. Grasso C, Jansen G, Giovannetti E. Drug resistance in pancreatic cancer: Impact of altered energy metabolism. Crit Rev Oncol Hematol (2017) 114:139–52. doi: 10.1016/j.critrevonc.2017.03.026
19. Halbrook CJ, Lyssiotis CA. Employing metabolism to improve the diagnosis and treatment of pancreatic cancer. Cancer Cell (2017) 31(1):5–19. doi: 10.1016/j.ccell.2016.12.006
20. Chang CH, Qiu J, O'Sullivan D, Buck MD, Noguchi T, Curtis JD, et al. Metabolic competition in the tumor microenvironment is a driver of cancer progression. Cell (2015) 162(6):1229–41. doi: 10.1016/j.cell.2015.08.016
21. Gunda V, Souchek J, Abrego J, Shukla SK, Goode GD, Vernucci E, et al. MUC1-mediated metabolic alterations regulate response to radiotherapy in pancreatic cancer. Clin Cancer Res (2017) 23(19):5881–91. doi: 10.1158/1078-0432.CCR-17-1151
22. Bader JE, Voss K, Rathmell JC. Targeting metabolism to improve the tumor microenvironment for cancer immunotherapy. Mol Cell (2020) 78(6):1019–33. doi: 10.1016/j.molcel.2020.05.034
23. Hoy AJ, Nagarajan SR, Butler LM. Tumour fatty acid metabolism in the context of therapy resistance and obesity. Nat Rev Cancer (2021) 21(12):753–66. doi: 10.1038/s41568-021-00388-4
24. Goto A, Yamaji T, Sawada N, Momozawa Y, Kamatani Y, Kubo M, et al. Diabetes and cancer risk: A mendelian randomization study. Int J Cancer (2020) 146(3):712–9. doi: 10.1002/ijc.32310
25. Kim DS, Scherer PE. Obesity, diabetes, and increased cancer progression. Diabetes Metab J (2021) 45(6):799–812. doi: 10.4093/dmj.2021.0077
26. Kang C, LeRoith D, Gallagher EJ. Diabetes, obesity, and breast cancer. Endocrinology (2018) 159(11):3801–12. doi: 10.1210/en.2018-00574
27. Giovannucci E, Harlan DM, Archer MC, Bergenstal RM, Gapstur SM, Habel LA, et al. Diabetes and cancer: A consensus report. Diabetes Care (2010) 33(7):1674–85. doi: 10.2337/dc10-0666
28. Boyd AD, Riba M. Depression and pancreatic cancer. J Natl Compr Canc Netw (2007) 5(1):113–6. doi: 10.6004/jnccn.2007.0012
29. Yuan C, Babic A, Khalaf N, Nowak JA, Brais LK, Rubinson DA, et al. Diabetes, weight change, and pancreatic cancer risk. JAMA Oncol (2020) 6(10):e202948. doi: 10.1001/jamaoncol.2020.2948
30. Duan X, Wang W, Pan Q, Guo L. Type 2 diabetes mellitus intersects with pancreatic cancer diagnosis and development. Front Oncol (2021) 11:730038. doi: 10.3389/fonc.2021.730038
31. Klein AP. Pancreatic cancer epidemiology: understanding the role of lifestyle and inherited risk factors. Nat Rev Gastroenterol Hepatol (2021) 18(7):493–502. doi: 10.1038/s41575-021-00457-x
32. Rebours V, Gaujoux S, d'Assignies G, Sauvanet A, Ruszniewski P, Lévy P, et al. Obesity and fatty pancreatic infiltration are risk factors for pancreatic precancerous lesions (PanIN). Clin Cancer Res (2015) 21(15):3522–8. doi: 10.1158/1078-0432
33. Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics (2012) 28(6):882–3. doi: 10.1093/bioinformatics/bts034
34. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. STRING v11: Protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res (2019) 47(D1):D607–13. doi: 10.1093/nar/gky1131
35. Hu Y, Qi Q, Zheng Y, Wang H, Zhou J, Hao Z, et al. Nomogram for predicting the overall survival of patients with early-onset prostate cancer: A population-based retrospective study. Cancer Med (2022) 11(17):3260–71. doi: 10.1002/cam4.4694
36. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612
37. Li T, Fu J, Zeng Z, Cohen D, Li J, Chen Q, et al. TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res (2020) 48(W1):W509–14. doi: 10.1093/nar/gkaa407
38. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med (2018) 24(10):1550–8. doi: 10.1038/s41591-018-0136-1
39. Chalmers ZR, Connelly CF, Fabrizio D, Gay L, Ali SM, Ennis R, et al. Analysis of 100,000 human cancer genomes reveals the landscape of tumor mutational burden. Genome Med (2017) 9(1):34. doi: 10.1186/s13073-017-0424-2
40. Lv J, Zhu Y, Ji A, Zhang Q, Liao G. Mining TCGA database for tumor mutation burden and their clinical significance in bladder cancer. Biosci Rep (2020) 40(4):BSR20194337. doi: 10.1042/BSR20194337
41. Geeleher P, Cox N, Huang RS. pRRophetic: an r package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One (2014) 9(9):e107468. doi: 10.1371/journal.pone.0107468
42. Shankavaram UT, Varma S, Kane D, Sunshine M, Chary KK, Reinhold WC, et al. CellMiner: a relational database and query tool for the NCI-60 cancer cell lines. BMC Genomics (2009) 10:277. doi: 10.1186/1471-2164-10-277
43. Pontén F, Schwenk JM, Asplund A, Edqvist PH. The human protein atlas as a proteomic resource for biomarker discovery. J Intern Med (2011) 270(5):428–46. doi: 10.1111/j.1365-2796.2011.02427.x
44. Zhao Z, Zhang D, Yang F, Xu M, Zhao S, Pan T, et al. Evolutionarily conservative and non-conservative regulatory networks during primate interneuron development revealed by single-cell RNA and ATAC sequencing. Cell Res (2022) 32(5):425–36. doi: 10.1038/s41422-022-00635-9
45. Sun D, Wang J, Han Y, Dong X, Ge J, Zheng R, et al. TISCH: a comprehensive web resource enabling interactive single-cell transcriptome visualization of tumor microenvironment. Nucleic Acids Res (2021) 49(D1):D1420–30. doi: 10.1093/nar/gkaa1020
46. Steele CW, Karim SA, Leach JDG, Bailey P, Upstill-Goddard R, Rishi L, et al. CXCR2 inhibition profoundly suppresses metastases and augments immunotherapy in pancreatic ductal adenocarcinoma. Cancer Cell (2016) 29(6):832–45. doi: 10.1016/j.ccell.2016.04.014
47. Jones S, Zhang X, Parsons DW, Lin JC, Leary RJ, Angenendt P, et al. Core signaling pathways in human pancreatic cancers revealed by global genomic analyses. Science (2008) 321(5897):1801–6. doi: 10.1126/science.1164368
48. Chen Y, Meng J, Lu X, Li X, Wang C. Clustering analysis revealed the autophagy classification and potential autophagy regulators' sensitivity of pancreatic cancer based on multi-omics data. Cancer Med (2022). doi: 10.1002/cam4.4932
49. Chen D, Huang H, Zang L, Gao W, Zhu H, Yu X. Development and verification of the hypoxia- and immune-associated prognostic signature for pancreatic ductal adenocarcinoma. Front Immunol (2021) 12:728062. doi: 10.3389/fimmu.2021.728062
50. Ye Y, Sun X, Lu Y. Obesity-related fatty acid and cholesterol metabolism in cancer-associated host cells. Front Cell Dev Biol (2020) 8:600350. doi: 10.3389/fcell.2020.600350
51. Qin C, Yang G, Yang J, Ren B, Wang H, Chen G, et al. Metabolism of pancreatic cancer: paving the way to better anticancer strategies. Mol Cancer (2020) 19(1):50. doi: 10.1186/s12943-020-01169-7
52. Deng T, Lyon CJ, Bergin S, Caligiuri MA, Hsueh WA. Obesity, inflammation, and cancer. Annu Rev Pathol (2016) 11:421–49. doi: 10.1146/annurev-pathol-012615-044359
53. Wang YH, Li JQ, Shi JF, Que JY, Liu JJ, Lappin JM, et al. Depression and anxiety in relation to cancer incidence and mortality: a systematic review and meta-analysis of cohort studies. Mol Psychiatry (2020) 25(7):1487–99. doi: 10.1038/s41380-019-0595-x
54. Lega IC, Lipscombe LL. Review: Diabetes, obesity, and cancer-pathophysiology and clinical implications. Endocr Rev (2020) 41(1):bnz014. doi: 10.1210/endrev/bnz014
55. Kettunen P, Larsson S, Holmgren S, Olsson S, Minthon L, Zetterberg H, et al. Genetic variants of GSK3B are associated with biomarkers for alzheimer's disease and cognitive function. J Alzheimers Dis (2015) 44(4):1313–22. doi: 10.3233/JAD-142025
56. Chen J, Wang M, Waheed Khan RA, He K, Wang Q, Li Z, et al. The GSK3B gene confers risk for both major depressive disorder and schizophrenia in the han Chinese population. J Affect Disord (2015) 185:149–55. doi: 10.1016/j.jad.2015.06.040
57. Taylan E, Zayou F, Murali R, Karlan BY, Pandol SJ, Edderkaoui M, et al. Dual targeting of GSK3B and HDACs reduces tumor growth and improves survival in an ovarian cancer mouse model. Gynecol Oncol (2020) 159(1):277–84. doi: 10.1016/j.ygyno.2020.07.005
58. Diniz BS, Talib LL, Joaquim HP, de Paula VR, Gattaz WF, Forlenza OV. Platelet GSK3B activity in patients with late-life depression: Marker of depressive episode severity and cognitive impairment? World J Biol Psychiatry (2011) 12(3):216–22. doi: 10.3109/15622975.2010.551408
59. Lin J, Song T, Li C, Mao W. GSK-3β in DNA repair, apoptosis, and resistance of chemotherapy, radiotherapy of cancer. Biochim Biophys Acta Mol Cell Res (2020) 1867(5):118659. doi: 10.1016/j.bbamcr
60. Darrington RS, Campa VM, Walker MM, Bengoa-Vergniory N, Gorrono-Etxebarria I, Uysal-Onganer P, et al. Distinct expression and activity of GSK-3α and GSK-3β in prostate cancer. Int J Cancer (2012) 131(6):E872–83. doi: 10.1002/ijc.27620
61. Mamaghani S, Patel S, Hedley DW. Glycogen synthase kinase-3 inhibition disrupts nuclear factor-kappaB activity in pancreatic cancer, but fails to sensitize to gemcitabine chemotherapy. BMC Cancer (2009) 9:132. doi: 10.1186/1471-2407-9-132
62. Ougolkov AV, Fernandez-Zapico ME, Bilim VN, Smyrk TC, Chari ST, Billadeau DD. Aberrant nuclear accumulation of glycogen synthase kinase-3beta in human pancreatic cancer: Association with kinase activity and tumor dedifferentiation. Clin Cancer Res (2006) 12(17):5074–81. doi: 10.1158/1078-0432.CCR-06-0196
63. Edderkaoui M, Chheda C, Soufi B, Zayou F, Hu RW, Ramanujan VK, et al. An inhibitor of GSK3B and HDACs kills pancreatic cancer cells and slows pancreatic tumor growth and metastasis in mice. Gastroenterology (2018) 155(6):1985–98.e5. doi: 10.1053/j.gastro.2018.08.028
64. Aghdassi A, Sendler M, Guenther A, Mayerle J, Behn CO, Heidecke CD, et al. Recruitment of histone deacetylases HDAC1 and HDAC2 by the transcriptional repressor ZEB1 downregulates e-cadherin expression in pancreatic cancer. Gut (2012) 61(3):439–48. doi: 10.1136/gutjnl-2011-300060
65. Fabbi M, Carbotti G, Ferrini S. Context-dependent role of IL-18 in cancer biology and counter-regulation by IL-18BP. J Leukoc Biol (2015) 97(4):665–75. doi: 10.1189/jlb.5RU0714-360RR
66. Deswaerte V, Nguyen P, West A, Browning AF, Yu L, Ruwanpura SM, et al. Inflammasome adaptor ASC suppresses apoptosis of gastric cancer cells by an IL18-mediated inflammation-independent mechanism. Cancer Res (2018) 78(5):1293–307. doi: 10.1158/0008-5472.CAN-17-1887
67. Liu W, Hu M, Wang Y, Sun B, Guo Y, Xu Z, et al. Overexpression of interleukin-18 protein reduces viability and induces apoptosis of tongue squamous cell carcinoma cells by activation of glycogen synthase kinase-3β signaling. Oncol Rep (2015) 33(3):1049–56. doi: 10.3892/or.2015.3724
68. Li B, Wang F, Ma C, Hao T, Geng L, Jiang H. Predictive value of IL-18 and IL-10 in the prognosis of patients with colorectal cancer. Oncol Lett (2019) 18(1):713–9. doi: 10.3892/ol.2019.10338
69. Kim KE, Song H, Kim TS, Yoon D, Kim CW, Bang SI, et al. Interleukin-18 is a critical factor for vascular endothelial growth factor-enhanced migration in human gastric cancer cell lines. Oncogene (2007) 26(10):1468–76. doi: 10.1038/sj.onc.1209926
70. Guo X, Zheng L, Jiang J, Zhao Y, Wang X, Shen M, et al. Blocking NF-κB is essential for the immunotherapeutic effect of recombinant IL18 in pancreatic cancer. Clin Cancer Res (2016) 22(23):5939–50. doi: 10.1158/1078-0432.CCR-15-1144
71. Tamma R, Annese T, Ruggieri S, Marzullo A, Nico B, Ribatti D. VEGFA and VEGFR2 RNAscope determination in gastric cancer. J Mol Histol (2018) 49(4):429–35. doi: 10.1007/s10735-018-9777-0
72. Zhou P, Xiong T, Chen J, Li F, Qi T, Yuan J. Clinical significance of melanoma cell adhesion molecule CD146 and VEGFA expression in epithelial ovarian cancer. Oncol Lett (2019) 17(2):2418–24. doi: 10.3892/ol.2018.9840
73. Claesson-Welsh L, Welsh M. VEGFA and tumour angiogenesis. J Intern Med (2013) 273(2):114–27. doi: 10.1111/joim.12019
74. Do HTT, Lee CH, Cho J. Chemokines and their receptors: Multifaceted roles in cancer progression and potential value as cancer prognostic markers. Cancers (Basel) (2020) 12(2):287. doi: 10.3390/cancers12020287
75. Durgeau A, Virk Y, Corgnac S, Mami-Chouaib F. Recent advances in targeting CD8 T-cell immunity for more effective cancer immunotherapy. Front Immunol (2018) 9:14. doi: 10.3389/fimmu.2018.00014
76. Anderson NM, Simon MC. The tumor microenvironment. Curr Biol (2020) 30(16):R921–5. doi: 10.1016/j.cub.2020.06.081
77. Vivier E, Tomasello E, Baratin M, Walzer T, Ugolini S. Functions of natural killer cells. Nat Immunol (2008) 9(5):503–10. doi: 10.1038/ni1582
78. Pan X, Zhou J, Xiao Q, Fujiwara K, Zhang M, Mo G, et al. Cancer-associated fibroblast heterogeneity is associated with organ-specific metastasis in pancreatic ductal adenocarcinoma. J Hematol Oncol (2021) 14(1):184. doi: 10.1186/s13045-021-01203-1
79. Kim J, Bae JS. Tumor-associated macrophages and neutrophils in tumor microenvironment. Mediators Inflammation (2016) 2016:6058147. doi: 10.1155/2016/6058147
80. Cui C, Chakraborty K, Tang XA, Zhou G, Schoenfelt KQ, Becker KM, et al. Neutrophil elastase selectively kills cancer cells and attenuates tumorigenesis. Cell (2021) 184(12):3163–77.e21. doi: 10.1016/j.cell.2021.04.016
81. Jin L, Kim HS, Shi J. Neutrophil in the pancreatic tumor microenvironment. Biomolecules (2021) 11(8):1170. doi: 10.3390/biom11081170
82. Nywening TM, Belt BA, Cullinan DR, Panni RZ, Han BJ, Sanford DE, et al. Targeting both tumour-associated CXCR2+ neutrophils and CCR2+ macrophages disrupts myeloid recruitment and improves chemotherapeutic responses in pancreatic ductal adenocarcinoma. Gut (2018) 67(6):1112–23. doi: 10.1136/gutjnl-2017-313738
83. Sangro B, Sarobe P, Hervás-Stubbs S, Melero I. Advances in immunotherapy for hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol (2021) 18(8):525–43. doi: 10.1038/s41575-021-00438-0
84. Quhal F, Mori K, Bruchbacher A, Resch I, Mostafaei H, Pradere B, et al. First-line immunotherapy-based combinations for metastatic renal cell carcinoma: A systematic review and network meta-analysis. Eur Urol Oncol (2021) 4(5):755–65. doi: 10.1016/j.euo.2021.03.001
85. Zhao L, Cao Y. PD-L1 expression level displays a positive correlation with immune response in pancreatic cancer. Dis Markers (2020) 2020:8843146. doi: 10.1155/2020/8843146
86. Samstein RM, Lee CH, Shoushtari AN, Hellmann MD, Shen R, Janjigian YY, et al. Tumor mutational load predicts survival after immunotherapy across multiple cancer types. Nat Genet (2019) 51(2):202–6. doi: 10.1038/s41588-018-0312-8
87. Liu L, Bai X, Wang J, Tang XR, Wu DH, Du SS, et al. Combination of TMB and CNA stratifies prognostic and predictive responses to immunotherapy across metastatic cancer. Clin Cancer Res (2019) 25(24):7413–23. doi: 10.1158/1078-0432.CCR-19-0558
88. Goldstein D, El-Maraghi RH, Hammel P, Heinemann V, Kunzmann V, Sastre J, et al. Nab-paclitaxel plus gemcitabine for metastatic pancreatic cancer: Long-term survival from a phase III trial. J Natl Cancer Inst (2015) 107(2):dju413. doi: 10.1093/jnci/dju413
Keywords: pancreatic cancer, neuroendocrine regulation, metabolism, prognosis, immune, mutation, drug sensitivity
Citation: Zhang B, Yuan Q, Zhang B, Li S, Wang Z, Liu H, Meng F, Chen X and Shang D (2023) Characterization of neuroendocrine regulation- and metabolism-associated molecular features and prognostic indicators with aid to clinical chemotherapy and immunotherapy of patients with pancreatic cancer. Front. Endocrinol. 13:1078424. doi: 10.3389/fendo.2022.1078424
Received: 24 October 2022; Accepted: 23 December 2022;
Published: 20 January 2023.
Edited by:
Ruiqin Han, State Key Laboratory of Medical Molecular Biology, Chinese Academy of Medical Sciences, ChinaReviewed by:
Qingxiang Liu, Mayo Clinic, United StatesDan Zhang, Institute of Semiconductors, Chinese Academy of Sciences (CAS), China
Jialin Meng, University of Science and Technology of China, China
Copyright © 2023 Zhang, Yuan, Zhang, Li, Wang, Liu, Meng, Chen and Shang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Xu Chen, MTE2OTc3NzkyOEBxcS5jb20=; Dong Shang, c2hhbmdkb25nQGRtdS5lZHUuY24=
†These authors have contributed equally to this work