Patterns of gene expression characterize T1 and T3 clear cell renal cell carcinoma subtypes

Renal carcinoma is the 20th most common cancer worldwide. Clear cell renal cell carcinoma is the most frequent type of renal cancer. Even in patients diagnosed at an early stage, characteristics of disease progression remain heterogeneous. Up-to-date molecular classifications stratify the ccRCC samples into two clusters. We analyzed gene expression in 23 T1 or T3 ccRCC samples. Unsupervised clustering divided this group into three clusters, two of them contained pure T1 or T3 samples while one contained a mixed group. We defined a group of 36 genes that discriminate the mixed cluster. This gene set could be associated with tumor classification into a higher stage and it contained significant number of genes coding for molecular transporters, channel and transmembrane proteins. External data from TCGA used to test our findings confirmed that the expression levels of those 36 genes varied significantly between T1 and T3 tumors. In conclusion, we found a clustering pattern of gene expression, informative for heterogeneity among T1 and T3 tumors of clear cell renal cell carcinoma.


Introduction
Renal tumors are classified as the 20th most common malignancy worldwide, both based on incidence and death rates [1]. Clear cell renal cell carcinoma (ccRCC) is the most frequent renal tumor (80-90% of cases) [2,3]. Multiple morphotypes have been described within RCC [4,5] and a growing body of evidence suggests that those morphotypes represent different molecular entities [6][7][8].
There are several classification systems used to describe renal tumors. Grading is performed by Fuhrman system, based on the nuclear and nucleolar features, and recently modernized by International Society of Urologic Pathology [4]. The most important for prognosis is the stage of the tumor which is evaluated by American Joint Committee on Cancer / The Union for International Cancer Control (AJCC/UICC) TNM system [9]. Although ccRCC cases are usually diagnosed at early stages (in TCGA database, T1 stage represents 48% of all ccRCC cases), a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 clinical outcomes remain heterogeneous within each staging group, suggesting the existence of molecular features unaccounted for by pathology assessment [6,7]. A significant challenge is that metastatic potential and clinical outcome are not well correlated with tumor size and stage [6,7].
In the up-to-date molecular classifications, ccRCC samples are classified into two groups [10]. The authors annotate those clusters ccA and ccB and state that ccA tumors have markedly improved disease-specific survival compared to ccB. Their analyses suggests that the proposed classification was independently associated with survival. However, the heterogeneity within described clusters is significant. An important step in progression of cancer is extension of the tumor beyond the natural limits of the affected organ. In current classification, T1 and T2 tumors differ by their size only, and both are confined to the kidney, while both T3 and T4 tumors extend beyond this organ. Therefore, we decided to select T1 and T3 samples for our study. We aim to verify whether gene expression patterns reflect stage of the disease and to investigate the heterogeneity based on gene expression within the current classification systems in T1 and T3 tumors. Gene expression in ccRCC was studied extensively in the past (exemplary papers: [11][12][13]). Our study provides additional information on heterogeneity within the samples from various tumor stages as well as points out towards potential mechanisms of transition between these stages. 23 ccRCC tumor samples were collected during radical nephrectomy at the Department of Urology, JUMC. Samples were fixed with formalin and embedded in paraffin at the Department of Pathology for microscopic evaluation and transferred to the Center for Medical Genomics OMICRON for gene expression studies. The study was approved by the Bioethics Committee of the Jagiellonian University.

Sample collection
All patients signed written informed consent forms. Experiments conform to the provisions of the Declaration of Helsinki in 1995 (as revised in Edinburgh 2000). Patient tumors were classified into T1 (13) or T3 (10) stages by a pathologist and independently re-evaluated. Selection of T1 and T3 tumors, as a basis of sample collection for our study, gave prospect to investigate clinically most frequent specimens. In addition, each study group remains homogeneous and sample selection parallels kidney restriction of the tumors in T1 group and extension beyond the kidney in T3 group. Additional clinical data were collected, along with immunohistochemical information summarized in S1 Table. RNA isolation RNA was isolated from 10 x 5 μm slides from Formaldehyde Fixed-Paraffin Embedded (FFPE) block, using Maxwell 16 FFPE Tissue LEV DNA Purification Kit (Promega). Briefly, 300 μl of Mineral Oil and 250 μl of lysis master mix were added per sample and incubated in 56˚C for 15 min and subsequently at 80˚C for 1 hour. DNA was degraded by DNase I treatment (15 min, RT). The aqueous phase was transferred to Maxwell FFPE Cartridge and RNA was isolated according to Promega RNA-FFPE protocol. 50 μl of Nuclease-Free Water was used for RNA elution. The RNA quantity was measured using NanoDrop 1000 (Thermo Scientific) device and quality was assessed on 2200 TapeStation System (Agilent, RNA ScreenTape), according to manufacturer instructions. DV200 parameter, describing percentage of RNA fragments longer than 200 bases was used for sample classification (S1 Table). Samples with DV200 > 30% were classifies as suitable for further analysis.

Whole genome DASL assay
The Illumina Whole Genome-DASL assay was performed using 200 ng of RNA following the manufacturer's instructions. Briefly, RNA was reverse transcribed to cDNA using biotinylated primers, followed by immobilization to streptavidin-conjugated paramagnetic particles. Biotinylated cDNAs were then simultaneously annealed to a set of assay-specific oligonucleotides. Extension and ligation of the annealed oligonucleotides generated PCR templates that were amplified using Titanium Taq DNA Polymerase (Clontech). Labeled PCR products were washed and denatured to yield single-stranded fluorescent molecules, which were hybridized to the HumanHT12 v4.0 Whole Genome Gene Expression BeadChips for 16 h at 58˚C. The Illumina HiScan was used to scan the arrays.

Data analysis
Microarray data in � .IDAT format were uploaded and pre-processed in R environment. The 'beadarray' package was used to upload the data and 'lumi' for normalization and filtration of the data.

Differential expression analysis
The differentially expressed probes were detected via the Generalized Linear Model framework implemented in the package 'limma'.
For the comparison between T1 and T3 groups as well as the groups defined via hierarchical clustering (A1, A2, A3) the functions 'contrast.fit' and 'eBayes' were used. For analysis of differential expression in the TCGA cohort the framework implemented in the package 'edgeR' was utilized. Gene counts were normalized with the default options and subsequently filtered to reduce the number of hypotheses tested. After estimating the dispersion parameter, the Generalized Linear Model was fitted and tests for coefficients were performed. Since this was used as a replication cohort, we have only recorded the number of genes differentially expressed between the two study groups with the standard level of statistical significance 0.05.

Hierarchical clustering
The 23 T1 and T3 samples were clustered based on expression of 543 probes. To this aim the function 'hclust' with complete linkage as implemented in the 'heatmap.2' procedure was used. The noticeable pattern where the dendrogram is divided into three main groups was further confirmed with the use of the 'cutreeDynamic' function in the 'dynamicTreeCut' package. The faithfulness of clustering was evaluated using the cophenetic correlation coefficient.
Both the T1 and the T3 samples were clustered based on normalized gene expression values (pseudocounts) generated with 'edgeR' package. To overcome the issue of the Euclidean metric being driven by highly expressed genes, the Renyi divergence function was used as the measure of similarity. Renyi divergence was previously used by the authors of [14] in the context of liver cancer. Once the similarity matrix was estimated, hierarchical clustering was performed as implemented in the function 'hclust'. The optimal number of clusters on each dendrogram, was established via analysis of gap statistics as implemented in the function 'clusGap'.
The ROC analysis was based on logistic models with the indicator of the event that the sample is T3 used as the dependent variable. For each of the probes used for hierarchical clustering 300 random training and testing sets were selected (each time the testing set was of size 7) and ROC as well as AUC was calculated as implemented in packages 'ROCR' and "OptimalCutpoints'. Subsequently, for each probe the median AUC was calculated for each sample (taken as the median AUC over all testing sets which contained a given sample). For each sample, the 'goodness' of classification was quantified as the median of these median AUC values over all probes.

Results
We analyzed 23 ccRCC samples on a microarray platform. Our samples belonged to T1 and T3 stages, as the T2 and T4 stages are rarely diagnosed (only 69 (13%) T2 and 11 (2%) T4 samples in TCGA database). Our main interests were to (1) test the hypothesis if gene expression reflects the histological classification of the JUMC samples (in particular, does the gene expression pattern allow for discrimination between T1 and T3 cases via unsupervised clustering) and (2) whether we will be able to find molecular features that reflect the observed diversity of disease progression.
Pathway Enrichment performed on the differentially expressed gene set (adjusted pvalue < 0.1, 481 genes) was narrowed down to those in level 3 in the Genome Ontology (GO) hierarchy. This returned a list of enriched terms, presented in S2 Fig. Further narrowing the results with 'Use GO Term Fusion' option reduced the list to 9 GO biological processes terms ( Fig 1A) including 'kidney development' with corrected p-value 1.31x10-3 (18 associated genes). Interestingly, genes associated with this term were downregulated in T3 samples vs. T1 samples. Analysis of genes with higher log-fold-change values and more stringent adjusted pvalue cut-off (0.05) ( Table 1) revealed one enriched pathway ( Fig 1B)-'response to copper ion' with three downregulated genes: aquaporin 1 (Colton blood group, AQP1), amine oxidase, copper containing 1 (AOC1) and aldolase B, Fructose-Bisphosphate (ALDOB).

Sample clustering
Unsupervised hierarchical clustering, based on expression of 481 genes, divided 23 T1 and T3 samples into three distinct clusters: A1, A2, and A3. Two of the clusters contain populations of T1 (A1) or T3 (A3) samples only, whereas the third cluster includes samples from both groups (A2). This three-cluster pattern (two 'pure' and one mixed) is not present when all (~34K) probes are used for analysis. Therefore, it is unlikely that it is due to a batch effects.
As the distances between the clusters suggests that the A2 cluster is more closely related with A3, despite containing samples from both T1 and T3, we aimed to investigate which expression profiles characterize the A2 group. A heatmap presenting relative gene expressions is shown in Fig 1C.

Dimension reduction by t-SNE algorithm in the context of sample clustering
To further test whether the pre-selection of features (based on differential expression) allows for faithful sample classification between T1 and T3, an additional machine-learning approach has been adapted. Three sets of probes were used in this analysis: (1) the probes used for hierarchical clustering (aligned to 481 genes); (2) top 40 differentially expressed probes, and (3) all 34476 probes. Subsequently, using these sets of features, samples were projected, using the t-SNE algorithm (see Methods section), on a 3-dimensional space. For the unbiased case (all probes) no association between tumor size and the three components is present. Interestingly, for the two sets of pre-selected features, not only do we see a separation between T1 and T3 samples in the 3D space, but also a separation between the three clusters defined in the previous section. The results are presented in S3 Fig. Additionally, to further test the three-cluster pattern, we applied the UMAP algorithm [16] to project the entire dataset (~33000 probes) onto a 10-dimensional space. Subsequently, we selected three dimensions for which the projection has the strongest association with the clinical diagnosis (T1 vs T3) and visualize the projected data. Interestingly, even in this agnostic approach (with features not being pre-selected) we see a further support for the 'intermediate cluster' to appear (see S4 Fig).

ROC-based classification of T1 and T3 samples in the context of sample clustering
To further test whether there are indeed samples more difficult to correctly classify as T1 or T3 (i.e. samples in the 'mixed cluster'), a ROC-based analysis was performed. For each of the probes aligned to 481 genes, the AUC was calculated for 300 random test subsets of size 7 for a (logistic) model fitted on the remaining 16 training samples. Subsequently, the median for each sample/probe was calculated and the median of these 500 number was assigned to each sample as a measure of 'goodness' of classification. Fig 2 includes violin plots for the 23 samples divided according to the three-cluster pattern. It is clear that the AUC in the 'mixed cluster' is lower than for the two remaining 'pure' clusters.

Differential variability and clustering faithfulness
In the current study, we use a relaxed threshold (p<0.1) in the process of selection of probes for sample clustering. We wish to support this choice by demonstrating that probes, which are differentially variable between the study groups are more informative about the clustering of samples than the ones with similar variances. To this aim we first compare the variances between T1 and T3 tumors using Levene's test and detect six probes (ILMN_1762410 (SLC22A2), ILMN_ 1716246 (FRZB), ILMN_1677851 (RARRES1), ILMN_1746128 (ACSM2B), ILMN_3311035 (miR-1251) and ILMN_1793309 (BEND4)) with FDR below the standard 0.05 significance threshold. Secondly, we compare the cophenetic correlation coefficient for two different clusterings: (1) based on differentially expressed probes (p<0.1) with p-value in the Levene's test above the median, and (2) based on differentially expressed probes (p<0.1) with p-value in the Levene's test below the median. We note that for the first set of probes the coefficient equals 0.72 and the second 0.87. Note, that in both of the above clusterings, we use the same number of probes for analysis.

Characterization of Intermediate Cluster
A2 vs A1. First the A2 cluster was compared to A1. In total 13 genes with adj. pvalue < 0.05 were found, with the largest log-fold-change = -1.98 achieved by interleukin 6.
A2 vs A3. Secondly, A2 and A3 clusters were compared. 22 differentially expressed genes (adj. p-value < 0.05) were identified and the top 15 (with |log-fold-change|>1.5) of them are presented in Table 2. In ClueGO analysis no enriched pathways with at least three genes were found.
A list of all differentially expressed probes between A2 and A3 is presented in S3 Table. Main groups/families of genes represented in the results are (trans)membrane proteins, ion- channel proteins or carrier proteins, suggesting a role of regulatory genes and modulation of signal transduction in the observed outcome heterogeneity. A3 vs A1. Differential expression analysis of A3 and A1 clusters revealed a larger set of differentially expressed genes than A1 vs A2 and A3 vs A2. A list of 58 down-regulated and 101 up-regulated probes with |logFoldChange| > |1.5| is presented in S4 Table. ClueGO-based analysis resulted in network depicted in Fig 3A and 3B. Interestingly, genes with lower expression in A3 are associated with morphogenesis and stress response related GO's and those that are overexpressed with metabolic processes.
Validation of results with TCGA data. Our sample size was relatively small, therefore we used TCGA RNA-seq data as a larger replication cohort. Of 481 genes, differentially expressed between T1 and T3 groups, 394 had expression levels available in the TCGA database. Illumina Probe IDs were converted to ENSG# using BioMart. A Gene was considered for further analysis if it was expressed in at least 80% of samples and the median read count exceeded 10.

T3 vs T1 comparison
Validation of our primary analysis revealed that almost 67% of differentially expressed genes (264; non-adj. p-value < 0.05) were also differentially expressed in the TCGA RNA-seq data. We additionally note that the correlation coefficient for the logFC's between the two cohorts equals 0.78, as presented in the

Hierarchical clustering with Renyi divergence
The TCGA cohort was further used to test the observations made with the use of unsupervised classification. T1 and T3 samples from TCGA were clustered based on each of the 394 differentially expressed genes in the UJ CM cohort. These 394 clusters were then evaluated, using Renyi Divergence measures, for heterogeneity with the expectation that those genes driving the clustering observed in the UJ CM cohort will show evidence of heterogeneity. To this aim, differential expression analysis was performed (between samples in a given cluster versus the largest, reference cluster). The results of this analysis were compared to the set of 36 genes which discriminate between A2 from A1 and A2 from A3.

T1
Using above-described procedure, T1 samples were divided into 8 clusters (C1-C8), where C1 was the largest and was further used as reference. Results of this analysis are presented in Table 3. Clusters 7 and 8 were excluded from further analysis due to sample size (i.e. the disproportion in the sample size in the case-control design versus the largest cluster).

T3
T3 samples were also divided into 8 clusters. Clusters 5 to 8 were excluded from further analysis based on sample size. Results of the analysis are presented in Table 3.

Discussion
Clear cell renal cell carcinoma is the most frequent kidney neoplasm in adults, comprising 80-90% cases of renal tumors [2]. A characteristic feature of ccRCC is large heterogeneity of individual survival times and disease outcomes, even within the same TNM classification groups. Existing pathological classifications do not reflect the molecular basis of the disease [10]. The inability to predict treatment outcome and metastasis in ccRCC could be attributed to the molecular heterogeneity of tumor cells [6,7,17]. Since the high molecular heterogeneity within staging groups could implicitly account for treatment outcome and disease recurrence, we investigated the molecular landscape of ccRCC. We characterized differences in gene expression patterns between T1 and T3 stages in search of genes associated with the molecular heterogeneity of tumors. This approach aimed to identify genes which would be altered between the pure and mixed group (A1 vs A2 and A3 vs A2). The detected genes are involved in regulatory processes and signal transduction. Therefore we hypothesize that the sample heterogeneity can be accounted for by accumulation of subtle deviations in metabolic processes caused by changes in gene expression. We repeated this analysis on the TCGA ccRCC cohort and confirmed 67% of our results. We verified the usability of the gene set to depict the molecular heterogeneity of ccRCC samples.
One of the main goals of the study was to emphasize heterogeneity of expression patterns in the context of discrimination between study groups. Therefore, as noted in the Results section (see Differential Variability and Clustering faithfulness) we choose to relax the statistical significance threshold (from 0.05 to 0.1) to include in further analysis genes which have more heterogeneous expression profiles in our cohort and thus higher chance of falling above the standard significance level. Patterns of gene expression in clear cell renal cell carcinoma Table 3

Cluster 5 (n = 42)
Cluster 6 (n = 11) Pathway enrichment analysis emphasized the role of copper metabolism, which is an important process in renal tissue in general, and has a role in cancer development. However, presented genes do not take direct part in pathways regarding those issues.
Other differentially expressed genes include molecular transporters (SLC22A12, SLC22A6, SLC22A2, SLC5A10), (trans)membrane proteins (AOC1, TMEM27, FLRT3, STEAP3, GPRC5A, TMEM145) and other channel proteins (TRPM3, AQP1) involved in regulation and signal transduction in cell metabolism and response to external stimuli. This suggests that dysregulation of signal transduction maybe important in defining the observed diversity of ccRCC outcomes.

Sample clustering
Clustering of 23 samples, based on all significantly differentially expressed genes, revealed partition of T1 and T3 samples into 3 distinct clusters (Fig 1C). Two of them (A1 and A3) contained only T1 and T3 samples respectively, whereas A2 contained samples from T1 and T3. Interestingly, the gene expression profiles in A2 show no clear pattern of up or downregulation, in contrast to the other two clusters. Therefore we aimed to identify genes involved in molecular heterogeneity-i.e. differentially expressed between A1 vs A2 or A3 vs A2.
Comparison of a A2 with A1 cluster revealed a role for IL-6. Overexpression of IL-6 is associated with enhanced invasiveness and epithelial-mesenchymal transition (EMT) and IL-6 is involved in a JAK/STAT signaling pathway [36]. Although there has been reported lack of correlation between expression of this protein and tumor size or grade [37] our analysis suggests another evidence on regulative role of IL-6 in clear cell renal cell carcinoma.  [48]) or kidney metabolism: TMEM171, CYS1 [49,50]. One of the two down-regulated genes-IGF2BP3 is not expressed in normal adult tissues and is known to promote tumor invasion and metastasis [45,51,52]. Some of these same genes were identified as differentially expressed in T3 vs. T1 comparison: HAO2, AOC1, SLC22A2, GBA3, TMEM27, SLC5A10, NPR3, PAX2, IGF2BP3.
The differences shown here lead us to postulate that the isolated intermediate cluster reflects the tumors that are less metastatic prone/aggressive. Several statistically significantly disturbed genes (IL6, GBA3, TMEM27) show contradictory expression change trend to expression changes described in the literature and associated with tumor progression and metastasis [32,38,43].
The 36 genes obtained from A2 vs A1 and A2 vs A3 comparisons code for proteins associated with intracellular signaling and metabolic processes, but lack driver genes or commonly known cancer master regulators, yet these modulators account for the observed sample heterogeneity. This is in line with the previous results of T3 vs T1 comparison and underlies the significance of regulatory/modulatory genes in the progression of the disease.

Validation of results in TCGA ccRCC cohort
Use of TCGA ccRCC cohort confirmed almost 70% of our results. We tested whether genes differentially expressed between A1/A3 and A2 can be used to measure heterogeneity in a larger dataset. For that purpose, we clustered TCGA ccRCC samples separately T1 and T3 and used the A2-specific genes as an input for differential expression. We found that expression changes in gradual fashion for T3 clusters (Fig 5) suggesting growing dysregulation. Interestingly, the largest T1 clusters (cluster 2 and 5) show contradictory changes in expression suggesting opposite directions of regulatory processes in these samples.
In conclusion we propose that expression of certain RNAs can be used to study the molecular basics of the heterogeneity of ccRCC.
We have found a clustering pattern reflecting heterogeneity of samples. Furthermore, we detected genes associated with diversity of ccRCC samples. We postulate that genes associated with regulatory or signal transduction modulation roles are related to diverse representations of ccRCC occurring regardless of the histological classifications. Further functional research is needed to test these observations. Supporting information S1 Table. Clinical parameters of analyzed samples. T, N, M-classification of samples, T_expanded T classification, diameter-measured in the widest dimension, Grade-ISUP modified Fuhrman grade, survival time-calculated as the number of days between collection date and date of death (calculated when applicable), mdm2 -result of histochemical staining of mdm2 protein, p53-result of histochemical staining of p53 protein, procedure-name of the procedure at which the sample was obtained, necrosis-was the tumor tissue necrotic, DV 200 -Illumina proposed parameter for description of quality of FFPE derived RNA samples (over 30% qualifies sample as sufficient for further analysis). (DOCX) S2 Table. List of all differentially expressed probes between T3/T1 comparison with adjusted p value under 0.01. ILMN ID-Illumina Probe ID, logFC-logFoldChange of probe expression, AveExpr-average expression of the given probe, P.Value-p value, adj.P.Val-p value adjusted for multiple testing. (DOCX) S3 Table. List of differentially expressed genes in A2 vs A1 and A2 vs A3 comparisons. All probes that reached adj. p. value < 0.05 cut-off value. ILMN ID-Illumina probe ID, logFC-log Fold Change, AveExpr-average probe expression value, P.Value-p value, adj.P.Val-p value adjusted for multiple testing. (DOCX) S4 Table. List of differentially expressed genes in A3 vs A1 comparison. All probes that reached adj. p. value < 0.05 and logFC > 1.5 cut-off values. ILMN ID-Illumina probe ID, logFC-log Fold Change, AveExpr-average probe expression value, P.Value-p value, adj.P. Val-p value adjusted for multiple testing. For the unbiased case (all probes) no association between tumor size and the three components is present. Interestingly, for the two sets of pre-selected features, not only do we see a separation between T1 and T3 samples in the 3D space, but also a separation between the three clusters defined in the previous section. (TIF) S4 Fig. Results of analysis using the UMAP algorithm. Entire dataset (~33000 probes) was projected onto a 10-dimensional space. Three dimensions for which the projection has the strongest association with the clinical diagnosis (T1 vs T3) was selected and projected data visualized. Interestingly, even in this agnostic approach (with features not being pre-selected) we see a further support for the '