Molecular Portrait of Hypoxia in Breast Cancer: A Prognostic Signature and Novel HIF-Regulated Genes

Increased cell proliferation and oxygen consumption result in lower oxygen availability in solid tumors as compared with normal tissue (1, 2). Intratumoral hypoxia has been associated with invasion, metastasis, treatment failure, and patient mortality (3). Although mechanistic data from preclinical studies and survival data from clinical studies indicate that intratumoral hypoxia and increased HIF-1α expression are associated with aggressive cancers (4–6), the ability to identify individual patients at increased risk is still limited. This is likely a reflection of the heterogeneous nature of the disease (7), both across patients with the same diagnosis and within individual tumors as a function of time, treatment, and spatial localization of hypoxia within a tumor.

Cancer cells survive and adapt to hypoxic conditions, in part, through the activation of hypoxia-inducible factor 1 (HIF-1) and HIF-2, which induce the expression of gene products involved in angiogenesis, glucose utilization, invasion, and metastasis (8). HIF-1 (-2) is a heterodimeric protein composed of a constitutively expressed HIF-1β subunit and an O₂-regulated HIF-1α (-2α) subunit that act as transcription factors (9, 10). HIF-1 or HIF-2 regulate the expression of potentially thousands of genes. How HIF response differs between cell types of different organs or between cells in the same organ that have undergone transformation still remains elusive.

Several studies have used gene expression microarrays (see review; ref. 11) or chromatin immunoprecipitation (ChIP) of HIF-1α or HIF-2α followed by microarray or RNA sequencing (RNA-seq; refs. 12–14) to attempt to identify global transcriptional responses to hypoxia and/or HIFs. This has typically been carried out in cell lines exposed to normal tissue culture conditions (20% O₂) as compared with cells exposed to hypoxic conditions (0.5%–2.0% O₂) and has resulted in a variety of “hypoxia gene signatures” (11). A recent review aimed at comparing hypoxia signatures identified 32 published gene signatures and concluded that no single gene was present among all signatures (11). The authors suggest that this could reflect a difference between cell lines and clinical samples from different tumor types. Twenty-seven of the studies used at most 7 cell lines to derive a hypoxia signature. The remaining 7 studies used classic hypoxia-regulated genes in cluster-type analyses to identify additional hypoxia-regulated genes. Currently, there is no consensus as to the best way to derive a hypoxia signature. However, it is clear that such a signature would be useful in the stratification of patients who may benefit from therapies currently under development that are aimed at targeting the hypoxic response. For example, a recent study highlighted the utility of a 26-gene hypoxia signature to predict benefit from hypoxia-modifying treatment in laryngeal cancer (15).

Seminal studies by the Gray group, Perou group, The Cancer Genome Atlas (TCGA), and METABRIC (16–21) have uncovered the diverse nature of breast cancer. An analysis of a large panel of breast cancer cell lines has demonstrated that the same heterogeneity in copy number and expression abnormalities found in cell lines are also found in primary tumors (16). The most common cell line panels used in breast cancer research represent almost all of the recurrent genomic abnormalities associated with clinical outcome in primary tumors (16). Early studies using unsupervised clustering of mRNA expression data from cell lines and patient samples led to the identification of “intrinsic subtypes”: Luminal A and B, basal-like, HER2⁺, and normal-like (17, 22), although subsequent groups have further stratified these disease entities (23–25). The intrinsic molecular subtypes complement, but do not fully overlap, pathologic classification by estrogen receptor (ER), PR, and HER2 status (26). Given that in breast cancer tissue, the measured partial pressure of oxygen (PO₂) ranges between 2.5 and 28 mm Hg, with a median value of 10 mm Hg (∼1% O₂), as compared with 65 mm Hg in normal breast tissue (27), we sought to determine whether hypoxia could alter the intrinsic subtype classification of the cell lines and to identify a conserved hypoxia-inducible gene signature in breast cancer.

To assess hypoxia-induced transcriptional profile changes specifically in breast cancer cells, we exposed a panel of 31 breast cancer cell lines or normal human mammary epithelial cells to 20% or 1% O₂. Our results show that over 1,000 genes are induced or repressed in each cell line in response to hypoxia; however, only 42 genes shared a conserved response to hypoxia. We demonstrate that cells classified as normal-like have an overall stronger response to hypoxia both in terms of number of genes induced by hypoxia as well as level of expression. An extensive literature search as well as loss-of-function assays for HIF expression demonstrate that HIF-1α or HIF-2α are required for the induction of the 42-gene signature under hypoxic conditions. This hypoxia gene signature is also enriched in basal-like cell lines and basal-like tumors over and above other subtypes. The signature was prognostic in patients with basal-like breast cancer from two METABRIC cohorts (19, 20) as well as the 2017 KMplotter (28) cohort and demonstrated statistically significant differences in survival upon stratification. Of the 42 genes identified, we show that caspase-14 (CASP14) is a novel hypoxia-regulated gene bound by HIF-1 but not HIF-2. We also show for the first time that Fucosyltranferase 11 (FUT11), Dynein axonemal heavy chain 11 (DNAH11), and TRPM8 channel–associated factor 2 (TCAF2) are induced under hypoxic conditions in an HIF-dependent manner. The 42-gene hypoxia signature is prognostic for breast cancer patient outcome, and may provide a robust set of hypoxia biomarkers for clinical use.

All cell lines except SUM159 and SUM149 were obtained directly from the ATCC or Asterand (SUM lines). Cells were cultured as described by the ATCC or Asterand, and experiments were conducted on early passage cells within 2 to 3 passages after receipt. The SUM149 and SUM159 cells were gifts from the Sukumar laboratory and were authenticated by STR sequencing and confirmed to be Mycoplasma free. We classified each cell line according to the PAM50 subtype classification as recently published by Marcotte and colleagues (25). In brief, the authors (25) utilized a centroid-based classification method. The expression of each gene in the classifier was median-centered across 78 cell lines classified by Marcotte and colleagues. Pearson correlation was used to match each cell line to an intrinsic subtype, defined as the subtype with the highest associated correlation coefficient. We used the reported subtype values for each of the cell lines from their analysis as the subtype calls for the cell lines in our current study. Human telomerase immortalized mammary epithelial cells were not considered in the study above and thus classified in the current study as htert-HME. Hypoxic cells were maintained at 37°C in a modular incubator chamber (Billups–Rothenberg) flushed with a gas mixture containing 1% O₂, 5% CO₂, and 94% N₂ for 24 hours.

Total RNA was extracted from cells using TRIzol (Invitrogen) and treated with DNase I (Ambion). All samples had an RIN value of >9.0 when measured on an Agilent Bioanalyzer. One microgram of total RNA was used for first-strand DNA synthesis with the iScript cDNA Synthesis system (Bio-Rad). qPCR was performed using human-specific primers and iTaq SYBR Green Universal Master Mix (Bio-Rad). The expression of each target mRNA relative to 18S rRNA was calculated on the basis of the threshold cycle (C_t) as 2^−Δ(ΔC_t⁾, where |$\Delta {C_{\rm{t}}} = {\rm{ }}{C_{\rm{t}}}_{_{{\rm{target}}}} - {\rm{ }}{C_{\rm{t}}}_{_{_{{\rm{18S}}}}}$| and |$\Delta ( {\Delta {C_{\rm{t}}}} ){\rm{ }} = {\rm{ }}\Delta {C_{\rm{t}}}_{_{{\rm{test}}}} - {\rm{ }}\Delta {C_{\rm{t}}}_{_{{\rm{control}}}}$|⁠. Primer sequences are shown in Supplementary Table S1.

Libraries for RNA-seq were prepared with KAPA Stranded RNA-Seq Kit. The workflow consisted of mRNA enrichment, cDNA generation, end repair to generate blunt ends, A-tailing, adaptor ligation, and 14 cycles of PCR amplification. Unique adaptors were used for each sample in order to multiplex samples into several lanes. Sequencing was performed on Illumina Hiseq 3000/4000 with a 150bp pair-end run by Quick Biology. A data quality check was done on Illumina SAV. Demultiplexing was performed with Illumina Bcl2fastq2 v 2.17 program.

Transcript level counts for the cell line RNA-seq data were quantified from reads with salmon version 0.8.2 (29) using human genome build h19. Gene-level counts were computed with Tximport version 1.2.0 (30) using Gencode v26 annotations. Normalized and raw data were uploaded to GEO (GSE111653). TCGA level 3 RNA-seq RSEM V2 data for breast cancer tumors were used for analyses, with subtypes defined in the consortium marker article (21). Normalized gene expression data for genes in the hypoxia signature defined in this article from both METABRIC cohorts (19, 20) were obtained using the CGDS-R package to access data from the cBioPortal database (31).

Differential expression statistics were computed on read counts with DESeq2 version 1.14.1 (32). Heatmaps were generated with ComplexHeatmap version 1.12.0 (33) with default hierarchical clustering on read counts normalized with the rlog function in DESeq2. We perform paired tests that retain cell-line pairing using a multi-factor design. Genes with FDR-adjusted P values below 0.05 and log fold change greater than 1 were called statistically significant. The calculation of the number of genes with fold changes in different ranges were computed based upon differences between hypoxic and nonhypoxic rlog output values and converted from log₂ differences to report values as fold changes. Gene set statistics were computed with Hallmark gene sets from MSigDB version 6.0 (34) and HIF gene target sets curated from the literature (Supplementary Table S2; ref. 35) using a one-sided Wilcoxon gene set test in LIMMA (36). Gene sets with FDR-adjusted P values below 0.05 were considered significantly enriched. These analyses were performed in R using scripts for normalization and analysis available as Supplementary Files S1 and S2, respectively.

LentiCRISPR v2 plasmid used for generating a CRISPR-Cas9 endonuclease was a gift from Feng Zhang [Broad Institute, Massachusetts Institute of Technology, Cambridge, MA, obtained via Addgene (Addgene plasmid #52961)]. HIF-1α/2α knockout by CRISPR/Cas9 was performed as described previously with slight modifications (37). Insert oligonucleotides that include a guide RNA sequence were annealed and inserted into the BsmBI cloning site. After bacterial transformation and DNA purification, all plasmid constructs were confirmed by Sanger sequencing. The LentiCRISPR v2 HIF-1α or HIF-2α or a nontargeting gDNA control (NTC) were transfected along with helper vectors into 293T cells using PolyJet transfection reagent (SignaGen Laboratories) according to the manufacturer's instructions. Media were refreshed 16 to 24 hours following transfection, and filtered viral supernatant is collected 24 to 48 hours later. Viral supernatant was added to MCF-7 or MDA-MB-231 cells, and puromycin (0.5 μg/mL) was added to select for cells expressing the lentiviral vector.

Cells were lysed in IGEPAL CA-630 buffer (150 mmol/L NaCl, 1% IGEPAL CA-630, 50 mmol/L Tris-HCl, pH 8.0 and protease inhibitors) for 10 minutes on ice, centrifuged for 10 minutes at 13,000 rpm at 4°C, and the insoluble debris were discarded. Whole-cell lysates were fractionated by 12% SDS-PAGE and transferred to a nitrocellulose membrane (Bio-Rad). The membrane was incubated for 1 hour with 5% milk in TBS-T (Tris-buffered saline and 0.1% Tween-20) and then incubated overnight with primary antibodies diluted in blocking buffer. Antibodies against the following proteins were used: HIF-1α (BD Biosciences), HIF-2α (Novus Biologicals), Actin (ProteinTech). The membrane was then washed and incubated with the corresponding HRP-conjugated secondary antibody (Azure Biosystems) for 2 hours. Chemiluminescence signal was detected on an AZURE C300 (Azure Biosystems) using ECL (PerkinElmer).

BT474 cells were exposed to hypoxia for 4 hours and immediately cross-linked with formaldehyde, quenched with glycine, and lysed with SDS lysis buffer (1% SDS, 10 mmol/L EDTA, and 50 mmol/L Tris, pH 8.1). Chromatin was sheared by sonication using a Covaris sonicator [settings: power (W), 150; duty factor, 5%; cycles, 200; treatment time, 420 seconds (7 minutes)], and size distribution of DNA fragments ranging from 300 to 800 bp was confirmed by running samples on an agarose gel. Following sonication, lysates were precleared with salmon sperm DNA/protein A-agarose slurry (Millipore) and incubated with antibody against HIF-1α (Santa Cruz Biotechnology), HIF-1β (Novus Biologicals), or HIF-2α (Novus Biologicals), or with IgG (Santa Cruz Biotechnology), as described previously (38). Following elution, DNA enrichment of the CASP14 promoter was determined using primers sequences listed in Supplementary Table S1.

We performed an RNA-seq analysis of 31 well-characterized breast cancer cell lines or normal mammary epithelial cells exposed to control tissue culture conditions (20% O₂) or hypoxic (1% O₂) conditions for 24 hours (in triplicate). The cell lines in this study (see Supplementary Table S3) include basal-like, Luminal A, Luminal B, HER2, and normal-like cells as classified using PAM50 analysis by Marcotte and colleagues (25). Prior to performing RNA-seq, each sample was tested for the induction of three established HIF-1α inducible RNA transcripts using real-time PCR (RT-PCR): prolyl 4-hydroxylase 1 (P4HA1), BCL2/adenovirus E1B 19 kDa protein-interacting protein 3 (BNIP3), and VEGFA in order to verify that the cells respond to hypoxia for validation purposes (39–41). The expression of P4HA1, BNIP3, and VEGFA was significantly increased in each of the 31 cell lines analyzed (Fig. 1A–C).

Following the experimental verification of hypoxia exposure, RNA libraries were constructed from RNA lysates followed by RNA-seq to identify gene expression alterations among all cell lines. We performed an unsupervised hierarchical clustering analysis of the 1,000 most variable genes across all samples to delineate the main sources of variability in the whole group of samples (Supplementary Fig. S1). The analysis revealed Luminal A, Luminal B, and HER2 cells clustered together, whereas basal-like cells clustered with cells classified as normal-like. This is in keeping with the finding that noncancerous MCF10A and MCF12A cells clustered with Basal B in studies by Neve and colleagues (16). We also found that hTERT-immortalized HMEs clustered with basal-like cell lines. The analysis demonstrated that hypoxia did not alter clustering based on breast cancer subtype.

We performed a differential expression analysis to characterize global differences in RNA transcript levels induced under hypoxic conditions among all cell lines. Specifically, we performed a paired test between 1% O₂ and 20% O₂ conditions retaining the cell line pairing using a multi-factor design with DESeq2 (32). Volcano plots highlight genes (in red) with a mean log fold change across all cell lines (LFC) greater than 1 and an FDR-adjusted P value less than 0.05 that are significantly increased in cells in response to hypoxia (Fig. 1D). We identified 42 genes (including noncoding RNAs and pseudogenes) that met these criteria and used them to define the “conserved 42-gene hypoxia signature” (Supplementary Table S4). Interestingly, no genes with an FDR-adjusted P value below 0.05 had a LFC < −1, suggesting that the conserved response to hypoxia when grouping all cell lines together involves the induction rather than suppression of genes.

The 42 genes that were significantly induced by hypoxia include many previously reported HIF-dependent hypoxia-responsive genes including genes encoding proteins with central roles in glucose metabolism [glucose transporter type 1 (SLC2A1), fructose-2,6-bisphosphate synthesis (PFKFB4), lactate dehydrogenase A (LDHA), aldolases (ALDOA, ALDOC)], angiogenesis [angiopoietin-like 4 (ANGPTL4, VEGFA)], cell proliferation and apoptosis [BNIP3L, BNIP3, neuregulin 1 (NDRG1)], and extracellular matrix synthesis [P4HA1, lysyl oxidase (LOX)]. Four novel genes (CASP14, FUT11, DNAH11, and TCAF2) and 5 unreported transcripts coding for long noncoding RNAs or pseudogenes were also identified as part of the 42-gene hypoxia signature (Supplementary Table S4).

Next, we performed a gene set analysis using differential expression statistics with Hallmark gene sets from MSigDB (34). We identified 22 pathways that are upregulated under hypoxic conditions with an FDR-adjusted P value below 0.05 (Supplementary File S2). Not surprisingly, the top pathway identified was the hypoxia response pathway (P < 10⁻⁶⁰). In addition, many well-known hypoxia-regulated pathways such as the epithelial-to-mesenchymal transition, glycolysis, TGFβ signaling, and NOTCH signaling pathways were among the top 22 altered pathways.

We clustered gene expression data for each cell line by the normalized RNA gene expression levels of the 42-gene signature and examined the cell line clustering in cells that were not exposed to hypoxia. The analysis revealed that half of the basal-like cell lines (SUM1315MO, HS578T, BT549, MDA-MB-231, SUM159, SUM149, HBL100, HCC1569, and DU4475) have the highest level of expression of the 42-gene hypoxia signature under normal tissue culture conditions (20% O₂; Fig. 2A). Gene set analysis confirms that the set of 42 genes are significantly differentially expressed between basal-like cell lines relative to Luminal A (P = 7 × 10⁻⁷) and basal-like cell lines relative to Luminal B (P = 1.0 × 10⁻⁷), but not between basal-like and HER2 cells (P = 0.45) under nonhypoxic conditions.

Next, we compared the level of induction of the 42-gene hypoxia signature after exposure to hypoxia by performing a hierarchical clustering analysis on the fold changes (Supplementary Fig. S2) and gene expression values (Fig. 2). Cell lines exposed to hypoxia did not cluster based on PAM50 but rather clustered based on culture conditions (20% O₂ vs. 1% O₂ conditions; Fig. 2B). MCF12A, MCF10A, SUM149, SUM159, HCC1569, and MCF-7 cells have the greatest induction of the 42-gene signature in response to hypoxia. No obvious commonalities have been identified among these cell lines. For example, the cell lines do not share a similar ER, BRCA, p53, or PI3K status, suggesting that mechanisms outside of a specific genetic aberration may be involved.

In addition to the 42-gene signature, we also assessed the total number of genes in each cell line by fold change induction or repression between hypoxic and nonhypoxic conditions (Fig. 3A and B, respectively). Each of the cell lines had greater than 1,000 (out of 15,945 genes that were detected) genes induced under hypoxic conditions with a fold change of 2 or greater (Fig. 3A). Likewise, over 800 genes were repressed under hypoxic conditions with a fold change of 2 or greater in each of the cell lines analyzed (Fig. 3B).

Given the heterogeneity in the genes that are induced by hypoxia among the cell lines tested, we sought to determine whether a subset of genes was preferentially expressed on the basis of breast cancer cell subtype. We performed a paired test between hypoxic and nonhypoxic cells of each subtype and generated subtype-specific volcano plots (Fig. 4A–D). According to PAM50 classification, only SKBR3 and SUM185 cells were considered as an HER2 subtype. Differential expression analysis did not reveal differences between HER2-subtyped cells cultured under 20% or 1% O₂ conditions. On the other hand, cells classified as normal-like differentially express 139 genes that are upregulated under hypoxic conditions and 15 genes that are downregulated under hypoxic conditions and that are exclusive to normal-like cells (Fig. 4E). The normal-like cell lines also had the highest degree of overlap with other subtypes. Several subtype exclusive genes were also identified.

Given that in our analysis, many basal-like cell lines have high expression of the 42-gene signature in the absence of hypoxic conditions, we hypothesized that a similar stratification occurs in patient samples. We performed hierarchical clustering of RNA-seq data from patients with breast cancer in the TCGA database (21) for genes in our hypoxia signature. The analysis revealed subtype-specific clusters within the 42-gene hypoxia signature (Fig. 5). A one-sided gene set analysis using our hypoxia signature as the gene set performed on differential expression statistics for basal-like tumors relative to other subtypes confirms that genes in this hypoxia signature are significantly overexpressed in basal-like tumors relative to tumors of other subtypes (P = 0.029).

Given that a subset of the basal-like cell lines showed high levels of the hypoxia signature and still had a significant induction under hypoxic conditions, we hypothesized that the hypoxia gene signature could be used to stratify patients based upon survival. We defined clustered gene expression data for basal-like tumors in the two METABRIC sample cohorts (19, 20). Dividing the samples into two clusters identifies sets of basal-like tumor samples with low expression among genes in the hypoxia signature and basal-like tumor samples with high expression among genes in the hypoxia signature (Supplementary File S2). In both cohorts, the tumor samples from the basal subgroup with lower expression among genes in the hypoxia signature had significantly better prognosis than those in the higher expression group (P value of 0.031 for the 2012 cohort and P value of 0.0037 for the 2016 cohort; Supplementary Fig. S3A and S3B, respectively). Further analysis of basal tumors from available datasets last updated in 2017 in KMplotter (28) also confirmed that basal patients that had expression of the 42-gene hypoxia signature lower than median level of expression had a better prognosis than those that had high expression (P = 0.022; Supplementary Fig. S3C).

Hypoxia stimulates a variety of adaptive responses, many mediated via the HIF family of transcriptional complexes. However, the transcript levels of HIF-1 and HIF-2 are not significantly induced by hypoxic conditions in the cell lines assayed (Supplementary Fig. S4). This is not surprising given that the regulation of HIF-1 and HIF-2 is known to occur at the protein level via oxygen-dependent hydroxylation of proline and asparagine residues followed by ubiquitation and subsequent degradation through the proteasome (42).

HIF-1 and -2 directly upregulate gene transcription by binding to putative hypoxia response elements (5‘-RCGTG-3′) within target genes (35). To determine whether the genes in our signature were likely the result of HIF-1 or HIF-2 activation, we performed a gene set enrichment analysis using HIF gene target sets curated from the literature (Supplementary Table S2; ref. 35) and we demonstrated that these HIF targets are significantly enriched under 1% O₂ versus 20% O₂ conditions (one sided P < 1.1 × 10⁻¹⁶; Supplementary Fig. S5).

To further explore the relationship between hypoxia and HIF expression in our gene signature, we also compared our hypoxia gene list with 3 independent ChIP microarray or sequencing studies (Table 1). Thirty-three of the genes on our list were enriched by ChIP of either HIF-1α, HIF-2α, or both in at least one of three of the ChIP followed by microarray or DNA sequencing studies. Twenty-five genes on the list were validated in transactivation assays. Twenty-two of the genes were further validated as HIF-1α or HIF-2α target genes in standard DNA-binding assays such as ChIP or EMSA. If only a single assay was performed to implicate a role for HIFs in the induction of the specified gene under hypoxic conditions, we also queried the literature to identify a study using gain-or-loss of function or loss-of-function of HIF-1, HIF-2, or both to assess HIF dependence. Taken together, 31 of the 42 genes identified in the study were validated as HIF-regulated in two or more independent assays and/or publications, suggesting that HIFs are responsible for the conserved response to hypoxia in breast cancer cells.

Our study revealed CASP14 as a novel hypoxia-regulated gene product involved in antiapoptotic responses. To confirm this result, we utilized a subset of 12 breast cancer cell lines to confirm the induction of CASP14 expression upon exposure to 24 hours of hypoxia (Fig. 6A). We also confirmed the induction of FUT11, DNAH11, and TCAF2 under hypoxic conditions (Fig. 6A). Next, we generated lentiviral vectors expressing CAS9 and a unique gDNA to target HIF-1α or HIF-2α, in MDA-MB-231 and MCF-7 cells, for CAS9-mediated destruction. Western blot analysis demonstrates the suppression of HIF-1α or HIF-2α in the CRISPR-mediated knockout cell lines under hypoxic conditions (Supplementary Fig. S6A and S6B). RT-PCR of cDNA lysate confirms that blocking HIF-1α inhibits P4HA1 induction as we have previously shown (Supplementary Fig. S6C and S6D; ref. 39). RT-PCR analyses of the knockout cells demonstrate that FUT11, DNAH11, CASP14, and TCAF2 all require HIF-1α (but not HIF-2α) for induction under hypoxic conditions (Fig. 6B and C).

In order to determine whether CASP14 is directly regulated by HIF-1 binding to the promoter region, the human CASP14 gene sequence was examined for matches to the putative hypoxia response element 5′-(A/G)CGTG-3′. Three candidate sites were identified in the 5′UTR of CASP14 (Fig. 6D). Two binding sites preceded the first noncoding exon and the third was found in intron one. We performed ChIP assays of BT474 cells exposed to 20% or 1% O₂ for 4 hours to assess potential HIF-binding sites. To validate our samples and ChIP procedures, we tested for enriched binding of HIF-1α and HIF-1β at the previously identified PDK1 HIF-binding site (43) under 1% compared with 20% O₂ conditions (Supplementary Fig. S6E). We determined that all three candidate binding sites were enriched for HIF-1α and HIF-1β binding under hypoxic conditions, whereas HIF-1α and HIF-1β did not bind to RPL13A promoter region (Fig. 6E; Supplementary Fig. S6E). Taken together, the data indicated that HIF-1 activates CASP14 transcription under hypoxic conditions by binding to the CASP14 promoter.

The prevalence of hypoxia and its contribution to poor prognosis and treatment failure has provided a compelling rationale for clinical-stage drug development, but the recent failure of two major clinical trial efforts (tirapazamine and evofosfamide) calls for reflection (44). Translating these compounds into the clinic has been confounded by the complex biology and heterogeneity of hypoxia and failure to clearly define clinical indications in which hypoxia is treatment-limiting. Patient selection criteria used to select patients that are predicted to benefit from hypoxia-targeted agents are lacking in the design of clinical trials aimed at targeting hypoxic cells. Therefore, there is a need to identify a robust technique that can be used in the clinic to guide treatment decisions for the potential application of this class of drugs. Historically, polarographic oxygen electrodes have been utilized in clinical trials to measure microregional PO₂ in multiple locations within a tumor. Results from clinical studies in cervix, head and neck, and breast carcinomas indicate that hypoxic tumors have a worse prognosis (45–48). However, the invasive nature of the technique has limited its widespread use. Given the emerging use of microarray technologies and the development of various gene panels that attempt to predict clinical outcome and potentially the response to therapy, an appropriate panel that may serve as a proxy for hypoxia may prove useful. Although there have been many published “hypoxia signatures” that contain genes that have been validated in terms of prognostication in retrospective studies of patient outcome (11), there remains no consensus on the genes that should be included in such a hypoxic signature and no accepted validation procedure for signature development.

Our goal was to provide the first comprehensive study of the transcriptional response of 31 breast cancer cell lines to hypoxia. Our study reveals that although more than 1,000 genes may have altered expression patterns under hypoxic conditions in any one breast cancer cell line, less than 50 genes could be considered as having a conserved response to hypoxia. This highlights the complex heterogeneity in the hypoxic response and suggests a need for further studies to identify what leads to these differences. Cells classified as normal-like had the largest congruence with the transcriptional profile changes of all other subtypes under hypoxia (Fig. 4E), suggesting that as cells undergo transformation, evolve, and acquire unique mutations, they become more divergent in their response. We did identify 9 genes with differential expression in luminal (and ER⁺) cells compared with other subtypes. Three of the genes (STC2, ILVBL, and FAM13A) contain both ER and HIF-binding sites in their promoter, suggesting a potential cross-talk between ER-HIF as a potential mechanism for a unique luminal cancer cell response to hypoxia.

In 2012, the TCGA project analyzed breast cancer tissue samples with transcriptional profiling (21). The group utilized a PARADIGM analysis to look for pathways that were highly active in basal compared with luminal subtype tumors and reported that the HIF-1α/ARNT pathway was preferentially active in samples from patients with basal breast cancer. Increased HIF-1α pathway activation in this setting could reflect a high level of intratumoral hypoxia in basal breast tumors compared with other subtypes. However, our results show that the hypoxia gene signature was enhanced in subconfluent cultures of basal subtype breast cancer cells without the exposure of the cells to hypoxic conditions. Our findings are corroborated in a study that molecularly profiled a large panel of breast cancer cell lines and identified high HIF pathway activation in triple-negative breast cancer (TNBC) cells by gene set enrichment analysis (49). Increased HIF-1α protein levels have also been noted in basal subtype cell lines compared with luminal cells grown under well-oxygenated conditions (50). We note that we did not see increased HIF1α or HIF2α mRNA levels in basal compared with luminal breast cancer cells, suggesting a transcription-independent mechanism.

This highlights the need to determine how HIF pathway activation is maintained even in the presence of oxygen. One recent report suggests that TNBC cells have higher endogenous levels of HIF-1α due to a paracrine signaling mechanism that involves glutamate secretion (50). Glutamate production in turn inhibits the glutamate-cystine antiporter, xCT, lowering intracellular cysteine levels. Decreased intracellular cysteine levels inhibit the HIF prolyl-hydroxylases, leading to HIF1α stabilization (50). Even though many basal cell lines do have a high endogenous level of HIF-inducible genes under well-oxygenated conditions, it is important to note that the expression of hypoxia-inducible genes can be further increased upon exposure to hypoxia. Future studies are warranted to identify whether and how a “hypoxic memory” may be bestowed upon basal-like cells and whether HIFs are required to maintain the transcriptional profile. Our previously reported data (38, 39, 51) and data from Fig. 6 show that the knockout of HIF-1 causes a reduction in the genes that were studied in unexposed cells, suggesting a role for HIFs in the maintenance of hypoxia-induced transcription patterns.

Although, many HIF-dependent and independent pathways regulate the response to hypoxic conditions, our study suggests that HIFs mediate the most highly conserved response. For example, of the 42 most conserved transcripts identified among the 31 cell lines, 37 have been validated as HIF-1α or HIF-2α dependent, and the remaining 5 transcripts correspond to long-coding RNAs or pseudogenes that have yet to be studied. In addition, we did not find a single gene with a conserved repression upon exposure to hypoxia among the cell lines tested. This is in keeping with two recent genome-wide analyses that have suggested that both HIF-1α and HIF-2α mainly function as transcriptional activators and not as repressors (12, 13). Not surprisingly, many of the genes that are induced by hypoxia but that are unique to a handful of cell lines are also controlled by HIF regulation. Any number of mechanisms may alter the specific genes that are regulated by hypoxia in a given cell line. For example, HIFs interact with an extensive list of proteins and their expression level and association with HIFs may impact the level of response of some transcripts. Furthermore, tumor hypoxia also causes DNA hypermethylation. A recent study shows that loss of oxygen-dependent ten-eleven translocation (TET) enzymes increases hypermethylation at gene promoters in vitro (52), suggesting that promoter methylation may likely result in differences in gene expression under hypoxia between cell lines This is consistent with our finding that cell lines classified as normal-like had the strongest response to hypoxic conditions and the most overlap with each subtype. In addition, HIF-1α has been reported to alter nucleosome structure in an HIF-dependent manner and this appears to vary by promoter (53). Moreover, expression of genes encoding for histone demethylases, including JMJD1A [KDM3A, demethylase for Lys9 of histone H3 (H3K9)], which are induced by HIFs, may also dictate cell-specific responses to hypoxia.

Our study revealed four novel HIF-regulated genes: CASP14, TCAF2, FUT11, and DNAH11. We identified three HIF-1α–binding sites in the 5′ UTR of the CASP14 gene. Interestingly, CASP14 expression has been associated with triple-negative phenotypes and was recently shown to be a breast cancer stem cell marker (54). CASP14 has also been implicated as a gene found in brain metastatic compared with nonbrain metastatic breast cancer (55). FUT11 has been shown to be a potential biomarker of clear cell renal cell carcinoma progression based on meta-analysis of gene expression data (56), which may reflect VHL-dependent HIF activation. Further investigation into the functional consequences of the upregulation of these novel HIF-1–regulated genes may reveal additional insight into the role of hypoxia in breast cancer progression and metastasis.

No potential conflicts of interest were disclosed.

Conception and design: I.C. Ye, E.J. Fertig, D.M. Gilkes

Development of methodology: I.C. Ye, D.M. Gilkes

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): I.C. Ye, J.W. DiGiacomo

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): I.C. Ye, E.J. Fertig, J.W. DiGiacomo, M. Considine

Writing, review, and/or revision of the manuscript: I.C. Ye, E.J. Fertig, J.W. DiGiacomo, I. Godet, D.M. Gilkes

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): I.C. Ye, E.J. Fertig, I. Godet, D.M. Gilkes

Study supervision: E.J. Fertig, D.M. Gilkes

The authors would like to thank their funding sources for their support: Work in the Gilkes lab is supported by U54-CA210173, R00-CA181352, V Scholar Foundation, Susan G. Komen Foundation (CCR17483484), The Jayne Koskinas Ted Giovanis Foundation for Health and Policy, and the Breast Cancer Research Foundation. Work in the Fertig lab is supported by the Experimental and Computational Genomics Core (ECGC) of the Sidney Kimmel Comprehensive Cancer Center (P30CA006973). We also would like to thank João Brás, Seung Yeon Choi, Grace Kim, Jisu Shin, and Daniel Shade for support in the lab that aided in this project. We thank Dr. Ludmila Danilova for providing the TCGA data for analysis and Drs. Gregg Semenza, Ben Ho Park, and Sarawati Sukumar for careful reading and suggestions for this manuscript.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.