Identification of the acute transcriptional response to hypoxia

In order to identify rapid changes in transcriptional activity upon hypoxia, we performed PRO-seq in HCT116 wild-type (WT) and HIF1A^−/− cells under normoxia (21% O₂) and after 90 min of exposure to hypoxia (1% O₂) (Fig. 1a). This short treatment time was sufficient to induce both intracellular hypoxia, as detected in real-time using the fluorogenic hypoxia reporter Image-iT Red (Supplementary Fig. 1a), and stabilization of HIF1A and HIF2A (Fig. 1b). In addition, we measured HIF1A chromatin binding by ChIP-seq and changes in the steady-state transcriptome using poly(A) + RNA-seq after longer-term exposure to hypoxia (24 h) (Fig. 1a). This later time point was chosen to capture persistent changes in steady-state mRNA and because HIF1A chromatin binding patterns have been shown to be stable over time¹⁴.

a Genome browser view of PRO-seq, ChIP-seq, and RNA-seq signals across the ENO1 locus for HCT116 wild-type (WT) and HIF1A^−/− cells exposed to normoxia (21% O₂, blue) or hypoxia (1% O₂, red). b Western blot analysis of HIF1A and HIF2A levels in HCT116 wild-type cells exposed to 21% or 1% O₂ for the indicated times. Images are representative of at least two replicates with similar results. Source data are provided as a Source Data file. c DESeq2 differential expression analysis of PRO-seq signal within gene body regions for HCT116 WT (blue) and HIF1A^−/− (green) cells. Horizontal dashed lines indicate an FDR threshold of 10% for negative binomial Wald test; numbers and colored points indicate significant genes at this threshold. d Comparison of fold changes in PRO-seq gene body signal induced by 90 min hypoxia in HCT116 WT versus HIF1A^−/− cells. Red/green points denote genes with significant up/down regulation; points representing all other genes are colored by density. e Gene set enrichment analysis (GSEA) of “Hallmark” gene sets against genes ranked by gene body fold changes induced by 90 min hypoxia in HCT116 WT cells. Red/green points denote gene sets with significant (FDR 10%) positive/negative enrichment. f Distributions of fold changes in PRO-seq signal at gene bodies and transcription start sites (TSS) and pausing index (PI) for upregulated genes in HCT116 WT (blue) and HIF1A^−/− (green) cells. Horizontal spread of data points is proportional to density; boxes indicate medians and upper and lower quartiles; numbers above the HIF1A^−/− plots indicate effect size estimates and p-values for two-sided Mann–Whitney U tests against the corresponding measure in WT cells. g Metagene showing typical PRO-seq signal profiles across upregulated genes in HCT116 WT (blue) and HIF1A^−/− (green). Data are represented as splined linear fit lines with 95% confidence intervals in gray. h Upregulated genes ranked by WT/HIF1A^−/− fold-change in normoxia, comparing the effect of HIF1A in normoxia (green/black/red points) and hypoxia (gray points). i Heatmap of gene body RPKM Z-scores for select example genes in HCT116 WT and HIF1A^−/− cells in normoxia and hypoxia. See also Supplementary Figs. 1, 2 and Supplementary Data 1, 2.

At most human protein-coding genes, shortly after initiation, elongating RNAPII pauses at ~20–60 nucleotides downstream of transcription start sites (TSS), due to the action of negative elongation factors and before recruitment of positive elongation factors enables further elongation^17,18. This regulated rate-limiting step results in typical PRO-seq profiles with peaks of high read density near the TSS, representing paused RNAPII, and lower coverage throughout gene bodies, corresponding to elongating RNAPII. Thus, to identify changes in productive transcription in response to acute hypoxia, we first quantified PRO-seq signals within gene body regions and tested for differential transcriptional activity using DESeq2 (ref. ¹⁹) on data from biological replicates (see Methods, Supplementary Data 1 and Supplementary Fig. 1b, c for replicate comparison). Rapid transcriptional activation and repression were apparent after 90 min of hypoxia in WT cells (Fig. 1c). This early transcriptional response was largely dependent on the presence of HIF1A, with both activation and repression reduced in HIF1A^−/− cells (Fig. 1c, d), and included many well-characterized HIF1A targets such as ENO1, PDK1, HK1, and SLC2A1 (Supplementary Fig. 1d). Genome browser views for example downregulated and unchanged genes are shown in Supplementary Fig. 1e. Thus, while other O₂-sensitive factors such as lysine demethylases may rapidly modulate the chromatin environment^8,9, our data indicate that HIF1A plays a critical functional role in driving the acute transcriptional response to hypoxia. The residual changes observed in HIF1A^−/− cells could be driven by HIF2A, which is expressed in HCT116 cells and stabilized upon hypoxia (Fig. 1b).

Gene set enrichment analysis (GSEA) of the acute transcriptional response indicates positive enrichment of known hypoxia- and glycolysis-related genes (Fig. 1e, Supplementary Fig. 1f), and the upstream regulator analysis module of Ingenuity Pathway Analysis (IPA) suite correctly predicts inactivation of the prolyl hydroxylases that repress HIF activity (EGLNs) concurrently with activation of HIF1A (Supplementary Fig. 1g). Interestingly, GSEA reveals early repression of MYC targets (Fig. 1e, Supplementary Fig. 1f), which could be explained by HIF1A-dependent induction of the MYC repressor MXI1 (refs. ^20,21,22,23), a gene that we found is indeed induced at the early time point by PRO-seq (Fig. 1c, d, Supplementary Fig. 1d, Supplementary Data 1). Importantly, a comparison of the acute transcriptional response identified by PRO-seq with 22 published studies characterizing the hypoxic response using more traditional approaches revealed that many genes in our data set have not previously been linked to hypoxic signaling (Supplementary Data 2), thus providing a resource to further characterize hypoxia-regulated pathways.

In many transcriptional networks, stimulus-responsive genes are associated with the release of pre-loaded, paused RNAPII¹⁷. Indeed, relative to their gene bodies, acute hypoxia-inducible genes display relatively modest changes in transcriptional activity at TSSs and consequent decreases in pausing index (PI), with an impaired response in HIF1A^−/− cells (Fig. 1f, g). By contrast, acute hypoxia-repressed genes (defined by decreases in gene body activity) display mildly decreased signals at TSS along with increased PI values (Supplementary Fig. 2a–c). These observations are consistent with the ability of HIF1A to induce recruitment of transcription elongation factors to sites of active chromatin^15,24,25.

Notably, depletion of HIF1A does not have an obvious global impact on transcription during normoxia or hypoxia, with WT and HIF1A^−/− cells displaying highly correlated transcriptional activity at both TSSs and gene bodies (Supplementary Fig. 2d, e). However, the metagene profiles of HIF1A^−/− cells show mildly decreased signals at TSS (Fig. 1g, Supplementary Fig. 2b, c). Using DESeq2 analysis, we identified 386 TSS regions and 1979 gene body regions with significant changes in HIF1A^−/− cells in normoxia (Supplementary Fig. 2f, g, Supplementary Data 1). These observations led us to test if basal levels of HIF1A contribute to the expression of hypoxia-inducible genes during normoxia. Indeed, we found that a subset of 107 acute hypoxia-inducible genes showed lower basal expression in HIF1A^−/− cells (Supplementary Fig. 2h), including well-characterized HIF1A targets such as ETS1, ELF3, and DDIT4 (Fig. 1h-i, Supplementary Data 1). Other key HIF1A targets did not show this basal HIF1A-dependence in normoxia, with some displaying lower transcription in WT cells, such as the glucose transporters SLC2A1 and SLC2A3 (Fig. 1h, i). These observations are consistent with low but measurable HIF1A activity in normoxia, which is further supported by analysis of CRISPR genetic screen data under normoxic conditions (see later).

Acute transactivation associates with strong HIF1A binding

We next investigated the relationship between acute activation and repression, and HIF1A chromatin binding as measured by ChIP-seq (Supplementary Data 3). Expectedly, HIF1A ChIP-seq peaks identified in this data set are enriched for hypoxia response element (HRE) motifs (Supplementary Fig. 3a, b) and essentially all display increased signal in hypoxia (Fig. 2a). Although ChIP-seq enrichment signal increases at TSSs for all classes of genes upon hypoxia, only those that are acutely upregulated display stronger binding around their TSS during hypoxia relative to genes with non-significant (n.s.) differences in transcription, with HIF1A peaks being most frequent within 2.5 kbp of the TSS (Fig. 2b, Supplementary Fig. 3c). Although hypoxia induces HIF1A binding near hundreds of promoter regions, peaks associated with acutely transactivated genes tend to display a stronger enrichment signal relative to either repressed or unaffected genes (Fig. 2c). Furthermore, HIF1A enrichment signal at TSS or at peaks within 50 kbp tends to be higher for genes that depend on HIF1A for their upregulation during hypoxia than for genes that display HIF1A-independent upregulation (Supplementary Fig. 3d, e). When HIF1A peaks are classified as “proximal” and “distal” according to enrichment signal and distance from TSS (Supplementary Fig. 3f), the high-confidence proximal peaks are significantly overrepresented near upregulated (~36%) but not downregulated (~10%) genes (Supplementary Fig. 3g). Earlier studies suggest that HIF1A can also operate over larger distances by acting on preformed enhancer-promoter interactions^26,27. Our data support these findings, with distal HIF1A binding sites also being significantly overrepresented near acute upregulated (~21%) but not downregulated (~15%) genes or genes with no significant change during acute hypoxia (Supplementary Fig. 3h).

a Comparison of enrichment signal in normoxia and hypoxia for HIF1A peaks called in hypoxic HCT116 wild-type (WT) cells (n = 3996 peaks). Points are colored by density. b Meta profile showing typical HIF1A occupancy profile over transcription start sites (TSS) of upregulated (red), downregulated (green), and not significantly regulated (n.s., gray) genes during acute hypoxia. Data are represented as splined linear fit lines with 95% confidence intervals in gray. c Enrichment signal for HIF1A peaks within 50 kbp of TSS, separated by gene body differential expression class: not significantly regulated (n.s., gray), upregulated (red), downregulated (green). Horizontal spread of data points is proportional to density; boxes indicate medians and upper and lower quartiles; horizontal bars with numbers indicate FDR-adjusted p-values for two-sided Mann–Whitney U tests. d Meta profiles showing typical PRO-seq signal profile across distal HIF1A peak regions associated with genes upregulated (Productive) or not upregulated (Unproductive) during acute hypoxia in HCT116 WT (blue) and HIF1A^−/− (green) cells. Data are represented as splined linear fit lines with 95% confidence intervals in lighter shades, with positive and negative values reflecting signal density on + and − genomic DNA strands, respectively. e Distributions of fold changes in PRO-seq signal (+/− strands combined) within ±250 bp of distal HIF1A peaks, separated by association with genes upregulated (Productive) or not upregulated (Unproductive) during acute hypoxia in HCT116 WT (blue) and HIF1A^−/− (green) cells. Horizontal spread of data points is proportional to density; boxes indicate medians and upper and lower quartiles; numbers indicate FDR-adjusted p-values for two-sided Mann–Whitney U tests. f Distributions of DNaseI accessibility or ChIP-seq signal enrichment for various chromatin marks within regions corresponding to distal HIF1A peaks associated with genes upregulated (Productive, pink) or not upregulated (Unproductive, teal) during acute hypoxia. Horizontal spread of data points is proportional to density; boxes indicate medians and upper and lower quartiles; numbers indicate p-values for two-sided Mann–Whitney U tests. g, h Comparison of PRO-seq (gene body) fold-change with ChIP-seq peak signal for proximal (g) and distal (h) HIF1A peaks. Gene class indicates PRO-seq gene body differential expression: not significantly regulated (n.s., gray), upregulated (red), downregulated (green). Solid lines in corresponding colors denote linear model fits to the data, with 95% confidence intervals in gray. Numbers in upper left are Pearson correlation coefficients in corresponding colors; * denotes significant correlations (p-value < 0.05). i Proportions of genes in each class associated with proximal and/or distal HIF1A peaks. j Distributions of fold changes in gene body activity for upregulated genes associated with different classes of HIF1A peaks: proximal and distal (Both, pink), proximal only (green), distal only (teal), or no peak (purple). Horizontal spread of data points is proportional to density; boxes indicate medians and upper and lower quartiles; numbers indicate FDR-adjusted p-values for two-sided Mann–Whitney U tests. k Heatmap of activation Z-scores for upstream regulator predictions by the Ingenuity Pathway Analysis suite for acutely upregulated genes lacking associated HIF1A binding. l Heatmap showing relative expression of putative indirect regulators of the acute hypoxic response. Data are represented as row-wise Z-scores calculated from RPKM values. See also Supplementary Fig. 3 and Supplementary Data 3.

Closer examination of distal HIF1A peaks revealed that those associated with upregulated genes, referred herein as “productive” binding sites, exhibit bidirectional transcriptional activity even in normoxia, prior to upregulation of their associated genes (Fig. 2d). Bidirectional transcription at intergenic regions is a recognized hallmark of active enhancers^28,29,30. Transcription at these sites increases upon exposure to hypoxia and is dampened in HIF1A^−/− cells (Fig. 2d, e). By contrast, distal peaks not associated with upregulated genes, referred herein as “unproductive” binding sites, display much lower transcriptional activity that is mostly unaffected by hypoxia or the absence of HIF1A (Fig. 2d, e). When we examined available ENCODE data for HCT116 cells, we found that productive distal sites display increased DNaseI accessibility and increased enrichment of histone modifications associated with enhancers (e.g. H3K27ac, H3K4me1) relative to unproductive sites (Fig. 2f).

HIF1A enrichment signal for both proximal and distal peaks is positively correlated with fold changes in gene body transcription at upregulated genes, with a stronger correlation being observed for proximal peaks (Fig. 2g, h). Overall, genes associated with proximal and/or distal HIF1A peaks account for ~52% of acute transactivation events (Fig. 2i) and genes with both peak types tend to display larger increases in productive transcription (Fig. 2j). Together, these data indicate that HIF1A binding is associated with acute transactivation, and do not support a direct role for HIF1A in acute transcriptional repression, which may be explained in part by indirect repression of the MYC transcriptional program via MXI1 induction²² (Fig. 1d, e, Supplementary Fig. 1d, f). Thus, from this point forth, we focused on the acute transactivation response.

Curiously, nearly half of the acute transactivation events exhibit dependence on HIF1A despite having no associated HIF1A peaks within 50 kbp (Fig. 2i, j, “no peak” category), suggesting either that some indirect transactivation events can occur by 90 min of hypoxia, or that these genes are regulated by HIF1A enhancers located at much greater distances. In support of the former hypothesis, IPA upstream regulator analysis predicts that some of these genes are regulated by signaling pathways, transcription factors, and chromatin regulators that are themselves directly induced by HIF1A, such as TGFB signaling (which could be explained by the strong induction of the SMAD family members SMAD3 and SMAD9), FOXO3, and KDM3A (Fig. 2k, l, Supplementary Fig. 3i). Therefore, the bulk of the acute transactivation response dependent on HIF1A can be attributed to strong HIF1A binding near direct target genes and the indirect action of early targets of HIF1A.

CDK8 kinase activity is required for acute transactivation

We previously identified the Mediator-associated kinase CDK8 as a transcriptional cofactor of HIF1A³¹. However, our studies employed measurements of steady-state mRNA, which could not properly define direct versus indirect contributions of CDK8 to the transcriptional program during acute hypoxia. Furthermore, our previous studies did not test the role of CDK8 catalytic activity versus structural or scaffolding functions in HIF1A-driven transactivation. Importantly, in different settings, the kinase activity of CKD8 has been involved in both transcriptional repression and activation^32,33, and Mediator-associated kinases have been shown to have kinase-independent effects in some transcriptional networks^33,34. Therefore, we employed PRO-seq analysis of HCT116 cells engineered to express an “analog-sensitive” allele of CDK8 (CDK8^as/as) that can be specifically inhibited by bulky ATP analogs. The generation and characterization of this cell line were previously reported³⁵. Notably, the CDK8^as/as alleles behave as hypomorphs, showing decreased kinase activity even in the absence of bulky ATP analogs, while the remaining kinase activity can be fully blocked by the ATP analog 3MB-PP1 (ref. ³⁵). We thus repeated our PRO-seq analysis in HCT116 WT and CDK8^as/as cells treated with vehicle (DMSO) or 3MB-PP1 and exposed to 21% or 1% O₂ for 90 min (Supplementary Data 4, see Supplementary Fig. 4a, b for evaluation of biological replicates). Analysis of gene body transcription activity in DMSO- and 3MB-PP1-treated HCT116 cells revealed a widespread requirement for CDK8 kinase activity for hypoxia-driven transactivation, with a more obvious impact in the +3MB-PP1 conditions (Fig. 3a, Supplementary Fig. 4c, d). Among 1547 genes significantly transactivated upon hypoxia in this experiment, 405 were expressed at significantly lower levels upon full CDK8 inhibition with 3MB-PP1, while 99 genes showed increased expression (Fig. 3b, Supplementary Fig. 4c). Thus, among hypoxia-inducible genes, CDK8 behaves mostly as a positive regulator of transcription. This requirement for CDK8 kinase activity was also evident, albeit with milder quantitative effects, in DMSO-treated cells expressing the hypomorph CDK8^as/as alleles (Supplementary Fig. 4c, d).

a Heatmap showing relative gene body signal for genes with increased transcription activity after 90 min hypoxia in wild-type (WT) and CDK8^as/as HCT116 cells treated with 3MB-PP1 during normoxia (21% O₂) and hypoxia (1% O₂). Data are represented as row-wise Z-scores calculated from RPKM values. b Upregulated genes ranked by CDK8^as/as/WT fold-change in hypoxia + 3MB-PP1. Genes with significant differences in gene body activity during CDK8 inhibition are highlighted in red (increased activity) and green (decreased activity). c Comparison of top 10 pathway/function clusters enriched among CDK8 dependent and CDK8 repressed genes, as identified by Metascape. Heatmap color represents −log₁₀(p-value) from hypergeometric enrichment test. d Comparison of transcription regulators with enriched targets among CDK8-dependent and CDK8-repressed genes, as identified by the TRRUST module of Metascape. Heatmap color represents −log₁₀(p-value) from hypergeometric enrichment test. e Bubble plots showing relative transcription activity in WT and CDK8^as/as HCT116 cells exposed to normoxia (pink) or hypoxia (blue) for CDK8-dependent glycolytic genes. Circle areas correspond to gene-wise Z-scores calculated from gene body RPKM values. f Metagene showing typical PRO-seq signal profiles across CDK8-dependent upregulated genes in HCT116 WT (blue) and CDK8^as/as (orange) cells treated with 3MB-PP1. Data are represented as splined linear fit lines with 95% confidence intervals in gray. g PRO-seq fold-change distributions for gene body, transcription start site (TSS) and pausing index (PI) of CDK8-dependent upregulated genes in HCT116 WT (blue) and CDK8^as/as (orange) cells treated with 3MB-PP1. Horizontal spread of data points is proportional to density; boxes indicate medians and upper and lower quartiles; numbers below the CDK8^as/as plots indicate effect size estimates and p-values for two-sided Mann–Whitney U tests against the corresponding measure in WT cells. See also Supplementary Fig. 4 and Supplementary Data 4.

Pathway analysis reveals stronger enrichment of gene signatures associated with hypoxia, glycolysis, and extracellular matrix remodeling among hypoxia-inducible genes that require CDK8 kinase activity (Fig. 3c). Upstream regulator analysis identifies enrichment of genes regulated by HIF1A among those that required CDK8 for transactivation, whereas no clear upstream regulators are identified among “CDK8 repressed genes” (Fig. 3d). CDK8-dependent genes include prominent HIF1A targets involved in glycolysis (Fig. 3e), consistent with our previous finding that CDK8 inhibition impairs glycolysis³⁵. Among hypoxia-inducible genes that require CDK8 kinase activity for transactivation, CDK8 inhibition significantly reduced productive elongation, with lesser impacts on transcription at TSSs (Fig. 3f, g, Supplementary Fig. 4e, f). Thus, upon CDK8 inhibition, pausing indices were not decreased at these genes during hypoxia, reinforcing the notion of HIF1A and CDK8 working coordinately to stimulate RNAPII elongation (Fig. 3g).

Interplay between acute transactivation and mRNA expression

We next examined the relationship between the acute transactivation response and subsequent changes in steady-state mRNA levels at 24 h hypoxia (Supplementary Data 5). A majority (~59%) of acute transactivation events correspond to significant increases in steady-state polyadenylated mRNAs at the later time point (460 genes, Fig. 4a). We define this gene set henceforth as acute response genes. Expectedly, hypoxia induced many more changes at the mRNA level at the later time point that did not correspond to early transactivation events, and we refer to these genes as late response genes (1672 genes, Fig. 4a). We then asked how HIF1A binding related to these two groups. When probing only for proximal HIF1A binding, ~45% (209 out of 460) of acute response genes display nearby HIF1A binding, compared to only ~13% (217 out of 1672) of late response genes (Supplementary Fig. 5a). When extending the analysis to also include distal peaks, the percentages increase to ~58% (268 out of 460) for acute response genes and ~27% (457 out of 1672) for late response genes. Thus, genes that are hypoxia-inducible as observed by both PRO-seq and RNA-seq are more likely to be bound by HIF1A. Furthermore, HIF1A peaks associated with acute response genes tend to display stronger levels of HIF1A binding relative to peaks associated with late response genes (Fig. 4b). This exercise reveals that simply cross-referencing steady-state mRNA data with HIF1A ChIP-seq data would predict the existence of hundreds of putative direct HIF1A targets that are not immediately induced by hypoxia (i.e. 457 genes when using proximal + distal sites, Supplementary Fig. 5b). Additionally, this approach would miss dozens of acute response genes identified by PRO-seq for which nearby HIF1A is not evident (Supplementary Data 2).

a Qualitative overlap between genes called as upregulated in PRO-seq (90 min hypoxia) and RNA-seq data sets (24 h hypoxia). Odds ratio (OR) and p-value (p) are for two-sided Fisher’s exact test. b Distributions of HIF1A ChIP-seq enrichment signal for HIF1A peaks associated with genes upregulated in both PRO-seq and RNA-seq data sets (Acute response, pink) or RNA-seq only (Late response, teal). Horizontal spread of data points is proportional to density; boxes indicate medians and upper and lower quartiles. Number below plot indicates p-value for two-sided Mann–Whitney U test between the two distributions. c Comparison of transcription regulators with enriched targets among Acute response and Late response genes, as identified by the TRRUST module of Metascape. Heatmap color represents −log₁₀(p-value) from hypergeometric enrichment test. d Comparison of transcription activity at 90 min hypoxia (gene body RPKM) with mRNA levels after 24 h hypoxia (RNA RPKM) for Acute response genes. Points are colored by density and selected example genes are labeled. Pearson correlation coefficient is shown in upper left, * denotes significant correlation (p < 2.2e−16). e Bubble plots showing relative transcription activity in wild-type (WT) and CDK8^as/as HCT116 cells exposed to normoxia (pink) or hypoxia (blue) for CDK8-dependent glycolytic genes. Circle area corresponds to log₂(RPKM) values, scaled within each data set. f Comparison of transcription activity at 90 min hypoxia (gene body RPKM) with enrichment signal (mean normalized tag counts) at TSS regions (for ENCODE project chromatin and RNAPII data) or peak regions (for HIF1A ChIP-seq data) for Acute response genes. Points are colored by density. Blue lines denote linear model fits to the data, with 95% confidence intervals in gray. Numbers in upper left are Pearson correlation coefficients, * denotes significant correlations (10% FDR). FDR-corrected p-values are as follows: H3K27ac 4.07e−29, H3K79me2 1.29e−27, H3K9me3 5.36e−2, DNaseI 3.93e−6, RNAPII 3.70e−15, H3K4me3 1.26e−7. See also Supplementary Figs. 5, 6 and Supplementary Data 5, 6.

Interestingly, the subset of transcripts that were detected as hypoxia-inducible only by PRO-seq (234 genes, Fig. 4a) included many antisense RNAs, miRNAs, and long non-coding RNAs (lncRNAs) (Supplementary Fig. 5c,d, Supplementary Data 6). Given the evolving annotation of lncRNAs, many of which have been reclassified as eRNAs upon further examination³⁶, we hope our PRO-seq data set will enable a more detailed characterization of these transcripts in future studies.

Upstream regulator analysis demonstrates strong enrichment of HIF1A targets among acute response genes, whereas targets of a much larger set of transcriptional regulators are enriched among the late response genes (Fig. 4c), involving many other pathways (Supplementary Fig. 5e). From this point forth, we focused our analyses on acute response genes.

When comparing gene-body PRO-seq signals to mRNA RNA-seq signals for acute response genes, we observed that both RNA output levels and fold-change are positively correlated between the two measurements (Fig. 4d, Supplementary Fig. 5f), with many changes in transcriptional activity being amplified in terms of steady-state changes (Supplementary Fig. 5f). However, fold changes in gene body activity have a weaker correlation with output levels (Supplementary Fig. 5g). This suggests that early transcriptional output is an important driver of steady-state mRNA levels at much later time points, with hypoxia-inducible genes displaying expression levels covering several orders of magnitude (Fig. 4d). For example, HIF1A targets such as ENO1 and DDIT4 produce >100-fold more nascent RNA early on and steady-state mRNA later on relative to lowly transcribed genes such as PTPRR and NEDD9 (Fig. 4e). To investigate the mechanisms modulating the strength of transcriptional output, we analyzed variations in chromatin environment, RNAPII occupancy, and HIF1A binding. Transcriptional output at 90 min of hypoxia correlates positively with the pre-existing levels of histone marks associated with gene activity, such as H3K27 acetylation (H3K27ac) and H3K79 dimethylation (H3K79me2), but not so with H3K9 trimethylation (H3K9me3), a mark of repressed chromatin (Fig. 4f, Supplementary Fig. 6a). Highly transcribed acute response genes also show higher levels of DNAse I accessibility, RNAPII occupancy, and H3K4me3 (a mark of transcriptional initiation) around their promoters in normoxia (Fig. 4f, Supplementary Fig. 6a). By contrast, these chromatin features are not correlated with the fold-change in transcriptional activity (Supplementary Fig. 6b). Lastly, HIF1A peaks associated with highly transcribed genes tend to show higher enrichment signal in hypoxia (Fig. 4f). Thus, the strength of HIF1A binding correlates with both absolute transcriptional output (Fig. 4f) and fold-change in transcriptional output (Fig. 2g) among its target genes.

Altogether, these observations indicate that an open chromatin environment during normoxia is positively associated with the strength of transcriptional output (but not fold-change) during the acute response to hypoxia and with subsequent steady-state mRNA levels.

Conservation of the acute transcriptional response to hypoxia

To investigate the conservation of the acute transcriptional response to hypoxia, we repeated RNA-seq and HIF1A ChIP-seq analyses for three additional cancer cell types: RKO (colorectal carcinoma), A549 (lung cancer), and H460 (lung cancer) (Supplementary Data 3,5). Although we observed high overall diversity among hypoxia-induced mRNAs detected by RNA-seq after 24 h hypoxia (Supplementary Fig. 7a), ~90% of the acute response genes identified in HCT116 cells were induced in one or more additional cell lines (Fig. 5a, b). Genes upregulated in all four cell lines (core) tend to display greater increases in both nascent RNA and steady-state mRNA relative to less-conserved (shared 3, shared 2) or HCT116-specific genes (Fig. 5c, d).

a Proportions of genes upregulated at both mRNA and PRO-seq levels (Acute response) or at the mRNA level (All mRNAs) in HCT116 cells with mRNA upregulation in one or more cell types. b Heatmap showing relative mRNA levels for acute response genes, grouped by conservation as in (a), across all four cell types during normoxia (21% O₂) or hypoxia (1% O₂). Data are represented as Z-scores calculated from RNA-seq RPKM values. c, d Fold-change distributions for gene body transcription activity (c, PRO-seq) and mRNA levels (d, RNA-seq) in HCT116 cells, separated by conservation group as in (a): core (pink), shared 3 (green), shared 2 (blue), and HCT116 only (purple). Horizontal spread of data points is proportional to density; boxes indicate medians and upper and lower quartiles; horizontal bars with numbers indicate FDR-adjusted p-values for two-sided Mann–Whitney U tests. e Heatmap of enrichment signal (±2.5 kbp, 100 bp bins) in HCT116, RKO, A549, and H460 cells exposed to hypoxia (1% O₂ for 24 h) for HIF1A peaks called in all four cell lines (Common) and peaks called only in HCT116 cells. Rows are sorted by HCT116 enrichment signal. f Comparison of enrichment signal in hypoxic HCT116 vs. H460 cells for HIF1A peaks called in one or both cell types. Colors highlight common (gray) and cell type-specific peaks (green, red) as well as peaks associated with genes with conserved upregulation in all four cell lines (core genes, blue). g Distributions of enrichment signal in hypoxic HCT116 cells for proximal HIF1A peaks associated with acute response genes, separated by conservation group as in (a): core (pink), shared 3 (green), shared 2 (blue), and HCT116 only (purple). Horizontal spread of data points is proportional to density; boxes indicate medians and upper and lower quartiles; horizontal bars with numbers indicate FDR-adjusted p-values for two-sided Mann–Whitney U tests. See also Supplementary Fig. 7.

Comparison of HIF1A peaks identified in each cell type indicates that, although HIF1A binding is consistently induced by hypoxia in all four cell types (Supplementary Fig. 7b) and consistently associated with underlying HREs (Supplementary Fig. 7c), there is vast cell type-specificity in HIF1A peak calls (Supplementary Fig. 7d). However, within this massive diversity, common HIF1A peaks show stronger HIF1A binding (Fig. 5e) and those associated with core genes are mostly conserved across cell types (Fig. 5f) and display significantly higher HIF1A signals than peaks associated with non-core genes (Fig. 5g). Furthermore, HIF1A peaks are associated with greater fractions of core upregulated genes in each cell type than less-conserved or cell type-specific genes (Supplementary Fig. 7e).

Overall, these results reveal that acute response genes are more conserved across cancer cell types than late response genes, commonly associated with high-occupancy HIF1A binding sites, and thus more likely to play a role in the cellular adaptation to hypoxia across multiple cancer types.

A tumor-suppressive role for HIF1A linked to mTOR suppression

Having identified the acute transcriptional response to hypoxia, which is more conserved than the late response across several cancer cell lines and also more clearly associated with strong and conserved HIF1A binding events, we then set out to investigate the contribution of these acute response genes to HIF1A-dependent processes in cancer cells. First, we analyzed genome-wide CRISPR-mediated knockout data for 625 cancer cell lines from the Cancer Dependency Map project³⁷. In this data set, depletion of tumor suppressor genes leads to improved cell fitness and positive gene effect scores (e.g. RB1), while depletion of oncogenes has the opposite effect (e.g. KRAS) (Fig. 6a). Gene effect scores for HIF1A are positive in most cell types with a distribution similar to that of RB1, whereas the HIF repressors HIF1AN, EGLN1, and VHL clearly promote cell fitness (Fig. 6a). Therefore, in the context of in vitro cell culture under normoxia and nutrient-rich conditions, HIF1A behaves as a strong tumor suppressor through cell-autonomous mechanisms. Furthermore, this reveals a clear cellular role for HIF1A under normoxia, likely explained by the existence of additional HIF1A-activating stimuli such as growth factor signaling and genetic alterations in cancer cells^2,3.

a Gene effect score distributions for HIF1A and select known regulators, alongside the tumor suppressor RB1 and oncogene KRAS for comparison, Horizontal spread of data points is proportional to density; boxes indicate medians and upper and lower quartiles; dashed line indicates threshold below which a gene is considered essential. b Ranked Spearman correlation coefficients of gene effect scores for each gene against HIF1A, with selected genes labeled, and acute response genes highlighted in red. * denotes significant correlations (FDR < 10%) — see Supplementary Data 7 for exact FDR values. c Comparison of gene effect scores with select genes for HIF1A. Points representing cell lines are colored by density and dashed lines represent linear fits to the data with 95% confidence intervals in gray. d Ranked Spearman correlation coefficients of gene effect scores for each gene against VHL, with selected genes labeled, and acute response genes highlighted in red. * denotes significant correlations (FDR < 10%) — see Supplementary Data 7 for exact FDR values. e Comparison of gene effect scores with select genes for VHL. Points representing cell lines are colored by density and dashed lines represent linear fits to the data with 95% confidence intervals in gray. f Ranked Spearman correlation coefficients of gene effect scores for each gene against DDIT4, with selected genes labeled, and acute response genes highlighted in red. * denotes significant correlations (FDR < 10%) — see Supplementary Data 7 for exact FDR values. g Comparison of gene effect scores with select genes for DDIT4. Points representing cell lines are colored by density and dashed lines represent linear fits to the data with 95% confidence intervals in gray. See also Supplementary Fig. 8 and Supplementary Data 7.

Next, to assess the contribution of acute response genes to HIF1A suppressive effects in vitro, we analyzed co-dependency relationships in this large data set (Supplementary Table 7). Reassuringly, the HIF1A cofactor ARNT shows a strong positive co-dependency with HIF1A, while the HIF1A repressors EGLN1-3, VHL, and HIF1AN/FIH exhibit strong negative relationships (Fig. 6b, c, Supplementary Fig. 8a), confirming the validity of this approach for identifying bona fide functional relationships. Conversely, for VHL the strongest negative correlation is with HIF1A, while the HIF1A repressors EGLN1 and EGLN3 as well as the VHL cofactors ELOB and CUL2 exhibit positive correlations. (Fig. 6d, e, Supplementary Fig. 8b).

Notably, the strongest positive correlation for HIF1A is DDIT4 (REDD1) (Fig. 6b, c), a core acute response gene (Supplementary Fig. 8c) previously characterized as a direct HIF1A target³⁸. DDIT4 was previously shown to repress mTOR signaling through the TSC1-2 complex³⁹, a notion that is recapitulated by our co-dependency analysis, whereby DDIT4 shows positive correlations with TSC1-2 and negative correlations with MTOR and LAMTOR1 (Fig. 6f, g, Supplementary Fig. 8d). There is also a strong positive relationship between HIF1A and FOXO3 (Fig. 6b, c, Supplementary Fig. 8a), another acute response gene that encodes a transcription factor with central roles in multiple stress responses⁴⁰, and a known repressor of mTORC1 (ref. ⁴¹).

Altogether, these results point to cell-autonomous anti-proliferative effects for HIF1A before the onset of hypoxia, likely through suppression of mTOR signaling via induction of DDIT4, a target gene sensitive to basal levels of HIF1A (Fig. 1h, i).

Acute response genes associated with cancer progression

Having shown that acute transcriptional output and early fold changes in transcriptional activity have proportional impacts on steady-state mRNA levels at late time points of hypoxia, we then set out to investigate the prognostic value of acute response genes through analysis of ~11,000 human tumors across 27 cancer types for which gene expression are available via The Cancer Genome Atlas (TCGA) and for which curated patient survival data have been derived⁴². We first tested for association between progression-free interval (PFI) and an aggregate expression score across all acute response genes present in this data set using an iterative log-rank approach to find the optimal split between samples with high and low scores (see Methods, Supplementary Table 8). Importantly, high expression scores were predominantly associated with decreased survival, reaching statistical significance (FDR < 10%) in six cancer types: stomach adenocarcinoma (STAD), adenoid cystic carcinoma (ACC), cervical squamous cell carcinoma (CESC), low-grade gliomas (LGG), kidney renal papillary cell carcinoma (KIRP), and urothelial bladder carcinoma (BLCA) (Fig. 7a, b, Supplementary Fig. 9a), indicating that high expression of hypoxia-inducible genes as a group is preferentially associated with worse outcome. Only one tumor type, skin cutaneous melanoma (SKCM), presented better prognosis when expressing high levels of acute response genes (Fig. 7a, Supplementary Fig. 9a). We then analyzed acute response genes individually (Supplementary Table 8), focusing on those significantly associated with PFI in multiple cancer types (FDR < 5%). High expression of individual acute response genes was more often associated with unfavorable than favorable prognosis (161 versus 107, respectively), with some genes exhibiting divergent associations in different cancer types (83 genes) (Supplementary Fig. 9b). Of note, MXI1, the acute response gene known to repress the MYC transcriptional program, was consistently associated with favorable prognosis, as would be expected for a suppressor of MYC signaling (Supplementary Fig. 9c, d).

a Iterative log-rank analysis of Acute Response Gene Score vs. progression-free interval (PFI): FDR-adjusted p-value vs. PFI ratio (mean survival time in high group/mean survival time in low group) for the indicated cancer types; red/green denotes significant association of high score with PFI (FDR < 10%). b Kaplan–Meier plots for indicated cancer types, separated by Acute Response Gene Score; split indicates the proportions of samples placed into high and low scoring groups, with initial numbers at risk per arm indicated at lower left. P-values are from log-rank analysis with the indicated split. c Top 10 functional enrichment clusters from Metascape pathway analysis among individual unfavorable acute response genes, limited by significance in increasing numbers of cancer types (left to right). Heatmap color represents −log₁₀(p-value) from hypergeometric enrichment test. d Heatmap showing significant adjusted p-values across cancer types for ECM-related genes associated with lower PFS. Heatmap color represents −log₁₀(p-value) from hypergeometric enrichment test. e Kaplan–Meier plots for selected individual genes in the indicated cancer types, separated by expression level; split indicates the proportions of samples placed into high and low expression groups, with initial numbers at risk per arm indicated at lower left or right. P-values are from log-rank analysis with the indicated split. f Selected ECM-related genes with significant (FDR < 10%) association with PFI using Cox regression analysis. Blue boxes indicate estimated hazard ratio, with box size proportional to FDR-adjusted p-value; bars indicate 95% confidence intervals. See also Supplementary Fig. 9 and Supplementary Data 8.

Importantly, when we analyzed the pathways enriched among acute response genes significantly associated with decreased PFI across increasing numbers of cancer types, we observed increased enrichment only of genes related to ECM remodeling (Fig. 7c, d). This subset of acute response genes included key enzymes involved in collagen remodeling (e.g. PLOD2, PLOD1, P4HA1, LOXL2), collagen subunits (COL1A1, COL13A1), and proteins previously involved in cell migration and invasion, such as SERPINE1 (ref. ⁴³), CD109 (ref. ⁴⁴), and calpastatin (CAST)⁴⁵. Some of these genes have previously been linked to hypoxia-induced ECM remodeling, local invasion, and/or metastasis, such as PLOD1/2, P4HA1, and LOXL2 (refs. ^3,46). When we repeated this analysis using Overall Survival (OS) as an endpoint instead of PFI, high expression of these ECM-related genes was also associated with shorter time to death (Supplementary Fig. 9e). Consistently, high expression of these genes in diverse tumor types is associated with poor prognosis (Fig. 7e, Supplementary Fig. 9d). Cox regression analysis with additional variables (see Methods) further confirmed association of these genes with adverse outcome (Fig. 7f, Supplementary Table 8). Notably, PFI and OS data may not be sufficiently mature for some cancer types in the TCGA data set, which prompts cautious interpretation of negative results in this analysis.

Overall, these results highlight the pleiotropic nature of the acute hypoxic response, with tumor-suppressive and oncogenic roles in different contexts through the action of diverse HIF1A targets, while also revealing a prominent role for hypoxia-induced ECM remodeling in cancer progression.