RNA-seq analysis of PHD and VHL inhibitors reveals differences and similarities to the hypoxia response.

Background: Hypoxia-inducible factor (HIF) transcription factors are well known to control the transcriptional response to hypoxia. Given the importance of cellular response to hypoxia, a number of pharmacological agents to interfere with this pathway have been developed and entered pre-clinical or clinical trial phases. However, how similar or divergent the transcriptional response elicited by different points of interference in cells is currently unknown. Methods: We performed RNA-sequencing to analyse the similarities and differences of transcriptional response in HeLa cells treated with hypoxia or chemical agents that stabilise HIF by inhibiting components of the hypoxia signalling pathway – prolyl hydroxylase (PHD) inhibitor or von Hippel–Lindau (VHL) inhibitor. Results: This analysis revealed that hypoxia produces the highest changes in gene transcription, with activation and repression of genes being in large numbers. Treatment with the PHD inhibitor IOX2 or the VHL inhibitor VH032 led mostly to gene activation, majorly via a HIF-dependent manner. These results were also confirmed by qRT-PCR using more specific and/or efficient inhibitors, FG-4592 (PHDs) and VH298 (VHL). Conclusion: PHD inhibition and VHL inhibition mimic gene activation promoted by hypoxia via a HIF-dependent manner. However, gene repression is mostly associated with the hypoxia response and not common to the response elicited by inhibitors of the pathway.


Introduction
Hypoxia, or reduced oxygen availability, is associated with many physiological processes, such as embryonic development and high altitude living but also pathological processes such as stroke and cancer (Rocha, 2007). A major regulator of oxygen sensing and response is the family of transcription factors called hypoxia-inducible factors (HIFs). HIFs are activated in response to hypoxia to initiate a transcriptional program, and ultimately restore oxygen homeostasis and promote cell survival (Kenneth & Rocha, 2008). HIFs are heterodimeric transcription factors composed of a constitutively stable β-subunit (HIF-1β) and an oxygen-labile α-subunit (HIF-α) (Kenneth & Rocha, 2008). HIF-α is encoded by three different genes: HIF-1α, HIF-2α and HIF-3α, and each could function differently depending on tissue localisation (Kenneth & Rocha, 2008). HIF-α is rapidly degraded by the proteasome under normal oxygen levels as prolyl hydroxylase domain (PHD) enzymes and factor inhibiting HIF (FIH) utilise molecular oxygen as a co-factor, in addition to Fe 2+ and 2-oxoglutarate, to hydroxylate HIF-α proteins (Kenneth & Rocha, 2008). Hydroxylated prolines on HIF-α create a recognition site with a substantial increase in affinity over the parent protein containing unmodified proline, for the E3 ubiquitin ligase, von Hippel-Lindau (VHL) tumour suppressor that poly-ubiquitinates HIF-α, targeting the protein for proteasomal degradation (Hon et al., 2002). In hypoxia, however, HIF-α evades degradation and is stabilised as a result of insufficient oxygen molecules for PHDs to function. The accumulated HIF-α dimerises with HIF-1β and binds to the consensus motif hypoxia response elements (HREs) of HIF target genes to activate the expression of a wide range of genes associated with key biological processes including metabolism, angiogenesis, cell differentiation, apoptosis and autophagy, for adaptation to hypoxia (Liu et al., 2012).
Considering the pharmacological use and therapeutic potential of PHD inhibitors and the newly emerging VHL inhibitors, it is important to identify gene expression responses elicited by such agents. As such, we employed RNA-sequencing (RNA-seq) to determine the similarities and differences of the transcriptional response under hypoxia, the inhibitor of PHD, IOX2 (Chowdhury et al., 2013), as well as the VHL inhibitor VH032 (Frost et al., 2016;Galdeano et al., 2014). We show that IOX2 and VH032 mimic the hypoxia response and that these predominantly induce a HIF-dependent gene signature. On the other hand, hypoxia produces the broader transcriptional response amongst all the inducers used, with significant numbers of genes being induced and repressed.

Treatments
For hypoxia induction, cells were incubated at 1% O 2 in an InVIVO 300 hypoxia workstation (Ruskin Technologies). To prevent reoxygenation, cells were lysed for protein or RNA extraction in the hypoxia workstation. DMSO was used as vehicle control for compound treatment. PHD inhibitors IOX2 and FG-4592 were purchased from from Sigma (SML0652) or Selleckchem (S2919) and Selleckchem (S1007), respectively. Drugs were added to cells for the indicated length of time. VHL inhibitors VH032 and VH298 were synthesised by Pedro Soares (Ciulli lab, University of Dundee) as previously described: VH032 (ligand 7 in Galdeano et al. (2014); compound 1 in (Soares et al., 2018) andVH298 (Frost et al., 2016). VH298 was also purchased from Sigma (SML1896).
RNA preparation for RNA-seq HeLa cells were seeded in 35 mm dishes one day prior to treatments with 0.05% DMSO, hypoxia (1% O 2 ), 250 µM IOX2, or 250 µM VH032 for 16 h. RNA was extracted using the RNeasy Mini Kit (Qiagen; 74104) according to manufacturer's instruction. Genomic DNA was removed from RNA samples using RNase-free DNase from Qiagen (79254) as per manufacturer's protocol. Experiments were performed in triplicates.

RNA-seq library preparation and sequencing
Library preparation and sequencing were performed by the University of Dundee Genome Sequencing Unit. The library was prepared using TruSeq Stranded Total RNA Library Preparation Kit with Ribo-Zero TM Human/Mouse/Rat kit (Illumina; RS-122 2201) to remove ribosomal RNA (rRNA). RNA ERCC ExFold RNA Spike-In Mix (Mix1 and Mix2) was distributed throughout the RNA-seq experiment according to manufacturer's protocol (4456739, Thermo Scientific). Paired-end Illumina sequencing was performed on the NextSeq 500 platform.

RNA-seq data analysis
The raw sequence reads from each replicate were aligned to the Ensembl human genome GRCh37 and ERCC sequence with STAR version 2.4.2a. The aligned reads were combined and number of reads for each gene was counted with subread-featureCounts pipeline version 1.4.6-p4 (Liao et al., 2014). The files were found to contain ribosomal DNA (rDNA) contaminations, the majority of which were the following two mitochondrial DNA: ENSG00000211459 and ENSG00000210082 -which were removed manually. Differential gene expression analysis was performed by the R package edgeR (v3.24.1) according to its user guide (Robinson et al., 2010), and differentially expressed genes were identified at FDR of <0.05 and log2 fold change > 0.58.
Integrative analysis was performed manually in R (v1.1.453) to obtain lists of genes that overlap to publicly available datasets. Briefly, a list of differentially expressed genes upregulated in hypoxia, IOX2, VH032 or in all three conditions was compared to publicly available data and overlapping genes were exported into excel sheet. Enrichment analysis of transcription factors and chromatin binding proteins on our datasets was carried out using TFEA.ChIP (v1.2.2) according to its user guide (Puente-Santamaria & Del Peso, 2019). Gene set enrichment analysis was performed using GSEA MSigDB online tool (Liberzon et al., 2011;Subramanian et al., 2005) for hallmark genes with FDR < 0.05 and p value < 0.05.
Sequence data are available from Gene Expression Omnibus GSE120675.

IOX2 and VH032 induce a similar transcriptional response profile, while hypoxia induces a broader response in cells.
The hypoxia inducible factors (HIFs) can be induced in a variety of different ways, from the physiological inducer of low oxygen, to the pharmacological inhibition of proteins involved in the HIF pathway, as well as by changes in alternative pathways such as transcription and translation (Ferreira et al., 2015;Moniz et al., 2015;Semenza, 2003).
To understand the similarities and differences between the transcriptional responses to several HIF inducers, an unbiased high-throughput RNA-sequencing (RNA-seq) was performed. Human cervical cancer HeLa cells were exposed to 0.05% DMSO (vehicle control), hypoxia (1% O 2 ), PHD inhibitor IOX2 or VHL inhibitor VH032 for 16 hours prior to profiling for global transcriptomic analysis using RNA seq. Differential expression analysis of data collected from RNAseq showed that DMSO samples cluster together with weak correlation to the other treatments, whilst hypoxia, IOX2 and VH032 conditions were grouped close to one another ( Figure 1A-B). Furthermore, heatmaps generated using the top 100 most differentially expressed (DE) genes in each experimental condition comparing to DMSO control showed that hypoxia ( Figure 1C), VH032 ( Figure 1D) and IOX2 ( Figure 1E) displayed The Pearson correlation was used to compute distances between genes and samples, and the clustering was performed using average linkage. Each column corresponds to a sample and each row corresponds to a specific gene. (F) Heatmaps of Pearson correlations between replicates of the same conditions. Each data had been normalised to their total counts. (G) Each dot represents a differentially expressed gene comparing the condition stated in the heading legend to DMSO vehicle control. Blue dots represent genes with increased expression (logFC > 0.58; to the right) or decreased expression (logFC < -0.58; to the left) at false discovery rate (FDR)<0.05. Blue triangles (present at -log10 FDR of 30) represent genes with logFC > 0.58 or < -0.58 and -log10 FDR > 30. similar transcriptional profiles to each other and were noticeably distinct from the vehicle control DMSO. These observations are likely due to the activation of HIF pathway as the three treatments activate HIF. Correlation analysis heatmaps for each condition showed a strong correlation between 0.95 and 1 across replicates of the same experimental condition ( Figure 1F), and replicates of each treatment cluster together in heatmaps of top 100 most DE genes ( Figure 1C-E), demonstrating that replicates of each condition were similar and statistically close to each other.
To investigate the nature of the transcriptional data we obtained, we performed integrative analysis of our hypoxia dataset with publicly available hypoxia-inducible gene sets (Table 1, Dataset 1 (Frost, 2019)). Overlap analysis showed that 36% (410 out the 1,133) genes upregulated in hypoxia were present in at least one of the reported datasets, with 115 genes found upregulated in 16 cell lines (Ortiz-Barahona et al., 2010), 129 genes in HeLa dataset (Mense et al., 2006), as well as 75 and 307 genes in in MCF7 cells (Chan et al., 2016;Elvidge et al., 2006). This analysis confirmed the cell-specific and time-dependent transcriptional responses elicited by hypoxia exposure in cells. We also compared our IOX2 and VH032 datasets with the previously reported gene sets to assess the extent to which genes upregulated in IOX2 or VH032 also hypoxia-inducible. We identified a large portion of genes upregulated by IOX2 (39%; 325 out of 827) or VH032 (56%; 200 out of 362) to be present in at least one of these reported hypoxia datasets (Table 1, Dataset 1 (Frost, 2019)). This analysis showed that VH032 is predominantly regulating hypoxia-inducible genes, consistent with specific on-target effects on VHL (Frost et al., 2016).
Genes upregulated in hypoxia, IOX2 and VH032 are HIFdependent Comparative analysis of upregulated genes distinctly showed that the majority of VH032-induced genes (~87%; 315 out of 362) are also upregulated by hypoxia (Figure 2A). On the other hand, IOX2-induced genes are only partially shared with the hypoxia signature (~68%; 559 out of 827). Notably, nearly all of the VH032-induced genes (93%) are shared with IOX2. Overall, 306 genes are upregulated in all of the three experimental conditions (Figure 2A; Dataset 1 (Frost, 2019)).
Given that hypoxia, IOX2, and VH032 all induce HIF transcriptional activity, we next investigated the extent to which these 306 overlapped genes upregulated in all three conditions were regulated by HIF transcription factors. We performed integrative analysis on the overlapped genes with reported datasets of validated  (Frost, 2019)). This analysis revealed that 132 out of these 306 (43%) genes were HIF-dependent (Table 2, Dataset 1 (Frost, 2019)). A total of 33 out of the 306 shared genes was present in the list of 93 validated HIF-1 target genes (Dataset 1 (Frost, 2019)). Analysis using MCF7 ChIP-seq dataset showed that 62 (20%) and 33 (11%) of the 306 upregulated genes contained HIF-1α and HIF-2α binding sites, respectively (Dataset 1 (Frost, 2019)). A higher level of overlap was observed when we analysed the HepG2 ChIP-seq dataset, revealing that 90 out of the 306 genes (29%) contained either HIF-1α or HIF-2α binding sites (Dataset 1 (Frost, 2019)).
We next utilised TFEA.ChIP that exploits publicly available ChIP-seq datasets to perform enrichment analysis of transcription   Figure 2B).
Taken together, these comparative analyses demonstrate the level of HIF dependency for genes upregulated by both hypoxia, IOX2 and VH032.
To investigate the cellular processes induced by hypoxia, IOX2 or VH032, gene set enrichment analysis (GSEA) was performed according to the molecular signature database (MSigDB) (Subramanian et al., 2005;Liberzon et al., 2011). All treatments induced a similar set of enrichment for genes involved in the "cellular response to hypoxia", "glycolysis", "epithelial-mesenchymal transition", "mTORC1 signalling" and "NF-κB signalling" ( Figure 2C). However, genes repressed by the different treatments mapped to quite diverse cellular pathways and responses ( Figure 2D). Furthermore, the group of 306 commonly upregulated genes in all three conditions was enriched with genes found in these same hallmarks ( Figure 2B) Figure 3A). We replaced IOX2 with the PHD inhibitor FG-4592 that is currently in clinical trial phase III (Provenzano et al., 2016). Furthermore, VH032 was replaced with the more potent VHL inhibitor VH298 (Frost et al., 2016). After exposure to 16 hours of hypoxia, FG-4592 (50 µM) or VH298 (100 µM), mRNA levels of these genes increased significantly with similar induction profiles in both HeLa ( Figure 3B) and HFF cells ( Figure 3C). Hypoxia showed the strongest induction profiles in nearly all genes examined in both cell lines (Figure 3). Moreover, the changes in transcript levels were also reflected at the protein level ( Figure 4). Accumulation of the products of these HIF target genes, as well as GLUT1 protein, another HIF target which we had previously characterised at mRNA level (Frost et al., 2016), was detected following 24 hour treatment of hypoxia, VH298 or FG-4592 in HeLa and HFF cells ( Figure 4A). In both cell lines, the three treatments induced similar levels of protein expression for the majority of HIF targets assessed, with NDRG1 being the most prominent in hypoxia. FG-4592 was the strongest inducer of CA9 in HeLa cells and BNIP3 in both cell lines ( Figure 4A). The difference in the levels these genes and therefore proteins were increased could be due to the single time point used; particularly since the three conditions act on the HIF pathway at different stages (1% oxygen level limits the activities of PHD and FIH, FG 4592 inhibits PHDs and VH298 inhibits VHL downstream of hydroxylation by PHD). To address this question, we performed a time course analysis for the three inducers and investigated protein levels of the different HIF-target genes ( Figure 4B). This revealed that hypoxia is the strongest inducer of all the proteins we have analysed at the 24 h post-treatment time point.

RNA-seq validation, genes solely upregulated in hypoxia and IOX2.
Hypoxia and IOX2 share the larger overlap of 252 upregulated genes that are not found in VH032 (Dataset 1 (Frost, 2019)). On the other hand, there are 30 upregulated common genes between IOX2 and VH032, but not hypoxia, as well as the 9 common upregulated genes in hypoxia and VH032, but not IOX2 ( Figure 2A). As recent studies have revealed additional targets of PHD  enzymes, we analysed several of these 252 genes to determine whether PHDs induce transcriptomic changes independent of HIF activity. We selected four genes, including IDH2, RNF187, FAM117B and JMJD6 from the list of 252 genes upregulated solely in hypoxia and IOX2 for validation by qRT-PCR. The results, however, show that mRNA levels of these genes increased significantly in all the three conditions, including the VHL inhibitor VH298 ( Figure 5A-B). Analysis of the RNA-seq data revealed an increase in each of the four genes in VH032 treatment (Dataset 1 (Frost, 2019)); however, this level was insufficient to reach the threshold of log2FC of 0.58 ( Figure 5C). As VH298 is more potent than VH032 (Frost et al., 2016), VH298 is predicted to induce a more pronounced effect on gene expression of target genes. It is likely that the 252 upregulated genes were found to be enriched solely in hypoxia and IOX2, are also induced by the more potent VHL inhibitor VH298, indicative of a common regulator between these treatments. Furthermore, these 252 genes showed significant enrichment of genes involved in pathways similar to commonly upregulated genes ( Figure 5D, Figure 2B), as well as enriched for HIF binding sites ( Figure 5E).

Discussion
Here, we used high-throughput RNA-sequencing to investigate the similarity and differences in the transcriptional response towards hypoxia, the PHD inhibitor IOX2 and the VHL inhibitor VH032. Although genome-wide expression profiling comparing hypoxia and IOX2 has previously been reported (Chan et al., 2016), to our knowledge this is the first report of gene expression profiling comparing side-by-side responses of hypoxia and PHD inhibitors to VHL inhibitors. These three treatments activate the HIF transcription factors, but via limiting or inhibiting different components of the hypoxia signalling pathway.
Our results provide insights into the effects of inhibiting PHD or VHL on HIF target genes, and unique responses in each condition. While hypoxia induced the broadest transcriptional changes, IOX2 and VH032 possessed similar transcriptional responses. The three conditions upregulated a common group of 306 genes (Dataset 1 (Frost, 2019)), the majority of which are regulated by HIF transcription factors ( Figure 2B). From this list, we were able to validate a number of known HIF targets in HeLa and HFF cells (Figure 3, Figure 4). Furthermore, we also found that 132 of these 306 genes were either validated HIF targets or possess HIF-1α/2α binding sites (Dataset 1 (Frost,  2019)). This suggest that while the 132 genes are likely HIF targets, the remaining 174 genes (Dataset 1 (Frost, 2019)) could also be potential novel HIF targets.
As hypoxia, VH032 and IOX2 activate HIF, our datasets of genes induced in these conditions are predominantly enriched for HIF transcription factors ( Figure 6A-C). Beside gene activation, hypoxia also promotes gene repression. Our results show that hypoxia downregulated a significantly larger number of genes  compared to IOX2 and VH032 ( Figure 1G). Various mechanisms of transcriptional repression in hypoxia are known (Batie et al., 2018) and one mechanism includes the activity of SIN3A. A recent bioinformatics study showed that SIN3A and a number of its co-repressors including HDAC1 were overrepresented in the proximity of genes transcriptionally repressed by hypoxia (Tiana et al., 2018). Consistent to the reported roles of SIN3A and HDAC1 in hypoxia signalling (Batie et al., 2018), we found that our datasets of downregulated genes in response to hypoxia were enriched for SIN3A and HDAC1 ( Figure 6D). The transcription factor REST was also enriched in genes repressed in hypoxia ( Figure 6D) and this is consistent to a recent finding that REST transcriptionally repressed genes in hypoxia (Cavadas et al., 2016).
Geneset enrichment analysis suggests that hypoxia, VH032 and IOX2 commonly upregulated genes enriched for NF-κB signalling ( Figure 2B). NF-κB is a transcription factor that has been shown to respond to cellular stress, including hypoxia and PHD inhibition (Cummins et al., 2006). Under hypoxia, NF-κB is activated and induces increased angiogenesis and decreased apoptosis (D'Ignazio & Rocha, 2016).
Overall, we reveal for the first time a comparison of genomewide gene expression analysis of HIF activators, including the physiological inducer hypoxia, and small molecule inhibitors of PHD enzymes and VHL. Understanding the differential regulation of genes in response to the three conditions will help to determine the functions of PHD and VHL in hypoxia signalling, as well as revealing novel HIF-dependent and -independent genes. Furthermore, considering the potential use of PHD inhibitors that are currently in clinical trials and the potential of VHL inhibitors for therapeutic benefits, our report contributes to the further understanding of the pharmacological effects of these inhibitors. In this manuscript, Frost have used high-throughput RNA sequencing to compare the cellular et al transcriptional response to hypoxia with that of pharmacologic HIF hydroxylase or pVHL inhibition. Of interest, hypoxia was a more potent stimulus of gene expression and repression than either hydroxylase or pVHL inhibition. Given the recent introduction of pharmacologic targeting of the HIF hydroxylases into clinical medicine, understanding the cellular consequences of exposure to these agents in terms of global gene expression (while already studied to some degree) is of continuing interest. The manuscript is well written, addresses a clear and interesting hypothesis and delivers clear and important data. The experiments are controlled and appear to be well carried out. The introduction includes a nice up to date review of the pharmacology of the HIF pathway.

Specific Comments:
The authors should discuss the choice of IOX2 as a pharmacologic hydroxylase inhibitor for the transcriptomic aspect of the study. What was the rationale of choosing this inhibitor over some of the more commonly used and well characterized agents such as DMOG. Furthermore, the authors should explain why the use of IOX2/VH032 was changed for FG-4592/VH298 for the RT-PCR validation studies (would it not have made more sense to stick with the same compound).
It should be made clear why Hela cells and not a more physiologically relevant cell line were used for the transcriptomic experiments. As these studies involve pharmacologic and hypoxic treatments, could primary cells not have been used? Hela cells carry many chromosomal abnormalities and have limited relevance to normal cell function. This limitation should be discussed.
In Figure 1, it would be good to demonstrate (e.g. by western blot or HRE-liuciferase reporter assay or both) the distinct temporal profile of HIF-1 and HIF-2 expression/activity elicited by the three treatments in Hela cells over the 16 hour time course of exposure. It is possible that these treatments (due to their differential mechanisms) elicit quite different temporal changes in HIF expression which may be of relevance to the genes expressed.

Is the work clearly and accurately presented and does it cite the current literature? Yes
Is the study design appropriate and is the work technically sound?