Extensive regulation of the non-coding transcriptome by hypoxia: role of HIF in releasing paused RNApol2

Hypoxia is central to both ischaemic and neoplastic diseases. However, the non-coding transcriptional response to hypoxia is largely uncharacterized. We undertook integrated genomic analyses of both non-coding and coding transcripts using massively parallel sequencing and interfaced this data with pan-genomic analyses of hypoxia-inducible factor (HIF) and RNApol2 binding in hypoxic cells. These analyses revealed that all classes of RNA are profoundly regulated by hypoxia and implicated HIF as a major direct regulator of both the non-coding and coding transcriptome, acting predominantly through release of pre-bound promoter-paused RNApol2. These findings indicate that the transcriptional response to hypoxia is substantially more extensive than previously considered.


Introduction
Cells respond to physiological insults by altering their transcriptional output. Hypoxia (low oxygen level) is an important environmental stress with a central role in many physiological responses including adaptation to altitude, exercise, growth and development, as well as major pathophysiological processes such as ischaemic vascular disease, inflammation, wound healing and cancer [1]. The transcription factor, hypoxia-inducible factor (HIF) orchestrates many of these responses to limit the disturbance of oxygen homeostasis or promote repair processes, through the transactivation of protein-coding genes with key roles in pathways such as apoptosis, differentiation, proliferation, energy metabolism and growth factor production [2]. Given the fundamental role of disturbed oxygen homeostasis in human disease and the potential for therapeutic manipulation of the HIF pathway, much interest has focused on understanding the extent and architecture of these pathways.
To date, pan-genomic analyses of the transcriptional response to hypoxia using microarrays have focussed on protein-coding RNAs (mRNAs) and some microRNAs (miRNAs) [3,4]. However, recent genomic analyses have transformed our perspective of regulatory transcriptional networks. It is now recognised that mRNA forms only part of a much broader transcriptional output that includes even greater numbers of RNAs that are not translated to protein [5]. These 'non-coding' RNAs have important regulatory roles that further shape the transcriptional output of the cell. They include short non-coding RNA (<200 nucleotides) such as micro (miRNA), small nuclear/nucleolar (snRNA), piwi-interacting RNAs (piwiRNAs), and transfer RNAs (tRNAs) as well as long non-coding RNAs (lncRNA) (>200 nucleotides). Recent work has demonstrated altered noncoding RNA expression in many types of cancer and revealed functions in cell cycle regulation, apoptosis, carcinogenesis and metastasis [5,6]. However, little is known about the mechanisms regulating the expression of the non-coding transcriptome.
Here, we provide the first pan-genomic analysis of both non-coding and coding transcriptional responses to hypoxia, including lncRNAs, miRNAs, piwiRNAs, snRNAs, and tRNAs. The analysis reveals marked bidirectional hypoxia-inducible changes in transcript abundance across all classes of RNA and demonstrates the importance of HIF in regulating non-coding as well as coding transcriptional responses to hypoxia. Correlation with changes in histone marks and RNA polymerase 2 (RNApol2) positioning, indicates that activation of gene expression by HIF occurs commonly through long-range interactions that cause the release of pre-bound, promoter-paused RNApol2.

Results and Discussion
Extensive regulation of the coding and non-coding transcriptome by hypoxia To gain a global overview of the transcriptional response to hypoxia, we integrated expression profile and chromatin maps of human breast cancer MCF-7 cells incubated in either 21% ambient oxygen (normoxia) or 1% ambient oxygen (hypoxia) (Fig 1A). Independent directional and nondirectional cDNA libraries were generated from polyadenylated RNA (polyA+) and total RNA depleted of ribosomal RNAs (ribo-) and analysed by high-throughput sequencing. For transcripts present in both data sets, expression levels determined by the two methods showed a high degree of correlation (supplementary Fig S1), although a subset of transcripts showed incomplete adenylation [7].
To determine transcript abundance across all classes of RNA, the ribosome depleted, directional RNA-seq dataset was mapped sequentially to databases of different RNA classes (Fig 1B). 42989 non-ribosomal transcripts were detected in normoxia and 43169 in hypoxia. The raw data for each class is presented in Fig 1C, whilst Fig 1D shows the median fold-change in hypoxia, together with the interquartile and full ranges, after exclusion of low-abundance transcripts where this ratio might be unreliable. This demonstrates that all classes of RNA are regulated by hypoxia, but reveals major classspecific differences. Several classes of transcript (snRNAs, tRNAs and piwiRNAs) demonstrated overall downregulation, whereas mRNAs, lncRNAs and miRNAs demonstrated overall upregulation. To confirm these findings, representative samples of snRNA, mRNA and lncRNA transcripts exhibiting regulation close to the median values were tested by quantitative (q)PCR. This analysis recapitulated differences between the classes, consistent with the finding of overall class-specific regulation (supplementary Fig S2). In addition to these global changes in hypoxia, a number of transcripts in each category show strong up-or downregulation, ranging from 15-fold up to 8-fold down (piwiRNA), >100-fold up to 30-fold down (miRNA) and 4-fold up to 12-fold down (lncRNA) (supplementary Table S1).
Several microarray-based pan-genomic analyses have described heterogeneous responses in the miRNA response to hypoxia depending on cell type, hypoxic stimulus and profiling platformreviewed in McCormick et al [4]. A number of these miRNAs (including mir-210) were also increased in our study (supplementary Table S2). In addition, upregulated miRNAs from our analysis were also induced across a panel of breast cancer cell lines (supplementary Fig S3).
As few lncRNAs have been included in previous pan-genomic analyses our data identifies lncRNAs as a new class of transcript that are commonly regulated by hypoxia. The majority of these transcripts remain functionally uncharacterized. However, two lncRNAs (NEAT1 and MALAT1) included in a previous microarray-based analysis [8] were the most strongly upregulated lncRNAs in our dataset and were also upregulated in a panel of breast cancer cell lines (supplementary Fig S4). NEAT1 is a component of interchromatin-paraspeckles, which regulate nuclear retention of Adenosine to Inosine (A-to-I) edited RNAs [9]. MALAT1 promotes cellular proliferation, is associated with tumour growth, metastasis and recurrence and is frequently upregulated or mutated in solid tumours [10]. In addition, H19, another hypoxia-inducible oncogenic lncRNA [11] was also induced in our analysis but did not have a HIF-binding site close by.

Hypoxic induction of non-annotated intergenic and antisense transcripts
A proportion of reads did not map to an annotated transcript of any class. To define whether these represent previously non-annotated A RNA-seq was performed following 24 h culture in 21% (normoxia) or 1% (hypoxia) ambient oxygen. B Outline of the pipeline for mapping ribosomal depleted directional RNA-seq reads. C Raw data plotted as log2 fold-change by hypoxia on the vertical axis versus expression level on the horizontal axis. Vertical lines denote thresholds (higher for classes of longer RNAs) to select transcripts for further analysis. D Box-and-whisker plots of log2 fold-change by hypoxia for filtered transcripts in each class of RNA. The vertical dotted line denotes no foldchange.
transcripts, they were first examined for contiguity, using the Cufflinks de novo transcript assembler. This identified 91 novel RNA transcripts with a median length of 9988 bp (range 269-38 2549 bp) many of which were upregulated in hypoxia (supplementary Tables  S3 and S4) and were confirmed by qPCR (supplementary Fig S5).
In keeping with assignment as bona fide non-annotated transcripts, reference to UCSC data (http://genome.ucsc.edu/cgi-bin/ hgTrackUi?hgsid=336849981&g=cpgIslandExt) revealed that 50 (55%) have a CpG Island within 1000 bp of the putative promoter. All 91 transcripts were also detected in the independently derived polyA+ RNA-seq dataset. Of the 91 non-annotated transcripts, 37 did not overlap with any annotated RefSeq gene, and were therefore classified as non-annotated intergenic transcripts (supplementary Table S3). All had a low coding potential [12], indicating that they belong to the lncRNA class of RNA. The remaining 54 overlapped with a previously documented RNA, but were expressed from the opposite DNA strand indicating that they are non-annotated antisense transcripts (NATS) (supplementary Table S4). Anti-sense transcripts may act in cis to regulate expression of the overlapping sense transcript [13]. Interestingly, the majority were upregulated by hypoxia, induction being commonly associated either with counterregulation (e.g. HIF-1a, TBX2) or with co-regulation of the overlapping sense transcript (e.g. SPAG4, CELSR2) (supplementary Fig S6). Thus, hypoxia-induced antisense transcripts are allied to both induction and repression in cis adding further tiers to regulation of the transcriptome by hypoxia.
Transcriptional regulation of long non-coding RNAs by HIF As HIF is a major transcriptional regulator of the mRNA response to hypoxia [14], we next examined the role of HIF in the regulation of other classes of RNA. As a first step we examined for spatial association between each class of transcript and HIF-binding sites, using well validated, previously published HIF-1-and HIF-2-binding sites in MCF-7 cells [14]. That work considered only protein coding genes and reported a non-random distribution of HIF binding across the genome. Although HIF-binding sites were strongly enriched at the promoters of coding genes, many HIF-binding sites were observed to be remote from these promoters. In the light of extensive regulation of non-coding RNAs by hypoxia, we reanalysed the distribution of HIF-binding in relation to the location of all transcripts. The number of HIF-binding transcripts (defined as the closest promoter to each HIF-binding site) in each category broadly mirrored the total number in each class. Although the numbers in many classes were too small to permit firm conclusions, no enrichment of HIF binding was observed in the vicinity of classes of RNA showing global downregulation (piwiRNAs, snRNAs and tRNAs) suggesting that they are unlikely to be directly regulated by HIF. In contrast, HIF binding was enriched in the vicinity of both mRNAs and lncRNAs (Fig 2A) with HIF-2 binding a slightly higher proportion of lncRNAs than HIF-1. Nevertheless, despite promoter enrichment, the majority of HIF-1-and HIF-2-binding sites (50% and 70% respectively) were still found to lie in excess of 2.5-kb from the promoter of an expressed transcript of any class (Fig 2B and C).
To assess the functional role of HIF in the regulation of different classes of RNA by hypoxia, we next examined for association of HIF binding with regulation by hypoxia across all classes of RNA using gene set enrichment analysis (GSEA) [15]. All expressed transcripts were pooled and ranked according to their hypoxic induction. HIF binding was strongly enriched amongst transcripts upregulated (but not downregulated) by hypoxia (Fig 2D). Only mRNA, lncRNA, miRNA and non-annotated transcripts were present in the coreenrichment group of HIF-binding hypoxia-inducible transcripts. GSEA analysis by subgroup confirmed enrichment for mRNAs ( Fig 2E) and lncRNAs (Fig 2F), but the number of miRNAs and non-annotated RNAs was too small to permit statistical analysis. Taken together, this establishes a non-random spatial association between HIF-binding sites and hypoxia-inducible transcripts for both mRNAs and lncRNAs, pointing to a major new role for HIF in the transcriptional regulation of lncRNAs.
To pursue this further, we undertook polyA+ RNA-seq analysis of mRNA and lncRNA in hypoxic MCF-7 cells transfected with siRNAs targeting HIF-1a, HIF-2a or both isoforms simultaneously. Approximately 25% of all lncRNA transcripts detected in the ribosomal-depleted RNA-seq were detected in this analysis. Fig 3A  portrays the results of HIF-a siRNA interventions on all HIF-binding lncRNA transcripts detected in this dataset. As expected, the heatmap confirms that many of these RNAs were upregulated by hypoxia and reveals a strong correlation with down regulation by siRNAs targeting one or both HIF-a isoforms, a result that was confirmed by GSEA (supplementary Fig S7). Individual datasets are illustrated for MALAT1 (downregulated by both HIF-1a and HIF-2a siRNA, Fig 3B) and NEAT1 (downregulated by HIF-2a but not by HIF-1a siRNA, Fig 3C). Thus, the presence of HIF binding, hypoxic upregulation and subsequent downregulation by the suppression of HIF-a defines lncRNAs as a new group of HIF target genes.

HIF-dependent recruitment and release of promoter-paused RNApol2 in hypoxia
To pursue the mechanism of transcriptional activation by hypoxia further, we next examined RNApol2 binding and histone H3K4me3 modification in normoxia and hypoxia. As expected, the promoterassociated abundance of each mark correlated with the basal transcript level. Strong peaks of RNApol2 and H3K4me3 signal were also seen at the putative transcriptional start sites (TSS) of the previously non-annotated transcripts consistent with assignment as bona fide genes (supplementary Fig S8) [16,17].
Consideration of changes in RNApol2 binding in hypoxia revealed that for the majority of hypoxia-inducible genes, RNApol2 was already bound at the promoter in normoxic cells. Although there were exceptions, at the majority of promoter sites hypoxia did not increase promoter-bound RNApol2 substantially. This is illustrated in Fig 4A which depicts changes in RNApol2 ChIP signal in the vicinity of the promoters of the hundred most strongly hypoxia upregulated transcripts of any class. In hypoxia, little increase in loading of RNApol2 at the transcriptional start was seen; instead an increase in RNApol2 is observed across the body of the gene. Similar patterns of change in hypoxia were seen for both mRNAs and lncRNAs (supplementary Fig S9). Actively transcribing RNApol2 can recruit histone methyl transferase activity that trimethylates H3K4 either directly [18] or through p300 [19]. The distribution of H3K4me3 showed a very similar response to that of RNApol2, with little change in signal at the promoter in hypoxia, but an increase downstream of the transcriptional start site, consistent with the hypothesis that RNApol2 recruits histone methyl transferases that trimethylate H3K4 (Fig 4B). In contrast, consideration of DNAse hypersensitivity analyses in normoxic and hypoxic MCF7 cells reported by the Encyclopedia of DNA Elements (ENCODE) consortium [17] revealed the expected hypersensitivity signals at these promoter regions, but did not reveal changes in hypoxia (Fig 4C).
On close inspection the maximum RNApol2 signal was observed in the promoter-proximal region, 40-45 bp downstream of the transcriptional start. This suggests that promoter-proximal pausing is an important rate-limiting step for transcription of hypoxia-inducible genes and that reversal of pausing is a key step through which hypoxia regulates transcription [20]. To determine whether this is a feature of direct transcriptional regulation by HIF, as opposed to indirect effects of HIF or hypoxia itself on transcription, we performed similar analysis of RNApol2 distribution for genes that were downregulated by hypoxia and which (consistent with indirect regulation) do not show promoter enrichment of HIF binding. We also compared hypoxia-inducible changes in patterns of RNApol2 binding and H3K4 trimethylation amongst groups of genes that were matched for the extent of upregulation by hypoxia, but were selected as being either the closest gene to a HIF-binding site (likely to represent direct transcriptional regulation) or were located >1 Mb from the nearest HIF-binding site (likely to represent indirect actions of HIF or hypoxia). Whilst a clear change in the ratio of promoter paused RNApol2 to that distributed along the gene body (travelling ratio) was observed for genes predicted to be direct HIF transcriptional targets this was not observed for genes predicted to be upregulated or downregulated indirectly by HIF or by hypoxia itself (supplementary Fig S10). Finally, ChIP-qPCR at a selection of these genes confirmed that the hypoxic induction of RNApol2 within the body of the gene was dependent on the presence of HIF-a subunits (Fig 4D-F).
Amongst genes predicted to be direct transcriptional targets of HIF we observed a clear change in the distribution of travelling ratios in hypoxia, reflecting release of promoter-paused RNApol2 at most loci ( Supplementary Fig S9E, F and Fig 4G). However, we also observed infrequent HIF target genes where RNApol2 was recruited de novo in hypoxia (Fig 4H). No association was detected between this type of behaviour and biological function. However, the promoters of these transcripts also show constitutive DNAse1 HSS signal indicating that other factors must limit normoxic RNApol2 binding at these sites. As DNA-seq cannot distinguish individual nucleosome positions, this might include mechanisms such as nucleosome positioning as well as DNA/histone methylation [21]. Taken together, our findings indicate that the regulation of gene expression by hypoxia is substantially more complex than has previously been considered. All RNA classes were regulated bi-directionally, our analysis identifying both general responses characteristic of particular types of RNA and high-amplitude effects on individual species. Both the non-random spatial association between HIF-binding sites and hypoxia-inducible genes and direct functional intervention by siRNA implicated HIF in the regulation of extensive networks of non-coding RNAs. RNApol2 profiling indicates that this direct regulation by HIF is mediated through activation of pre-bound RNApol2. These findings confirm and extend the role of HIF in the regulation of miRNAs and identify both intergenic and antisense lncRNAs as a new target class for the HIF transcriptional response. The ability of non-coding RNAs to regulate gene expression both in trans and in cis suggests the hitherto unrecognised involvement of HIF in a far broader network of gene regulation. Given the fundamental role of hypoxia in human pathology and emerging data on the importance of non-coding RNAs in cancer and other diseases, these findings open new avenues to the better understanding of human disease processes.
High throughput sequencing All libraries were sequenced using the GAII or HiSeq platforms (Illumina) according to the manufacturers protocol. Raw and mapped data for RNA-seq and ChIP-seq is available at EMBL-EBI Array Express (E-MTAB-1994, E-MTAB-1995).

Statistical analysis
All assays were performed in duplicate or triplicate. Statistical analyses using two-tailed t-tests were performed in SPSS (IBM Corp, New York, NY, USA) and R (http://www.R-project.org).