Transcriptional response to hypoxic stress in melanoma and prognostic potential of GBE1 and BNIP3

INTRODUCTION

Cancer immunotherapy has recently emerged as an important treatment modality. FDA approval of sipuleucel-T (Provenge^*), ipilimumab (Yervoy^*), nivolumab (Opdivo^*) and pembrolizumab (Keytruda^*) has started to deliver on the long awaited promise of cancer immunotherapy [1]. Many new modalities of immunotherapies targeting cytotoxic T lymphocytes (CTLs) responses, such as adoptive cell therapies and vaccines, are in advanced clinical trials. The ultimate goal of most cancer immunotherapy strategies is to induce a strong cytotoxic T lymphocyte (CTL) response [2]. The prevailing view is that induced or boosted CTLs will eradicate tumor cells. However, this view has been seriously challenged by clinical observations [3]. Immunotherapy effectiveness is dependent on the qualitative and/or quantitative features of the killer cells and the complexity of the genomic aberrations harbored by neoplastic cells, but is also regulated by numerous dynamic properties of the tumor microenvironment [4]. Besides tumor cells, the tumor microenvironment harbors a variety of host-derived cells. It is a complex system playing an important role in tumor development and progression [5]. It involves soluble factors and metabolic changes. Among the metabolic changes, hypoxia plays a key role in sculpting tumor microenvironment [6]. Hypoxia arises due to a combination of excessive oxygen consumption by growing tumor cells and the disorganized tumour-associated vasculature [7]. Accumulating evidence indicates that hypoxia plays an important role in tumor progression, affecting both metastatic spread and selection of cells with more aggressive phenotypes [8]. It is well established that hypoxic stress is a feature of most solid tumors and is associated with poor prognosis in several cancer types [9, 10]. In the context of tumor microenvironment, tumors impose several limitations to dampen T cell immunity as T cells, experiencing the metabolic framework of growing tumors, fail to activate distinct pathways to accomplish their functional requirements. Tumor microenvironmental hypoxia is in this regard a relevant example demonstrating how the tumor microenvironment can paralyze and neutralize T cell functions [11]. It is a negative prognostic and predictive factor owing to its multiple contributions to chemoresistance, radioresistance, angiogenesis, resistance to cell death, altered metabolism and genomic instability [12].

The master regulator of the hypoxic response is the Hypoxia-inducible factor 1 (HIF-1). Several reports have identified links between cancer outcomes and the level of HIF-1α protein [6, 12]. While incipient angiogenesis in small tumors may occur independently of hypoxia, growing tumors will at some point inevitably experience inadequate nutrient and oxygen supply [12, 13]. This deprivation triggers an angiogenic switch, which is associated with reduced sensitivity to cytotoxic and genotoxic treatment and more aggressive metastatic behavior [13]. Therefore, precise knowledge of the hypoxic state of a tumor not only provides a valuable entry point to understanding tumor progression, but also holds considerable prognostic value.

Despite the ubiquity and importance of HIF-1 response and our knowledge about the variation in the global transcriptional response to hypoxia among different cell types, little is known about gene expression signatures that might relate to melanoma development and response to treatments. It is becoming clear that HIF-1α expression alone is not a reliable marker of tumor response to hypoxia. Although treatment efficacy has been improved for patients with melanoma using checkpoints inhibitors [14], the overall 5-year survival rate is still about 50%. In fact, tumors respond very heterogeneously to this treatment and biomarkers are needed. Here, we have focused on determining transcriptional response to hypoxic stress in melanoma. RNA microarray analysis was performed to examine the effect of hypoxia on gene expression. We identified a gene expression signature, and discuss the putative prognostic and predictive potential of BNIP3 and GBE1 genes in the clinical outcome of melanoma patients treated with anti-PD1 (pembrolizumab).

RESULTS

Transcriptional changes associated with hypoxia in human primary melanoma cell lines

We analyzed the changes in global transcript levels in response to hypoxic stress. For this purpose, we used DNA microarrays to examine the gene expression program in response to hypoxia (1% O2) in different melanoma cell lines established from melanoma patients. We profiled global mRNA levels at the 24h time point selected on the basis of our previous finding [15]. The different mRNA samples were analyzed by hybridization to DNA microarrays. Analysis of gene expression profiles of human primary cell lines of melanoma cultured at 1% oxygen vs 21% oxygen, allowed to establish a signature of 35 genes (Figure 1A and Table 1). 26 genes up-regulated (fold-change ≥ 2.5) and 9 genes down-regulated under hypoxic conditions (fold change ≤-2) and an adjusted p-value (FDR) < 0.05 (Figure 1B). Clustering analysis clearly showed seperation of the samples into two groups according to their hypoxic status (Figure 1A). To investigate putative interactions between these genes, we then used STRING (Search Tool for Recurring Instances of Neighbouring Genes). The analysis revealed the existence of a functional interaction between 15 genes of the hypomel signature (13 genes overexpressed: BNIP3, AK4, SLC2A1, ADM, PFKFB4, ENO2, VEGFA, DDIT4, PGK1, GBE1, ALDOC, CCL28 and 2 genes underexpressed: EPRS and DDX21) (Figure 1C). A further examination of the hypomel signature using Gene Set Enrichment Analysis (GSEA) showed that several pathways were involved including hypoxia, metabolic and glucose catabolic processes (Figure 1D).

We next performed RT-qPCR analysis and confirmed the expression levels of hypoxia-associated genes identified by microarray assay (Figure 1E). An analysis of the correlation between expression level fold changes derived from microarray and RT-qPCR experiments is depicted in Figure 1F and revealed a significant correlation (Pearson’s r = 0.76; p<0.0001). Thus, the microarray data highly corroborate with those obtained by RT-qPCR for the 35 gene expression in each of the 10 tested human cell lines. Data depicted in Supplementary Figure 1 indicate that the genes associated with hypoxic stress are specifically induced in melanoma and not in peripheral blood mononuclear cells (PBMC) cultured under hypoxic conditions.

Comparative analysis of hypoxia-associated gene expression in primary and metastatic melanoma

We next examined the gene expression, by RT-qPCR, in response to hypoxic stress on a second set of melanoma cell lines freshly established from primary and metastatic tumors derived from 3 melanoma patients (Figure 2A, 2B and 2C). Further analysis indicates a correlation between the fold change expression levels obtained by microarray in the previous assay and the RT-qPCR analysis (Figure 2D) (Pearson’s r = 0.61; p<0.0001). Furthermore, comparative analysis of gene expression levels, as assessed by RT-qPCR, in primary vs metastatic cell lines indicates a strong correlation of gene expression between the two cell types (Figure 2E) (Pearson’s r = 0.84; p<0.0001). Using the same experimental system, Western blot analysis on two highly expressed genes (ANGPTL4 and BNIP3) confirmed their induced expression concomitant with HIF-1α induction in the cells when cultured under hypoxic conditions at 24h and 48h (Figure 2F, 2G and 2H). These results, obtained from materials freshly isolated from patient tumors, substantiate the deregulation of genes of the hypomel signature in melanoma cells. It further indicates that primary and their metastatic counterparts will respond similarly, with respect to these genes under hypoxic stress.

Hypoxia-associated gene expression in hypoxic zones detected by anti-HIF-1α immunohistochemistry (IHC) in 8 melanoma patient tissues

We then investigated tissue specimens from patients with melanoma. Using IHC and HIF-1α on FFPE serial sections, we could identify (Figure 3) non hypoxic (A1) and hypoxic zones (A2, A3, A4, A5) in primary melanoma from 4 patient tissues. A1 and A2 are from the same primary melanoma. In this survey, It should be noted that HIF-1α positivity was considered specific only when the staining was found in cancer cell nuclei (arrows) while cytoplasmic staining was considered as non specific. Analysis of metastatic tissues illustrates hypoxic (A7 and A9) and non hypoxic (A6 and A8) zones, in 2 metastatic lymph nodes and 2 cutaneous metastases. As depicted in Figure 3, the number of HIF-1α positive cells as well as the intensity of staining can vary in the tumor zone. When the HIF-1α positive or negative zones had a limited surface, laser dissection of hypoxic or non hypoxic zone of the patient specimen was required for subsequent transcriptional analysis. To this end, we dissected the hypoxic and non hypoxic areas in melanoma samples from 8 patients who developed primary melanoma, metastatic lymph node, or cutaneous metastasis. RNAs were extracted then amplified by RT-qPCR for the 26 up-regulated genes of interest in tumor zones delimited after HIF-1α immunostaining (Figure 4A, 4B and 4C) and compared to microarray fold changes obtained in the hypomel signature derived from cell lines.

In Figure 4A, RT-qPCR data obtained were comparable for the 4 hypoxic zones of primary melanoma when normalized to a non hypoxic zone. We noted that the gene expression fold changes measured by qRT-PCR for the zones of interest (melanoma primary, metastatic lymph node and cutaneous metastatic) were similar in most instances, but often slightly reduced as compared with those obtained in the hypomel signature derived from primary cell lines.

This difference could be explained in part by the fact that in the different hypoxic zones, only a fraction of cells were actually imposed to hypoxic conditions (as shown by HIF-1α immunostaining), whereas in our cell line-based assays, virtually all the cells were exposed to hypoxic conditions. Moreover, by comparing HIF-1α immunostainings in patients with primary melanoma (Figure 5A and 5C) or metastatic lymph nodes (Figure 5B), we observed that tumor hypoxic zones with the strongest staining of HIF-1α were associated with a higher expression of most of the 26 up-regulated genes of the hypomel signature, suggesting a link between hypoxia, HIF-1 expression enrichment of the hypomel signature.

The recent application of the NanoString as a reliable gene expressiosn analysis prompted us to test whether the expression of our selected genes was correlated to treatment response in 19 melanoma patients treated with anti-PD1 (Supplementary Table 1). The NanoString approach offers a valuable alternative to RT-qPCR, is more accurate, and needs less materials because of direct quantification of gene copy number without the need for enzymatic amplification. We first validated this approach using 8 samples from patient tissues and cell lines with known hypoxic status. As shown in Supplementary Figure 2, a correlation exists between Nanostring data and the data previously obtained by microarray and RT-qPCR on cell lines (Pearson’s r = 0.506 ; p=0.00835), reinforcing the robustness of the NanoString analysis.

We next applied this approach for the 26 up-regulated genes within the hypomel signature using the 19 melanoma patient samples treated with anti-PD1 (9 responders (R) and 10 non responders (NR)). Finally, we investigated whether some genes within the signature could have a predictive value for the clinical outcome.

To compensate for processes unrelated directly to hypoxia, we used the “biomarker pair” approach under the rationale that such pairs would normally share common biological properties yielding similar expression levels, however could differ in their response to treatment [16]. We identified 18 pairs of highly correlated pairs (correlation level > 0.80) based on the raw expression levels of the 26 studied genes. Then, we applied Nonnegative matrix factorization analysis (NMF) to the matrix consisting of samples in (rows) and expression ratios between each selected pair of genes as variable. NMF seeks matrix factors containing only nonnegative elements. The resulting factorization often leads to substantial improvements in interpretability of the factors and in prediction of the outcome [17]. Prototypical expression patterns are estimated, in effect yielding a bi-ordering of samples and genes within clusters. In case of association with responder status, a majority of samples with a given status should be grouped together on top of a given cluster. In order to assess the significance of such grouping, a permutation test is performed, as described further in the methods section. The NMF-based ordering of samples appeared significantly associated with the responder status (p = 0.041 based on 100000 permutations). The dual ordering of genes pointed to the BNIP3/GBE1 differential pair, which had the largest leverage on the second cluster. Responder status appeared significantly associated to BNIP3/GBE1 differential expression level (p=0.0138) following a one-way analysis of the differential pair (Figure 6A and 6B).

DISCUSSION

Tumor hypoxia is one of the most important features of the tumor microenvironment, exerting an adverse effect on tumor aggressiveness and patient prognosis. In the course of these studies, we have identified a gene-expression profile of the hypoxia response in freshly established melanoma cell lines. We established a signature (called hypomel for HYPOxia MELanoma) of 35 genes: 26 genes up-regulated (fold-change ≥ 2.5) and 9 genes down-regulated under hypoxic conditions (fold change ≤ -2). The hypomel signature was validated in tumors FFPE from patients who developed primary melanoma, metastatic lymph node, or cutaneous metastasis. We have analyzed transcriptional expression by RT-qPCR of the hypomel signature in hypoxic zones delimited after IHC staining with anti-HIF-1α. RNAs were extracted from the hypoxic zones which were levied by microdissection laser. We found that a high staining with anti-HIF-1α in tumor hypoxic zone associated with a high expression of genes belonging to the hypomel signature. Based on the evidence currently available, it now appears that both the adaptive and innate immune systems can recognize and eliminate tumors. The problem we face, however, is that the tumor microenvironment is able to neutralize and paralyze both responses. One challenge for tumor immunologists in the future, is identifying patients for who will respond favorably to immunotherapy. Accordingly, in the course of these studies, we attempted to examine whether this identified signature correlates with clinical outcome of melanoma patients treated with anti-PD1 monoclonal antibody pembrolizumab and explored the quantitative NanoString technology to investigate whether some genes within the signature could have a predictive value for response following pembrolizumab treatment of patients with advanced melanoma. More importantly, we looked for pairs of highly correlated genes based on their raw expression levels (correlation level > 0.80). Non-Negative Matrix Factorization, NMF, was applied to the split matrix [17], yielding a dual clustering of samples and genes into two clusters. Samples and genes were ordered by decreasing leverage on their respective cluster. Finally, the association between the ordering of samples and responder status was assessed through a permutation test. Specifically, responder status was permuted among patients, the association score re-calculated and compared to the original score.

We identified 18 pairs of highly correlated pairs based on the raw expression levels of the 32 studied genes corresponding to 26 up-regulated genes and 6 housekeeping genes. The dual ordering of genes pointed to the BNIP3/GBE1 differential pair, which had largest leverage. Responder status appeared significant following a one-way analysis of the differential pair. We could demonstrate that the levels of BNIP3 and GBE1 correlate with the clinical response in melanoma patients treated with anti-PD1.

Several studies have shown that autophagy constitutes a potential target for cancer therapy and that the induction of autophagy in response to therapeutics can be viewed as having a prodeath or a prosurvival role [18], contributing to drug resistance. Our very recent studies highlighted a new hypoxia-induced pathway in which NANOG activates BNIP3L expression, contributing to autophagy induction in hypoxic tumor cells and their resistance to killing by CTL [19] further suggesting a link between hypoxia-induced resistance and autophagy-related stemness. BNIP3 is a member of the BCL-2 family of proteins with reported pro-death as well as pro-autophagic and cytoprotective functions, depending on the type of stress and cellular context [20]. In line with this, increased BNIP3 levels in melanoma patients appear to be linked with poor prognosis [20]. The induction of autophagy in response to metabolic and therapeutic stresses can have a prodeath or a prosurvival role. It should be noted that autophagy can contribute to fulfilling acute metabolic needs under starvation conditions by degrading and recycling the cargos [18]. Moreover, accumulating evidence indicates that among the various metabolic adaptations used by cancer cells to adjust to the conditions imposed by the tumor microenvironment, changes in glycogen metabolism are emerging as an essential response [21]. Of note, the induction of the melanogenic pathway may lead to robust upregulation of HIF-1-dependent and independent pathways in cultured melanoma cells, suggesting a key role for melanogenesis in regulation of cellular metabolism [22, 23]. In addition several reports indicate that melanogenesis can affect disease or therapeutic outcomes [24-26], metabolism [22, 27] or immune functions [28].

Hypoxia is known to induce adaptive changes in cell metabolism that include a switch from oxidative phosphorylation to glycolysis and increased glycogen synthesis. In this context, glycogen provides a convenient glucose reservoir during energy stress, glucose deprivation or senescence [29]. In addition, cancer genomics data indicate that elevated levels of the glycogenic enzyme GBE1 are associated with poor survival in AML [30]. Consistent with these reports, our studies point to a role of elevated glycogenic flux that correlates with a poor clinical response in melanoma treated with checkpoint inhibitor anti-PD1. These studies further suggest that channeling of glucose through glycogen may promote the survival of melanoma cells under hypoxia. Through elevated glycogenic flux and induction of autophagy, hypoxia appears to be a potential critical molecular program that could be considered as a prognostic factor for melanoma.

Nordsmark and colleagues found that pretreatment tumor oxygenation is a highly significant prognostic factor for survival after primary radiotherapy applied alone or combined with chemotherapy, surgery, or radiation sensitizers in patients with locally advanced HNSCC [31]. Increased HIF-1α expression has been reported by Mouriaux F et al. to correlate with cell proliferation and vascular markers CD31 and VEGF-A in uveal melanoma [32]. However, there was no correlation found between high HIF-1α expression and patient survival However very recently, the association of HIF-1α expression with clinicopathological characteristics and overall survival (OS) of patients with OSCC was evaluated by Zhou J et al. who found an association of HIF-1α overexpression with tumor size, tumor stage, lymph node metastasis, and overall survival [33].

The current studies establish a link between hypoxia, autophagy, glucose metabolism and clinical outcome. Indeed, much remains to be learned to further dissect this relationship. Next experiments will provide more insights into the role of hypoxic stress in shaping the anti-tumor response and its impact in cancer immunotherapy.

MATERIALS AND METHODS

Clinical samples

Clinical samples were collected from 19 patients, 12 females and 7 males who were treated with anti-PD1 (pembrolizumab) for a primary or metastatic melanoma. All patients had given their written informed consent in accordance with the declaration of Helsinki to participate. The main characteristics of the patients and their clinical responses according to the classical RECIST1.1 evaluation criteria are reported in the Supplementary Table 1.

Human tumor melanoma cell lines and peripheral blood mononuclear cell (PBMC)

Human melanoma cell lines ME204 AI/ER, ME260 LN/DG, T921 DUF/ALE, ME300PB, ME290mH and NA8 derived from the primary lesion, were provided by Dr Pedro Romero (Ludwig Center for Cancer Research, Lausanne, Switzerland). M74 derived from the primary lesion was established by Pr Jotereau (CRCNA, Inserm UMR892, Nantes, France). RIOUP2 was derived from the primary lesion, by the laboratory. The 3 pairs of human melanoma cell lines : T1 and G1, M4T and M4T2, RAYI2 and M1 were derived from the primary lesion and the metastatic lymph node, respectively, of 3 patients in the laboratory. All the melanoma cells were cultured in RPMI 1640 with glutamax supplemented with 10% FCS, 1% penicillin-streptomycin, 1% sodium pyruvate at 37°C in a humidified atmosphere containing 5% CO2. All culture reagents were from ThermoFisher Scientific (Waltham, MA, USA). PBMC were generated from one healthy donor and cultured in RPMI 1640 with glutamax supplemented with IL-2 (25 U/ml; Roussel-Uclaf, Romainville, France), 10% Human AB serum (Institut Jacques Boy, Reims, France), 1% penicillin-streptomycin, 1% sodium pyruvate at 37°C in a humidified atmosphere containing 5% CO2.

Hypoxic conditioning of tumor cells

Hypoxic treatment was conducted in a hypoxia workstation (Invivo2 400, Ruskinn, UK) in a humidified atmosphere containing 5% Co2, 1% O2 and 94% N2 at 37°C (24h and 48h). Melanoma cells for RNA and protein analysis were lysed directly in the hypoxia workstation without reoxygenation.

Microarray assay

Gene expression analysis were performed with Agilent^® SurePrint G3 Human GE 8x60K Microarray (Agilent Technologies, Santa Clara, CA, USA) with the following dual-color design: the test samples (Hypoxic samples) were labeled with Cy5 whereas the control samples (normoxic samples) were labeled in Cy3 using the two-color Agilent labeling kit (Low Input Quick Amp Labeling Kit 5190-2306) adapted for small amount of total RNA (100 ng total RNA per reaction). Hybridization was then performed following the manufacturer instructions. Microarray images were analysed by using Feature Extraction software version (10.7.3.1) from Agilent technologies. Defaults settings were used.

Microarray data processing and analysis

Raw data files from Feature Extraction were imported into R with LIMMA [34], an R package from the Bioconductor project, and processed as follows: gMedianSignal and rMedianSignal data were imported, controls probes were systematically removed, and flagged probes (gIsSaturated, gIsFeatpopnOL, gIsFeatNonUnifOL, rIsSaturated, rIsFeatpopnOL, rIsFeatNonUnifOL) were set to NA. Intra-array normalization was performed by a loess normalization, followed by a quantile normalization of both Cy3 and Cy5 channel. Then inter-array normalization was performed by quantile normalization on M values. To get a single value for each transcript, taking the mean of each replicated probes summarized data. Missing values were inferred using KNN algorithm from the package ‘impute’ from R bioconductor. Normalized data were then analyzed. To assess differentially expressed genes between two groups, we start by fitting a linear model to the data. Then we used an empirical Bayes method to moderate the standard errors of the estimated log-fold changes. The top-ranked genes were selected with the following criteria: an absolute fold-change ≥ 2.5 and ≤-2 and an adjusted p-value (FDR) < 0.005. To interogate interactions between hypomel genes and pathways, we used STRING (Search Tool for Recurring Instances of Neighbouring Genes). Gene Set Enrichment Analysis (GSEA) used the Molecular Signatures Database (MSigDB).

RNA isolation and real-time quantitative polymerase chain reaction (RT-qPCR)

Total RNAs were extracted from cell samples using TRIzol solution (Invitrogen). The quality of RNAs was assessed using a Bioanalyzer instrument (Agilent) and then quantified using a Biospecnano (Shimadzu, Kyoto, Japan). cDNA synthesis was prepared from 1 μg of total RNA with random hexamers using Applied Biosytems Reverse Transcription kit according to the supplied protocols. Gene expression was quantified by SYBR Green qPCR method using the Maxima™ SYBR Green/ ROX qPCR Master Mix on an StepOnePlus Real Time PCR system (ThermoFisher Scientific). Relative expression was calculated by using the comparative Ct method (2-ΔΔCt). Primer sequences for the quantification of 35 genes were purchased from Sigma and are available upon request. Transcript levels of HPRT for PBMC, or 18S for melanoma cells lines were used as endogenous control.

Western blot

Melanoma cells lines from 3 patients were grown in two different conditions normoxia (21% PO2) and hypoxia (1% PO2) for 24h and 48h at 37°C. Cancer cells were washed twice in phosphate-buffered saline and lysed in plates with lysis buffer (62.5 mM Tris-HCl [pH 6.8], 2% weight/volume sodium dodecyl sulfate, 10% glycerol, 1 mM orthovanadate, 2 mM phenylmethylsulfonyl fluoride, 25 μM leupeptin, 5 mM benzamidine, 1μM pepstatin, and 25 μM aprotinin). Lysates were sonicated on ice, resolved by sodium dodecyl sulfate-polyacrylamide gel electrophoresis (30μg/well), and transferred onto nitrocellulose membranes. The membranes were blocked in blocking buffer then probed overnight at 4°C with the indicated primary Abs. Primary antibodies (Abs) against HIF-1α (mouse Ab, Clone 54/HIF-1α 610959), BNIP3 (mouse Ab, clone ANa40, 10433), ANGPTL4 (Rabbit polyclonal Ab, SAB1410901) and β-Actin (mouse Ab, clone AC-15) were purchased respectively from BD Biosciences (San Jose, CA, USA), Abcam (Cambridge, UK) and Sigma-Aldrich (St Louis, MO, USA). The labeling was performed following incubation with horseradish peroxidase (HRP)-conjugated secondary Abs: a HRP goat anti-mouse for HIF-1α and BNIP3, and a HRP goat anti-rabbit for ANGPTL4, and detection with an enhanced chemiluminescence kit (GE Healthcare, Chicago, IL, USA). Blots were scanned and processed by Adobe Photoshop 7.0 software.

Immunohistochemistry staining for HIF-1α

The tissue collection was composed of 4 human melanoma primary tumors, 2 human cutaneous metastasis tumors and 2 human melanoma lymph node metastases.

For each patient, four micrometer sections of Fixed-formaldehyde paraffin embedded (FFPE) melanoma human tumors were prepared and stained with HES (Hematoxylin Eosin Safran). Deparaffinized tissue sections were treated with Antigen Retrieval Solution (citrate buffer, pH 8.0, concentrated 10×, T0010 (Diapath, Martinengo, Italy) in water bath at 95°C. Tissue sections were then incubated with H2O2 3% for 10 min and solution PowerVision IHC/ISH Super Blocking PV6122 (MM France, Brignais, France) for 20 min. Histological slides were incubated over night at 4°C with a polyclonal rabbit anti-human HIF-1α antibody NB100-479 (Novus Biologicals, Littleton, CO, USA) or a rabbit polyclonal IgG Ab 27472 (Abcam). For signal amplification, slides were then incubated with rabbit alkaline phosphatase conjugated secondary antibody PowerVision poly-AP anti-Rabbit IgG PV3133 (MM France). The signal was revealed with the Liquid Permanent Red K0640 (Dako, Les Ulis, France) and Mayer’s hemalun solution counterstain (Merck Millipore, Billerica, MA, USA). In this analysis, was considered as positive nuclear immunostaining of HIF-1α in cancer cells, whereas cytoplasmic immunostaining was considered as non-specific.

Laser microdissection and pressure catapulting of glomerules

Laser microdissection was performed with a PALM^® RoboSoftware 4.6 MicroBeam system (PALM Microlaser Technologies, Zeiss Micro-Imaging, Munich, Germany) coupled to an inverted microscope Axio Observer.Z1. Serial 20 μm-thick sections from 3 human melanoma FFPE fragments (primary melanoma 1, 2 and 4, Figure 3) were spread onto polyethylene naphthalate (PEN) membrane-coated slides (Carl Zeiss Micro Imaging, Munich, Germany). After sectioning, the slides were incubated for 10 min in a toluene solution followed by 10 min in absolute alcohol to completely remove the paraffin embedding. Staining for 1 min in a Mayer’s Hematoxylin solution, 30 s in a saturated lithium carbonate solution and 30 s in a solution of erythrosine were successively done followed by 30 s each of absolute alcohol and toluene, respectively.

A total surface of 4 areas (between 6.6 × 106 μm² and 18 × 106 μm²) was collected from several sections representing the 3 identified hypoxic tumoral areas (zone A2, A3 and A5, Figure 3) and 1 non hypoxic tumoral area (zone A1, Figure 3). The total of surface was determined by the computer during microdissection process. All the surfaces were computed and the total surface was determined accordingly. Each sample was recovered in lysis buffer from the RNeasy FFPE kit (Qiagen, Hilden, Germany) for further molecular analyses.

Macrodissection with scalpel

Serial 20 μm-thick sections from 5 human melanoma FFPE tumors (primary melanoma 4, metastatic lymph node 1 and 2, cutaneous metastasis 1 and 2, Figure 3) were prepared as samples obtained by laser microdissection. A total of 5 areas was collected with a scalpel from several sections representing the 3 identified hypoxic tumoral areas (zone A4, A7, and A9, Figure 3) and 2 non hypoxic tumoral areas (zone A6 and A8, Figure 3). Each sample was recovered in lysis buffer from the RNeasy FFPE kit (Qiagen,) for further molecular analyses.

Pathway specific gene expression profiling of hypoxic and non hypoxic tumoral FFPE samples

Total RNAs were extracted from 9 area melanoma FFPE tissue sections obtained by laser microdissection or macrodissection with scalpel (zone A1 to A9, Figure 3), using RNeasy FFPE kit (Qiagen, Hilden, Germany). RNA quantity and quality was assessed using the Nanodrop-ND-1000 (Nanodrop Technologies, Wilmington, USA). First-strand cDNA was synthesized using a High-Capacity cDNA Reverse Transcription Kit (ThermoFisher Scientific) according to the manufacturer’s protocol. Signaling pathways transcript analyses were conducted in duplicates using a personalised Human qPCR SignArrays^® 384 system (gene profiling analysis Human qPCR SignArrays^® 384 kit for 26 genes of interest; and Perfect MasterMix SYBR Green (AnyGenes, Paris, France)) on a LightCycler 480 (Roche, Rosny-sous-Bois, France) as described by the manufacturer, in 9 FFPE samples. Quality control of qPCR data for consequent analysis was based on positive and negative PCR controls. Briefly, a total volume of 20μl of PCR mix, which included 10μl of Perfect MasterMix SYBR Green^®, 8μl of PCR grade water and 2μl of cDNA was loaded into each well of the qPCR array. PCR amplification was conducted at 95°C for 10 min, followed by 40 cycles of 95°C for 10 sec and 60°C for 30 sec. The mRNA expression for each gene was normalized using the average expression of 3 housekeeping genes: peptidylprolyl isomerase A (cyclophilin A, PPIA), b-actin (ACTB), and Glyceraldehyde-3-phosphate deshydrogenase (GAPDH). Data analysis was conducted using AnyGenes^® Excel analysis tools based on ΔΔCt method by calculating fold changes for each gene as the difference in gene expression between positive HIF-1α and negative HIF-1α samples.

Nanostring analysis

Total RNAs were extracted from a cohort of 19 tissues patients (Supplementary Table 1). 8 RNA samples already analyzed by microarray and/or by RT-qPCR were included as controls in the analysis : 6 RNA samples from 3 cell lines (Mel-1, Mel-6 and Mel-10) cultured under normoxia and hypoxia, 1 RNA sample extracted from hypoxia FFPE zone positive of a primary melanoma (zone A4, Figure 3) and 1 RNA sample extracted from hypoxia FFPE zone negative of a primary melanoma (zone A1, Figure 3). Total RNA concentration and purity (Ratio 260/280 and ratio 260/230nm) were calculated using a Nanodrop ND8000 spectrophotometer (Ozyme, Saint-Quentin en Yvelines, France). Total RNA integrity was assessed using a micro electropheresys (RNA6000 LabChip, Agilent technologies), and RIN or percentage of fragment longer than 300 nt were calculated, upon a total RNA migration. Direct quantification of mRNA was achieved according a Nanostring Custom Elements approach. 50 ng of total RNA were used as template to detect 32 targets corresponding to 26 mRNA of interest and 6 housekeeping genes. The Nanostring nSolver software was used to control raw data, and to normalize data based on geometric mean of positive controls, and water signals to deduce unspecific counts.

With the goal of identifing biomarker of response to therapy, we looked for pairs of highly correlated genes based on their raw expression levels (correlation level > 0.80), under the rationale that such correlated genes would normally share common biological properties yielding similar expression levels, however could differ in their response to treatment [16]. Within each pair of correlated genes, the difference between genes, called “differential pair”, was calculated, thereby ensuring that the variation in expression, not caused by the response to treatment itself, will be optimally filtered out. The matrix of differential pairs was further split into two parts corresponding to positive and negative differences. Within first part, negative differences were replaced by 0. Within second part, positive differences were replaced by 0, whereas negative differences were replaced by their absolute value. Non-Negative Matrix Factorization (NMF) was applied to the split matrix [17]. Specifically, two expression profiles among genes, each one linked to a particular expression pattern, were simultaneously estimated. A dual clustering of samples and genes into two clusters was achieved by grouping together samples and genes that loaded most on a particular pattern or expression profile. Finally, within each cluster, samples and genes were ordered by descending leverage on their respective cluster. The association between the ordering of samples and responder status was assessed through a permutation test. Responder status was permuted among patients and the rank-sum score was re-calculated and compared to the original score.

Abbreviations

HIF, hypoxia inducible factor

Author contributions

Conception and design: S Chouaib; Acquisition of the data: S Buart, S Terry, L Legrès, MZ Noman, P Dessen; Technical support: S Buart, S Terry, MZ Noman, V Marty, L Legrès, N Kamsu-Kom, M El Bouchtaoui; Analysis and interpretation of the data: S Chouaib, S Buart, P Dessen, S Terry, P Vielh, S Ferlicot, Y Gaston-Mathe, P Fogel; Drafting the manuscript: S Chouaib; Critical revision of the manuscript for important intellectual content: A Eggermont, C Robert, S Terry, C Boutros; Statistical analysis: P Dessen, Y Gaston-Mathe, P Fogel, E Lanoy, G Meurice; Obtaining funding: S Chouaib, JC Soria, C Robert; Clinical support: E Deutsch, J Muret, S Roy, C Boutros, E Routier; Supervision: S Chouaib.

ACKNOWLEDGMENTS

This work was supported OSEO, la Ligue Contre le Cancer (EL2015.LNCC/SaC). We would like to acknowledge Lisa Bain for carefully editing this manuscript. We thank Olivia Bawa from Gustave Roussy, Villejuif, France, for her assistance in immunohistochemistry staining. We thank David Gentien and Benoit Albaud from Institut Curie, Paris, France, for generating raw data of the Nanostring analysis.

CONFLICTS OF INTEREST

We declare no competing financial interests.

REFERENCES

1. Sharma P, Wagner K, Wolchok JD, Allison JP. Novel cancer immunotherapy agents with survival benefit: recent successes and next steps. Nat Rev Cancer. 2011; 11:805-812.

2. Boon T, Cerottini JC, Van den Eynde B, van der Bruggen P, Van Pel A. Tumor antigens recognized by T lymphocytes. Annu Rev Immunol. 1994; 12:337-365.

3. Chouaib S. Integrating the quality of the cytotoxic response and tumor susceptibility into the design of protective vaccines in tumor immunotherapy. J Clin Invest. 2003; 111:595-597.

4. Zippelius A, Batard P, Rubio-Godoy V, Bioley G, Lienard D, Lejeune F, Rimoldi D, Guillaume P, Meidenbauer N, Mackensen A, Rufer N, Lubenow N, Speiser D, et al. Effector function of human tumor-specific CD8 T cells in melanoma lesions: a state of local functional tolerance. Cancer Res. 2004; 64:2865-2873.

5. Hanahan D, Folkman J. Patterns and emerging mechanisms of the angiogenic switch during tumorigenesis. Cell. 1996; 86:353-364.

6. Semenza GL. Defining the role of hypoxia-inducible factor 1 in cancer biology and therapeutics. Oncogene. 2010; 29:625-634.

7. Mazzone M, Dettori D, de Oliveira RL, Loges S, Schmidt T, Jonckx B, Tian YM, Lanahan AA, Pollard P, de Almodovar CR, De Smet F, Vinckier S, Aragonés J, et al. Heterozygous deficiency of PHD2 restores tumor oxygenation and inhibits metastasis via endothelial normalization. Cell. 2009; 136:839-851.

8. Miroshnikova YA, Mouw JK, Barnes JM, Pickup MW, Lakins JN, Kim Y, Lobo K, Persson AI, Reis GF, McKnight TR, Holland EC, Phillips JJ, Weaver VM. Tissue mechanics promote IDH1-dependent HIF1alpha-tenascin C feedback to regulate glioblastoma aggression. Nat Cell Biol. 2016; 18:1336-1345.

9. Crespo J, Sun H, Welling TH, Tian Z, Zou W. T cell anergy, exhaustion, senescence, and stemness in the tumor microenvironment. Curr Opin Immunol. 2013; 25:214-221.

10. Brizel DM, Scully SP, Harrelson JM, Layfield LJ, Bean JM, Prosnitz LR, Dewhirst MW. Tumor oxygenation predicts for the likelihood of distant metastases in human soft tissue sarcoma. Cancer Res. 1996; 56:941-943.

11. Noman MZ, Hasmim M, Messai Y, Terry S, Kieda C, Janji B, Chouaib S. Hypoxia: a key player in antitumor immune response. A Review in the Theme: Cellular Responses to Hypoxia. Am J Physiol Cell Physiol. 2015; 309:C569-79.

12. Semenza GL. Targeting HIF-1 for cancer therapy. Nat Rev Cancer. 2003; 3:721-732.

13. Hendriksen EM, Span PN, Schuuring J, Peters JP, Sweep FC, van der Kogel AJ, Bussink J. Angiogenesis, hypoxia and VEGF expression during tumour growth in a human xenograft tumour model. Microvasc Res. 2009; 77:96-103.

14. Robert C, Schachter J, Long GV, Arance A, Grob JJ, Mortier L, Daud A, Carlino MS, McNeil C, Lotem M, Larkin J, Lorigan P, Neyns B, et al, and KEYNOTE-006 investigators. Pembrolizumab versus Ipilimumab in Advanced Melanoma. N Engl J Med. 2015; 372:2521-2532.

15. Noman MZ, Buart S, Romero P, Ketari S, Janji B, Mari B, Mami-Chouaib F, Chouaib S. Hypoxia-inducible miR-210 regulates the susceptibility of tumor cells to lysis by cytotoxic T cells. Cancer Res. 2012; 72:4629-4641.

16. Sheinerman KS, Tsivinsky VG, Abdullah L, Crawford F, Umansky SR. Plasma microRNA biomarkers for detection of mild cognitive impairment: biomarker validation study. Aging (Albany NY). 2013; 5:925-38. https://doi.org/10.18632/aging.100624.

17. Fogel P, Gaston-Mathe Y, Hawkins D, Fogel F, Luta G, Young SS. Applications of a novel clustering approach using non-negative matrix factorization to environmental research in public health. Int J Environ Res Public Health. 2016; 13. 509, https://doi.org/10.3390/ijerph13050509

18. Jacobs KM, Bhave SR, Ferraro DJ, Jaboin JJ, Hallahan DE, Thotala D. GSK-3beta: a bifunctional role in cell death pathways. Int J Cell Biol. 2012; 2012:930710.

19. Hasmim M, Janji B, Khaled M, Noman MZ, Louache F, Bordereaux D, Abderamane A, Baud V, Mami-Chouaib F, Chouaib S. Cutting edge: NANOG activates autophagy under hypoxic stress by binding to BNIP3L promoter. J Immunol. 2017; 198:1423-1428.

20. Maes H, Van Eygen S, Krysko DV, Vandenabeele P, Nys K, Rillaerts K, Garg AD, Verfaillie T, Agostinis P. BNIP3 supports melanoma cell migration and vasculogenic mimicry by orchestrating the actin cytoskeleton. Cell Death Dis. 2014; 5:e1127.

21. Brahimi-Horn MC, Bellot G, Pouyssegur J. Hypoxia and energetic tumour metabolism. Curr Opin Genet Dev. 2011; 21:67-72.

22. Slominski A, Kim TK, Brozyna AA, Janjetovic Z, Brooks DL, Schwab LP, Skobowiat C, Jozwicki W, Seagroves TN. The role of melanogenesis in regulation of melanoma behavior: melanogenesis leads to stimulation of HIF-1alpha expression and HIF-dependent attendant pathways. Arch Biochem Biophys. 2014; 563:79-93.

23. Slominski A, Tobin DJ, Shibahara S, Wortsman J. Melanin pigmentation in mammalian skin and its hormonal regulation. Physiol Rev. 2004; 84:1155-1228.

24. Brozyna AA, Jozwicki W, Carlson JA, Slominski AT. Melanogenesis affects overall and disease-free survival in patients with stage III, IV melanoma. Hum Pathol. 2013; 44:2071-2074.

25. Brozyna AA, Jozwicki W, Roszkowski K, Filipiak J, Slominski AT. Melanin content in melanoma metastases affects the outcome of radiotherapy. Oncotarget. 2016; 7:17844-17853. https://doi.org/10.18632/oncotarget.7528.

26. Slominski RM, Zmijewski MA, Slominski AT. The role of melanin pigment in melanoma. Exp Dermatol. 2015; 24:258-259.

27. Li W, Slominski R, Slominski AT. High-resolution magic angle spinning nuclear magnetic resonance analysis of metabolic changes in melanoma cells after induction of melanogenesis. Anal Biochem. 2009; 386:282-284.

28. Slominski A, Zbytek B, Slominski R. Inhibitors of melanogenesis increase toxicity of cyclophosphamide and lymphocytes against melanoma cells. Int J Cancer. 2009; 124:1470-1477.

29. Favaro E, Bensaad K, Chong MG, Tennant DA, Ferguson DJ, Snell C, Steers G, Turley H, Li JL, Gunther UL, Buffa FM, McIntyre A, Harris AL. Glucose utilization via glycogen phosphorylase sustains proliferation and prevents premature senescence in cancer cells. Cell Metab. 2012; 16:751-764.

30. Bhanot H, Reddy MM, Nonami A, Weisberg EL, Bonal D, Kirschmeier PT, Salgia S, Podar K, Galinsky I, Chowdary TK, Neuberg D, Tonon G, Stone RM, et al. Pathological glycogenesis through glycogen synthase 1 and suppression of excessive AMP kinase activity in myeloid leukemia cells. Leukemia. 2015; 29:1555-1563.

31. Nordsmark M, Bentzen SM, Rudat V, Brizel D, Lartigau E, Stadler P, Becker A, Adam M, Molls M, Dunst J, Terris DJ, Overgaard J. Prognostic value of tumor oxygenation in 397 head and neck tumors after primary radiation therapy. An international multi-center study. Radiother Oncol. 2005; 77:18-24.

32. Mouriaux F, Sanschagrin F, Diorio C, Landreville S, Comoz F, Petit E, Bernaudin M, Rousseau AP, Bergeron D, Morcos M. Increased HIF-1alpha expression correlates with cell proliferation and vascular markers CD31 and VEGF-A in uveal melanoma. Invest Ophthalmol Vis Sci. 2014; 55:1277-1283.

33. Zhou J, Huang S, Wang L, Yuan X, Dong Q, Zhang D, Wang X. Clinical and prognostic significance of HIF-1alpha overexpression in oral squamous cell carcinoma: a meta-analysis. World J Surg Oncol. 2017; 15:104.

34. Smyth GK. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol. 2004; 3:Article3.