Transcriptome Sequencing (RNAseq) Enables Utilization of Formalin-Fixed, Paraffin-Embedded Biopsies with Clear Cell Renal Cell Carcinoma for Exploration of Disease Biology and Biomarker Development

Formalin-fixed, paraffin-embedded (FFPE) tissues are an underused resource for molecular analyses. This proof of concept study aimed to compare RNAseq results from FFPE biopsies with the corresponding RNAlater® (Qiagen, Germany) stored samples from clear cell renal cell carcinoma (ccRCC) patients to investigate feasibility of RNAseq in archival tissue. From each of 16 patients undergoing partial or full nephrectomy, four core biopsies, such as two specimens with ccRCC and two specimens of adjacent normal tissue, were obtained with a 16g needle. One normal and one ccRCC tissue specimen per patient was stored either in FFPE or RNAlater®. RNA sequencing libraries were generated applying the new Illumina TruSeq® Access library preparation protocol. Comparative analysis was done using voom/Limma R-package. The analysis of the FFPE and RNAlater® datasets yielded similar numbers of detected genes, differentially expressed transcripts and affected pathways. The FFPE and RNAlater datasets shared 80% (n = 1106) differentially expressed genes. The average expression and the log2 fold changes of these transcripts correlated with R2 = 0.97, and R2 = 0.96, respectively. Among transcripts with the highest fold changes in both datasets were carbonic anhydrase 9 (CA9), neuronal pentraxin-2 (NPTX2) and uromodulin (UMOD) that were confirmed by immunohistochemistry. IPA revealed the presence of gene signatures of cancer and nephrotoxicity, renal damage and immune response. To simulate the feasibility of clinical biomarker studies with FFPE samples, a classifier model was developed for the FFPE dataset: expression data for CA9 alone had an accuracy, specificity and sensitivity of 94%, respectively, and achieved similar performance in the RNAlater dataset. Transforming growth factor-ß1 (TGFB1)-regulated genes, epithelial to mesenchymal transition (EMT) and NOTCH signaling cascade may support novel therapeutic strategies. In conclusion, in this proof of concept study, RNAseq data obtained from FFPE kidney biopsies are comparable to data obtained from fresh stored material, thereby expanding the utility of archival tissue specimens.


Introduction
Clear cell renal cell carcinoma (ccRCC) makes up the majority of primary renal neoplasms with increasing incidence and considerable morbidity and mortality. Metastasis reflects a major cause of patient death [1,2]. Renal cell cancer ranks among the ten most frequent cancers in women and men accounting for up to 2-3% of all adult cancers or malignancies [2][3][4][5][6].
The ccRCC is only curable by early surgical tumor removal. Thus, efforts to unravel molecular mechanisms of this disease for the search of prognostic markers and novel drug targets are important, e.g. by applying gene expression detection technologies to develop molecular signatures of disease progression.
In this study, we applied RNA sequencing (RNAseq), a method for measuring mRNA abundance based on next generation sequencing (NGS) technology. NGS can identify transcripts even at a low expression level and provides an increased dynamic range for gene expression measurements compared to microarrays [7,8].
Current technologies for whole genome gene expression analyses are largely dependent on "high quality" RNA with low level of degradation. We wanted to test whether lower quality, partially degraded RNA obtained from archival formalin-fixed and paraffin-embedded (FFPE) renal tissues could serve as appropriate source of information.
The quality of RNA extracted from FFPE samples can vary widely among different specimens, or within different samples from the same specimen. RNA undergoes substantial chemical modification during formalin fixation, nucleic acids are cross-linked to proteins and RNA transcripts are degraded to smaller fragments [9]. Differences in formalin fixation methods and age of archival tissue samples add further variation to RNA quality. The Illumina TruSeq RNA Access Kit 1 holds promise to overcome these challenges for RNA sequencing applications by isolating mRNA through a sequence-specific capture protocol resulting in reduced ribosomal RNA and enriched exonic RNA sequences. The TruSeq RNA Access library preparation kit was designed to ensure high quality RNA sequencing data from degraded FFPE samples and to allow comparison across samples that vary in quality.
Transcriptome sequencing of RNA from concurrently harvested FFPE and fresh stored kidney biopsies with subsequent analysis of transcripts and pathways underlying ccRCC in our patient group served as indication of the comparability of the two sources of RNA. The comparison to published data helped to estimate the biological and clinical plausibility of our results.

Study design
This study includes 16 adult patients from Haukeland University Hospital with ccRCC undergoing partial (n = 10) or full (n = 6) nephrectomy between November 2013 and August 2014 (Table 1). Each patient donated four core biopsies, including two with ccRCC and two from adjacent non-affected tissue ("normal"). One pair of ccRCC and normal tissue per patient was then stored in FFPE, the other pair in RNAlater 1 . This paired design allows comparison of mRNA abundance level differences between ccRCC and normal in FFPE and in RNAlater 1 , and to evaluate the impact of storage condition on expression profiles using RNAseq. decision to publish, or preparation of the manuscript. The library preparation, sequencing and partly the bioinformatics analysis were provided by the Genomics Core Facility (GCF), which is funded by the Faculty of Medicine at NTNU and Central Norway Regional Health Authority. The specific roles of this author is articulated in the 'author contributions' section." Note, O.E. was an essential part in performing experiments and manuscript writing. Andreas Scherer (A. S.) is the sole owner/employee of Spheromics (http://spheromics.com/) and provided help in RNA sequencing data analysis and preparation of respective manuscript parts. Illumina Inc. was not involved in the present study at all; accordingly, no funding was obtained. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing Interests: The affiliation with Spheromics does not alter the authors' adherence to PLOS ONE policies on sharing data and materials. Illumina Inc. was not involved in the study.

Quality of Extracted RNA
To assess the quality of the 64 samples of extracted RNA we determined the Agilent RNA integrity number (RIN). Currently, the RIN is the most commonly used measure to determine RNA quality for gene expression analysis [10]. However, RIN values from FFPE samples are not a sensitive measure of RNA quality nor are they a reliable predictor of successful library preparation. Accordingly, previous investigators have used mean RNA fragment size as a determinant of RNA quality for the RNA sequencing library preparation (Illumina TruSeq RNA Access Kit 1 ) when working with RNA obtained from FFPE tissues [11][12][13].
We have therefore also used the DV 200 metric, the percentage of RNA fragments >200 nucleotides to evaluate the RNA quality according to the recommendation of the manufacturer and as described [11][12][13]. Using DV 200 to accurately assess FFPE RNA quality, and by adjusting RNA input amounts, high-quality libraries can be prepared from poor-quality FFPE samples. In this respect, a sufficient DV 200 value of as low as 30% was reported [13].

Gene Expression (mRNA Abundance)
The number of detected genes, which passed an expression filter of more than 15 cpm in at least 8 samples per dataset, for FFPE was n = 9164 and for RNAlater 1 n = 9205. Notably, about 94% of the genes in each dataset (n = 8893) were common to both FFPE and RNAlater 1 datasets; correlation of the logarithmic fold change was R 2 = 0.93, and correlation of the average expression R 2 = 0.97, as shown in S1 Fig. To find sources of similarity in the dataset consisting of all 64 samples and the expression values of expression-filtered 8893 genes, we applied multidimensional scaling (MDS). Samples segregate into two large groups along the leading log-fold change in the dimension 1 of the MDS plot. The leading log-fold change is the average (root-mean-square) of the largest  9  70  Male  24  Partial  >60 pT1AcN0cM0  18  2  I   10  69  Male  34  Partial  >60 pT3AcN0cM0  15  2  III   11  37  Male  27  Partial  >60 pT1AcN0cM0  19  2  I   13  63  Male  24  Full  40 pT3AcN0cM0  69  4  III   15  68  Male  28  Partial  >60 pT1AcN0cM0  21  2  I   16  53  Male  33  Full  56  pT3bN0M1  100  2  IV   18  78  Male  27  Full  47  T3AcN0cM0  60  2  absolute log-fold change between each pair of samples. As deducible from sample annotation in Fig 1A, the major known factor explaining the similarity of biopsy samples was attributed to "Diagnosis" (i.e. tumor and normal). Storage condition (FFPE or RNAlater 1 ) did not appear to cause sample segregation (Fig 1B). In a next step, we identified for each dataset the genes with differential expression changes between ccRCC and normal, and compared the two sets. The FFPE dataset demonstrated 1367 differentially regulated genes and the RNAlater 1 dataset 1418 genes (Benjamini-Hochberg adjusted p value 0.05, and abs FC 2); comparison of the non-tumorous, normal FFPE tissues versus the corresponding normal tissues from the RNAlater 1 group revealed a very high concordance with only 37 differentially expressed genes (data not shown).
In the MDS analysis, plotting values for differentially expressed genes indicates less withingroup variance compared to the analysis of all detected genes, and the shrinkage of log-fold changes indicates that some non-differentially expressed genes can have quite large fold changes (Fig 1C and 1D).
Each of these two datasets shared 1106 (about 80%) of differentially expressed genes with each other. The correlation of the average expression of these 1106 genes was R 2 = 0.97 ( Fig  2A). The log2 fold changes of these differentially expressed genes correlated by R 2 = 0.96 ( Fig  2B). All those genes in both datasets had the same direction of change. Table 2 shows the 20 most significantly affected genes with largest absolute fold changes in the FFPE dataset and the corresponding values of the RNAlater 1 dataset; 17 of these 20 genes were differentially expressed in both datasets, 3 did not pass the expression filter in the RNAlater 1 dataset. Amongst the 17 genes, 14 were among the top 20 ranking differentially expressed genes in the RNAlater 1 dataset. Vice versa, all top 20 differentially expressed genes of the RNAlater 1 dataset were differentially expressed in the FFPE dataset, 14 of which ranking among the top 20 in both datasets (not shown).

Immunohistochemistry
Immunohistochemistry of the three most regulated genes according to Table 2 confirmed strong overrepresentation of neuronal pentraxin-2 (NPTX2) and carbonic anhydrase 9 (CA9) as well as the underrepresentation of uromodulin (UMOD) in ccRCC [14][15][16]. The results are depicted in Fig 3, which also presents respective mRNA abundance plots.

Pathway Analyses
To test whether disease-relevant pathways have been captured in our experiment, we performed Ingenuity Pathway Analyses (IPA) of differentially expressed genes. 91 canonical pathways in the FFPE dataset and 109 pathways in the RNAlater 1 dataset were affected (adjusted p-value0.05) with an overlap of 75%. The most affected pathways to a good extent reflect humoral and adaptive immune responses (Table 3). Sorting the pathways by smallest adjusted  Table 2. Gene expression analyses. The 20 most up-or down-regulated genes in the FFPE data set with corresponding RNAlater 1 values (upper panel), and the 20 most up-or down regulated genes in the RNAlater 1 dataset with corresponding FFPE values (lower panel), filtered by adjusted p-value0.05. Rank indicates the rank of the gene within the list of differentially genes sorted by largest to smallest absolute fold change. 14 genes are shared between the two lists. TU: tumour, NO: normal, FC: fold change, ND: not detected, did not pass the expression filter.

Comparison with Published Data
We compared our ccRCC gene expression changes with findings described in a recently published meta-analysis of ccRCC datasets [17]. All 10 most up-regulated genes and 7 of the 10 most down-regulated genes from Zaravinos et al. [17] were found in the present study and are differentially expressed in FFPE and RNAlater 1 datasets ( Table 4). The remaining genes did not pass our expression filter. The direction of fold changes was identical for all listed genes. We further compared the findings from the FFPE and the RNAlater 1 datasets in relation to the known involvement of vascular endothelial growth factor (VEGF) in ccRCC [18,19]. As demonstrated in Fig 4, many genes of the VEGF and NOTCH signaling cascades were retrieved in the FFPE and the RNAlater 1 datasets with very similar fold changes and agreement in direction of changes. We can also confirm a link to epithelial to mesenchymal transition (EMT) by the overrepresentation of mesenchymal markers, e.g. vimentin (VIM), endothelin 1 (EDN1), fibronectin 1 (FN1), or transforming growth factor-β (TGFB1), and  Table 3. Pathway analysis. The 20 most affected canonical pathways in each NGS dataset with the corresponding values and ranks. Rank indicates the place of the pathway within the list of pathways sorted by largest to smallest -log(adjusted p-value). 12 of 20 pathways are shared between both datasets. TU: tumour, NO: normal, FC: fold change, ND: not detected, did not pass the expression filter. underrepresentation of epithelial markers such as epithelial cell adhesion molecule (EPCAM) or E-cadherin (CDH1). The transcription factor grainyhead-like 2 (GRHL2), which inhibits EMT, is about 10 fold underrepresented [20]. IPA revealed TGFB1 as one the most important regulator of gene expression in our ccRCC datasets, as shown in Fig 5. Of the 1367 differentially expressed genes in the FFPE dataset, the expression levels of 237 genes (17%) are influenced by TGFB1 in the FFPE dataset (Fig 5A), and 253 of the 1418 (18%) differentially affected genes in the RNAlater dataset (Fig 5B). TGFB1 itself was overrepresented 2.3 fold and 2.8 fold in the FFPE and the RNAlater dataset, respectively (Fig 4).

Classifier Analysis
We further wanted to test whether the RNAseq data from the FFPE dataset could be used to develop a molecular classifier for ccRCC. Hence, in a proof of concept approach, we first selected 100 genes with the largest absolute fold change and smallest adjusted p-value among the group of differentially expressed genes in the FFPE dataset. To avoid overfitting, we initially tested the performance of classifier models with 15 or fewer genes, where we preferred those Table 4. Comparison of our gene expression data with data from literature [17]. Twenty genes with smallest p-values and largest absolute fold changes in a meta-analysis of five microarray studies are compared to the corresponding genes and their fold changes and p-values of the NGS datasets. The median fold changes and standard deviations for the meta-analysis are presented. All shown genes were differentially expressed in only 2 or 3 microarray datasets. Large standard deviations indicate a large spread of values in the individual microarray studies. 17 of the 20 genes were found differentially expressed in both NGS datasets, 13 of these with fold changes within the fold change range of the microarray meta-analysis. ND: not detected, did not pass initial expression filter.
Zaravinos et al. [17] Eikrem  Fig 6A-6C. One misclassified sample was a normal sample classified as tumor. However, importantly, this specimen contained some admixture of tumor tissue detected at a second look. The other misclassified sample from a different patient was a tumor sample with some adjacent tissue that was judged to be normal.
In the RNAlater 1 dataset, the single gene classifier model assigned one sample with the histological classification "normal" to the group of tumor samples, yielding an accuracy ACC = 96.8%, AUC = 1.0, and a specificity of 93.8% and a sensitivity of 100%. We then tested the single gene classifier model in an external dataset on a different technology platform. The publically available Gene Expression Omnibus (www.ncbi.nlm.gov/geo/) dataset GSE53757 contains Affymetrix HG-U133 microarray gene expression data from 72 human renal biopsies with four stages of ccRCC, and 72 matched normal samples [14]. The CA9-model correctly classified 139 of 144 samples independent of cancer stage (ACC = 96.5%, ROC AUC = 0.98). Results of this CA9 validation are shown in Fig 6D and 6E.

Serum Analyses of CA9 Levels
Optimally, biomarkers such as the gene panel classifiers are further developed into clinically applicable tests. In our simulation study, we wanted to examine, whether CA9-assisted detection of ccRCC could be translated into a less-invasive clinical application going beyond the information obtainable from tissue samples. To that end, we measured CA9 protein in the serum of our patients with early T1a tumor stage and compared the results of these subjects with patient groups suffering from a more advanced disease, because a strong association between serum levels of CA9 with tumor stage has recently been reported [15].
The between group analyses showed significant differences between patients with T1a tumor stage and either with high tumor load (p = 0.0031) or with metastases (p = 0.0158). The comparison between the latter two groups showed no significant difference.
Additional potential novel classifiers have been found, but await further examination and validation. For example, expression values of the highly up-regulated TNFAIP6 (tumor necrosis factor, alpha-induced protein 6; Fig 4) showed similar performance as CA9 in the FFPE, RNAlater 1 , and the microarray dataset (ACC = 96.9%, 96.7%, 94.4%, respectively). We are presently collecting more material and data to expand and confirm these findings.

Discussion
Our proof of concept study compares transcriptome sequencing of RNA extracted from human renal biopsies of ccRCC and matched adjacent non-tumorous tissue; samples were preserved in two different storage conditions (FFPE and RNAlater 1 ). High similarity of the two datasets indicates that archival FFPE-samples can be utilized in respective studies. We chose RNAlater 1 storage as the comparator. RNAlater 1 is considered to be an excellent RNA stabiliser [21] and many studies show that RNA yields and gene RNA abundance with RNAlater 1 are comparable to those obtained using frozen tissues [22]. Furthermore, the utilization of RNAlater 1 is more practical allowing also decentralized tissue harvesting without special equipment [22,23].
To the best of our knowledge, there has been no in depth report yet comparing matched RNAlater 1 and FFPE storage conditions for parallel RNA sequencing and we are among the first to demonstrate the usability of the new Access kit (Illumina) also allowing low FFPE RNA amounts to generate RNA sequencing libraries. A related study has also demonstrated good concordance of RNA sequencing between the two storage conditions but has used different technology for only two renal cancers [24]. Obviously, the TruSeq Access kit is focused on studying mature mRNA levels in biological samples. A recent study has shown that other approaches, such as DSN (Duplex-specific nuclease)-seq and Ribo-zero-Seq can be used to investigate intergenic and intronic RNA species, reportedly giving information on slightly more mRNA species than polyA-enrichment methods, but at the expense of requiring more sequencing effort [25]. Where it is sufficient to study the human transcriptome coding regions, the TruSeq Access kit provides a cost-effective, highly reliable method, as our study shows.
Recent publications have studied the effect of storage time (up to 10 years) in FFPE on RNA quality and quantity, and the usability in mRNA expression experiments, both microarrays and RNAseq [26][27][28]. In concordance with our own unpublished data where we measured RNA quality and quantity from up to 30 year-old FFPE samples indicating their suitability for RNA sequencing, the publications agree that, RNA is still usable for RNAseq transcriptome studies although the RNA quality suffers with increasing time of FFPE-preservation.
Our approach is further supported by a recent publication showing that a newly developed exon capture RNAseq library preparation protocol for highly degraded RNA provided accurate estimates of RNA abundance, uniform transcript coverage and broad dynamic range investigating FFPE and flash frozen cancer tissues [29].
However, for the genome-wide detection of novel transcripts, whole exome enrichment of RNA might be a necessary additional step [30].
We detected a high degree of similarity between the gene expression results for the two datasets: 94% of the transcripts passing the initial expression filter were shared between the FFPE and RNAlater 1 sample groups, 80% of differentially expressed genes were in common, and 75% of the differentially affected pathways were found in both datasets. The differences in gene expression can probably be mostly explained by the cell-composition variation of the respective biopsies. This well described intra-tumor heterogeneity precluded the detection of an even higher number of common, differentially regulated genes and pathways [31]. Also, the capture process during library preparation could be different depending on the RNA quality. However, the very high concordance between FFPE non-tumor, normal tissue vs. normal tissue stored in RNAlater 1 further emphasizes the high similarity of the two data sets.
Despite some limitations, we have shown a striking similarity between the FFPE and the RNAlater 1 datasets, maintaining biologically relevant information at large. Immunohistochemistry confirmed the three most regulated genes of both data sets. CA9 is essentially not expressed in the normal nephron but specifically in ccRCC [5]. Thus, CA9 is an extensively investigated biomarker of ccRCC and also a predictor of outcome following anti-VEGF therapy [19,32]. In a microarray study with nine patients, UMOD was the gene with the strongest under-representation in RCC [16]. The over-representation of NPTX2 is in accordance with the literature [14].
We also show good concordance with microarray gene expression profiling studies of ccRCC ( Table 4). Directions of gene expression changes between ccRCC and normal samples were identical for a set of differentially expressed genes in the microarray studies (14) and in the NGS studies. 17 of the 20 genes with largest absolute fold changes in the microarray metaanalysis were also differentially expressed in the NGS datasets (Table 4), and most fold changes were within the same range across the studies.
However, limitations and uncertainties in this comparison come from the large discrepancy in the fold changes detected in the microarray studies, and from the fact that all genes in the Table 4 were differentially expressed in only 2 or 3 of five microarray studies used in the metaanalysis. Different amplitudes in fold changes between the microarray dataset and the NGS dataset have been reported before [33]. The authors believe, one reason is that microarray probes might hit some, but not all, isoforms of a gene, and as a result the reported fold change of the probe set does not necessarily represent the expression change of the entire gene [33]. Furthermore, NGS is more sensitive in measurement of abundance differences of lowly or highly expressed genes. Microarrays reach a saturation level in the case of highly expressed genes, but NGS technology with its wider dynamic range of detection is more likely to detect fold changes. This may explain some of the fold change differences observed in the comparison of microarray and NGS data. Nevertheless, our dataset confirmed the trend of expression changes observed in microarray studies.
Our data also support and in part confirm novel therapeutic avenues, such as targeted at activated VEGF /NOTCH /DLL4 signaling cascades [18,[34][35][36][37]. The up-regulated NOTCH ligand Delta 4 (DLL4) is stimulated by VEGF and plays a role in tumor progression also predicting bad outcome [36,38,39]. EMT is augmented in our cancer data and is known to be a relevant feature in ccRCC [40]. Up-regulated TGFB1 was the most significantly affected gene regulator in our study. Accordingly, TGFB1 inhibition was shown to attenuate the invasive capacity of ccRCC cells [34]. However, potential cancer therapy targeted at TGFB1 remains to be developed.
Classifier models consisting of features such as gene expression data in combination with a decision algorithm are powerful tools to support diagnostic and prognostic evaluation of patient data. Gene expression data for CA9-supplemented by CA9 serum protein datashowed an excellent performance both in our datasets and in an independent ccRCC microarray dataset. Thus, our data expand previous reports, which promote CA9 as a diagnostic tool in ccRCC [5,19,41,42].
Taken together, we show that in our hands RNAseq FFPE data are comparable to matched RNAlater 1 data. We used the proof of concept data to explore and to confirm published biological findings, and findings which may be worth following up in larger cohorts, leading to possible novel therapeutic strategies, e.g. based on TGFB1-regulated genes, the NOTCH signaling cascade, and EMT. Also of note, FFPE tissues have the distinctive advantage that material designated for RNA sequencing can be concurrently investigated by light microscopy.
Conclusions: Our study opens the door to transcriptome analyses of the archival, FFPE stored tissues from patients with ccRCC and supports CA9 as a potential marker for ccRCC.

Materials and Methods Patients
Adult patients (n = 16) from Haukeland University Hospital with ccRCC undergoing partial (n = 10) or full (n = 6) nephrectomy and with the possibility to undergo biopsies for this project were included consecutively from November 2013 until August 2014. Patients had a mean age of 58.2±6.8 years (3 females and 13 males). Patients had pT tumor stages T1a (n = 10), T2a or b (n = 2) and T3a or b (n = 4) [43]; additional patient characteristics can be found in Table 1. The regional ethics committee of Western Norway has approved our studies (REC West no. 78/05). All participants provided written consent as requested by our ethics committee.

Kidney Biopsies
Core biopsies have been obtained by O.E., L.L. and T.S. with a 16g needle from 16 patients undergoing (partial) nephrectomy in the operating room itself exactly at the time of surgery. Four paired biopsies from each patient with histologically-confirmed clear cell renal cell carcinoma (ccRCC) and adjacent non-tumorous ("normal") tissue were either stored as FFPE tissue or in an RNA-stabilizing agent (RNAlater 1 , Qiagen, Germany). Total RNA was extracted with miRNeasy FFPE kit or miRNeasy micro kit (Qiagen), respectively.

RNA Library Preparation and Sequencing
RNA sequencing libraries were prepared using TruSeq RNA Access library kit (Illumina, Inc., San Diego, CA, USA) according to the manufacturer`s protocol.
Initially total RNA concentration was measured using Qubit 1 RNA HS Assay Kit on a Qubit 1 2.0 Fluorometer (Thermo Fisher Scientific Inc., Waltham, MA, USA). Integrity was assessed using Agilent RNA 6000 Nano Kit on a 2100 Bioanalyzer instrument (Agilent Technologies, Santa Clara, CA, USA) and the percentages of fragments larger than 200 nucleotides were calculated.
Thereafter, RNA samples (100 ng total RNA) were fragmented at 94°C for 8 minutes on a thermal cycler. First strand cDNA syntheses were performed at 25°C for 10 minutes, 42°C for 15 minutes and 70°C for 15 minutes, using random hexameres and SuperScript II Reverse Transcriptase (Thermo Fisher Scientific Inc., Waltham, MA, USA). In a second strand cDNA synthesis the RNA templates were removed and a second replacement strand was generated by incorporation dUTP (in place of dTTP, to keep strand information) to generate ds cDNA. AMPure XP beads (Beckman Coulter, Inc., Indianapolis, IN, USA) were used to clean up the blunt-ended cDNA from the second strand reaction mix. The 3`ends of the cDNA were then adenylated to facilitate adaptor ligation in the next step. After ligation of indexing adaptors, AMPure XP beads were used to clean up the libraries. In a first PCR amplification step, PCR (15 cycles of 98°C for 10 seconds, 60°C for 30 seconds and 72°C for 30 seconds) were used to selectively enrich those DNA fragments that have adapter molecules on both ends and to amplify the amount of DNA in the library. After validation of the libraries, using Agilent DNA 1000 kit on a 2100 Bioanalyzer instrument, the first hybridization step were performed using exome capture probes. Before hybridization a 4-plex pool of libraries were made, by combining 200 ng of each DNA library. The hybridization was performed by 18 cycles of 1 minute incubation, starting at 94°C, and then decreasing 2°C per cycle. Then streptavidin coated magnetic beads were used to capture probes hybridized to the target regions. The enriched libraries were then eluted from the beads and prepared for a second round of hybridization. This second hybridization (18 cycles of 1 minute incubation, starting at 94°C, and then decreasing 2°C per cycle) were required to ensure high specificity of the capture regions. A second capture with streptavidin coated beads were performed, followed by two heated wash procedures to remove non-specific binding form the beads. The enriched libraries where then eluted from the beads and cleaned up by AMPure XP beads prior to a second PCR amplification. The amplification step were performed by 10 cycles (98°C for 10 seconds, 60°C for 30 seconds and 72°C for 30 seconds) followed by a second PCR clean up using AMPure XP beads. Finally, the libraries were quantitated by qPCR using KAPA Library Quantification Kit-Illumina/ABI Prism 1 (Kapa Biosystems, Inc., Wilmington, MA, USA) and validated using Agilent High Sensitivity DNA Kit on a Bioanalyzer. The size range of the DNA fragments were measured to be in the range of 200-650 bp and peaked around 270 bp.
Libraries were normalized to 22 pM and subjected to cluster and single read sequencing was performed for 50 cycles on a HiSeq2500 instrument (Illumina, Inc. San Diego, CA, USA), according to the manufacturer's instructions. Base calling were done on the HiSeq instrument by RTA 1.17.21.3. FASTQ files were generated using CASAVA 1.8.2 (Illumina, Inc. San Diego, CA, USA). Data are available in the repository Gene Expression Omnibus, http://www.ncbi. nlm.nih.gov/geo/query/acc.cgi?acc=GSE76207.

Statistics and NGS Data Processing
We have a sample size of 64 samples, which is equivalent to 32 paired samples (tumor samples vs. normal samples). Within both the FFPE and in the RNAlater dataset, we have 16 sample pairs (tumors vs. normals). This sample size is sufficient to achieve a power of 0.85, where we apply a standard deviation of 0.7 of the expressed genes, an effect size of 2, and an alpha of 0.05 (R package RNASeqPower in https://www.bioconductor.org).
Assembly of reads and alignment of the contigs to the Human genome assembly GRCh38 was guided by Tophat and Bowtie. An empirical expression filter was applied, which left genes with more than 15 counts per million (cpm) in more than 8 samples per dataset. Comparative analysis was done using voom/Limma R-package. Differential gene expression was defined as Benjamini-Hochberg adjusted p-value 0.05, and an absolute fold change of 2. Pathway analysis was performed with Ingenuity Pathway Analysis (Qiagen, USA; version 24718999). The Ingenuity Knowledge Base information was used as reference set. Canonical pathways were sorted by smallest Benjamini-Hochberg-adjusted p-value.
Classifier analysis was performed with the KNNX Validation package in GenePattern (http://www.broadinstitute.org/cancer/software/genepattern). Leave-one-out method was used as internal cross validation method. Euclidean distance was used as distance measure, where three neighbors were considered. Data visualization was performed with JMP Pro 11 (www.sas. com), and Graphpad (www.graphpad.com).
ELISA for CA9 Serum Levels CA9 serum concentrations of 38 patients was measured using the Quantikine Human Carbonic Anhydrase IX Immunoassay (R&D Systems, Minneapolis, USA, catalogue number DCA900) according to instructions of the manufacturer, but with an overnight incubation at 4°C after having added the serum. Results were assessed with the Kruskal-Wallis and Dunn's test [44].