The test cohort was composed of 60 freshly frozen samples derived from renal resections of 30 ccRCC patients provided by the Tumour Bank of the Department of Pathology, Albert Szent-Györgyi Medical School, University of Szeged. The sample size of the cohort was not determined statistically prior to experimentation. All diagnoses were histologically proven and classified according to the 2016 WHO classification of renal neoplasms [2]. Matched healthy tissue parts of the resected renal portions were also pathologically examined therefore considered as “normal” samples. Patient data are summarised in Table 1. All patients signed written informed consent forms. The study was performed with the approval of the Scientific and Research Ethical Committee of Hungarian Scientific Council (ETT TUKEB, #IV/5376-2/2020 EKU), and the experiments conform to the Declaration of Helsinki in 1995 (revised in Edinburgh in 2000).

Total RNA isolates were prepared from fast-frozen tumour sections using Promega ReliaPrep Total RNA Miniprep System (Promega, Madison, WI), according to the manufacturer’s protocol. The amounts and purity of RNA preparations were analysed by NanoDrop Spectrophotometer (Thermo Fischer Scientific, Waltham, MA, United States). First-strand cDNA synthesis reactions were performed by SuperScript VILO cDNA Synthesis Kit (Invitrogen, Thermo Fischer Scientific, Waltham, MA, United States) using 400 ng of each RNA sample, according to the producer’s protocol.

Genomic DNA isolates were prepared from fast-frozen tumour sections using Macherey Nagel NucleoSpin Tissue Kit (Macherey Nagel, Düren, Germany), according to the manufacturer’s protocol. The amounts and purity of the DNA preparations were determined by NanoDrop Spectrophotometer (Thermo Fischer Scientific, Waltham, MA, United States), and 20 ng of each sample were used for further application.

The qPCR reactions were performed on Thermo Scientific PIKO 96-well Thermal Cycler (Thermo Fischer Scientific, Waltham, MA, United States) platform using SYBR Green chemistry (GoTaq qPCR Master Mix, Promega, Madison, WI, United States). Sequences of gene-specific primer pairs are detailed in Table 2. Sample absolute gene expression and copy number variant (CNV) Cq values were normalised to the Cq average of 18S rRNA and cyclophilin B internal control genes or an intergenic region, respectively. These values are referred to as normalised gene expression and normalised CNV in the text, respectively. Relative gene expressions and CNV values were calculated by the −ΔΔCq method referred to as relative gene expression and relative CNV in the text, respectively [16]. For the mRNA and CNV measurement of VHL, the same oligo pair was used in DNase treated RNA and RNA free gDNA samples, respectively.

Statistical comparisons of normal vs. tumour datasets were performed using SigmaPlot 12.5 software (Systat Software Inc., San Jose, CA, United States). Regression between the gene doses and gene expression was determined by using a Polynomial linear equation. The distribution of datasets was examined by Shapiro-Wilk normality and Equal Variance tests. The values of variances were determined by ANOVA or Kruskal-Wallis ANOVA on Ranks methods followed by multiple comparisons using Holm-Sidak test or Tukey-test, respectively. Diagnostic abilities and cutoff values were defined by the receiver operating characteristic curve (ROC curve) analysis tool of SigmaPlot 12.5 software using a 99% confidential interval.

Tumorigenesis of ccRCC is highly linked to alterations of somatic drivers of the 3p21 and 3p25 chromosome regions, namely VHL, SETD2, PBRM1 and BAP1 genes. These alterations include single nucleotide variations (SNVs), small insertions or deletions (INDELs), dinucleotide substitutions (DNVs) and somatic copy number alterations (SCNAs) like the deletion of chromosomal arm 3p as well as the alteration of DNA methylation patterns [8]. The loss of large genomic parts can be easily traced by determining the gene doses. Analysing the relative CNV results, we found double doses of VHL in two, and half gene doses in four ccRCC samples out of the 30 tumour samples tested (Table 4, column 3). Single copies of SETD2 and BAP1 were detected in three and six cases, respectively (Table 4 column 5 and 9, respectively). Finally, double doses of PBRM1 were observed in four, and half in two samples. Most of these alterations afflicted only one gene out of the four; multiple gene hits were found exclusively in 2 cases (Table 4, column 7). However, we could not find any significant difference between the two categories when comparing the normalised CNVs or the medians of the four somatic drivers in ccRCC tumour tissues and their normal counterparts (Figures 1A–E). Our findings reflect that in most of the examined cases, the genetic alteration of chromosome 3p arm was other than large range deletion. Indeed, in some patients, point- and frameshift mutations, small deletions or insertions of VHL were proven by targeted sequencing (Table 3).

Manifestation of genetic or epigenetic alterations of a given locus in transcript or protein levels can be regulated in multiple manners. In this sense, such dysregulation on genomic modifications could be the primary cause of impaired final cellular processes resulting incomplete cellular reprogramming. Thus, to better understand tumour behaviour, it would be necessary to identify genetic and even transcriptional changes occurring during cancer formation and progression. To carry out such a multilevel investigation, we measured the mRNA levels of VHL, SETD2, PBRM1 and BAP1 genes by qPCR on the same tumour-normal paired sample set used for our CNV characterisation. At first, we carried out a pairwise comparison of gene expressions in each patient sample. The relative mRNA levels of VHL showed a marked reduction in 24 out of the 30 tumours, which were not suspected by the data of gene copy numbers. In five patients among the six outliers, we did not detect a significant difference between (from 0.73- to 1.55-fold change) the tumorous and normal parts of the kidney. However, in one case, we measured a dramatic, 5.75-fold relative mRNA level increase of VHL in the ccRCC tissue sections despite that this patient had a normal copy number of VHL (Table 4, column 2 vs. 3). Similarly, the transcription of PBRM1 was not altered, or just low level changed in 8 patients (between 0.73- and 1.59-fold change) and showed significant downregulation in 22 patients. As we found in the case of VHL, the relative mRNA levels of PBRM1 did not correlate with the detected gene copy numbers (Table 4, column 6 vs. 7). A far more BAP1 gene expression activity reduction was observed in 26 patient samples (Table 4, column 8 vs. 9). We also detected the parallel decrease of gene doses and transcription levels in only a few cases. The reduction of SETD2 mRNA level was the most dramatic, with relative values between 0.04 and 0.49, which again did not reflect the detected copy numbers of the gene. Unaltered transcription was detected only in two patients (Table 4, column 4 vs. 5). Finally, we concluded that almost all the tested tumour samples showed reduced gene expression of at least one out of the four genes compared to their normal tissue pairs. Only one exceptional patient sample (patient #11) identified with unchanged and either upregulated transcription activity of SETD2, PBRM1, BAP1 and VHL, respectively, in his ccRCC sections. Nevertheless, none of the general pathological tests showed any particular characteristic features that this ccRCC sample or the indicated patient possessed.

We performed regression analyses on each gene and the medians of the detected values to evaluate the statistical correlation between the detected relative gene expression and -CNV results. We found that none of the genes nor their medians showed a linear correlation between the relative transcription activity and gene doses (Figures 2A–E, respectively) which further verify that in most of the examined cases, the genetic alteration was other than large range deletion, such as epigenetic dysregulation or point mutation [17–19].

To have a more specific sight of the manifested transcription alterations of the four genes, we also compared the normalised absolute gene expression levels detected in the normal and tumorous tissue samples. The distribution of the values indicated highly significant differences (p < 0.001) of each of the examined genes in control versus tumour comparisons (Figure 3). The mRNA levels of PBRM1 showed a week overlap between the normal and tumorous samples (Figure 3C), while those of BAP1 and VHL grouped more discretely (Figures 3D,A, respectively). Additionally, the transcription of SETD2 was sharply discriminated between the two groups as the middle half of the gene expression values aggregated in distances, and only some extreme values overlayed (Figure 3B). The normalised mRNA level medians of the four genes were also significantly different in the tumorous samples compared to the normal ones with negligible overlapping (Figure 3E).

The loss or alteration of the 9p21.3 chromosome region, mainly affecting the gene of cyclin-dependent kinase inhibitor 2A (CDKN2A), is highly associated with the formation of ccRCC originated metastatic lesions [9]. To test whether the gene expression of one of the several transcription variants of CDKN2A, p14 ^ARF reflects ccRCC-associated changes, we determined and compared the normalised mRNA levels of each normal and tumour sample. Similarly to SETD2, we observed almost non-overlapping sharply different transcription levels between the two groups (p < 0.001) (Figure 4A). Nonetheless, dislike the genes of chromosome 3p, the relative gene expression of p14 ^ARF showed dramatic overexpression in most of the tumorous samples (Table 5), indicating that almost all the patients were susceptible for metastases formation; however, the pathological records did not confirm this. The normalised gene expression medians of tumorous sections completed with the data of p14 ^ARF were still significantly distinct from those of the normal tissues (p < 0.001) (Figure 4B). However, the values of the mRNA level medians of the four genes of 3p regions grouped noticeable more precisely in the two categories (Figure 3E vs. Figure 4B).

Our data indicate that the gene expression pattern derived from the transcription alteration of 3p21, 3p25 and 9p21.3 chromosome regions is different in normal kidney and ccRCC tumour tissues. Particularly, examining mRNA levels of VHL, SETD2, PBRM1, and BAP1 serve as discrete signatures for identifying ccRCC patients.

Recently in tumour biology, many efforts have been made to identify valuable RNA signatures applied in diagnostics. Our results raised the further question of whether the transcription signature of VHL, SETD2, PBRM1 and BAP1 as a molecular marker has a valuable potential in ccRCC molecular diagnostics. We performed receiver operating characteristic (ROC) curves to calculate the sensitivity and specificity of the normalised gene expression data obtained from our cohort study (Figures 5A,B). Although SETD2 was likely the most promising of the four genes, based on its data distribution (Figure 3B), the transcription pattern of BAP1 showed the highest sensitivity and specificity, possessing a 0.92 area under the ROC curve. According to the statistical evaluation, measuring the mRNA level of BAP1 is suitable for distinguishing ccRCC tumour specimens from normal kidney tissues with only an 8 per cent false rate (in 99% CI) (Figure 5A, green line). On the other hand, SETD2 showed a modest utility in diagnosis with an 87% likelihood of categorising a given patient sample (Figure 5A, red line). Using PBRM1 and VHL revealed the most inaccurate 22% and 21% false rates, respectively (Figure 5A, black and yellow lines, respectively). Regarding that, the distribution of the normalised gene expression data (box plots, see the previous section) is not in concert with the result of ROC analyses, using only the best candidate, BAP1, as a diagnostic marker may be oversimplified. As the four genes mapped on two distinct chromosomal regions of chromosome 3, applying the median normalised mRNA levels of VHL, SETD2, PBRM1 and BAP1 as transcription signatures is reasonable. Performing a ROC curve analysis on these datasets, 82% of sensitivity and specificity was revealed (Figure 5A, blue line).

As we proved that the gene expression of CDKN2A transcript variant, p14 ^ARF , could be the indicator of alteration of an independent genomic locus involved in the progression of tumorigenesis of ccRCC, we investigated the potential diagnostical value of it as well. The area under the ROC curve indicated that statistically, using the normalised mRNA level of p14 ^ARF, 95% of the cases could be identified correctly (Figure 5B). Nevertheless, if the median values of gene expression were complemented with the data of p14 ^ARF , the sensitivity and specificity of the test dropped down to 79% (data not shown). This surprising result could be that opposite to the significantly downregulated transcription of the four genes of the 3p chromosome, the mRNA levels of p14 ^ARF were dramatically elevated in ccRCC tumours compared to the normal specimens mitigating the effectiveness of the test.

Using ROC curve analysis, possible diagnostical cutoff values can be evaluated as well. These cutoffs could serve as thresholds in the differentiation of normal kidneys from ccRCC tumorous tissues. The optimal cutoff values of VHL, SETD2, PBRM1 and BAP1 and their medians were determined from the pre-test probability and cost ratio by SigmaPlot to minimise negative and positive false results (Figure 5C). According to our results, PBRM1 proved to be the less reliable marker. Only 57% of the samples could be sorted in the correct category, using the ideal threshold value of normalised absolute gene expression, −7.04 (Figure 5C first and second lanes). Surprisingly, the sensitivity and specificity of SETD2 testing capacity (77% and 83%, respectively) statistically were predicted lower (Figure 5C second and third lanes) than expected based on the distribution of its expression. The highest sensitivity, 0.93, was evaluated for VHL with a cutoff value of −7.39. Nonetheless, determination of normalised mRNA level of VHL classifies only 73% of the specimens correctly (Figure 5C, seventh and eighth lanes). The testing ability of the transcription measurement of BAP1 appeared to be the most reliable as it was also proved by its ROC curve. 83% of sensitivity can be achieved along with 97% of specificity, using the calculated cutoff, −9, of BAP1 normalised mRNA levels as differentiation value which is the highest one among the four genes we tested (Figure 5C, fifth and sixth lanes). We believe that examining the transcription signature of the four genes encoded in 3p25 and 3p21 gene regions could serve as a more accurate indicator than determining the gene expression of a single but very reliable gene. For distinguishing ccRCC specimens from normal renal tissues, we calculated the cutoff value of median data of normalised mRNA levels. We found that using the cutoff as a threshold value (−7.62) of the median of the normalised transcription vales could discriminate the normal and cancerous categories with 87% of sensitivity and 77% of specificity, indicating a valuable possible diagnostic value of this approach (Figure 5C, ninth and tenth lane).

In our study, the transcription alteration of p14 ^ARF was also a candidate for a diagnostic marker. Our findings showed that the mRNA level of p14 ^ARF is significantly different in ccRCC specimens compared to normal kidney tissues. Thus, we determined the cutoff value of the p14 ^ARF normalised gene expression dataset as well. We found that using the calculated cutoff (−9.5) as a testing threshold, we can reach 100% of sensitivity but only 65% of specificity statistically (Figure 5D). Furthermore, we showed that the transcription level of p14 ^ARF was elevated in most of the ccRCC tumour samples opposite those of the four genes of 3p25 and 3p21 chromosome regions. For this reason, we suppose that the transcription evaluation of p14 ^ARF might be used as an independent indicator as a supplementer for the gene expression signature of VHL, SETD2, PBRM1 and BAP1.

Our results indicate that besides the single genetic alterations, epigenetic changes also have to be considered, which can occur during tumorigenesis, and are worth being characterised in multiple ways. Furthermore, observation of the transcriptional level of specific genes in the affected genomic loci can indicate the manifestation of mutations and chromosomal arrangements, as well as epigenetic dysregulation, which has a high impact on the entire behaviour of a given tumour. Therefore, we believe that monitoring the gene expression signature of tumour driver genes gives a unique opportunity to use more valuable diagnostic markers for supplementing the existing genetic testing.