Development and validation of a hypoxia-stemness-based prognostic signature in pancreatic adenocarcinoma

Background: Pancreatic adenocarcinoma (PAAD) is one of the most aggressive and fatal gastrointestinal malignancies with high morbidity and mortality worldwide. Accumulating evidence has revealed the clinical significance of the interaction between the hypoxic microenvironment and cancer stemness in pancreatic cancer progression and therapies. This study aims to identify a hypoxia-stemness index-related gene signature for risk stratification and prognosis prediction in PAAD.

Methods: The mRNA expression-based stemness index (mRNAsi) data of PAAD samples from The Cancer Genome Atlas (TCGA) database were calculated based on the one-class logistic regression (OCLR) machine learning algorithm. Univariate Cox regression and LASSO regression analyses were then performed to establish a hypoxia-mRNAsi-related gene signature, and its prognostic performance was verified in both the TCGA-PAAD and GSE62452 corhorts by Kaplan-Meier and receiver operating characteristic (ROC) analyses. Additionally, we further validated the expression levels of signature genes using the TCGA, GTEx and HPA databases as well as qPCR experiments. Moreover, we constructed a prognostic nomogram incorporating the eight-gene signature and traditional clinical factors and analyzed the correlations of the risk score with immune infiltrates and immune checkpoint genes.

Results: The mRNAsi values of PAAD samples were significantly higher than those of normal samples (p < 0.001), and PAAD patients with high mRNAsi values exhibited worse overall survival (OS). A novel prognostic risk model was successfully constructed based on the eight-gene signature comprising JMJD6, NDST1, ENO3, LDHA, TES, ANKZF1, CITED, and SIAH2, which could accurately predict the 1-, 3-, and 5-year OS of PAAD patients in both the training and external validation datasets. Additionally, the eight-gene signature could distinguish PAAD samples from normal samples and stratify PAAD patients into low- and high-risk groups with distinct OS. The risk score was closely correlated with immune cell infiltration patterns and immune checkpoint molecules. Moreover, calibration analysis showed the excellent predictive ability of the nomogram incorporating the eight-gene signature and traditional clinical factors.

Conclusion: We developed a hypoxia-stemness-related prognostic signature that reliably predicts the OS of PAAD. Our findings may aid in the risk stratification and individual treatment of PAAD patients.

Introduction

Pancreatic adenocarcinoma (PAAD) is one of the most aggressive and fatal gastrointestinal malignancies and has become the fourth leading cause of cancer-related deaths worldwide, severely threatening people’s lives (Rawla et al., 2019). Due to its occult and atypical clinical symptoms, PAAD is difficult to diagnose early. Most patients are diagnosed at a locally advanced stage or metastatic stage and are not eligible for surgical resection, which is the only curative therapy (Skelton et al., 2017). Despite new advances in comprehensive treatment, targeted molecular therapy and immunotherapy, PAAD patients still experience a dismal prognosis with an average 5-year survival rate of less than 5% due to late detection, drug resistance and postoperative metastasis and recurrence (Siegel et al., 2019). For patients with an increased risk of poor outcomes, individualized systemic treatments may help prolong survival and improve quality of life. Therefore, there is still an urgent need to develop an effective predictive model to accurately evaluate the survival outcomes of PAAD patients and provide support for clinical decision making.

Increasing evidence has shown that tumor cell proliferation and growth are highly dependent on the existence of a functional subpopulation of cancer stem cells (CSCs) that play critical roles in tumor metastasis, recurrence and chemoresistance (Tanase et al., 2014; Alferez et al., 2018). Prior research suggests that CSCs promote the development and progression of pancreatic cancer, and high expression levels of CSC markers such as CD44 and CD133 are strongly correlated with neoplasm recurrence and poor prognosis in PAAD (Stoica et al., 2020). Intratumoral hypoxia is a prominent feature of the pancreatic tumor microenvironment facilitating tumor metastasis and invasion and is closely associated with disease progression and poor survival in patients with pancreatic cancer (Erkan et al., 2016). Recent work has established that hypoxia contributes to the induction and maintenance of CSC stemness by upregulating the expression of hypoxia-inducible factors (HIFs), while blockade of HIF-1 activity decreases CSC marker expression and weakens the CSC population (Hao, 2015; Tong et al., 2018; Zeng et al., 2019; Mortezaee, 2021). In addition, Zhang et al. (2018) proposed that hypoxia might synergistically potentiate chemotherapy resistance through stemness induction, highlighting the clinical significance of hypoxia-CSC interactions in the PAAD microenvironment.

Recently, Malta et al. (2018) derived a novel stemness index (mRNAsi) that could reflect the stemness features of cancer samples based on the theory of CSCs using a one-class logistic regression (OCLR) machine learning algorithm. mRNAsi was found to be closely correlated with overall survival (OS) in pancreatic cancer, and patients with high mRNAsi values showed poor prognosis (Tang et al., 2021). In addition, Huang et al. (2021) established an mRNAsi-related prognostic model to successfully predict patient survival in pancreatic cancer. Moreover, multigene prognostic signatures based on the hypoxic microenvironment have been verified to be prognostic markers for pancreatic cancer and could be applied for risk stratification and clinical treatment (Abou Khouzam et al., 2021; Chen et al., 2021; Ding et al., 2021). Hence, targeting the hypoxic microenvironment and establishing reliable prognostic signatures based on the combination of hypoxia and cancer stemness may provide new perspectives for personalized disease management and therapeutic strategies in PAAD.

In this study, we aimed to explore the prognostic value of a hypoxia-stemness-based gene signature in PAAD. Based on the mRNA profile data of PAAD patients from The Cancer Genome Atlas (TCGA) database, we analyzed the correlation between mRNAsi and prognosis. A total of 45 hypoxia-stemness-related genes (HSRGs) with prognostic value were then identified, and functional enrichment analysis was performed to reveal the potential functions of these genes in the pathogenesis and progression of PAAD. A novel prognostic risk model including eight genes (JMJD6, NDST1, ENO3, LDHA, TES, ANKZF1, CITED and SIAH2) was then established, and its predictive performance was verified in both the training dataset and external validation dataset. Besides, we further validated the expression levels of signature genes using the TCGA, GTEx and HPA databases as well as qPCR analysis. Finally, we constructed a nomogram for prognostic prediction in PAAD and analyzed the correlations among the risk score, immune infiltrates and immune checkpoint genes. Overall, our prognostic gene signature and nomogram might be useful tools for the risk stratification and prognosis prediction of PAAD patients.

Materials

Data collection and processing

The RNA-sequencing (RNA-seq) data and clinical information of 180 PAAD samples were downloaded from the TCGA database. The expression profile and clinical data of 130 PAAD samples in the GSE62452 microarray dataset were downloaded from the Gene Expression Omnibus (GEO) database. After the removal of the samples without survival information, a total of 162 PAAD tumor samples in the TCGA database were included in the further analysis as a training cohort, and a total of 65 PAAD tumor samples in GSE62452 were enrolled as an external validation cohort. Since the data used in the present study were obtained from the open database for free, the approval of the Ethics Committee was not needed.

mRNAsi in PAAD and its prognostic value

Through the “gelnet” R package (Version 1.2.1, https://CRAN.R-project.org/package=gelnet), the mRNAsi score of PAAD tumor samples and control samples in the TCGA dataset was calculated by using the OCLR machine learning algorithm (Wei et al., 2021). A significant difference in mRNAsi values between PAAD tumor tissues and normal samples was defined using the intergroup t test, and all PAAD tumor samples were then divided into two groups, namely, the low-mRNAsi and high-mRNAsi groups, according to the median value of mRNAsi. In addition, the association between mRNAsi and the OS of PAAD patients in the two groups was analyzed using the Kaplan-Meier curve method with the “survival” R package.

Identification of HSRGs

The hypoxia-related genes were downloaded from the hallmark gene set in the GSEA database (http://www.gseamsigdb.org/gsea/downloads.jsp), and their expression levels were extracted from TCGA PAAD samples. The Pearson correlation coefficient (PCC) between the expression level and mRNAsi value of each hypoxia gene was then calculated using the Cor function in R software (http://77.66.12.57/R-help/cor.test.html). Genes with thresholds of |PCC| > 0.4 and p < 0.05 were determined to be significantly associated with mRNAsi.

PPI network construction and functional enrichment analysis

The STRING database (http://string-db.org/) was utilized to analyze the functional interactions among the encoded proteins of HSRGs, and those interaction pairs with an interaction score >0.4 were selected to construct a protein-protein interaction (PPI) network, which was then visually presented through Cytoscape (version 3.6.1, http://www. cytoscape.org/). In addition, we performed enrichment analysis of the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways to reveal the underlying biological functions and signaling pathways of HSRGs. p < 0.05 and FDR <0.05 were considered significantly different.

Construction and validation of a prognostic gene signature model

Through the “survival” R package, we conducted univariate Cox regression analysis to identify HSRGs significantly associated with OS (p < 0.05) and ultimately screened out survival-related HSRGs. The least absolute shrinkage and selection operator (LASSO) regression algorithm in the “lars” R package was subsequently applied to acquire the optimal OS-related HSRGs. Then, the risk model was constructed based on the expression levels and LASSO coefficients of the eight-gene signature. Finally, the risk score of each sample was calculated as follows:

$\mathrm{Risk}\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{=}{\displaystyle \sum \mathrm{Coe}{\mathrm{f}}_{\mathrm{genes}}\mathrm{\times Ex}{\mathrm{p}}_{\mathrm{genes}}}$

In the formula, Coef_genes represents the LASSO coefficient of the target gene, and Exp_genes represents the gene expression level.

The risk scores of PAAD patients in TCGA training cohort and GSE62452 validation cohort were calculated, and all patients were then classified into high-risk and low-risk group according to the median value of risk score served as the cutoff value. Kaplan-Meier analysis was performed to analyze the correlation between the risk score and OS. Receiver operating characteristic (ROC) curves were then used to assess the predictive accuracy of the risk model. To further verify the predictive performance of the prognostic risk model, we also performed survival analysis and ROC analysis in the GSE62452 external validation dataset.

Exploration of the mRNA and protein expression levels of the eight signature genes

The expression levels of JMJD6, NDST1, ENO3, LDHA, TES, ANKZF1, CITED, and SIAH2 were compared between PAAD tumor tissues and normal tissues using the TCGA and Genotype-Tissue Expression (GTEx) databases, and the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/) was used to investigate the protein expression levels.

Cell culture and qRT-PCR analysis

The human pancreatic cancer cell lines (BxPC-3, SW1990 and PANC-1) and the normal human pancreatic ductal epithelial cell line (HPNE) were purchased from Cell Bank of Chinese Academy of Sciences (Shanghai, China) and Mingzhou Biotechnology Co., LTD. (Ningbo, China). BxPC-3 and PANC-1 cells were cultured in RPMI 1640 medium (HyClone, United States) with 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin at a humidified incubator at 37°C with 5% CO₂, while SW1990 cells was 90% L-15 medium with 10% fetal bovine serum (FBS) at a humidified incubator at 37°C with 100% air. Meanwhile, HPNE cells were cultured in 70.5% glycoprival DMEM medium with 5% fetal bovine serum (FBS), 23.5% M3 medium, 10 ng/ml human recombinant EGF, 750 mg/ml puromycin, 2.5 g/L D-glucose and penicillin-streptomycin at a humidified incubator at 37°C with 95% air and 5% CO₂.

Total RNA from the normal human pancreatic ductal epithelial cells and human pancreatic cancer cells was extracted using TRIzol reagent (Thermo Fisher Scientific, MA, United States), and 1 μg of total RNA for reverse transcription was prepared using the PrimeScript RT reagent Kit with gDNA Eraser (Toyobo, Japan). Reverse transcription quantitative PCR was performed under the following cycling conditions: 95°C for 3 min, followed by 40 cycles of 95°C for 10 s and 60°C for 30 s. The relative mRNA expression levels of JMJD6, NDST1, ENO3, LDHA, TES, ANKZF1, CITED, and SIAH2 were normalized to GAPDH expression, and calculated by the 2^−ΔΔCt method. The primer sequences are presented in Supplementary Table S1.

Establishment and assessment of a prognostic nomogram

The clinical characteristics of PAAD patients, including age, sex, tumor size, lymph node metastasis, distant metastasis, tumor stage, chronic pancreatitis history, diabetes history, alcohol history, tobacco history, radiotherapy and recurrence, were extracted from the TCGA dataset, and Fisher’s exact test was used to compare the differences between the low- and high-risk groups. Subsequently, univariate and multivariate Cox regression analyses were conducted to identify independent prognostic factors (log rank p < 0.05). Through the “rms” R package, we established a nomogram integrating the risk score, age and targeted molecular therapy to predict the 1-, 3- and 5-year OS rates of PAAD patients. A calibration curve was plotted to evaluate the agreement between the predicted and observed OS probabilities. The concordance index (C-index) was calculated with the “survcomp” R package to evaluate the predictive accuracy of the nomogram, and a C-index >0.70 indicated a good predictive model.

Gene set variation analysis

To further investigate the variation in biological pathways in PAAD patients between the low- and high-risk groups in the training dataset, GSVA enrichment analysis was carried out using the “GSVA” R package. The KEGG and genetic data were downloaded from the GSEA database for GSVA. Adjusted p < 0.05 was considered statistically significant.

Correlation between the risk score and the immune microenvironment

Next, we explored the association between the risk score and the tumor immune microenvironment (TIME), including immune cell infiltration patterns and immune checkpoint molecules. The immune cell type in tumors was classified by utilizing the CIBERSORT method (Chen et al., 2018), and we obtained a total of 22 kinds of immune infiltrating cells, including T cells, B cells, dendritic cells, monocytes, neutrophils, natural killer (NK) cells, eosinophils and mast cells. The immune checkpoint genes included in this study were PDCD1/PD-1, PD-L1/CD274, CTLA4, CD278, CD366, LAG3/CD223, CD73, CD47, BTLA, TIGIT, MYD1, TNFRSF4, TNFRSF9, and VTCN1.

Statistical analysis

Statistical analyses were conducted using R software (version 3.6.1) and GraphPad Prism 8 software (GraphPad Software, CA, United States). All statistical tests with p < 0.05 (two-sided) were considered statistically significant.

Results

mRNAsi in PAAD and survival analysis

This study was performed according to the flow chart presented in Figure 1. The mRNAsi of each PAAD sample was calculated via the OCLR algorithm based on the mRNA profile data downloaded from the TCGA database. As shown in Figure 2A, the mRNAsi in PAAD tumor tissues was remarkably higher compared with that in normal tissues. All PAAD patients were then divided into low- and high-mRNAsi groups based on the median value of mRNAsi, and survival analysis revealed that PAAD patients with high mRNAsi values had a worse prognosis than those with low mRNAsi values (Figure 2B).

FIGURE 1

FIGURE 1. Overall flow chart of our current work.

FIGURE 2

FIGURE 2. mRNAsi and its prognostic value in PAAD. (A) Comparison of mRNAsi between PAAD tumor tissues and normal tissues in the TCGA dataset. (B) Kaplan-Meier curve analysis of patients in the low- and high-mRNAsi groups.

Identification of HSRGs

We downloaded the hallmark gene set from the GSEA database and obtained a total of 200 hypoxia-related genes, and a total of 108 HSRGs were identified according to the threshold criteria |PCC| > 0.4 and p < 0.05 (Supplementary Table S2). These genes were then arranged in order of PCC from lowest to highest, and correlation analysis of these genes was performed to determine the association between gene expression and mRNAsi. The top 3 hypoxia genes exhibiting negative and positive correlations with mRNAsi are represented in Supplementary Figure S1.

PPI network construction and functional enrichment analysis

At the protein level, the interactions among 108 HSRGs products were evaluated using STRING, and a PPI network consisting of 92 nodes and 297 edges was constructed (Figure 3A). In addition, GO and KEGG pathway enrichment analyses were carried out to reveal the biological functions of these genes and participating pathways. The results showed that they were significantly enriched in biological processes, including glycosaminoglycan metabolic process, glycolytic process, angiogenesis, response to hypoxia and negative regulation apoptotic process (Figure 3B). KEGG pathway analysis suggested that these genes were principally involved in the HIF-1 signaling pathway and glycolysis/gluconeogenesis (Figure 3C).

FIGURE 3

FIGURE 3. Construction of a PPI network and functional enrichment analysis of HSRGs. (A) A PPI network with 29 nodes and 297 edges was constructed to evaluate protein interactions. (B) Top 28 enriched biological processes. (C) Top 8 enriched KEGG pathways.

Construction of the eight-gene prognostic signature in PAAD

A total of 45 HSRGs strongly correlated with OS were identified using univariate Cox regression analysis (Supplementary Table S3). A novel prognostic signature composed of eight genes, including Jumonji domain containing 6 (IMJD6), N-deacetylase-N-sulfotransferase-1 (NDST1), enolase 3 (ENO3), lactate dehydrogenase A (LDHA), testin LIM domain protein (TES), ankyrin repeat and zinc finger peptidyl tRNA hydrolase 1 (ANZF1), CBP/p300-interacting transactivator with glutamic acid/aspartic acid-rich carboxyl-terminal domain 2 (CITED2) and seven in Absentia Homolog 2 (SIAH2), was successfully established by LASSO-penalized regression analysis (Supplementary Figure S2). Subsequently, the risk score of PAAD patients in the training cohort was calculated according to the following formula: Risk score = [(−0.003993857) × Expression value of JMJD6] + [(−0.010590801) × Expression value of NDST1] + [(−0.090216539) × Expression value of ENO3]+[0.00552938 × Expression value of LDHA] + [0.085804528 × Expression value of TES] + [(−0.021876003) × Expression value of ANKZF1] + [(−0.040066157) × Expression value of CITED2] + [0.114809344 × Expression value of SIAH2].

Internal and external validation of the predictive performance of the eight-gene signature

The risk scores of PAAD patients in the TCGA training dataset and the GSE62452 external validation dataset were calculated, and the optimal cutoff values of the risk score in the TCGA and GSE62452 datasets were 0.5 and 0.51, respectively (Figures 4E,F). PAAD patents in both datasets were then divided into high-risk and low-risk groups according to the median value. The survival time, survival status, and gene expression of eight genes in each PAAD patient in both datasets are shown in Figures 4E–J. Kaplan-Meier survival analysis revealed a notable correlation between high risk scores and the dismal prognosis of PAAD patients in the TCGA dataset (Figure 4A). The prognostic value of the eight-gene signature was assessed by utilizing time-dependent ROC, and the area under the curve (AUC) values for 1-, 3-, and 5-year OS were 0.936 (95% CI: 0.901–0.974), 0.836 (95% CI: 0.782–0.883) and 0.840 (95% CI: 0.790–0.889), respectively (Figure 4C), demonstrating the good predictive efficacy of the eight-gene signature in predicting the OS of PAAD patients. Corresponding with the results of the TCGA dataset, the survival outcomes of patients in the high-risk group were significantly worse than those of patients in the low-risk group in the GSE62452 dataset (Figure 4B). Time-dependent ROC analysis suggested that the AUCs for 1-, 3-, and 5-year OS were 0.814 (95% CI: 0.762–0.853), 0.784 (95% CI: 0.709–0.821), and 0.714 (95% CI: 0.685–0.767), respectively (Figure 4D), which confirmed the stability of the eight-gene signature in survival prediction for PAAD.

FIGURE 4

FIGURE 4. Internal evaluation and external validation of the prognostic performance of the eight-gene signature. (A,B) Time-dependent ROC analysis of the eight-gene signature for survival prediction in the TCGA training cohort and GSE62452 testing cohort. (C,D) Kaplan-Meier analysis of the correlation between the risk score and the OS of PAAD patients. (E,F) The distribution of the eight-gene risk scores of each PAAD patient. (G,H) Survival status of PAAD patients ranked by risk score. (I,J) The mRNA expression heatmap of the eight genes in the low- and high-risk groups. Red represents upregulation, and blue represents downregulation.

Validation of the expression and alteration of the eight genes in PAAD tissues

To explore the clinical significance of the eight genes in the model, their mRNA expression levels were validated using the TCGA and GTEx databases. As shown in Figure 5A, the mRNA expression levels of LDHA, TES and SIAH2 were obviously increased in PAAD tissues, while those of JMJD6, ANKZF1, ENO3 and CITED2 were significantly decreased in PAAD. HPA analysis showed that the protein levels of LDHA and TES were highly upregulated in cancer tissues compared to normal pancreatic tissue, while those of NDST1, ANKZF1, and SIAH2 were significantly decreased in cancer tissues. However, ENO3 was not detected in PAAD tumor tissues in the HPA database (Figure 5C). In addition, the Kaplan-Meier curve method was used to explore the effect of these eight genes on PAAD prognosis. We found that PAAD patients with overexpression of JMJD6, NDST1, ENO3, ANKZF1, and CITED2 had a relatively favorable prognosis, while patients with high levels of LDHA, TES and SIAH2 had poor survival outcomes (Figure 5B). Furthermore, the qPCR results indicated the mRNA expression levels of JMJD6, ANKZF1, CITED2, and ENO3 was obviously decreased in pancreatic cancer cells versus the normal pancreatic ductal epithelial cells, whereas those of LDHA, TES, and SIAH2 were oppisite (Figure 6).

FIGURE 5

FIGURE 5. Validation of the expression of the eight signature genes in PAAD. (A) The mRNA expression profile of the eight genes in tumor tissues from the TCGA database and normal pancreatic tissues from the TCGA and GTEx databases. (B) Kaplan-Meier curve of the association between the mRNA expression levels of the eight genes and the OS of PAAD patients. (C) The protein expression of the eight genes in pancreatic tumor tissues and normal tissues. The data were obtained from the HPA database. ENO3 was not found in the database.

FIGURE 6

FIGURE 6. Further verification of the mRNA expression levels of seven genes in human pancreatic cancer cell lines and human pancreatic ductal epithelial cell line by RT-qPCR analysis. *p < 0.05, **p < 0.01, ***p < 0.001.

Evaluation of prognostic factors in PAAD

Among the clinicopathological characteristics of 162 PAAD patients acquired from the TCGA dataset, three clinical parameters (pathologic T stage, smoking history and radiotherapy) were observed to be strongly associated with high risk scores (Table 1), and the distribution proportion of the three factors as well as mRNAsi variables between different risk groups was visualized (Figure 7). Univariate and multivariate Cox regression analyses revealed that age and targeted molecular therapy were significantly associated with OS. Moreover, multivariate analysis suggested that the eight-gene risk model was an independent prognostic factor for PAAD (HR = 2.503, p < 0.001) (Figure 8A and Table 2). In addition, the HR value of the risk score was higher than that of age and targeted molecular therapy, indicating that the risk score was more valuable for survival prediction in PAAD.

TABLE 1

TABLE 1. Clinical characteristics and pathological parameters of PAAD patients.

FIGURE 7

FIGURE 7. Associations between the risk score and clinical data as well as mRNAsi values. (A) The association between the risk score and pathologic T grading. (B) The association between the risk score and tobacco history. (C) The association between the risk score and radiotherapy. (D) The association between the risk score and mRNAsi values.

FIGURE 8

FIGURE 8. Construction of a nomogram for OS prediction in the TCGA PAAD dataset. (A) Forest plot of the multivariate Cox regression analysis of the risk score and clinicopathological parameters in PAAD. *p < 0.05, ***p < 0.001. (B) The nomogram incorporating risk score and clinical factors for survival prediction in PAAD. (C) The calibration curve of the nomogram for predicting the 1-, 3- and 5-year OS rates of PAAD patients in the training cohort. The X-axis represents the predicted OS rates, and the Y-axis represents the actual OS rates. The dashed line at 45° indicates the ideal performance, and the C-index was calculated to reflect the predictive accuracy of the nomogram.

TABLE 2

TABLE 2. Univariate and multivariate Cox regression analysis of clinical factors.

Construction and evaluation of a predictive nomogram for PAAD prognosis

The 162 PAAD patients with complete clinical information in the TCGA dataset were used to construct the prognostic nomogram. By integrating the risk score and traditional prognostic factors, including age and targeted molecular therapy, we established a prognostic nomogram to predict the 1-, 3-, and 5-year OS of PAAD patients (Figure 8B). The C-indexes of the nomogram for 1-, 3-, and 5-year OS were 0.745 (95% CI: 0.698–0.803), 0.709 (95% CI: 0.624–0.763) and 0.678 (95% CI: 0.602–0.727), respectively (Figure 8C), suggesting that the 1-year, 3-year and 5-year survival rates predicted by the nomogram were close to the actual survival rates, and the nomogram showed superiority in 1- and 3-year survival prediction but not in 5-year survival prediction for PAAD patients.

GSVA analysis between the low- and high-risk groups

To illustrate the potential molecular mechanisms of the eight-gene prognostic signature, 162 patients in the training cohort were classified into low- and high-risk groups. We further performed GSVA to estimate the difference in the pathway activation state of PAAD samples between the low- and high-risk groups, and a heatmap was used to visualize the top 15 distinct KEGG pathways arranged in order of p value from small to large. We found that several activation pathways of oncogenesis were remarkably enriched in the high-risk group, such as adherens junction, tight junction, sphingolipid metabolism, and pathways in cancer (Figure 9), indicating that abnormal expression of signature genes in PAAD is closely related to cancer-related pathways.

FIGURE 9

FIGURE 9. GSVA enrichment analysis of biological behaviors between the low- and high-risk groups in the training dataset. The heatmap was applied to visualize the top 15 distinct KEGG pathways arranged from small to large according to the p value; red indicates activated pathways, and blue indicates inhibited pathways.

Correlation between the eight-gene signature and immune status

To estimate the immunity relevance of our prognostic signature, the correlation of the risk score with infiltrating immune cells and immune checkpoint genes in PAAD tumor samples was analyzed. A remarkable difference in CD8⁺ T cells (p < 0.007), regulatory T cells (Tregs, p < 0.044) and neutrophils (p < 0.019) was found between the low- and high-risk groups (Figure 10A). In addition, the expression levels of immune checkpoint genes, including PDCD1/PD-1, CTLA4, LAG3, BTLA and TNFRSF4, were significantly decreased in the high-risk group, while the level of CD47 was dramatically increased in the high-risk group compared with the low-risk group (Figure 10B). Moreover, Pearson correlation analysis suggested that the prognostic risk signature was strongly inversely associated with CTLA4 (r = −0.17; p = 0.026), LAG3 (r = −0.17; p = 0.034) and TNFRSF4 (r = −0.34; p = 1.5e-05), whereas positively correlated with CD47 (r = 0.34; p = 1.4e-05; Figures 10C,D). These findings indicated that the eight-gene signature might play a critical role in regulating the immune response.

FIGURE 10

FIGURE 10. Correlation analysis between risk score and immune status. (A) The correlation between the risk score and infiltrating immune cells. (B) The correlation between risk score and immune checkpoint genes. (C) Correlation between risk score and immune checkpoint genes, including BTLA, CD274, CTLA4, LAG3, TNFRSF4 and PDCD1/PD-1. (D) Correlation analysis between risk score and the expression levels of immune checkpoint inhibitors. Red represents a positive correlation, and green represents a negative correlation.

Discussion

PAAD is one of the deadliest cancers in humans and is a hallmark of cellular and phenotypic heterogeneity (Melendez-Zajgla and Maldonado, 2021). The presence of CSCs may help explain the high mortality of pancreatic cancer, since the vital role of CSCs in the occurrence, development and recurrence of pancreatic cancer has been previously described in several studies (Hermann et al., 2007; Nassar and Blanpain, 2016; Stoica et al., 2020). As a key regulator of the cell response to hypoxia in the tumor microenvironment, HIF-1α is known to contribute to the maintenance of CSC stemness (Hao, 2015). Lv et al. (2019) reported that HIF-1α modulated the proliferation and differentiation of stem cells via the Wnt/β-catenin signaling pathway under hypoxic conditions, and activation of HIF-1α enhanced hypoxia-induced cancer stemness, thus facilitating cancer progression and resulting in unfavorable prognosis in PAAD patients (Colbert et al., 2015; Chen et al., 2019). To date, studies have clarified the prognostic value of gene signatures in the survival prediction of patients with pancreatic cancers, and prognostic signatures combined with conventional clinicopathological parameters such as TNM staging and tumor grade have been demonstrated to be superior to a single biomarker (Metzger et al., 2020; Abou Khouzam et al., 2021; Chen et al., 2021). However, few reports have well established reliable prognostic signatures for PAAD based on the combination of hypoxia and cancer stemness.

In this work, we identified a novel eight-gene signature for the survival prediction of PAAD patients based on comprehensive analyses of publicly available data. A prognostic risk model was established based on the eight-gene signature, in which NDST1, LDHA, TES and SIAH2 were significantly upregulated and correlated with poor prognosis, whereas the declined levels of JMJD6, ENO3, ANKZF1 and CITED was associated with adverse survival outcomes of PAAD patients. In addition, the expression levels of eight genes in pancreatic tumor cells was also verified by HPA database and qPCR analyses. Besides, the eight-gene signature was identified as an independent prognostic factor, and its excellent predictive efficiency was confirmed in both the training and external validation datasets. Furthermore, a prognostic nomogram incorporating the eight-gene signature and clinical factors was constructed and could reliably predict 1- and 3-year OS of PAAD patients. Overall, our findings highlighted the important clinical value of the eight-gene signature in the risk assessment and survival prediction of PAAD.

To better elucidate the molecular mechanisms underlying PAAD, we performed functional enrichment analysis of the identified HSRGs and found that they were closely correlated with angiogenesis, hypoxia response and the apoptotic process, and the HIF-1 signaling pathway and glycolysis/gluconeogenesis were the most enriched signaling pathways, suggesting that suggest these genes are tightly associated with cancer-related pathways. Additionally, GVSA analysis showed that adherens junction, tight junction and sphingolipid metabolism were abnormally developed in PAAD patients with high risk scores, revealing that the molecular alterations in the high-risk group were highly associated with the occurrence and development of PAAD. Collectively, these results provide important insights into the pathogenesis and development of PAAD.

Among the eight genes, the crucial roles of ENO3, LDHA and SIAH2 in the progression of pancreatic cancer have been demonstrated. ENO3, as a member of the human enolase (ENO) family catalyzing the transformation of 2-phosphoglycerate to phosphoenolpyruvate during glycolysis, was reported to inhibit the growth of cancer cells (Kong et al., 2016; Feng et al., 2021). Tan et al. (2020) reported that the expression level of ENO3 mRNA was remarkably downregulated in tumor tissues of pancreatic ductal adenocarcinoma (PDAC), and patients with decreased ENO3 levels had poor survival, suggesting that ENO3 is a promising biomarker for prognostic prediction in PADC patients. LDHA is an enzyme that promotes glycolytic processes by converting pyruvate to lactate, and its high expression is correlated with poor outcome in cancer patients (Brand et al., 2016). Knockdown of LDHA prevented tumor growth and metastasis by increasing the production of reactive oxygen species in several cancers (Rizwan et al., 2013; Jin et al., 2017). Previous studies demonstrated the strong correlation of LDHA overexpression with the poor prognosis of patients with pancreatic cancer, and the anticancer effect of LDHA inhibitors offered proof of LDHA as a potential therapeutic target for human cancers (Feng et al., 2018). SIAH2 functions as an E3 ubiquitin ligase that mediates the stabilization of HIF-1α and activates its downstream transcription, which reversely promotes SIAH2 expression through the PI3K/AKT pathway (Matsui-Hasumi et al., 2017). Inhibitors targeting SIAH2 and the SIAH2-HIF-1 axis altered the response of tumor cells to hypoxia and thus affected tumorigenesis and cancer progression (Xu and Li, 2021), highlighting the vital role of SIAH2 as a potential target for cancer therapy. Collectively, these results indicate that the stemness-related hypoxia genes ENO3, LDHA and SIAH2 are promising therapeutic targets in PAAD.

Nevertheless, the clinical significance of JMJD6, NDST1, TES, ANKZF1 and CITED2 has not been revealed in PAAD. Existing evidence suggests the crucial role of JMJD6 as a tumorigenic factor in several cancers, and the high expression of JMJD6 was found to promote the proliferation and survival of tumor cells and predicted the dismal prognosis of patients (Wang et al., 2014; Wong et al., 2019; Biswas et al., 2020). Previous studies have demonstrated the tumorigenic effects of NDST1 in primary glioblastoma and breast cancer, suggesting that it is a promising target for anticancer therapy (He et al., 2014; Xue et al., 2019). TES serves as a tumor suppressor in several cancers, such as gastric cancer and breast cancer, and overexpression of TES reduces the oncogenicity and metastasis of tumor cells (Okolicsanyi et al., 2014; Xue et al., 2019). Xu et al. (2020) proved the association of a glycolytic risk signature, including ANKZF1, with cancer progression and patient prognosis in renal cell carcinoma. Zhou et al. identified ANKZF1 as a prognostic biomarker in colon cancer, and overexpression of ANKZF1 could predict poor OS and recurrence-free survival (Zhou et al., 2019; Chen et al., 2020). CITED2 is highly expressed in many malignancies, including breast cancer, lung cancer, colon cancer and gastric cancer, and it plays critical roles in cancer metastasis and invasiveness (An et al., 2020). In addition, CITED2 is a direct product of HIF-mediated transcription and acts as a vital regulator of stem cell transcription factors in cancer stem-like cells (Rasti et al., 2018; Usui-Ouchi et al., 2020), suggesting that targeting CITED2 expression might be a novel therapeutic strategy for cancer treatment and relapse prevention. Therefore, the roles of JMJD6, NDST1, TES, ANKZF1 and CITED2 in PAAD should be further investigated.

The tumor microenvironment is an immune-related complex environment conducive to tumor cell survival and development. Existing research recognizes the significant effect of immune infiltration on tumorigenesis and metastasis and its association with patient survival in pancreatic cancer (Zhou et al., 2021). Activation of CD8⁺ T cells is crucial for the prevention of tumorigenesis and the induction of tumor regression and correlates with the long survival of patients with pancreatic cancer (Zhang et al., 2020). The prevalence of Treg cells showed a close relationship with tumor metastasis and poor survival of PDAC patients (Tang et al., 2014). Notably, Mota Reyes et al. (2020) revealed that depletion of Tregs could alter the TIME and accelerate pancreatic carcinogenesis and disease progression partially by inducing pathological CD4⁺ T-cell responses. By analyzing the infiltration of immune cells in PAAD samples, we found a remarkable decrease in CD8⁺ T cells (p < 0.007) and Treg cells (p < 0.044) in the high-risk group, which may partly explain the poor prognosis in the high-risk group from an immunological perspective. CTLA-4 and PD-1 are the most studied coinhibitory receptors of T cell receptor signaling, and blockade of CTLA-4 and PD-1 results in T-cell activation and enhanced immune responses (Lee et al., 2018; Henriksen et al., 2019). A preclinical study proved that a CTLA-4 antagonist combined with gemcitabine was safe and tolerated in patients with metastatic PDAC (Aglietta et al., 2014), and PD-1 inhibitors along with chemotherapy showed antitumor effects on advanced PDAC patients, although single-agent checkpoint inhibitors showed disappointing limited activity in pancreatic cancer (Weiss et al., 2017). We further analyzed the association between the risk score and immune checkpoint genes and found that patients in the low-risk group showed higher levels of CTLA4, PD-1, BTLA, LAG3 and TNFRSF4, suggesting that these patients might be more likely to benefit from the combined use of immune checkpoint inhibitors. Collectively, our findings shed light on the close correlation between the risk score and immune status in PAAD.

Our study comprehensively analyzed and identified the HSRGs associated with prognosis on the basis of public TCGA-PAAD database. More importantly, we developed a eight-gene prognostic signature with good capacity in predicting the survival outcomes of PAAD patients. However, our research presents several limitations. The construction and verification of the nomogram for predicting the OS of PAAD patients were conducted based on the TCGA dataset, and external dataset validation needs to be performed in the future. Additionally, we validated the gene expression of the eight signature genes at the cellular level, their expression in PAAD patient specimens also needs further experimental verification. Moreover, more researches are required to elucidate the underlying molecular mechanisms of the eight genes in the pathogenesis and progression of PAAD.

Conclusion

In summary, we developed a novel eight-gene signature and a prognostic nomogram that could reliably predict patient survival in PAAD. Our prognostic gene signature could be beneficial for the risk stratification and prognostic prediction of PAAD.

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

XT designed the study, analyzed the data and wrote the manuscript. JZ collected data and prepared figures and tables. WM, GL, SC, JD, YZ and SC participated in the data analysis. BS and JL critically commented and revised the manuscript. NW conceived of the study and critically revised the manuscript. All authors have read and approved the submitted version.

Funding

This work was supported by grants from the Natural Science Foundation of Zhejiang Province (No. LQ21H160039), National Natural Science Foundation of China (No. 82000774), Natural Science Foundation of Zhejiang Province (No. LQ19H200001), and Taizhou Science and Technology Planning Project (No. 20ywa07).

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/fphar.2022.939542/full#supplementary-material

Supplementary Figure S1 | Top 3 hypoxia genes correlated with mRNAsi. Green represents the top 3 genes negatively correlated with mRNAsi, while red represents the top 3 genes positively correlated with mRNAsi.

Supplementary Figure S2 | Identification of the optimal genes used to develop the prognostic signature by the LASSO regression algorithm. (A) Selection of the optimal parameter (lambda) in the LASSO model. The two vertical lines are lambda.min and lambda.lse. (B) LASSO coefficient profiles of eight OS-related genes with nonzero coefficients determined by the optimal lambda.

Supplementary Table S2 | The 108 hypoxia genes significantly associated with mRNAsi.

Supplementary Table S3 | The 45 HSRGs strongly correlated with OS of PAAD patients.

References

Abou Khouzam, R., Rao, S. P., Venkatesh, G. H., Zeinelabdin, N. A., Buart, S., Meylan, M., et al. (2021). An eight-gene hypoxia signature predicts survival in pancreatic cancer and is associated with an immunosuppressed tumor microenvironment. Front. Immunol. 12, 680435. doi:10.3389/fimmu.2021.680435

PubMed Abstract | CrossRef Full Text | Google Scholar

Aglietta, M., Barone, C., Sawyer, M. B., Moore, M. J., Miller, W. H., Bagala, C., et al. (2014). A phase I dose escalation trial of tremelimumab (CP-675, 206) in combination with gemcitabine in chemotherapy-naive patients with metastatic pancreatic cancer. Ann. Oncol. 25 (9), 1750–1755. doi:10.1093/annonc/mdu205

PubMed Abstract | CrossRef Full Text | Google Scholar

Alferez, D. G., Simoes, B. M., Howell, S. J., and Clarke, R. B. (2018). The role of steroid hormones in breast and effects on cancer stem cells. Curr. Stem Cell Rep. 4 (1), 81–94. doi:10.1007/s40778-018-0114-z

PubMed Abstract | CrossRef Full Text | Google Scholar

An, B., Ji, X., and Gong, Y. (2020). Role of CITED2 in stem cells and cancer. Oncol. Lett. 20 (4), 107. doi:10.3892/ol.2020.11968

PubMed Abstract | CrossRef Full Text | Google Scholar

Biswas, A., Mukherjee, G., Kondaiah, P., and Desai, K. V. (2020). Both EZH2 and JMJD6 regulate cell cycle genes in breast cancer. BMC Cancer 20 (1), 1159. doi:10.1186/s12885-020-07531-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Brand, A., Singer, K., Koehl, G. E., Kolitzus, M., Schoenhammer, G., Thiel, A., et al. (2016). LDHA-associated lactic acid production blunts tumor immunosurveillance by T and NK cells. Cell Metab. 24 (5), 657–671. doi:10.1016/j.cmet.2016.08.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, B., Khodadoust, M. S., Liu, C. L., Newman, A. M., and Alizadeh, A. A. (2018). Profiling tumor infiltrating immune cells with CIBERSORT. Methods Mol. Biol. 1711, 243–259. doi:10.1007/978-1-4939-7493-1_12

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, S., Zhang, J., Chen, J., Wang, Y., Zhou, S., Huang, L., et al. (2019). RER1 enhances carcinogenesis and stemness of pancreatic cancer under hypoxic environment. J. Exp. Clin. Cancer Res. 38 (1), 15. doi:10.1186/s13046-018-0986-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, S., Cao, G., Wu, W., Lu, Y., He, X., Yang, L., et al. (2020). Mining novel cell glycolysis related gene markers that can predict the survival of colon adenocarcinoma patients. Biosci. Rep. 40 (8), BSR20201427. doi:10.1042/BSR20201427

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, D., Huang, H., Zang, L., Gao, W., Zhu, H., Yu, X., et al. (2021). Development and verification of the hypoxia- and immune-associated prognostic signature for pancreatic ductal adenocarcinoma. Front. Immunol. 12, 728062. doi:10.3389/fimmu.2021.728062

PubMed Abstract | CrossRef Full Text | Google Scholar

Colbert, L. E., Fisher, S. B., Balci, S., Saka, B., Chen, Z., Kim, S., et al. (2015). High nuclear hypoxia-inducible factor 1 alpha expression is a predictor of distant recurrence in patients with resected pancreatic adenocarcinoma. Int. J. Radiat. Oncol. Biol. Phys. 91 (3), 631–639. doi:10.1016/j.ijrobp.2014.11.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Ding, J., He, X., Cheng, X., Cao, G., Chen, B., Chen, S., et al. (2021). A 4-gene-based hypoxia signature is associated with tumor immune microenvironment and predicts the prognosis of pancreatic cancer patients. World J. Surg. Oncol. 19 (1), 123. doi:10.1186/s12957-021-02204-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Erkan, M., Kurtoglu, M., and Kleeff, J. (2016). The role of hypoxia in pancreatic cancer: a potential therapeutic target? Expert Rev. Gastroenterol. Hepatol. 10 (3), 301–316. doi:10.1586/17474124.2016.1117386

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, Y., Xiong, Y., Qiao, T., Li, X., Jia, L., Han, Y., et al. (2018). Lactate dehydrogenase A: a key player in carcinogenesis and potential target in cancer therapy. Cancer Med. 7 (12), 6124–6136. doi:10.1002/cam4.1820

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, Z., Li, K., Lou, J., Wu, Y., and Peng, C. (2021). An EMT-related gene signature for predicting response to adjuvant chemotherapy in pancreatic ductal adenocarcinoma. Front. Cell Dev. Biol. 9, 665161. doi:10.3389/fcell.2021.665161

PubMed Abstract | CrossRef Full Text | Google Scholar

Hao, J. (2015). HIF-1 is a critical target of pancreatic cancer. Oncoimmunology 4 (9), e1026535. doi:10.1080/2162402X.2015.1026535

PubMed Abstract | CrossRef Full Text | Google Scholar

He, D. X., Gu, X. T., Li, Y. R., Jiang, L., Jin, J., Ma, X., et al. (2014). Methylation-regulated miR-149 modulates chemoresistance by targeting GlcNAc N-deacetylase/N-sulfotransferase-1 in human breast cancer. FEBS J. 281 (20), 4718–4730. doi:10.1111/febs.13012

PubMed Abstract | CrossRef Full Text | Google Scholar

Henriksen, A., Dyhl-Polk, A., Chen, I., and Nielsen, D. (2019). Checkpoint inhibitors in pancreatic cancer. Cancer Treat. Rev. 78, 17–30. doi:10.1016/j.ctrv.2019.06.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Hermann, P. C., Huber, S. L., Herrler, T., Aicher, A., Ellwart, J. W., Guba, M., et al. (2007). Distinct populations of cancer stem cells determine tumor growth and metastatic activity in human pancreatic cancer. Cell Stem Cell 1 (3), 313–323. doi:10.1016/j.stem.2007.06.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, X. Y., Qin, W. T., Su, Q. S., Qiu, C. C., Liu, R. C., Xie, S. S., et al. (2021). A new stemness-related prognostic model for predicting the prognosis in pancreatic ductal adenocarcinoma. Biomed. Res. Int. 2021, 6669570. doi:10.1155/2021/6669570

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, L., Chun, J., Pan, C., Alesi, G. N., Li, D., Magliocca, K. R., et al. (2017). Phosphorylation-mediated activation of LDHA promotes cancer cell invasion and tumour metastasis. Oncogene 36 (27), 3797–3806. doi:10.1038/onc.2017.6

PubMed Abstract | CrossRef Full Text | Google Scholar

Kong, K. W., Abdul Aziz, A., Razali, N., Aminuddin, N., and Mat Junit, S. (2016). Antioxidant-rich leaf extract of Barringtonia racemosa significantly alters the in vitro expression of genes encoding enzymes that are involved in methylglyoxal degradation III. PeerJ 4, e2379. doi:10.7717/peerj.2379

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, B., Hutchinson, R., Wong, H. L., Tie, J., Putoczki, T., Tran, B., et al. (2018). Emerging biomarkers for immunomodulatory cancer treatment of upper gastrointestinal, pancreatic and hepatic cancers. Semin. Cancer Biol. 52 (Pt 2), 241–252. doi:10.1016/j.semcancer.2017.12.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Lv, Z., Liu, R. D., Chen, X. Q., Wang, B., Li, L. F., Guo, Y. S., et al. (2019). HIF-1α promotes the stemness of oesophageal squamous cell carcinoma by activating the Wnt/β-catenin pathway. Oncol. Rep. 42 (2), 726–734. doi:10.3892/or.2019.7203

PubMed Abstract | CrossRef Full Text | Google Scholar

Malta, T. M., Sokolov, A., Gentles, A. J., Burzykowski, T., Poisson, L., Weinstein, J. N., et al. (2018). Machine learning identifies stemness features associated with oncogenic dedifferentiation. Cell 173 (2), 338–354. doi:10.1016/j.cell.2018.03.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Matsui-Hasumi, A., Sato, Y., Uto-Konomi, A., Yamashita, S., Uehori, J., Yoshimura, A., et al. (2017). E3 ubiquitin ligases SIAH1/2 regulate hypoxia-inducible factor-1 (HIF-1)-mediated Th17 cell differentiation. Int. Immunol. 29 (3), 133–143. doi:10.1093/intimm/dxx014

PubMed Abstract | CrossRef Full Text | Google Scholar

Melendez-Zajgla, J., and Maldonado, V. (2021). The role of lncRNAs in the stem phenotype of pancreatic ductal adenocarcinoma. Int. J. Mol. Sci. 22 (12), 6374. doi:10.3390/ijms22126374

PubMed Abstract | CrossRef Full Text | Google Scholar

Metzger, P., Kirchleitner, S. V., Boehmer, D. F. R., Horth, C., Eisele, A., Ormanns, S., et al. (2020). Systemic but not MDSC-specific IRF4 deficiency promotes an immunosuppressed tumor microenvironment in a murine pancreatic cancer model. Cancer Immunol. Immunother. 69 (10), 2101–2112. doi:10.1007/s00262-020-02605-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Mortezaee, K. (2021). Enriched cancer stem cells, dense stroma, and cold immunity: interrelated events in pancreatic cancer. J. Biochem. Mol. Toxicol. 35 (4), e22708. doi:10.1002/jbt.22708

PubMed Abstract | CrossRef Full Text | Google Scholar

Mota Reyes, C., Teller, S., Muckenhuber, A., Konukiewitz, B., Safak, O., Weichert, W., et al. (2020). Neoadjuvant therapy remodels the pancreatic cancer microenvironment via depletion of protumorigenic immune cells. Clin. Cancer Res. 26 (1), 220–231. doi:10.1158/1078-0432.CCR-19-1864

PubMed Abstract | CrossRef Full Text | Google Scholar

Nassar, D., and Blanpain, C. (2016). Cancer stem cells: Basic concepts and therapeutic implications. Annu. Rev. Pathol. 11, 47–76. doi:10.1146/annurev-pathol-012615-044438

PubMed Abstract | CrossRef Full Text | Google Scholar

Okolicsanyi, R. K., van Wijnen, A. J., Cool, S. M., Stein, G. S., Griffiths, L. R., Haupt, L. M., et al. (2014). Heparan sulfate proteoglycans and human breast cancer epithelial cell tumorigenicity. J. Cell. Biochem. 115 (5), 967–976. doi:10.1002/jcb.24746

PubMed Abstract | CrossRef Full Text | Google Scholar

Rasti, A., Mehrazma, M., Madjd, Z., Abolhasani, M., Saeednejad Zanjani, L., Asgari, M., et al. (2018). Co-Expression of cancer stem cell markers OCT4 and NANOG predicts poor prognosis in renal cell carcinomas. Sci. Rep. 8 (1), 11739. doi:10.1038/s41598-018-30168-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Rawla, P., Sunkara, T., and Gaduputi, V. (2019). Epidemiology of pancreatic cancer: global trends, etiology and risk factors. World J. Oncol. 10 (1), 10–27. doi:10.14740/wjon1166

PubMed Abstract | CrossRef Full Text | Google Scholar

Rizwan, A., Serganova, I., Khanin, R., Karabeber, H., Ni, X., Thakur, S., et al. (2013). Relationships between LDH-A, lactate, and metastases in 4T1 breast tumors. Clin. Cancer Res. 19 (18), 5158–5169. doi:10.1158/1078-0432.CCR-12-3300

PubMed Abstract | CrossRef Full Text | Google Scholar

Siegel, R. L., Miller, K. D., and Jemal, A. (2019). Cancer statistics, 2019. CA Cancer J. Clin. 69 (1), 7–34. doi:10.3322/caac.21551

PubMed Abstract | CrossRef Full Text | Google Scholar

Skelton, R. A., Javed, A., Zheng, L., and He, J. (2017). Overcoming the resistance of pancreatic cancer to immune checkpoint inhibitors. J. Surg. Oncol. 116 (1), 55–62. doi:10.1002/jso.24642

PubMed Abstract | CrossRef Full Text | Google Scholar

Stoica, A. F., Chang, C. H., and Pauklin, S. (2020). Molecular therapeutics of pancreatic ductal adenocarcinoma: targeted pathways and the role of cancer stem cells. Trends Pharmacol. Sci. 41 (12), 977–993. doi:10.1016/j.tips.2020.09.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Tan, Z., Lei, Y., Xu, J., Shi, S., Hua, J., Zhang, B., et al. (2020). The value of a metabolic reprogramming-related gene signature for pancreatic adenocarcinoma prognosis prediction. Aging (Albany NY) 12 (23), 24228–24241. doi:10.18632/aging.104134

PubMed Abstract | CrossRef Full Text | Google Scholar

Tanase, C. P., Neagu, A. I., Necula, L. G., Mambet, C., Enciu, A. M., Calenic, B., et al. (2014). Cancer stem cells: involvement in pancreatic cancer pathogenesis and perspectives on cancer therapeutics. World J. Gastroenterol. 20 (31), 10790–10801. doi:10.3748/wjg.v20.i31.10790

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, Y., Xu, X., Guo, S., Zhang, C., Tang, Y., Tian, Y., et al. (2014). An increased abundance of tumor-infiltrating regulatory T cells is correlated with the progression and prognosis of pancreatic ductal adenocarcinoma. PLoS One 9 (3), e91551. doi:10.1371/journal.pone.0091551

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, R., Liu, X., Wang, W., Hua, J., Xu, J., Liang, C., et al. (2021). Identification of the roles of a stemness index based on mRNA expression in the prognosis and metabolic reprograming of pancreatic ductal adenocarcinoma. Front. Oncol. 11, 643465. doi:10.3389/fonc.2021.643465

PubMed Abstract | CrossRef Full Text | Google Scholar

Tong, W. W., Tong, G. H., and Liu, Y. (2018). Cancer stem cells and hypoxia-inducible factors (Review). Int. J. Oncol. 53 (2), 469–476. doi:10.3892/ijo.2018.4417

PubMed Abstract | CrossRef Full Text | Google Scholar

Usui-Ouchi, A., Aguilar, E., Murinello, S., Prins, M., Gantner, M. L., Wright, P. E., et al. (2020). An allosteric peptide inhibitor of HIF-1α regulates hypoxia-induced retinal neovascularization. Proc. Natl. Acad. Sci. U. S. A. 117 (45), 28297–28306. doi:10.1073/pnas.2017234117

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, F., He, L., Huangyang, P., Liang, J., Si, W., Yan, R., et al. (2014). JMJD6 promotes colon carcinogenesis through negative regulation of p53 by hydroxylation. PLoS Biol. 12 (3), e1001819. doi:10.1371/journal.pbio.1001819

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei, R., Quan, J., Li, S., Liu, H., Guan, X., Jiang, Z., et al. (2021). Integrative analysis of biomarkers through machine learning identifies stemness features in colorectal cancer. Front. Cell Dev. Biol. 9, 724860. doi:10.3389/fcell.2021.724860

PubMed Abstract | CrossRef Full Text | Google Scholar

Weiss, G. J., Waypa, J., Blaydorn, L., Coats, J., McGahey, K., Sangal, A., et al. (2017). A phase Ib study of pembrolizumab plus chemotherapy in patients with advanced cancer (PembroPlus). Br. J. Cancer 117 (1), 33–40. doi:10.1038/bjc.2017.145

PubMed Abstract | CrossRef Full Text | Google Scholar

Wong, M., Sun, Y., Xi, Z., Milazzo, G., Poulos, R. C., Bartenhagen, C., et al. (2019). JMJD6 is a tumorigenic factor and therapeutic target in neuroblastoma. Nat. Commun. 10 (1), 3319. doi:10.1038/s41467-019-11132-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, D., and Li, C. (2021). Regulation of the SIAH2-HIF-1 Axis by protein kinases and its implication in cancer therapy. Front. Cell Dev. Biol. 9, 646687. doi:10.3389/fcell.2021.646687

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, F., Guan, Y., Xue, L., Huang, S., Gao, K., Yang, Z., et al. (2020). The effect of a novel glycolysis-related gene signature on progression, prognosis and immune microenvironment of renal cell carcinoma. BMC Cancer 20 (1), 1207. doi:10.1186/s12885-020-07702-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Xue, J., Yang, M., Hua, L. H., and Wang, Z. P. (2019). MiRNA-191 functions as an oncogene in primary glioblastoma by directly targeting NDST1. Eur. Rev. Med. Pharmacol. Sci. 23 (14), 6242–6249. doi:10.26355/eurrev_201907_18443

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, S., Pottler, M., Lan, B., Grutzmann, R., Pilarsky, C., Yang, H., et al. (2019). Chemoresistance in pancreatic cancer. Int. J. Mol. Sci. 20 (18), E4504. doi:10.3390/ijms20184504

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Z., Han, H., Rong, Y., Zhu, K., Zhu, Z., Tang, Z., et al. (2018). Hypoxia potentiates gemcitabine-induced stemness in pancreatic cancer cells through AKT/Notch1 signaling. J. Exp. Clin. Cancer Res. 37 (1), 291. doi:10.1186/s13046-018-0972-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Lazarus, J., Steele, N. G., Yan, W., Lee, H. J., Nwosu, Z. C., et al. (2020). Regulatory T-cell depletion alters the tumor microenvironment and accelerates pancreatic carcinogenesis. Cancer Discov. 10 (3), 422–439. doi:10.1158/2159-8290.CD-19-0958

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, X., Shang, Y. N., Lu, R., Fan, C. W., and Mo, X. M. (2019). High ANKZF1 expression is associated with poor overall survival and recurrence-free survival in colon cancer. Future Oncol. 15 (18), 2093–2106. doi:10.2217/fon-2018-0920

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, X., Zhu, X., Yao, J., Wang, X., and Wang, N. (2021). Comprehensive analysis of clinical prognosis and molecular immune characterization of tropomyosin 4 in pancreatic cancer. Invest. New Drugs 39 (6), 1469–1483. doi:10.1007/s10637-021-01128-z

PubMed Abstract | CrossRef Full Text | Google Scholar