Non-coding RNAs predict recurrence-free survival of patients with hypoxic tumours

Hypoxia promotes tumour aggressiveness and reduces patient survival. A spectrum of poor outcome among patients with hypoxic tumours suggests that additional factors modulate how tumours respond to hypoxia. PIWI-interacting RNAs (piRNAs) are small non-coding RNAs with a pivotal role in genomic stability and epigenetic regulation of gene expression. We reported that cancer type-specific piRNA signatures vary among patients. However, remarkably homogenous piRNA profiles are detected across patients with renal cell carcinoma, a cancer characterized by constitutive upregulation of hypoxia-related signaling induced by common mutation or loss of von Hippel-Lindau factor (VHL). By investigating >3000 piRNA transcriptomes in hypoxic and non-hypoxic tumors from seven organs, we discovered 40 hypoxia-regulated piRNAs and validated this in cells cultured under hypoxia. Moreover, a subset of these hypoxia-regulated piRNAs are regulated by VHL/HIF signaling in vitro. A hypoxia-regulated piRNA-based score (PiSco) was associated with poor RFS for hypoxic tumours, particularly Stage I lung adenocarcinomas, suggesting that hypoxia-regulated piRNA expression can predict tumour recurrence even in early-stage tumours and thus may be of clinical utility.

Using our custom small RNA sequence analysis pipeline 18,39,40 we assessed 3,020 piRNA transcriptomes in hypoxic and non-hypoxic patient tumours derived from seven different organs (Fig. 1). We discovered that hypoxia is associated with increased expression of a subgroup of piRNAs in patient tumours, and that most of these piRNAs are also overexpressed in tumour cell lines exposed to hypoxic conditions or VHL disruption in vitro. Remarkably, we have also found that that the expression of these hypoxia-regulated piRNAs (hypoxia-regulated piRNAs) can be used to predict recurrence-free survival (RFS) of patients with hypoxic tumours.

Experimental design and identification of hypoxic status in tumours.
To determine if hypoxia can induce changes in piRNA expression, a total of 3,020 tumour samples processed by The Cancer Genome Atlas (TCGA) Research Network were analyzed for hypoxic status based on expression of 143 genes included in the Winter hypoxia metagene signature 41 . Consensus clustering analyses were performed for each tumour type 42 and the resulting tumour clusters were then compared with the Winter hypoxia metagene signature. The cluster with the most concordant gene expression signature was classified as "hypoxic", while the most discordant was defined as "non-hypoxic". Samples that did not match the Winter signature were defined as "unclassified". Cell lines exposed to hypoxia in vitro. We used 9 tumour cell lines that resembled the TCGA tumour types analyzed above. Selected cell lines were: A549 (lung), FaDu (pharynx), HCC15 (lung), HCC4006 (lung), Mia-Paca-2 (pancreas), NCI-H841 (lung), OE33 (esophagus), PANC1 (pancreas), and SiHa (cervix). Cell lines were obtained from ATCC and passaged using 0.1% trypsin citrate buffer + EDTA. For hypoxia exposure, 70-80% confluent plates were transferred to a controlled-atmosphere hypoxia chamber and incubated for 16 hours at 1% O 2 . Cells were then rinsed with cold PBS, scraped off the plate in 1 mL of PBS into a sealed eppendorf tube, and taken out of the hypoxia chamber on ice. After spinning the cells at 1700 rpm for 5 min and aspirating PBS, cells were lysed in 1 mL of Trizol.
RNA extraction for Sequencing. Total RNA was extracted using TRIzol reagent according to the manufacturer's instructions, and eluted in RNase-free water. RNA concentration and quality were determined using a Experimental design for discovery and validation of hypoxia-regulated piRNAs. Hypoxic status from 3,020 human tumors was derived using an established metagene signature 41 . Expression of 23,440 human piRNAs was assessed using our previously published small RNA sequence data analysis pipeline. Hypoxiaregulated piRNAs were determined by comparing hypoxic and non-hypoxic tumors in a per tissue basis. Candidate hypoxia-regulated piRNAS derived from these analyses were further validated by in-vitro studies, using cell lines matching the tissues using during the discovery stage. These validated hypoxia-regulated piRNAs were considered for downstream analyses.
Scientific REPORtS | (2018) 8:152 | DOI:10.1038/s41598-017-18462-z NanoDrop ™ 2100 spectrophotometer, and samples were stored at −70 °C. RNA was cleaned up as needed using a Zymogen kit. The same sequencing and analysis protocols were applied to TCGA tumours and tumour cell lines 43 . Small RNA sequencing. Small RNA sequencing libraries were generated at Canada's Michael Smith Genome Sciences Centre and sequenced using the Illumina Genome Analyzer and HiSeq. 2000 platform, for both TCGA tumours and cell lines. We extracted unaligned FASTQ reads from BAM files of each tumour from the TCGA cohort, which were retrieved from the cgHub data repository (dbgap Project ID: 6208) using the PartekFlow ™ platform. To make data analysis protocols comparable, we applied the same procedure for the hypoxic/normoxic cell lines sequenced. Reads were subjected to quality control in order to exclude non-biological reads, such as 3′ adapter sequences. To avoid overlapping with miRNA-derived reads, a secondary filtering step was applied based on read trimming, retaining reads that were ≥23 bp with phred scores ≥20.
Generation of piRNA transcriptomes. For each tumour cell line analyzed, filtered reads were aligned to the human genome (GRCh37/hg19) using the STAR aligner 44 , with specific parameters for short reads: i) only reads 23 bp or longer matched to the genome, ii) number of mismatches < = 5% of mapped length, and iii) splicing switched off. Mapped reads were subsequently quantified within the PartekFlow ™ platform using a piRNA-specific annotation file generated by the piRNABank (http://pirnabank.ibab.ac.in/) 45 . This reference transcriptome considers widely accepted piRNA sequence features, such as sequence bias for uracil in the 1st and adenine in the 10th position, although known biases for 2-o-methylation in the 3′ end were not considered.
VHL knockdown by siRNA. A549 cells were treated with siRNA to VHL at 5 nM using Dharmafect for 48 hrs at 21% O 2 (Ambion, Inc. Cat#(Lot): non-targeting control #4390843, siRNA 1 VHL #4390824(s14789), siRNA 2 VHL #4392420(s14790), siRNA 3 VHL #4392420(s14791)). RNA was extracted using Quick-RNA MiniPrep Kit column purification (Zymo Research, #R1055), and resuspended in nuclease-free water. cDNA conversion was performed on 20 ng of isolated RNA using the TaqMan (R) MicroRNA Reverse Transcription Kit (Thermo Fisher, #4366596) and custom reverse transcription primers for DQ580854 (context sequence: TGAGGAG CCAATGGGGCGAAGCTACCATC, target sequence: UGAGGAGCCAAUGGGGCGAAGCUACCAUC), D Q 5 9 0 4 0 4 ( c o n t e x t s e q u e n c e : T G G T G TAT G T G C T T G G C T G A G G A G C C A AT G G , target sequence: UGGUGUAUGUGCUUGGCUGAGGAGCCAAUGG), and DQ596992 (context sequence: GCAATAACAGGTCTGTGATGCCCT TAGA, target sequence: GCAAUAACAGGUCUGUG AUGCCCUUAGA). All primers used FAM as a reporter dye and NFQ as a reporter quencher. RT-qPCR using paired primers was performed for 40 cycles, and RQ values were calculated with the ddCt method, with each treatment normalized to its respective endogenous U6 control. Experiments were performed in biological duplicates.
Survival Analysis. TCGA samples with complete recurrence-free survival and piRNA expression data (RPKM) from different tumour types (HNSC, LUAD, and LUSC) were included in the Cox proportional hazard (COXPH) model. Only hypoxia-regulated piRNAs that were validated in hypoxic vs normoxic cell lines of the same cell type were considered for the generation of predictive signatures. Age at diagnosis and clinical stage were included as know stratification factors in the COXPH model for each tumour type, in addition to smoking status and history for LUAD and LUSC patients, but removed in the score evaluation. PiRNAs were included in the piRNA score (PiSco) if they were found to contribute significantly to the classification of patient survival by the COXPH model (p-value ≤0.1), and the optimized model with the lowest p-value was used for further analysis. PiSco was generated by multiplying the expression value of a given piRNA by its hazard coefficient, then summing the transformed gene expression values per sample 39,46 . Risk scores were then ranked and divided into tertiles, selecting upper and lower tertiles to compare using the log-ranked method with a significance threshold of p ≤ 0.05. The recurrence-free survival of the high-score cohort was then compared to that of the low-score tertile.  Table S1) for a consensus clustering analysis to stratify TCGA tumour samples into "hypoxic", "non-hypoxic", or "unclassified" tumour groups ( Figure S1A-H) 41,42 . In total, 3,020 human tumours derived from cervix (CESC, n = 306), head and neck (HNSC, n = 522), esophagus (ESCA, n = 196), lung adenocarcinoma (LUAD, n = 517), pancreas (PAAD, n = 179), lung squamous (LUSC, n = 501), kidney (KIRC, n = 534), and ovary (OVCA, n = 265) were used for our analyses (Table 1). Each tissue type was analyzed individually. Most tumours were clearly classified as hypoxic or non-hypoxic in concordance with the metagene signature, and we therefore focused on these tumours for subsequent discovery of hypoxia-associated piRNA expression. Kidney and ovarian tumours could not be clearly classified as hypoxic or non-hypoxic based on the Winter signature ( Figure S1G-H), and were therefore not included for further analysis.

Identification of hypoxic status in human tumours.
piRNA expression is selectively deregulated by hypoxia in human tumours. We then sought to identify differences in piRNA transcriptomes between hypoxic and non-hypoxic tumours. Expression levels of 23,440 human piRNAs were determined on a per tumour basis ( Figure S2) using our previously published custom small RNA sequencing analysis pipeline 18 . To identify the most robust changes in piRNA expression between hypoxic and non-hypoxic groups, piRNAs were only included in analyses if they had a median expression ≥10 RPKM in at least one of the groups (hypoxic and/or normoxic) and a minimum of 2-fold change in median RPKM expression values. Comparisons were performed using the nonparametric Mann-Whitney U test, with a false discovery rate correction for multiple testing using the Benjamini-Hochberg method 47 .
In vitro tumour models recapitulate hypoxia-regulated piRNA expression patterns. To validate whether changes in hypoxia-regulated piRNA expression observed in clinical samples are driven by hypoxia, piRNA expression levels were examined in a panel of nine cell lines from the corresponding tumour types used in the previous analysis, following exposure to hypoxic (1% O 2 for 16 hours) or normoxic (21% O 2 ) conditions (Fig. 1). Forty out of the 71 hypoxia-regulated piRNAs identified in hypoxic clinical tumours showed hypoxic regulation in the associated tumour cell lines. Of these validated hypoxia-regulated piRNAs, 36 were upregulated and 4 downregulated in hypoxic cells (Fig. 3).
von Hippel-Lindau factor (VHL)/HIF signaling axis regulates the expression of hypoxia inducible piRNAs. We previously observed remarkable homogeneneity in piRNA expression in RCC tumors 18 , which are known to have a high frequency of loss-of-function mutations in the VHL gene that causes constitutive stabilization of HIF-1α. We hypothesized that the induction of those piRNAs upregulated by exposure to 1% O 2 would be dependent on HIF-1α. To address this hypothesis we selected A549 lung adenocarcinoma cells (with the highest number of hypoxia-regulated piRNAs detectable at the sequencing level) to test the effect of siRNA-mediated knockdown of VHL on expression of DQ590404 and DQ596992, which were two of the most highly overexpressed hypoxia-regulated piRNAs in A549 cells. We found that siRNA-mediated knockdown of VHL increased expression of both piRNAs (Fig. 4A,B), and that this induction was inhibited when VHL knockdown was combined with siRNA-mediated knockdown of HIF-1α, indicating that HIF-1α stabilization induced by VHL knockdown was responsible for the increase in piRNA expression.
A hypoxia-regulated piRNA signature delineates risk of recurrence in patients with hypoxic tumours. Considering the adverse effects of tumour hypoxia in tumour recurrence, we next investigated if a hypoxia-regulated piRNA expression signature was associated with recurrence-free survival (RFS). We first restricted our analyses to patients with hypoxic tumours as defined by the Winter hypoxia metagene signature. Using a Cox proportional hazard (COXPH) model, we developed a piRNA-based score (PiSco) for each tumour type by weighting the expression levels of hypoxia-regulated piRNAs in the tumour type using the obtained Cox coefficients and summing the resulting values, as previously described 48 . Weighted contribution of piRNAs to tissue-specific PiSco signature is shown in Table 3. After correction by known cancer-specific risk factors (e.g., age, stage), comparison of the upper and lower PiSco tertiles revealed that the combined expression of these hypoxia-regulated piRNAs could classify patients with hypoxic HNSC, LUSC, and LUAD tumours as low or high risk of recurrence-free survival (Fig. 5A-C).
To verify that the association between PiSco and RFS was not driven exclusively by late-stage tumours in the cohorts, we applied the same analysis strategy now only considering hypoxic Stage I tumours. Since the lung cancer cohort exhibited the largest number of Stage I cases (54 and 61 for LUSC and LUAD, respectively), we restricted these analyses to this tumour type. While results from LUSC showed a tendency to significance (p value = 0.0985), PiSco was able to appropriately classify early-stage LUAD patients into low or high risk groups for recurrence-free survival (p value = 0.0051, Fig. 5D).

Discussion
Hypoxic tumours are known to be aggressive and resistant to radiotherapy and chemotherapy 6 . However, there is a spectrum of poor responses to therapy in patients with hypoxic tumours, and there is no effective way to identify those who will respond the worst. Thus, the development of classification tools intended to identify patients harbouring tumours at high risk of recurrence would be extremely beneficial for monitoring patient outcome after receiving treatment. Our results indicate that hypoxia can regulate the expression of piRNAs which, in turn, may aid in the prediction of prognosis for patients with hypoxic tumours. Specifically, we found that: 1) piRNA expression patterns are different between hypoxic and non-hypoxic human tumours derived from a variety of anatomical sites; 2) hypoxia-mediated regulation of piRNAs occurs in vitro, resembling patterns observed in human tumours; and 3) a hypoxia-regulated piRNA expression signature can identify hypoxic tumours at higher risk of recurrence in HNSC, LUSC, and LUAD patients.
We deduced the oxygenation status of tumour samples processed by TCGA using a previously characterized and validated hypoxia gene expression signature from Winter et al. 41,42 . Using this method, tumours from all tissue types were classified into hypoxic or non-hypoxic groups, except for ovarian and kidney tissues. While this lack of clear classification between hypoxic and non-hypoxic ovarian cancer samples was not anticipated, it was expected in renal cell carcinomas because >75% of sporadic RCCs contain mutated or lost VHL 49 . Loss of VHL function mimics hypoxia-mediated stabilization of HIF-1α, which then dimerizes with HIF-1β to form HIF-1 and induce transcription. Constitutively activated HIF-1 in renal tumours, in part attributed to mutated VHL, upregulates genes involved in the hypoxic response regardless of local cellular oxygen tensions. Accordingly, we have previously found that piRNA expression is remarkably homogeneous across RCC patient samples 18 , suggesting piRNA expression patterns in RCC may be related to VHL status in these tumours. Thus, we sought to directly assess the ability of VHL to affect piRNA expression in a non-renal setting. To do this, we silenced VHL in a lung cancer cell line that had robust upregulation of piRNAs in hypoxia (Fig. 3). We observed increased expression of two different piRNAs when VHL was knocked down in A549 cells, and piRNA expression was decreased to control levels by simultaneous knockdown of HIF-1α (Fig. 4B,C). These data indicate that the VHL-mediated induction of DQ590404 and DQ596992 piRNAs is HIF-1α dependent, and suggest that the induction of piRNAs in hypoxic cells more generally may also be regulated by HIF-1.
In our analysis, we identified 40 piRNAs that are associated with hypoxia in patient tumours and regulated by hypoxia in cell models of the same tissue type. Interestingly, we found that most of the hypoxia-regulated piRNA species (36 out of 40; 90%) were upregulated in hypoxic compared to normoxic conditions, which is consistent with previous observations in a glioblastoma cell lines 50 . This phenomenon is opposite to what is observed in miRNA, where hypoxia mainly induces a reduction in expression levels by modulating their biogenesis 10 . Since  Table 2. Shared and common expression of hypoxia regulated piRNAs. The number of piRNAs for a single tumour type or combinations indicate those piRNAs that are hypoxia regulated exclusively in that or those tissues (and not the total number of hypoxia-regulated piRNAs in the tissue).
piRNAs are generated through a Dicer-independent pathway, it is expected that the effect of low-oxygen tension will be different on piRNAs compared to miRNAs. Hypoxia-mediated regulation of piRNAs has also been described in breast cancer cell lines in a HIF dependent manner, although the absence of HIF binding sites in the vicinity of encoding loci suggest that are unlikely to be directly regulated by HIF 51 . Overall, our results provide novel mechanistic insights involving a VHL-mediated stabilization of HIF1a regulating the expression of a subset of piRNAs, and also highlight the importance of considering hypoxia-mediated regulation of piRNA expression in a tissue-specific manner. Previous work has shown individual piRNAs that correlate with clinical features in different tumour types 28 . Also, multipiRNA signatures have been used to stratify patients by prognosis 39,46 , although the presence of some hypoxia-regulated piRNAs in these published piRNA signatures suggest potential identification of hypoxic tumours that are expected to have poorer prognosis from the bulk patient populations. Our PiSco signature enables stratification of patients with the most aggressive hypoxic tumours into low or high risk of recurrence, suggesting that analysis of hypoxia-regulated piRNA expression can be used as a predictive tool to prioritize outcome monitoring after treatment, even in patients with early-stage hypoxic lung adenocarcinoma (Fig. 5).   In conclusion, this proof-of-principle study reveals the influence of the tumour microenvironment on DNA-level regulatory mechanisms with important implications for predicting recurrence in patients with hypoxic tumours. Furthermore, we find that hypoxia-regulated piRNA expression can be altered by VHL status in vitro. Our data encourage further exploration of hypoxia-regulated piRNAs as clinical tools for evaluating the likelihood of tumour recurrence, and to identify patients that would most benefit from adjuvant therapies and/or therapies designed to target hypoxic tumour cells.