Transcriptome sequencing revealed differences in the response of renal cancer cells to hypoxia and CoCl 2 treatment

Human cancer cells are subjected to hypoxic conditions in many tumours. Hypoxia causes alterations in the glycolytic pathway activation through stabilization of hypoxia-inducible factor 1. Currently, two approaches are commonly used to model hypoxia: an alternative to generating low-oxygen conditions in an incubator, cells can be treated with CoCl 2. We performed RNA-seq experiments to study transcriptomes of human Caki-1 cells under real hypoxia and after CoCl 2 treatment. Despite causing transcriptional changes of a much higher order of magnitude for the genes in the hypoxia regulation pathway, CoCl 2 treatment fails to induce alterations in the glycolysis / gluconeogenesis pathway. Moreover, CoCl 2 caused aberrant activation of other oxidoreductases in glycine, serine and threonine metabolism pathways.


Introduction
Hypoxia is characterized by reduced oxygen supply and appears in multiple pathological conditions including tumours. However, hypoxia can also have a functional role during normal mammalian development and embryogenesis 1 . Cells respond to hypoxic conditions both on biochemical and gene expression levels by switching from aerobic metabolism to anaerobic glycolysis and by expression of stress-related genes involved in regulation of cell death, erythropoiesis, angiogenesis and survival 2-4 . The activation of many O 2 -regulated genes is mediated by hypoxia-inducible factor (Hif1a). Under normoxia, Hif1a is hydroxylated by specific prolyl hydroxylases (PHD1, PHD2 and PHD3). This reaction requires oxygen, 2-oxoglutarate and ascorbate 5,6 . When Hif1a is hydroxylated, it interacts with the von Hippel-Lindau tumor suppressor protein (pVHL). pVHL forms the substrate-recognition module of an E3 ubiquitin ligase complex, which directs Hif1a poly-ubiquitylation and proteasomal degradation 7,8 . Under hypoxia (less than 5% O 2 ), PHD activity is inhibited by cytoplasmic reactive-oxygen species (ROS) which alter the oxidation state of Fe 2+ (a cofactor for PHD activity) to Fe 3+ . This alteration inhibits PHD activity and Hif1a hydroxylation, thus Hif1a cannot interact with pVHL and promotes HIf1a stabilization 9,10 . This anaerobic condition and stabilization of Hif1a are characteristic of many tumors. The most common molecular abnormality in renal cell carcinoma is the loss of VHL, which is found in about 50-70% of sporadic cases. Consequently, renal carcinomas with mutations in VHL have high steady-state levels of Hif1a expression and are hypoxic 11 . Some divalent cations such as cobalt (Co 2+ ), nickel (Ni 2+ ), and the iron-chelator deferoxamine (DFX), have been applied to mimic hypoxic conditions in cultured cells as they activate hypoxic signals by stabilizing HIF1a 12 . Transition metal Co 2+ could induce hypoxic response by inhibiting PHD activity via iron replacement. Therefore, treatment of a cell culture with cobalt chloride (CoCl 2 ) is a common model of hypoxia 13 . The second classical setup to study hypoxia is hypoxia induction in a CO 2 incubator with a regulated level of oxygen (less than 1% O 2 ). In this work, we performed RNA sequencing of Caki-1 clear cell renal cancer cell lines treated with hypoxia and with CoCl 2 to understand how adequate CoCl 2 treatment was as a hypoxia model. We propose that CoCl 2 is not a completely correct model for hypoxia, as it aberrantly induces various hydroxylases not involved in hypoxia pathways and fails to induce downstream biochemical pathways normally induced by hypoxia.

RNA preparation and RNA sequencing
Total RNA was extracted from Caki-1 cells with Trisol reagent according to the manufacturer's instructions (Invitrogen). Quality was checked with BioAnalyser and RNA 6000 Nano Kit (Agilent). PolyA RNA was purified with Dynabeads ® mRNA Purification Kit (Ambion). An Illumina library was made from polyA RNA with NEBNext ® mRNA Library Prep Reagent Set (NEB) according to the manual. Sequencing was performed on HiSeq1500 with 50 bp read length. 10 million reads were generated for each sample.
We considered a gene to be differentially expressed if the adjusted p-value in DESeq test was lower than 0.05 and fold-change values were higher than 2 (or lower than 1 2 ). These sets of differentially expressed genes were further used for gene category enrichment analysis. We took the subset of genes which were found differential in both hypoxia against normoxia controls and after CoCl 2 treatment against normal control and only in one of each experiments. These 3 sets of genes were analyzed with DAVID web service (version 6.7) 16 to find KEGG 17 pathways enriched with the genes.
A simple transcriptome-based hypoxia signature was constructed as follows: for every sample being evaluated (e.g., TCGA cancer sample), we considered only the genes which were differentially expressed between hypoxia and untreated Caki-1 cell line (DESeq test adjusted p-value< 0.05). For these genes, we multiplied their logarithmic fold-change (hypoxia vs untreated) to their expression in the evaluated sample. The resulting values were then summed up over the genes under consideration. This yielded a per-sample hypoxia score which would be higher in samples with increased expression of hypoxia-induced genes and decreased expression of hypoxia-suppressed genes. The samples are labelled as untreated_1, untreated_2, untreated_3, CoCl 2 _1, CoCl 2 _2, CoCl 2 _3, hypoxia_1, hypoxia_2, hypoxia_3 which correspond to 3 untreated samples, 3 samples treated with CoCl 2 and 3 samples in hypoxic conditions.

Results
To compare transcriptional effects caused by hypoxia and by CoCl 2 exposure, we performed transcriptome sequencing (RNA-seq) of Caki-1 clear cell renal carcinoma cell line in three conditions: untreated, treated with CoCl 2 , and exposed to hypoxic conditions (1% O 2 ). We searched for differentially expressed genes in two comparisons: CoCl 2 -treated against untreated cells and cells under hypoxia against untreated cells.
We first asked how transcriptomic changes after both treatments fit into the established hypoxia gene signature 1,18,19 . Figure 1A shows relative expression of the genes from hypoxia signature in all sequenced samples. We observed that genes in hypoxia signature were upregulated both in hypoxia and after CoCl 2 treatment, though the effect was much stronger in CoCl 2 -treated samples. To check if a higher magnitude of transcriptomic changes in CoCl 2 compared to hypoxic conditions was characteristic to other differentially expressed genes, we compared the distributions of logarithmic fold-changes (absolute value) for gene expression in these two experiments ( Figure 1B). Surprisingly, for the majority of genes, hypoxia-induced changes were higher than the ones induced by CoCl 2 . In other words, hypoxia resulted in broader transcriptome response than CoCl 2 treatment even though specific hypoxia-related genes were more affected after CoCl 2 treatment.
These results suggested that CoCl 2 treatment could be an incomplete model of hypoxia capturing only upstream signalling events in hypoxia pathways and not reflecting broader downstream effects. To further investigate the differences between CoCl 2 model and real hypoxia, we compared the sets of differentially up-and downregulated genes in CoCl 2 -treatment (against untreated) and in hypoxia-treated (against untreated) cells. Figure 2A summarizes the overlap between those gene sets. To understand what regulatory and biochemical pathways were affected in each treatment, we performed gene category enrichment analysis over KEGG 17 pathways with DAVID web service (version 6.7) 16 for genes affected in both treatments or exclusively in one treatment.
The genes which were significantly upregulated in hypoxic conditions (1% O 2 ) but not after treatment with CoCl 2 were significantly enriched in the glycolysis/gluconeogenesis pathway which was known to be related to hypoxia 1,15,18 . Unexpectedly, we detected no enrichment in glycolysis/gluconeogenesis pathway with the genes

A B
differentially expressed after CoCl 2 treatment which was confirmed by a heatmap for glycolysis/gluconeogenesis-related genes ( Figure 2B).
The genes upregulated in both hypoxia and under CoCl 2 treatment were enriched in the MAPK pathway. As the MAPK pathway was known to activate hypoxic response, MAPK activation was expected in hypoxia. Even though CoCl 2 affected directly HIF1a pathway, MAPK was activated in CoCl 2 -treated samples as well as in real hypoxia. This result supported a previous observation of MAPK-dependent activation of hypoxia response under CoCl 2 treatment 12 . Surprisingly, we observed systemic lupus erythematosus-related pathway activation in both treatments. The set of genes upregulated in this pathway (histone proteins H2A, H2B, H3, H4 and MHCII antigen-presenting genes) could be unrelated to lupus erythematosus, but rather could indicate increased proliferation and inflammation.
We also explored the genes specific to CoCl 2 treatment but not affected by hypoxia. Surprisingly, the genes upregulated after CoCl 2 treatment but not changed in hypoxia were enriched in the glycine/serine/threonine biosynthesis pathway (Figure 2A and 2B).
We hypothesized that Co 2+ ion could substitute metal cofactors of several oxidoreductases in the pathway and subsequently impair their activity. This, in turn, could require greater amounts of enzyme to be synthesized.
Hypoxic conditions within a tumour have been shown to predict worse clinical outcome 20 . We used our data on wholetranscriptome profiling of kidney cancer cell lines in normal and hypoxic conditions to extract a wide hypoxia signature and validate it with TCGA data on transcriptomes of kidney tumours 21 . Major variation in TCGA transcriptomes ( Figure 3A) is generated by the difference between tumours and adjacent normal samples. We projected our sequenced samples to the principal components derived from TCGA samples and observed that the direction of transcriptome changes between untreated and hypoxic cell lines is slightly similar to the difference between normal and tumour samples, though hypoxia-related changes couldn't explain normal-tumour difference. We constructed a transcriptome-based hypoxia signature as described in the Methods section. To test if this hypoxia signature predicted clinical outcome for kidney cancer patients, we explored the distribution of our hypoxia scores between disease free patients and patients in which the disease had recurred or progressed. Hypoxia scores were significantly higher for recurred/ progressed patients (Wilcoxon test p-value=0.0009, Figure 3B).

Discussion
The current study for the first time provides RNA-seq data revealing hypoxia-induced transcriptomic changes, which allows broader understanding of the processes related to hypoxia. In our analysis, we explored the limits of applicability of CoCl 2 treatment as the model of hypoxia. Briefly, we observed that CoCl 2 strongly alters expression of few genes important for hypoxia signalling, but fails to influence the essential downstream consequences of hypoxia, particularly the glycolysis/gluconeogenesis pathway. This might suggest the existence of alternative regulation mechanisms which trigger the downstream events in hypoxia together with main VHL/HIF1a pathway. CoCl 2 treatment also abberantly induced pathways which did not respond to hypoxia. These included the glycine, serine and threonine metabolism pathways. We hypothesized that its aberrant activation might be caused by Co 2+ ion binding to the enzymes (other than PHD proteins) involved in Glycine, Serine and Threonine biosynthesis. RNA-seq data was deposited to NCBI SRA under SRP066934 study accession code. The study contained experiments under the following accession codes: untreated (SRX1459966, SRX1459967, SRX1459969), treated with CoCl 2 (SRX1459974, SRX1459977, SRX1459978), exposed to hypoxia (SRX1459979, SRX1459981, SRX1459984).

F1000Research
Open Peer Review This project uses mRNA sequencing to monitor changes in gene expression in a human kidney cell line in response to direct hypoxia (5% O ) as compared to treatment with 300mM CoCl , frequently used as a surrogate for oxygen depletion. The results are very clear-cut: monitoring several genes known to be key hypoxia regulators shows CoCl to be more effective than oxygen depletion but in contrast, monitoring a much broader gene response showed O depletion to have greater effects. This contrast led the authors to concentrate attention on several specific biochemical pathways: a) glycolysis/gluconeogenesis is strongly activated in hypoxia (as expected) but not following treatment with CoCl ; b) the reverse is the case for glycine, serine, threonine metabolism which is strongly induced by CoCl but not by hypoxia.
These last examples make clear the non-equivalence of CoCl treatment and real hypoxia and this article therefore issues a strong caveat regarding the use of the inorganic surrogate. It is most certainly worth publishing for this reason alone.
Since hypoxic conditions in tumours are taken to be predictive of a poor clinical outcome, the authors tested whether the direction of changes in the transcriptome of their untreated and hypoxic cell lines mirrors the difference between tumours and their normal cell counterparts. Hypoxia scores from mRNA expression data were somewhat (but significantly) elevated in recurrent/progressive patients.However, this difference does not (unfortunately) appear to be sufficient to use clinically (if I am understanding their intentions correctly). The upshot is that the authors do not mention this comparison in their Discussion.
On practical matters: they might expand the caption to Fig 1. In A) I assume the three tracks are biological replicates: if so, write it. In B) a fuller explanation of the X-axis 'Index' would be helpful. Check on the printing of '300mM CoCl ' (in Methods, bottom of p2).

I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.
No competing interests were disclosed.