Hassan Sadozai1*†
Animesh Acharjee2*†
Hateem Z. Kayani1
Thomas Gruber3
Reginald M. Gorczynski4
Bernard Burke1*- 1Centre for Health and Life Sciences, Coventry University, Coventry, United Kingdom
- 2Institute of Cancer and Genomic Sciences, University of Birmingham, Birmingham, United Kingdom
- 3Independent Scholar, National Coalition of Independent Scholars, Visp, Switzerland
- 4Institute of Medical Sciences, University of Toronto, Toronto, ON, Canada
Introduction: Pancreatic ductal adenocarcinoma (PDAC), the most common form of pancreatic cancer, is a particularly lethal disease that is often diagnosed late and is refractory to most forms of treatment. Tumour hypoxia is a key hallmark of PDAC and is purported to contribute to multiple facets of disease progression such as treatment resistance, increased invasiveness, metabolic reprogramming, and immunosuppression.
Methods: We used the Buffa gene signature as a hypoxia score to profile transcriptomics datasets from PDAC cases. We performed cell-type deconvolution and gene expression profiling approaches to compare the immunological phenotypes of cases with low and high hypoxia scores. We further supported our findings by qPCR analyses in PDAC cell lines cultured in hypoxic conditions.
Results: First, we demonstrated that this hypoxia score is associated with increased tumour grade and reduced survival suggesting that this score is correlated to disease progression. Subsequently, we compared the immune phenotypes of cases with high versus low hypoxia score expression (HypoxiaHI vs. HypoxiaLOW) to show that high hypoxia is associated with reduced levels of T cells, NK cells and dendritic cells (DC), including the crucial cDC1 subset. Concomitantly, immune-related gene expression profiling revealed that compared to HypoxiaLOW tumours, mRNA levels for multiple immunosuppressive molecules were notably elevated in HypoxiaHI cases. Using a Random Forest machine learning approach for variable selection, we identified LGALS3 (Galectin-3) as the top gene associated with high hypoxia status and confirmed its expression in hypoxic PDAC cell lines.
Discussion: In summary, we demonstrated novel associations between hypoxia and multiple immunosuppressive mediators in PDAC, highlighting avenues for improving PDAC immunotherapy by targeting these immune molecules in combination with hypoxia-targeted drugs.
Introduction
Pancreatic ductal adenocarcinoma (PDAC), which constitutes approximately 90% of all pancreatic cancer cases, is a highly aggressive malignancy (1, 2). Notably, 80% of cases are diagnosed at a late stage, which precludes surgical resection, the only potential cure (2, 3). Currently, the use of polychemotherapy regimens such as mFOLFIRINOX and gemcitabine/nab-paclitaxel leads to only modest improvements in outcome and the 5-year survival rate is roughly 10% (2, 4). Furthermore, PDAC incidence is increasing in many countries in Europe and in the USA (5, 6). PDAC is primarily driven by mutations in 4 genes (KRAS, TP53, CDKN2A and SMAD4) and displays significant therapeutic resistance to both chemotherapies and radiation (7, 8). In the past decade, cancer immunotherapy, particularly immune checkpoint blockade with antibodies targeted to the checkpoint receptors CTLA-4 and PD-1, have shown remarkable response rates in some cancer types such as melanoma, non-small cell lung cancer and mismatch-repair deficient (dMMR) colorectal cancer (9). Conversely, PD-1 blockade was found to exhibit efficacy only in PDAC cases that are mismatch-repair deficient (10), which constitute only approximately 1% of all PDAC cases (11). Nevertheless, preclinical studies in mice, and some clinical trials have provided evidence that certain immunomodulatory treatments targeted to myeloid cells (e.g. CCR2 inhibitor) can lead to enhanced anti-tumor immunity in PDAC (12). Thus, there remains the possibility that additional research on the immune phenotypes of PDAC could yield novel immunotherapy-based approaches for treating this lethal disease.
A major barrier to effective treatment of PDAC is an exceptionally atypical tumor microenvironment (TME) marked by high levels of desmoplasia (fibrosis), poor vascularization and the abundance of multiple subsets of immunosuppressive myeloid cells such as myeloid-derived suppressor cells (MDSC) and tumor-associated macrophages (TAM) (13, 14). Furthermore, tumor hypoxia, which arises due to abnormal, insufficient vasculature, and enhanced oxygen demand by tumor and stromal cells, is a major pathological hallmark of PDAC (14, 15). The cellular response to hypoxia is primarily driven by the hypoxia-inducible factor (HIF) family of transcription factors namely HIF-1, HIF-2 and HIF-3, comprising an oxygen-sensitive HIF-α subunit which forms a dimer with the constitutively expressed beta subunit (16, 17). An early report demonstrated that HIF-1α was overexpressed in PDAC tissues and absent from non-malignant tissue as well as being associated with advanced disease stage (18). Moreover, multiple studies in the literature have documented that hypoxia modulates multiple facets of disease progression in cancer including genomic instability, EMT (epithelial-mesenchymal transition), metabolic reprogramming and immunosuppression (15–17). Targeting hypoxia to abolish immunosuppression is an exciting prospect for developing more effective treatments for PDAC. In a PDAC murine model, genetic deletion of HIF-2α, but not HIF-1α, in cancer-associated fibroblasts (CAFs), led to significantly delayed tumor growth (19). Moreover, these authors also observed that deletion of HIF-2α in CAFs also led to reduced intra-tumoral infiltration of M2 macrophages and regulatory T cells. However, further studies are warranted to dissect the immune landscape of hypoxic pancreatic adenocarcinomas. Recently, gene signatures for hypoxia have been used to study patient outcomes and characterize the TME using transcriptomics data from a number of cancers including gastric cancer (20), bladder cancer (21), and PDAC (22–24). One of the most widely utilized gene signatures for hypoxia is a 51-gene score developed Buffa et al. (25). Recently, this score was also shown to have prognostic utility in RNA sequencing (RNA-seq) cohorts of lung adenocarcinoma from The Cancer Genome Atlas (TCGA) (26) and was utilized to study hypoxia-related proteins in metastases from 17 distinct carcinomas in TCGA including PDAC (27).
In the present study, we used the Buffa signature, hereafter referred to as the hypoxia score, to profile PDAC cases from TCGA and compare the tumor immune microenvironments of cases with high versus low hypoxia score status. We demonstrated that high hypoxia scores are associated with a distinctly immunosuppressive TME and identified immune genes with a strong correlation to high hypoxia and confirmed these findings in a secondary cohort of RNA-seq data from PDAC cases.
Materials and methods
Cell culture
Human pancreatic cancer cell lines (BxPC3; CRL-1687, PANC-1; CRL-1469) were purchased from the American Tissue Culture Collection (ATCC, USA). PANC-1 cells were maintained in DMEM with 10% FBS, 100 U/ml penicillin, and 100 µg/ml streptomycin. BxPC3 cells were maintained in RPMI-1640 with 10% FBS, 100 U/ml penicillin, and 100 µg/ml streptomycin. PANC-1 (3x105 cells per flask) or BxPC3 (1x106 cells per flask) were seeded in 75cm2 cell culture flasks in 10ml of medium and were incubated under normoxic conditions (20.9% O2, 5% CO2) for 6 days or 4 days, respectively. The medium was then changed, and the cells incubated in normoxia for a further 24h before being cultured for a further 24h under either normoxic or hypoxic (0.5% O2, 5% CO2) conditions.
RNA isolation and Real-time RT-PCR
RNA was isolated using Quick-RNA Miniprep Kit (ZymoResearch) according to the manufacturer’s instructions. Reverse transcription was carried out using a Tetro cDNA Synthesis kit (Bioline). Amplifications were carried out on a BioRad CFX Connect Real-time PCR machine (BioRad, England) using the following cycling parameters: initial denaturation at 95°C for 5 min then 40 cycles of 95°C for 10 seconds, 60°C for 15 seconds, 72°C for 30 seconds. PCR data were normalized to the relative amount of β2M housekeeping gene mRNA determined by separate PCR on each sample. The primer pairs used for each gene are shown in Supplementary Table 1.
TCGA dataset and hypoxia score assessment
Multiple transcriptomics datasets were utilized in this study. For all datasets utilized in this study, no exclusion criteria were applied on the basis of clinical variables. All patients with available gene expression profiles and curated survival data were included in the study. We downloaded uniformly processed gene expression data for TCGA, TARGET and GTEx studies from UCSC Xena Browser (Accessed on 28 September 2020). The analyses in this manuscript were limited to the samples from TCGA PAAD (PDAC) cohort. Sample metadata, including previously curated clinical data (28), were also obtained from Xena Browser. We computed hypoxia scores using the 51-gene signature described by Buffa and colleagues (25). We reviewed the signature to ensure correct mapping to the TCGA expression data gene annotation. The score was computed across all the samples in TCGA PAAD cohort using the rank-based, single sample scoring method implemented in the singscore package (29), with log2-transformed TPM values as input for expression measurement. For the TCGA PAAD cohort, we identified 176 primary tumor samples for which survival outcomes and expression values were available in both TPM (transcripts per million) and RSEM (RNA-seq by Expectation Maximization) expected count formats. The samples in the bottom and top quartiles of hypoxia score were designated as HypoxiaLOW and HypoxiaHI, respectively (n=44 each). These cohorts were used to contrast samples with high and low levels of hypoxia throughout this manuscript, unless otherwise stated.
Survival analysis
Association between overall survival (OS) and hypoxia in the PAAD cohort including sex, age at diagnosis, tumor stage and histological grade as covariates were evaluated within a Cox Proportional Hazard regression model using survival package in R (30, 31). For this analysis, 5 samples with missing stage and grade annotations were excluded leaving 171 samples for the analysis. Granular tumor stages were collapsed into their parent categories (e.g. Stage IA and IB into Stage I). Due to low representation of advanced stages in the dataset, Stages III and IV were collapsed into a single level Stage III/IV. Similarly, samples with histological grade G3 and G4 were combined into a single cohort of G3/4. For this analysis, hypoxia was modelled as a continuous hypoxia score.
To visualise differences in OS and progression-free survival (PFS) between high and low hypoxia cohorts, we generated Kaplan-Meier plots using survminer package (32), contrasting survival curves of patients in the top and bottom quartiles of the hypoxia score. The p values shown on the Kaplan-Meier plots were derived by Log Rank test comparing survival between these two cohorts (30, 31).
Differential gene expression analysis
Differential expression analysis of genes between high and low hypoxia samples was performed using the edgeR package (33). To obtain raw counts, the RSEM expected count matrix obtained from Xena Browser was anti-log2 transformed and the prior count of 1 was subtracted (Accessed on 19 March 2023). To remove lowly expressed genes, the average gene-level log2 counts per million (CPM) values were computed across all samples in the expression matrix. Genes with average log2 CPM value below zero were removed, leaving 17214 genes for subsequent differential expression analysis. We computed scaling normalisation factors using trimmed mean of M values (TMM) method (34). Post-normalisation, we subsetted the expression dataset to samples in HypoxiaLOW and HypoxiaHI groups based on the bottom and the top quartiles of the hypoxia score (n=44 each). The significance of expression difference between the two cohorts was assessed for all of 17214 genes using a quasi-likelihood F test framework following the edgeR package user guide. The gene-wise p value was adjusted for multiple testing using Benjamini-Hochberg method (35). All subsequent references herein regarding genes differentially expressed between high and low hypoxia samples in TCGA dataset refer to this adjusted p value.
Estimating TME cell-type abundance
Cell abundance estimate was performed using Microenvironment Cell Populations-counter (MCP-counter tool) (36), as implemented in immunedeconv package (37). Anti-log2 transformed TPM matrix was used as the gene expression input. To restrict expression table to uniquely annotated genes, we removed rows with duplicated gene symbols keeping the entry with the highest standard deviation.
Statistical comparisons in immune cell frequency between high and low hypoxia groups were performed using Mann-Whitney U tests.
Figures and text referencing differences in cell frequencies report p values adjusted using the Benjamini-Hochberg method (35). Signature scores for cDC1 were computed using the singscore package (29). Log2-transformed TPM matrix was used as gene expression input and the following cDC1 signature from a recent publication: CLEC9A, XCR1, CLNK and BATF3 (38) was selected. Association between the hypoxia and cDC1 activation score was assessed by contrasting samples within top and bottom quartiles of the hypoxia score using a Mann-Whitney U test (HypoxiaHI vs HypoxiaLOW).
Analysis of hypoxia between PDAC and normal pancreas tissue
To compare hypoxia levels between tumor and normal pancreas tissue, a previously published dataset was utilized (39). The expression matrix, sample metadata and gene annotation tables were downloaded from GEO (GSE28735). The expression matrix was quantile-normalised using the affyPLM package (40), and log2-transformed. The resulting matrix was used as input to the singscore package (29), to compute hypoxia scores. The difference in hypoxia scores between matched tumor and normal samples was assessed using a paired Mann-Whitney U test.
Random forest, feature elimination process and validation in additional cohort
Random Forest (41), a machine learning ensemble method was used to find the immune gene, the expression level of which could best identify the difference between HypoxiaHI vs. HypoxiaLOW tumors. Random Forest uses a bootstrapping method to pick random samples from the dataset. The samples are further split into training samples (two thirds of the set) and testing samples (one third of the set). The testing samples, which are also known as “out-of-bag” samples, are used to determine the accuracy of future predictions. This approach warrants that the operator set the number of trees (ntree) and the number of genes that are randomly chosen as options at each split for the Random Forest model to work. The ntree value was set to 500, and the mtry value was the square root of the factors.
To select the best genes that predict HypoxiaHI vs HypoxiaLOW tumors, iterative fitting to create random forests was utilized. With each iteration, a new forest was begun by removing 20% of the factors that were the least important. The set of variables that was chosen is used to predict how well the model will fit so that the “out-of-bag” error rate can be checked. The varSelRF function from the varSelRF package in R was used to perform this recursive feature elimination method process (41, 42).
The 5 genes retained in the final model were validated for their association with hypoxia in an independent RNA-seq dataset of 51 pancreatic adenocarcinoma samples (GSE79668). Pre-processed gene-level counts and samples metadata from GEO were downloaded. A filtering threshold was applied for lowly expressed genes by removing genes with average log2 count per million below value of 0, while ensuring that all genes from Buffa signature were retained. Normalisation factors were calculated using the TMM method (34), and computed log2-transformed count per million (CPM) matrix with a prior count of 1. The resulting log2 CPM matrix was used as input into the singscore method to calculate hypoxia score and as the measure of gene expression to perform correlation analysis between each of the genes and the hypoxia score.
Statistical analyses
Statistical analyses were performed using R v3.6.0 or GraphPad Prism v10. The statistical tests used to compare data are presented in the figure legends. Where data are shown as boxplots, the box represents the IQR, and the whiskers extend 1.5x IQR below and above bottom and top quartiles respectively. A value of p<0.05 was deemed to be statistically significant. Low p values were reported as p<0.001.
Results
Assessment of the Buffa hypoxia score for PDAC profiling
The investigations performed in this study are depicted as a flow chart in Supplementary Figure 1. First, we profiled the expression levels of the previously developed hypoxia score in pancreatic cancer (Figure 1). Using a rank-based gene signature scoring method for single samples (singscore) (29), we profiled 176 PDAC cases from the TCGA cohort for hypoxia score expression (Figure 1A). Next, we examined the expression of the hypoxia score across histological tumor grades. Higher tumor grades represent increased anaplasia and disordered structure and are generally associated with aggressive disease (43). As represented in Figure 1B, we observed a significant increase in the hypoxia score across tumor grades from Grade 1 to Grade 2 (p=0.003) and from Grade 2 to Grade 3/4 (p=0.003) (Grades 3 and 4 were collapsed into a single category due to there being so few cases of Grade 4 in our cohort). Here, we demonstrate that PDAC cases display varying levels of hypoxia score expression but also that the hypoxia score is increased in higher grade tumors suggesting that hypoxia is associated with aggressive disease. To associate the immune phenotypes of cases with high hypoxia, we divided our cohort into HypoxiaHI and HypoxiaLOW groups based on a top versus bottom quartile dichotomization (n=44 patients each). All subsequent analyses, unless stated otherwise, were performed using these groups. Clinical and pathological parameters for both groups are shown in Supplementary Table 2.
Figure 1 Assessing Buffa signature score in PDAC. (A) Box plot of sample-level hypoxia scores across the cohort of PDAC cases in TCGA. Samples with the hypoxia score in top and bottom quartiles were assigned into HypoxiaHI and HypoxiaLOW cohorts respectively. Each dot represents a single sample. (B) Box plots comparing hypoxia score distribution across reported histological grades. Statistical significance was determined using pair-wise Mann-Whitney U tests between groups and Holm-Sidak multiple testing correction was performed. ** p<0.01.
In order to further explore the utility of the hypoxia score for pancreatic cancer, we performed the following analyses. First, we used a microarray dataset of 45 PDAC tumors and paired adjacent non-malignant pancreatic tissue (GSE28735) (39), to show that the hypoxia score was significantly elevated (p<0.0001) in pancreatic tumors relative to the non-tumor tissue sample group (Supplementary Figure 2). This dataset was used as it contained transcriptomic profiles from PDAC biopsies and paired non-malignant pancreatic tissue, the latter of which are not available in the TCGA PDAC cohort. We also sought to study the expression of key genes from our signature in PDAC cells incubated in hypoxic conditions. As such, we examined the expression of 4 key hypoxia-inducible genes from the 10 “seed” genes that were used to define the Buffa hypoxia score (25), in PDAC cell lines incubated in hypoxic conditions (0.5% O2) for 24 hours. Using two well-known PDAC cell lines, PANC-1 and BxPC3 (44), we interrogated the expression of the genes Solute Carrier Family 2 Member 1 (SLC2A1), also known as GLUT1, Vascular Endothelial Growth Factor A (VEGFA), Carbonic Anhydrase IX (CA9), and Hexokinase 2 (HK2) as shown in Supplementary Figure 3. In PANC-1, all four genes were significantly induced in hypoxia at 24h compared to normoxic controls: SLC2A1 (6.5-fold, p=0.007), VEGFA (2.9-fold, p= 0.01), CA9 (7.6-fold, p= 0.002), and HK2 (34-fold, p=0.008). Comparable results were also seen in the BxPC3 cell line relative to normoxic controls: thus, SLC2A1 (17.8-fold, p=0.005), VEGFA (3.2-fold, p=0.03), CA9 (73.4-fold, p=0.01), and HK2 (9.9-fold, p=0.004) when compared with normoxia controls. Taken together, these results further support the validity of this hypoxia score in PDAC. The score was found on average to be significantly elevated in PDAC relative to non-malignant pancreas tissues, and that a subset of hypoxia-inducible genes that comprise the score show induction in PDAC cell lines in hypoxic conditions. Next, we investigated if the hypoxia score was associated with patient outcomes. We compared the survival of cases with low and high hypoxia (HypoxiaHI vs HypoxiaLOW, see Figure 1B). It was observed that HypoxiaHI patients demonstrated significantly lower OS (Figure 2A) and lower PFS (Figure 2B), compared to HypoxiaLOW patients. Finally, we used the multivariate Cox Proportional Hazards model to assess the association between overall survival in the entire TCGA PDAC cohort (n=171 with available clinical information) and the hypoxia score as well as, clinical variables (Age, Sex, Stage and Grade). When examined in this multivariate regression analysis, we noted that only the hypoxia score was independently associated to OS in the TCGA cohort (Supplementary Figure 4). The association between the Buffa hypoxia score and reduced survival has been previously shown in PDAC (45). However, we demonstrated in a multivariate analysis that this hypoxia score is an independent prognostic variable in PDAC. As such, investigating HypoxiaHI cases might reveal immunological mechanisms associated with disease progression in these patients.
Figure 2 Prognostic analysis of hypoxia score. Kaplan-Meier plots displaying the (A) overall survival (OS) and (B) progression-free survival (PFS) of PDAC patients in HypoxiaHI versus HypoxiaLOW groups. The number of cases in each group and number at risk are shown in the tables below each plot. p value shown in the plots was derived using log rank test.
High hypoxia status is associated with multiple features of immunosuppression in PDAC
It is now well-established that hypoxia is a crucial mediator of immune escape, as well as of immunosuppressive signaling in the TME of solid organ cancers (16, 46). Immunohistochemical studies in patient biopsies of colorectal cancer have shown that the hypoxia marker CAIX is strongly expressed in “immune cold” tumors (47), which are marked by low or absence of T cell infiltrates and dendritic cells (DC) (48–50). As such, we performed analyses to dissect the immune microenvironment of HypoxiaHI and HypoxiaLOW cases. First, we used the MCP-counter method developed by Becht et al. for estimating cell-types in bulk tissue transcriptomes (36), to infer the abundance of 8 immune and 2 non-immune cell populations in both hypoxia groups. After correcting for multiple comparisons statistically, we found that the HypoxiaHI group displayed notable differences in the abundance of key cell types in the TME, compared to the HypoxiaLOW group (Figure 3). No statistically significant differences were observed for B lineage cells, CAFs, monocyte lineage cells and neutrophils between both hypoxia groups. In contrast, HypoxiaHI tumors displayed markedly lower abundance scores for endothelial cells (p=0.02), myeloid dendritic cells (p=0.02), natural killer (NK) cells (p=0.02), total T cells (CD3+ T cells) (p=0.02), and CD8+ T cells (p=0.006) compared to HypoxiaLOW cases. Compared to HypoxiaLOW cases, the HypoxiaHI group also displayed lower abundance of “cytotoxic lymphocytes”, a functionally defined signature meant to score for mRNA expression from both NK cells and T cells (36). Given the crucial role of DC in cancer immunity (51) and noting the difference between hypoxia groups in terms of myeloid DC scores, we sought to further explore this difference. Over the past 6 years, studies have shown that the conventional DC1 (cDC1) subset of DC, are critical for antigen cross-presentation and priming anti-tumor immunity (51). Using a previously defined gene signature for cDC1, we used singscore to compute cDC1 scores (Supplementary Figure 5). We demonstrate that HypoxiaHI cases exhibit significantly lower cDC1 scores compared to HypoxiaLOW cases. It is currently well-established that determining a tumor’s immunological phenotype depends not only on the presence of T cells but also the presence of important leukocytes such as NK cells and DC (50). We found that HypoxiaHI cases displayed multiple features of a “cold” tumor with diminished anti-tumor immunity.
Figure 3 Profiling the TME of high versus low hypoxia scores. Box plots comparing 8 immune and 2 non-immune cell scores between HypoxiaHI and HypoxiaLOW groups. Cell abundance scores were computed using the MCP-counter algorithm. Statistical comparisons were performed through Mann-Whitney U test followed by Benjamini-Hochberg multiple testing correction. *p<0.05, ** p<0.01.
In order to investigate the molecular mechanisms that might account for a “cold” tumor status in cases with high hypoxia scores, we performed differential gene expression profiling between HypoxiaHI and HypoxiaLOW cases and performed FDR corrections. Based on the results of the MCP-counter analyses, we compared the expression of selected immune-related genes coding for molecules known to negatively regulate anti-tumor immunity. First, we profiled key immune checkpoint molecules (which can be expressed both on tumor and non-tumor cells), that regulate both adaptive and innate immunity (52, 53). These include the “don’t-eat-me” signaling molecules which inhibit macrophage phagocytosis of tumor cells, CD24 and CD47, the checkpoint molecules CD274 (PD-L1), PDCD1LG2 (PD-L2) and CD276, as well as the non-classical major histocompatibility complex I (MHC-I) molecule, HLA-G which inhibits NK cells (52, 53). As shown in Figure 4A, we observed a statistically significant increase in the expression of CD47, CD276 and HLA-G but not of the other checkpoint molecules in HypoxiaHI tumors relative to HypoxiaLOW. Second, we profiled genes encoding 3 of the best characterized members of the galectin family (Galectin-1, Galectin-3, and Galectin-9), which have roles in both cancer progression and immune modulation (54), as well as Galectin-4, which was recently shown to be associated with immune escape in PDAC and capable of inducing T cell apoptosis (55). Remarkably, expression of all four galectin genes LGALS1, LGALS2, LGALS3 and LGALS4 was observed to be significantly higher in HypoxiaHI as compared to HypoxiaLOW cases (Figure 4B). Third, we interrogated both groups for genes encoding key enzymes known to play critical roles in immunosuppression in the TME as shown in Figure 5. We analyzed expression of the ectonucleotidases CD39 (ENTPD1) and CD73 (NT5E) which mediate distinct steps in the conversion of extracellular ATP to adenosine, a potent inhibitory signal for immune cells (56), and the potential inflammatory mediator cyclooxygenase-2 (COX-2, PTGS2) (57). We also analyzed the expression of genes for 6 genetically unrelated amino acid catabolizing enzymes known to have immunosuppressive functions. As reviewed previously (58–60), the amino acids they catabolize are tryptophan (tryptophan 2,3-dioxygenase: TDO2 as well as indoleamine 2,3-dioxygenase 1 and 2: IDO1 and IDO2), arginine (inducible nitric oxide synthase: NOS2 as well as arginase 1 and 2: ARG1 and ARG2), and phenylalanine (Interleukin 4-induced 1: IL4I1) (58–60). It was observed that ARG1 and IDO2 mRNA expression levels in our TCGA PDAC cohort were below the abundance threshold applied in differential gene expression analysis, and therefore these genes were not included in subsequent analyses. When we compared both hypoxia groups, we noted that relative to HypoxiaLOW tumors, HypoxiaHI cases displayed lower expression for ENTPD1 (CD39) and ARG2 but significantly elevated expression of NT5E (CD73) and PTGS2 (COX-2) (Figure 5). No differences were observed between groups for the expression of IDO1, TDO2, NOS2 and IL4I1. Thus, these results demonstrate the activation of distinct metabolic pathways in cases with a high versus low hypoxia status. Taken together, our gene expression profiling demonstrated that HypoxiaHI tumors displayed the upregulation of multiple molecules associated with an immunosuppressive TME. These included immune checkpoints, galectins and key metabolic mediators indicating that distinct pathways underlie the formation of an immunosuppressive TME in PDAC.
Figure 4 Comparisons of immune checkpoints and galectins gene expression profiles. Box plots comparing gene expression of (A) adaptive and innate immune checkpoint molecules and (B) selected galectin molecules with roles in immune escape. Gene expression is displayed as log2-transformed, normalized CPM (counts per million) values. Statistical significance shown on the plot represents FDR values from transcriptome-wide differential gene expression analysis between HypoxiaHI and HypoxiaLOW groups. *FDR<0.05, **FDR<0.01,***FDR<0.001.
Figure 5 Comparisons of immunoregulatory metabolic mediators. Box plots comparing gene expression of enzymes purported to play a key role in promoting an immunosuppressive TME. Gene expression is displayed as log2-transformed, normalized CPM (counts per million) values. Statistical significance shown on the plot represents FDR values from transcriptome-wide differential gene expression analysis between HypoxiaHI and HypoxiaLOW groups. **FDR<0.01,***FDR<0.001.
Machine learning-based feature selection to predict high hypoxia status
In Figures 4 and 5, we investigated 18 genes with known immunosuppressive functions to show that multiple, distinct genes are elevated in HypoxiaHI tumors. Next, we sought to further explore which molecules were most strongly associated with a high hypoxia status. Using a random forest machine learning algorithm together with feature selection method (VarSelRF package using R) (61), to find a minimum set of genes that could predict HypoxiaHI status. This method revealed 5 genes that could predict HypoxiaHI and ROC (Receiver Operating Characteristic) analysis for these 5 genes displayed an AUC (Area Under ROC Curve) value of 0.961 (Figure 6A). The 5 genes, ranked by their AUC values are shown in Figure 6B, are LGALS3, NT5E, CD276, LGALS1, and ENTPD1 (CD39). It is pertinent to note that ENTPD1 displayed an inverse correlation with HypoxiaHI status. To further assess the association of these genes with high hypoxia scores (i.e. HypoxiaHI) in PDAC, we performed Spearman correlation between expression levels of these genes and the hypoxia score in an additional cohort of RNA-seq data from 51 PDAC patients (GSE79668) (62). The expression of LGALS3, CD276, NT5E was positively correlated to hypoxia score expression in this cohort. The correlation coefficients and p values are shown in Figure 7. The gene LGALS3 displayed the highest and most statistically significant positive correlation (r=0.68, p<0.0001) with hypoxia score expression. Similar to what was observed for the TCGA cohort, ENTPD1, also exhibited a negative significant correlation to the hypoxia score in this secondary cohort. LGALS1 did not show a statistically significant positive or negative correlation with hypoxia score expression. Collectively, these findings suggest that immunomodulatory genes such as LGALS3, CD276 and NT5E are the primary molecular signals that distinguish a highly hypoxic TME in PDAC cases. Given that LGALS3 expression was observed to be strongly and positively correlated with the hypoxia score in both the TCGA cohort and the additional cohort, we interrogated the expression of LGALS3 in PDAC cell lines in hypoxic conditions in vitro. As shown in Figure 8, LGALS3 was significantly induced by 24h of hypoxia (0.5% O2) in BxPC3 cells (4.2-fold, p=0.02). While we did note a trend towards increased expression in PANC-1 cells, these results did not achieve statistical significance (1.7-fold, p=0.13).
Figure 6 Selection of immune genes with high correlation to HypoxiaHI status in TCGA. (A) ROC curve and AUC value for the top 5 gene features for classifying a case as HypoxiaHI. (B) Individual gene names and AUC values for each of the 5 gene features are provided in the table. Note that ENTPD1 exhibits an inverse relationship with hypoxia scores.
Figure 7 Confirmation of HypoxiaHI correlated immune genes in additional cohort. Scatter plots depicting correlation between hypoxia score and the top 5 genes that could best classify a case as HypoxiaHI in validation dataset of 51 PDAC cases (GSE79668). Spearman correlation coefficients r, and p values are plotted on the graphs. Each dot represents a single tumor sample.
Figure 8 Expression of LGALS3 in hypoxic pancreatic cancer cells. Barplots showing normalized mRNA expression for LGALS3 in PANC-1 and BxPC3 pancreatic cancer cell lines incubated in hypoxia (0.5% O2) for 24 hours. mRNA expression was normalized to β2M. Data are from four independent experiments and are shown as means ± SEM. Paired t-test was used to assess the significance. *p<0.05. ns, not significant.
Discussion
Hypoxia phenotyping in PDAC using gene signatures
Currently, PDAC remains one of the most treatment-refractory and clinically challenging cancers in clinical practice (2). A growing body of scientific literature now supports the idea that targeting tumor hypoxia in solid organ cancers can sensitize them to radiation, chemotherapy, and immunotherapies (15, 63, 64). Currently, hypoxic-tumor targeted therapies fall into four major categories: 1) Hypoxia-activated prodrugs that can kill tumor cells in hypoxic zones, 2) Molecules that inhibit HIF thereby abrogating HIF-downstream signaling, 3) Therapies that increase tumor oxygenation and 4) Drugs that modulate hypoxia-associated TME remodeling such as acidosis (e.g. CAIX) (15, 63, 65). Given the remarkable successes observed in certain tumor types (e.g. melanoma, renal cell carcinoma) using immune checkpoint blockade (9), identifying novel strategies for immune modulation of hypoxic pancreatic cancers is a critical priority. In this report, we utilized a well-established hypoxia signature to demonstrate that high hypoxia status is associated with distinct hallmarks of immunosuppression. Moreover, we utilized a machine-learning approach to identify LGALS3 (Gal-3) as a potential therapeutic target for abolishing hypoxia associated suppression in PDAC (Figures 6–8). Using this innovative approach for further exploration of immunological profiles of PDAC tumors with high hypoxia gene expression, we identified potential therapeutic targets to improve immunotherapy responses in hypoxic cancers. However, in vivo models will be essential to further validate these observations.
Multiple reports have reported the derivation of a hypoxia signature as reviewed by Abou Khouzam et al. earlier this year (15). It is pertinent to note that the same group reported the derivation of an 8-gene hypoxia signature that showed prognostic value in PDAC as well as an association with “immune cold” TME (24). Direct measurements of oxygen partial pressure in limited clinical samples have revealed that pancreatic and prostate cancers are highly hypoxic (66). However, in the absence of direct measurements, hypoxia signatures have been developed and reported in the literature as useful scoring tools to assess hypoxia status in bulk tumor transcriptomes (67). As such, in published reports on hypoxia, cases are either stratified by the expression levels of a hypoxia gene signature or by a hypoxia risk score, where the expression of hypoxia-linked genes (selected for the score through survival regression analyses) is multiplied by their regression coefficients (15). In the past few months, another report by Ren et al. utilized a compact 15 gene version of Buffa signature to derive a 6 gene risk score in PDAC (22). It is pertinent to note that all published hypoxia gene signatures and risk scores have been shown to be correlated with poor prognosis in PDAC (15, 22). In multiple published reports, patients were classified as high or low based on a median expression cut-off (22–24). In our study, we sought to compare cases with notable differences in hypoxia status, and therefore used a cut-off of top versus bottom quartiles of expression for our hypoxia score, an approach that is used to classify patients in certain gene expression profiling studies (68, 69). We selected the Buffa hypoxia signature, a widely reported signature that was recently used to identify hypoxia-related proteins involved in metastasis in a study of 17 different carcinomas (27). We further expanded upon previous reports using this signature, to demonstrate its utility for profiling PDAC tissues. Using an additional dataset of PDAC transcriptomic profiles, we showed that the PDAC tumor group exhibited significantly higher hypoxia scores relative to paired non-malignant pancreatic tissues (Supplementary Figure 2). We also demonstrated that key hypoxia-inducible genes were elevated in hypoxic PDAC cell lines relative to normoxic controls (Supplementary Figure 3). As a final point, we observed that hypoxia scores were significantly elevated in higher grade tumors (Figure 1B). This is corroborated by previous reports in the literature. Using a 9 gene hypoxia risk score, Zhuang et al. recently also demonstrated that higher histological tumor grades displayed increased hypoxia risk scores (45). Finally, we demonstrated for the first time using both univariate and multivariate survival analyses, that the Buffa hypoxia score is an important prognostic indicator in PDAC. Hypoxia scores and hypoxia risk scores reported by others also confirm this independent prognostic association with patient outcomes (24, 45), thereby confirming the importance of studying hypoxia-related gene expression in pancreatic cancer.
Deciphering the immune landscape in hypoxic PDAC cases
The mechanisms through which hypoxia shapes an immunosuppressive TME are complex, and reviewed elsewhere (64, 70, 71). In our study, we sought to dissect the distinct immunological mechanisms that could account for the “cold” tumor phenotype exhibited by HypoxiaHI cases. We identified multiple features suggestive of an immunosuppressive TME, that were associated with high hypoxia status. We observed that HypoxiaHI and HypoxiaLOW could be discriminated by both cell-type deconvolution using MCP-counter (Figure 3) and immunosuppressive gene expression profiling (Figures 4, 5). Using MCP-counter, we demonstrated multiple features representing an “immune cold” TME marked by a paucity of T and B lymphocytes, NK cells and myeloid DC (49). Notably, the inverse relationship between CD8+ T cells and hypoxia has also been reported in immunohistochemistry studies of colorectal, breast and ovarian cancers where CAIX is used as a marker for hypoxia (47, 72, 73). In fact, in a mouse melanoma model, treatment with a small molecule inhibitor of CAIX increased frequency of effector TH1 skewed CD4+ T cells and increased the response of these tumors to immune checkpoint therapy with anti-PD-1 and anti-CTLA-4 antibodies (74). The significant differences in abundance of CD8+ T cells between cases with high versus low hypoxia might also account for why we observed increased CD73 (NT5E) but reduced CD39 (ENTPD1) gene expression in HypoxiaHI cases. It is important to note that both ectoenzymes are known to be expressed on a wide range of cells in the TME (56). CD39 is a marker for T cell exhaustion and even though it was recently reported to be induced in terminally exhausted CD8+ T cells by tumor hypoxia (75), there are other features in hypoxic TMEs such as abnormal tumor vasculature and increased fibrosis, which mediate T cell exclusion (70). Thus, the exclusion of T cells from highly hypoxic tumors is one possible explanation for the reduced gene expression of CD39 observed in HypoxiaHI versus HypoxiaLOW cases. This was confirmed by our MCP-counter analyses (Figure 3), which revealed that HypoxiaHI cases have significantly lower abundance of total T cells, cytotoxic lymphocytes and CD8+ T cells.
According to previous findings, hypoxic regions produce signals that recruit multiple populations of myeloid cells such as MDSC, TAM and neutrophils (16, 64). We did not see any significant increase in the abundance scores for monocyte lineage cells or neutrophils (Figure 3). There is notable heterogeneity in the published literature on hypoxia gene signatures and their association with myeloid cell types (71). In the recent report by Abou Khouzam et al. (24), using an 8 gene hypoxia signature and the CIBERSORTx (76) deconvolution method in the TCGA PDAC cohort and an additional cohort, the authors found high hypoxia scores to be associated with diminished M2 macrophages in the TCGA dataset, and with an elevated abundance of M0, which purportedly accounts for gene expression signatures of undifferentiated macrophages that were not polarized into M1 or M2 in vitro (77). One further point that warrants consideration is that TAM and other myeloid cells are recruited into hypoxic niches, which also harbor necrotic zones (16, 78). Given that RNA from highly hypoxic and necrotic niches might undergo significant degradation (79, 80), bulk tumor transcriptomics data might not capture the myeloid diversity in highly hypoxic tumors. The absence of differences between major monocyte lineage cells between HypoxiaHI and HypoxiaLOW groups might also account for why we failed to see differences in expression of a number of amino-acid catabolizing enzymes which can be expressed in both non-immune cells and myeloid cells (59).
One notable phenotype of immunologically suppressed tumors that we observed was the interplay between NK cells, DC and COX-2, where we observed reduced signature scores for cDC1, lower MCP-counter abundance scores for myeloid DC and NK cells, as well as the increased gene expression of COX-2 (PTGS2) in HypoxiaHI vs HypoxiaLOW cases (Figures 3, 5 and Supplementary Figure 5). The induction of COX-2 expression in hypoxic conditions has been reported in both colon cancer and ovarian cancer cells (81, 82). Moreover, while COX-2 is reported to have tumor-intrinsic roles in the development and progression of PDAC in murine models (83), recent findings also demonstrated a key role for COX-2 in hindering anti-tumor immunity. The interplay between COX-2 and the NK cell-cDC1 axis was revealed in a landmark study using murine melanoma models (38). The authors showed that NK cells mediate the recruitment of the cDC1 subset, which are critical for anti-tumor immunity, and COX-2 mediated production of Prostaglandin E2 interferes with both NK cell chemokine production and cDC1 chemokine receptor expression. It is pertinent to note that COX-2 is over-expressed in multiple cancer types and influences the function of MDSC and lymphocytes as well as DC (57). An additional phenotype that we noted with MCP-counter, was a reduced abundance for endothelial cells in cases with high hypoxia scores. These results are supported by in vivo findings, where a recent study utilizing intravital fluorescence microscopy in a murine orthotopic pancreatic cancer model, revealed an inverse relationship between vascular density and the fraction of hypoxic cells (84). Thus, when cell-type deconvolution results with MCP-counter (Figure 3), support the premise that hypoxia mediates immune cell exclusion (70).
Immunomodulatory gene expression correlated with increased hypoxia in PDAC
We sought to further explore immunosuppressive signaling in hypoxic PDAC tissues by examining differentially expressed genes related to immune function. We saw significantly increased expression of a number of important genes that have a critical role in immune escape such as CD276 and CD47 (Figure 4A). These findings are corroborated by previous reports, as CD276 was reported to be associated with high hypoxia scores in the report by Abou Khouzam et al. (24), and it was shown in breast cancer cells, that HIF-1 directly induced CD47 transcription (85). However, we identified a novel association between Gal-4 and our hypoxia score (Figure 4B). Gal-4 was previously shown to be highly expressed at the protein level in PDAC tissue, and Gal-4 knockout murine tumors exhibited notable elevation of CD4+ and CD8+ T cell infiltration (55). Indeed, we demonstrate that all examined members of the galectin family in our analyses were expressed at higher levels in HypoxiaHI tumors relative to HypoxiaLOW. As reviewed recently, galectins are β-galactoside binding proteins that display important intracellular and extracellular functions and are secreted from cells through a non-canonical pathway (86). Galectins are purported to be involved in fibrotic and inflammatory diseases as well as cancer, and multiple galectin-targeted therapies are currently in clinical trials (54, 86). Our primary interest in profiling and comparing galectins was due to their roles in immunosuppression via direct T cell apoptosis but also as immune checkpoints (e.g. Gal-9 and TIM-3) (54, 86). Using a machine-learning recursive feature selection approach (61), we identified LGALS3 as the top gene that could classify a PDAC tumor as HypoxiaHI and confirmed the correlation between high hypoxia scores and LGALS3 expression in an additional cohort of PDAC RNA-seq data. In PDAC, studies have shown that Gal-3 induces inflammatory cytokine expression in pancreatic stellate cells (activated PSCs constitute a portion of the CAF population in PDAC) (87). We demonstrated induction of LGALS3 gene expression in hypoxic BxPC3 cells and a trend towards increased expression in PANC-1 cells albeit not reaching statistical significance (Figure 8). Our results are partially corroborated by an earlier study which showed LGALS3 expression in PANC-1 cells at 24h and 48h (88). Moreover, in the same year, Gonnermann et al. demonstrated that the protein expression of Gal-3 was significantly higher in BxPC3 and other PDAC cell lines compared to PANC-1 in normoxic conditions (89). This report also showed that galectin-3 could inhibit the proliferation of Gamma-delta T cells, which are of interest due to their tumor-killing ability (89). Thus, our observations on elevated expression of LGALS3 in PDAC tumors and PDAC cell lines, as well as the aforementioned report by Gonnerman et al. (89), highlight Gal-3 as a potential therapeutic target to reverse hypoxia-mediated immunosuppression.
Our study also had some limitations. First, the PAAD cohort in TCGA is comparatively small and therefore, comparing patients from the top and bottom quartile only yielded 44 cases per group. Third, while cell-type deconvolution approaches in bulk transcriptomics are more effective in generating abundance scores for cell types such as T cells, estimating the fractions of more heterogeneous cell types such as myeloid cells are challenging in silico (90). Finally, while multiple signatures for hypoxia have been reported in the literature (67, 71), these signatures warrant benchmarking in prospective studies where the hypoxia probe pimonidazole can be given to patients prior to operating (15), and RNA-seq can be performed to assess the correlation of each signature to levels of pimonidazole labeling in situ. However, we demonstrate multiple features of immunosuppression that delineate cases with high hypoxia score expression. Moreover, we described 5 genes associated with hypoxia high status, of which the top 3 (LGALS3, CD276 and NT5E) have already been identified as targets of interest in cancer (56, 86, 91). Moreover, the novel association between Gal-4 and high hypoxia status also warrants further investigation to reveal the biological role of this protein in PDAC progression and its potential relevance as a drug target. We anticipate that future studies using newly emerging spatial transcriptomics technologies will be essential to decrypt the gene expression profiles in hypoxic versus normoxic regions of the TME and in order to identify potential drug targets for hypoxic cancers (92).
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.
Ethics statement
This study utilized publicly available gene expression data from existing repositories and experiments performed in human cell lines. This study and all experiments performed within were reviewed and approved by the University Research Ethics Committee at Coventry University, UK.
Author contributions
HS: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. AA: Conceptualization, Formal analysis, Investigation, Methodology, Validation, Visualization, Writing – original draft, Writing – review & editing. HK: Formal analysis, Investigation, Methodology, Validation, Visualization, Writing – original draft, Writing – review & editing. TG: Data curation, Formal analysis, Investigation, Methodology, Writing – original draft, Writing – review & editing. RG: Conceptualization, Methodology, Validation, Visualization, Writing – original draft, Writing – review & editing. BB: Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing.
Funding
The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. HS and BB have received funding from the Research Excellence Development Fund at Coventry University. AA acknowledges support from the, HYPERMARKER (Grant agreement ID 101095480), and the MRC Health Data Research UK (HDRUK/CFC/01) and HDRUK midlands regional community project [QQ2], initiatives funded by UK Research and Innovation, Department of Health and Social Care (England) and the devolved administrations, and leading medical research charities. The views expressed in this publication are those of the authors and not necessarily those of the NHS, the National Institute for Health Research, the Medical Research Council or the Department of Health.
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/fimmu.2024.1360629/full#supplementary-material
References
1. Kamisawa T, Wood LD, Itoi T, Takaori K. Pancreatic cancer. Lancet. (2016) 388:73–85. doi: 10.1016/S0140-6736(16)00141-0
2. Jiang Y, Sohal DPS. Pancreatic adenocarcinoma management. JCO Oncol Pract. (2022) 19:19–32. doi: 10.1200/OP.22.00328
3. Andersson R, Haglund C, Seppänen H, Ansari D. Pancreatic cancer – the past, the present, and the future. Scand J Gastroenterol. (2022) 57:1169–77. doi: 10.1080/00365521.2022.2067786
4. Principe DR, Underwood PW, Korc M, Trevino JG, Munshi HG, Rana A. The current treatment paradigm for pancreatic ductal adenocarcinoma and barriers to therapeutic efficacy [Internet]. Front Oncol. (2021) 11. doi: 10.3389/fonc.2021.688377.
5. Partyka O, Pajewska M, Kwaśniewska D, Czerw A, Deptała A, Budzik M, et al. Overview of pancreatic cancer epidemiology in europe and recommendations for screening in high-risk populations. Cancers. (2023) 15(14), 3634. doi: 10.3390/cancers15143634.
6. Rahib L, Wehner MR, Matrisian LM, Nead KT. Estimated projection of US cancer incidence and death to 2040. JAMA Netw Open. (2021) 4:e214708–e214708. doi: 10.1001/jamanetworkopen.2021.4708
7. Beatty GL, Werba G, Lyssiotis CA, Simeone DM. The biological underpinnings of therapeutic resistance in pancreatic cancer. Genes Dev. (2021) 35:940–62. doi: 10.1101/gad.348523.121.
8. Orth M, Metzger P, Gerum S, Mayerle J, Schneider G, Belka C, et al. Pancreatic ductal adenocarcinoma: biological hallmarks, current status, and future perspectives of combined modality treatment approaches. Radiat Oncol. (2019) 14:141. doi: 10.1186/s13014-019-1345-6
9. Korman AJ, Garrett-Thomson SC, Lonberg N. The foundations of immune checkpoint blockade and the ipilimumab approval decennial. Nat Rev Drug Discovery. (2022) 21:509–28. doi: 10.1038/s41573-021-00345-8
10. Le DT, Durham JN, Smith KN, Wang H, Bartlett BR, Aulakh LK, et al. Mismatch repair deficiency predicts response of solid tumors to PD-1 blockade. Science. (2017) 357:409–13. doi: 10.1126/science.aan6733
11. Hu ZI, Shia J, Stadler ZK, Varghese AM, Capanu M, Salo-Mullen E, et al. Evaluating mismatch repair deficiency in pancreatic adenocarcinoma: challenges and recommendations. Clin Cancer Res. (2018) 24:1326–36. doi: 10.1158/1078-0432.CCR-17-3099.
12. Hu ZI, O’Reilly EM. Therapeutic developments in pancreatic cancer. Nat Rev Gastroenterol Hepatol. (2024) 21, 7–24. doi: 10.1038/s41575-023-00840-w.
13. Murakami T, Hiroshima Y, Matsuyama R, Homma Y, Hoffman RM, Endo I. Role of the tumor microenvironment in pancreatic cancer. Ann Gastroenterol Surg. (2019) 3:130–7. doi: 10.1002/ags3.12225.
14. Sherman MH, Beatty GL. Tumor microenvironment in pancreatic cancer pathogenesis and therapeutic resistance. Annu Rev Pathology: Mech Disease. (2023) 18:123–48. doi: 10.1146/annurev-pathmechdis-031621-024600.
15. Abou Khouzam R, Lehn JM, Mayr H, Clavien PA, Wallace MB, Ducreux M, et al. Hypoxia, a targetable culprit to counter pancreatic cancer resistance to therapy. Cancers. MDPI. (2023) 15(4), 1235. doi: 10.3390/cancers15041235.
16. Petrova V, Annicchiarico-Petruzzelli M, Melino G, Amelio I. The hypoxic tumor microenvironment. Oncogenesis. (2018) 7:10. doi: 10.1038/s41389-017-0011-9
17. Tao J, Yang G, Zhou W, Qiu J, Chen G, Luo W, et al. Targeting hypoxic tumor microenvironment in pancreatic cancer. J Hematol Oncol. (2021) 14:14. doi: 10.1186/s13045-020-01030-w
18. Kitada T, Seki S, Sakaguchi H, Sawada T, Hirakawa K, Wakasa K. Clinicopathological significance of hypoxia-inducible factor-1α expression in human pancreatic carcinoma. Histopathology. (2003) 43:550–5. doi: 10.1111/j.1365-2559.2003.01733.x
19. Garcia Garcia CJ, Huang Y, Fuentes NR, Turner MC, Monberg ME, Lin D, et al. Stromal HIF2 regulates immune suppression in the pancreatic cancer microenvironment. Gastroenterol. (2022) 162:2018–31. doi: 10.1053/j.gastro.2022.02.024
20. Zhou K, Cai C, Ding G, He Y, Hu D. A signature of six-hypoxia-related genes to evaluate the tumor immune microenvironment and predict prognosis in gastric cancer. BMC Med Genomics. (2022) 15:1–12. doi: 10.1186/s12920-022-01411-9.
21. Cao R, Ma B, Wang G, Xiong Y, Tian Y, Yuan L. Characterization of hypoxia response patterns identified prognosis and immunotherapy response in bladder cancer. Mol Ther Oncolytics. (2021) 22:277–93. doi: 10.1016/j.omto.2021.06.011.
22. Ren M, Zhang J, Zong R, Sun H. A novel pancreatic cancer hypoxia status related gene signature for prognosis and therapeutic responses. Mol Biotechnol. (2023), 1–20. doi: 10.1007/s12033-023-00807-x
23. Ding J, He X, Cheng X, Cao G, Chen B, Chen S, et al. A 4-gene-based hypoxia signature is associated with tumor immune microenvironment and predicts the prognosis of pancreatic cancer patients. World J Surg Oncol. (2021) 19:123. doi: 10.1186/s12957-021-02204-7
24. Abou Khouzam R, Rao SP, Venkatesh GH, Zeinelabdin NA, Buart S, Meylan M, et al. An eight-gene hypoxia signature predicts survival in pancreatic cancer and is associated with an immunosuppressed tumor microenvironment. Front Immunol. (2021) 12:680435. doi: 10.3389/fimmu.2021.680435.
25. Buffa FM, Harris AL, West CM, Miller CJ. Large meta-analysis of multiple cancers reveals a common, compact and highly prognostic hypoxia metagene. Br J Cancer. (2010) 102:428–35. doi: 10.1038/sj.bjc.6605450
26. Lane B, Khan MT, Choudhury A, Salem A, West CML. Development and validation of a hypoxia-associated signature for lung adenocarcinoma. Sci Rep. (2022) 12:1290. doi: 10.1038/s41598-022-05385-7.
27. López-Cortés A, Prathap L, Ortiz-Prado E, Kyriakidis NC, León Cáceres Á, Armendáriz-Castillo I, et al. The close interaction between hypoxia-related proteins and metastasis in pancarcinomas. Sci Rep. (2022) 12:11100. doi: 10.1038/s41598-022-15246-y
28. Liu J, Lichtenberg T, Hoadley KA, Poisson LM, Lazar AJ, Cherniack AD, et al. An integrated TCGA pan-cancer clinical data resource to drive high-quality survival outcome analytics. Cell. (2018) 173:400–416.e11. doi: 10.1016/j.cell.2018.02.052
29. Foroutan M, Bhuva DD, Lyu R, Horan K, Cursons J, Davis MJ. Single sample scoring of molecular phenotypes. BMC Bioinf. (2018) 19:404. doi: 10.1186/s12859-018-2435-4
30. Therneau TM, Grambsch PM. Modeling survival data: extending the cox model. New York: Springer (2000).
31. Therneau TM. A package for survival analysis (2023). Available online at: https://cran.r-project.org/package=survival.
32. Kassambara A, Kosinski M, Biecek P, Scheipl F. survminer: drawing survival curves using ‘ggplot2’. (2021).
33. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. (2010) 26:139–40. doi: 10.1093/bioinformatics/btp616.
34. Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. (2010) 11:R25. doi: 10.1186/gb-2010-11-3-r25
35. Benjamini Y, Hochberg Y. Controlling the false discovery rate: A practical and powerful approach to multiple testing. J R Stat Soc Ser B (Methodological). (1995) 57:289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x.
36. Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. (2016) 17:1–20. doi: 10.1186/s13059-016-1070-5.
37. Sturm G, Finotello F, List M. Immunedeconv: an R package for unified access to computational methods for estimating immune cell fractions from bulk RNA-sequencing data. Methods Mol Biol. (2020) 2120:223–32. doi: 10.1007/978-1-0716-0327-7_16
38. Böttcher JP, Bonavita E, Chakravarty P, Blees H, Cabeza-Cabrerizo M, Sammicheli S, et al. NK Cells Stimulate Recruitment of cDC1 into the Tumor Microenvironment Promoting Cancer Immune Control. Cell. (2018) 172:1022–1037.e14. doi: 10.1016/j.cell.2018.01.004
39. Pei H, Li L, Fridley BL, Jenkins GD, Kalari KR, Lingle W, et al. FKBP51 affects cancer cell response to chemotherapy by negatively regulating akt. Cancer Cell. (2009) 16:259–66. doi: 10.1016/j.ccr.2009.07.016
40. Bolstad BM, Collin F, Brettschneider J, Simpson K, Cope L, Irizarry RA, et al. Quality assessment of affymetrix geneChip data. Bioinf Comput Biol Solutions Using R Bioconductor. (2005), 33–47.
41. Díaz-Uriarte R, Alvarez de Andrés S. Gene selection and classification of microarray data using random forest. BMC Bioinf. (2006) 7:3. doi: 10.1186/1471-2105-7-3.
42. Acharjee A, Ament Z, West JA, Stanley E, Griffin JL. Integration of metabolomics, lipidomics and clinical data using a machine learning method. BMC Bioinf. (2016) 17:440. doi: 10.1186/s12859-016-1292-2.
43. Telloni SM. Tumor staging and grading: A primer BT - molecular profiling: methods and protocols. Espina V, editor. New York, NY: Springer New York (2017) p. 1–17. doi: 10.1007/978-1-4939-6990-6_1
44. Deer EL, González-Hernández J, Coursen JD, Shea JE, Ngatia J, Scaife CL, et al. Phenotype and genotype of pancreatic cancer cell lines. Pancreas. (2010) 39:425–35. doi: 10.1097/MPA.0b013e3181c15963.
45. Zhuang H, Wang S, Chen B, Zhang Z, Ma Z, Li Z, et al. Prognostic stratification based on HIF-1 signaling for evaluating hypoxic status and immune infiltration in pancreatic ductal adenocarcinomas. Front Immunol. (2021) 12. doi: 10.3389/fimmu.2021.790661.
46. Abou Khouzam R, Brodaczewska K, Filipiak A, Zeinelabdin NA, Buart S, Szczylik C, et al. Tumor hypoxia regulates immune escape/invasion: influence on angiogenesis and potential impact of hypoxic biomarkers on cancer therapies. Front Immunol. (2021) 11. doi: 10.3389/fimmu.2020.613114.
47. Craig SG, Humphries MP, Alderdice M, Bingham V, Richman SD, Loughrey MB, et al. Immune status is prognostic for poor survival in colorectal cancer patients and is associated with tumor hypoxia. Br J Cancer. (2020) 123:1280–8. doi: 10.1038/s41416-020-0985-5
48. Galon J, Bruni D. Approaches to treat immune hot, altered and cold tumors with combination immunotherapies. Nat Rev Drug Discovery. (2019) 18:197–218. doi: 10.1038/s41573-018-0007-y
49. Wang L, Geng H, Liu Y, Liu L, Chen Y, Wu F, et al. Hot and cold tumors: Immunological features and the therapeutic strategies. MedComm (Beijing). (2023) 4:e343. doi: 10.1002/mco2.343
50. Peterson EE, Barry KC. The natural killer–dendritic cell immune axis in anti-cancer immunity and immunotherapy. Front Immunol. (2021) 11. doi: 10.3389/fimmu.2020.621254.
51. Del Prete A, Salvi V, Soriani A, Laffranchi M, Sozio F, Bosisio D, et al. Dendritic cell subsets in cancer immunity and tumor antigen sensing. Cell Mol Immunol. (2023) 20:432–47. doi: 10.1038/s41423-023-00990-6
52. Guo Z, Zhang R, Yang AG, Zheng G. Diversity of immune checkpoints in cancer immunotherapy. Front Immunol. (2023) 14. doi: 10.3389/fimmu.2023.1121285.
53. Deng H, Wang G, Zhao S, Tao Y, Zhang Z, Yang J, et al. New hope for tumor immunotherapy: the macrophage-related “do not eat me” signaling pathway. Front Pharmacol. (2023) 14. doi: 10.3389/fphar.2023.1228962.
54. Chou FC, Chen HY, Kuo CC, Sytwu HK. Role of galectins in tumors and in clinical immunotherapy. Int J Mol Sci. (2018) 19(2), 430. doi: 10.3390/ijms19020430.
55. Lidström T, Cumming J, Gaur R, Frängsmyr L, Pateras IS, Mickert MJ, et al. Extracellular galectin 4 drives immune evasion and promotes T-cell apoptosis in pancreatic cancer. Cancer Immunol Res. (2023) 11:72–92. doi: 10.1158/2326-6066.CIR-21-1088
56. Xia C, Yin S, To KKW, Fu L. CD39/CD73/A2AR pathway and cancer immunotherapy. Mol Cancer. (2023) 22:44. doi: 10.1186/s12943-023-01733-x
57. Jin K, Qian C, Lin J, Liu B. Cyclooxygenase-2-Prostaglandin E2 pathway: A key player in tumor-associated immune cells. Front Oncol. (2023) 13. doi: 10.3389/fonc.2023.1099811.
58. Niu F, Yu Y, Li Z, Ren Y, Li Z, Ye Q, et al. Arginase: An emerging and promising therapeutic target for cancer treatment. Biomedicine Pharmacotherapy. (2022) 149:112840. doi: 10.1016/j.biopha.2022.112840.
59. Molinier-Frenkel V, Castellano F. Immunosuppressive enzymes in the tumor microenvironment. FEBS Lett. (2017) 591:3135–57. doi: 10.1002/1873-3468.12784
60. Castellano F, Correale J, Molinier-Frenkel V. Editorial: immunosuppressive amino acid catabolizing enzymes in heallth and disease. Front Immunol Switzerland;. (2021) 12:689864. doi: 10.3389/fimmu.2021.689864.
61. Diaz-Uriarte R. GeneSrF and varSelRF: a web-based tool and R package for gene selection and classification using random forest. BMC Bioinf. (2007) 8:328. doi: 10.1186/1471-2105-8-328.
62. Kirby MK, Ramaker RC, Gertz J, Davis NS, Johnston BE, Oliver PG, et al. RNA sequencing of pancreatic adenocarcinoma tumors yields novel expression patterns associated with long-term survival and reveals a role for ANGPTL4. Mol Oncol. (2016) 10:1169–82. doi: 10.1016/j.molonc.2016.05.004.
63. Semenza GL. Targeting intratumoral hypoxia to enhance anti-tumor immunity. Semin Cancer Biol. (2023) 96:5–10. doi: 10.1016/j.semcancer.2023.09.002
64. Jayaprakash P, Vignali PDA, Delgoffe GM, Curran MA. Hypoxia reduction sensitizes refractory cancers to immunotherapy. Annu Rev Med. (2022) 73:251–65. doi: 10.1146/annurev-med-060619-022830.
65. Singleton DC, Macann A, Wilson WR. Therapeutic targeting of the hypoxic tumor microenvironment. Nat Rev Clin Oncol. (2021) 18:751–72. doi: 10.1038/s41571-021-00539-4.
66. McKeown SR. Defining normoxia, physoxia and hypoxia in tumors-implications for treatment response. Br J Radiol. (2014) 87:20130676. doi: 10.1259/bjr.20130676.
67. Yang L, West CM. Hypoxia gene expression signatures as predictive biomarkers for personalising radiotherapy. Br J Radiol. (2019) 92:20180036. doi: 10.1259/bjr.20180036.
68. Bolen CR, Mattiello F, Herold M, Hiddemann W, Huet S, Klapper W, et al. Treatment dependence of prognostic gene expression signatures in de novo follicular lymphoma. Blood. (2021) 137:2704–7. doi: 10.1182/blood.2020008119
69. Cheng WC, Tsui YC, Ragusa S, Koelzer VH, Mina M, Franco F, et al. Uncoupling protein 2 reprograms the tumor microenvironment to support the anti-tumor immune cycle. Nat Immunol. (2019) 20:206–17. doi: 10.1038/s41590-018-0290-0
70. Pietrobon V, Marincola FM. Hypoxia and the phenomenon of immune exclusion. J Transl Med. (2021) 19:9. doi: 10.1186/s12967-020-02667-4
71. Abou Khouzam R, Zaarour RF, Brodaczewska K, Azakir B, Venkatesh GH, Thiery J, et al. The effect of hypoxia and hypoxia-associated pathways in the regulation of antitumor response: friends or foes? Front Immunol. (2022) 13. doi: 10.3389/fimmu.2022.828875.
72. Juhász P, Hasulyó D, Bedekovics J, Beke L, Kacsala N, Török M, et al. Carbonic anhydrase IX (CAIX) expressing hypoxic micro-environment hampers CD8+ Immune cell infiltrate in breast carcinoma. Appl Immunohistochemistry Mol Morphology. (2023) 31:26–32. doi: 10.1097/PAI.0000000000001082.
73. Guo N, Yang A, Farooq FB, Kalaria S, Moss E, DeVorkin L, et al. CD8 + T cell infiltration is associated with improved survival and negatively correlates with hypoxia in clear cell ovarian cancer. Sci Rep. (2023) 13:6530. doi: 10.1038/s41598-023-30655-3
74. Chafe SC, McDonald PC, Saberi S, Nemirovsky O, Venkateswaran G, Burugu S, et al. Targeting hypoxia-induced carbonic anhydrase IX enhances immune-checkpoint blockade locally and systemically. Cancer Immunol Res. (2019) 7:1064–78. doi: 10.1158/2326-6066.CIR-18-0657.
75. Vignali PDA, DePeaux K, Watson MJ, Ye C, Ford BR, Lontos K, et al. Hypoxia drives CD39-dependent suppressor function in exhausted T cells to limit antitumor immunity. Nat Immunol. (2023) 24:267–79. doi: 10.1038/s41590-022-01379-9
76. Newman AM, Steen CB, Liu CL, Gentles AJ, Chaudhuri AA, Scherer F, et al. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat Biotechnol. (2019) 37:773–82. doi: 10.1038/s41587-019-0114-2
77. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. (2015) 12:453–7. doi: 10.1038/nmeth.3337.
78. Bai R, Li Y, Jian L, Yang Y, Zhao L, Wei M. The hypoxia-driven crosstalk between tumor and tumor-associated macrophages: mechanisms and clinical treatment strategies. Mol Cancer. (2022) 21:177. doi: 10.1186/s12943-022-01645-2
79. Gallego Romero I, Pai AA, Tung J, Gilad Y. RNA-seq: impact of RNA degradation on transcript quantification. BMC Biol. (2014) 12:42. doi: 10.1186/1741-7007-12-42
80. Millier MJ, Stamp LK, Hessian PA. Digital-PCR for gene expression: impact from inherent tissue RNA degradation. Sci Rep. (2017) 7:17235. doi: 10.1038/s41598-017-17619-0.
81. Xue X, Shah YM. Hypoxia-inducible factor-2α is essential in activating the COX2/mPGES-1/PGE 2 signaling axis in colon cancer. Carcinogenesis. (2012) 34:163–9. doi: 10.1093/carcin/bgs313
82. Ding Y, Zhuang S, Li Y, Yu X, Lu M, Ding N. Hypoxia-induced HIF1α dependent COX2 promotes ovarian cancer progress. J Bioenerg Biomembr. (2021) 53:441–8. doi: 10.1007/s10863-021-09900-9
83. Hill R, Li Y, Tran LM, Dry S, Calvopina JH, Garcia A, et al. Cell intrinsic role of COX-2 in pancreatic cancer development. Mol Cancer Ther. (2012) 11:2127–37. doi: 10.1158/1535-7163.MCT-12-0342
84. Samuel T, Rapic S, O’Brien C, Edson M, Zhong Y, DaCosta RS. Quantitative intravital imaging for real-time monitoring of pancreatic tumor cell hypoxia and stroma in an orthotopic mouse model. Sci Adv. (2023) 9:eade8672. doi: 10.1126/sciadv.ade8672
85. Zhang H, Lu H, Xiang L, Bullen JW, Zhang C, Samanta D, et al. HIF-1 regulates CD47 expression in breast cancer cells to promote evasion of phagocytosis and maintenance of cancer stem cells. Proc Natl Acad Sci U S A. (2015) 112:E6215–23. doi: 10.1073/pnas.1520032112.
86. Mariño KV, Cagnoni AJ, Croci DO, Rabinovich GA. Targeting galectin-driven regulatory circuits in cancer and fibrosis. Nat Rev Drug Discovery. (2023) 22:295–316. doi: 10.1038/s41573-023-00636-2.
87. Zhao W, Ajani JA, Sushovan G, Ochi N, Hwang R, Hafley M, et al. Galectin-3 mediates tumor cell–stroma interactions by activating pancreatic stellate cells to produce cytokines via integrin signaling. Gastroenterology. (2018) 154:1524–37. doi: 10.1053/j.gastro.2017.12.014.
88. da Silva Filho AF, Tavares LB, Pitta MGR, Beltrão EIC, Rêgo MJBM. Galectin-3 is modulated in pancreatic cancer cells under hypoxia and nutrient deprivation. Biol Chem. (2020) 401:1153–65. doi: 10.1515/hsz-2019-0413.
89. Gonnermann D, Oberg HH, Lettau M, Peipp M, Bauerschlag D, Sebens S, et al. Galectin-3 released by pancreatic ductal adenocarcinoma suppresses γδ T cell proliferation but not their cytotoxicity. Front Immunol. (2020) 11:1328. doi: 10.3389/fimmu.2020.01328
90. Finotello F, Rieder D, Hackl H, Trajanoski Z. Next-generation computational tools for interrogating cancer immunity. Nat Rev Genet. (2019) 20:724–46. doi: 10.1038/s41576-019-0166-7
91. Getu AA, Tigabu A, Zhou M, Lu J, Fodstad Ø, Tan M. New frontiers in immune checkpoint B7-H3 (CD276) research and drug development. Mol Cancer. (2023) 22:43. doi: 10.1186/s12943-023-01751-9
Keywords: hypoxia, tumor microenvironment (TME), pancreatic ductal adenocarcinoma (PDAC), immune checkpoint, galectins
Citation: Sadozai H, Acharjee A, Kayani HZ, Gruber T, Gorczynski RM and Burke B (2024) High hypoxia status in pancreatic cancer is associated with multiple hallmarks of an immunosuppressive tumor microenvironment. Front. Immunol. 15:1360629. doi: 10.3389/fimmu.2024.1360629
Received: 23 December 2023; Accepted: 12 February 2024;
Published: 06 March 2024.
Edited by:
Dmitry Aleksandrovich Zinovkin, Gomel State Medical University, BelarusReviewed by:
Eldar Nadyrov, Gomel State Medical University, BelarusJale Yuzugulen, Eastern Mediterranean University, Türkiye
Copyright © 2024 Sadozai, Acharjee, Kayani, Gruber, Gorczynski and Burke. 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: Hassan Sadozai, sadozaih@uni.coventry.ac.uk; Animesh Acharjee, a.acharjee@bham.ac.uk; Bernard Burke, ac2561@coventry.ac.uk
†These authors have contributed equally to this work