Mutational heterogeneity between different regional tumour grades of clear cell renal cell carcinoma.

Only a limited number of studies have explored the possible associations between tumour grade and mutated genes in clear cell renal cell carcinoma (ccRCC), and we set out to investigate this further using a multiple sampling and next generation sequencing (NGS) approach in a series of ccRCCs. Multiple regions were sampled from formalin-fixated paraffin-embedded ccRCC tumour blocks from seven patients. In 27 samples from six patients, we performed targeted NGS using a custom 42-gene panel based on the most frequently mutated genes in ccRCC reported in public databases. In four samples from the seventh patient, we performed whole exome sequencing (WES) and array comparative genomic hybridisation for detection of copy number variants (CNVs). Mutated genes and the tumour grades of the samples in which they had been identified were compared both within and between all individual tumours. CNVs were compared across all samples from patient 7. We identified clear genetic heterogeneity within and across tumours, but VHL mutations were seen in all patients. Looking across all samples, we identified eleven genes that were only mutated in samples with one particular tumour grade. However, these genes were never mutated in all samples with that tumour grade. Increasing chromosomal instability corresponded with increasing tumour grade, but we observed minimal association between tumour grade and total mutational load in the WES data. Our study confirms the genetic heterogeneity and tumour grade heterogeneity of ccRCC. Although a relatively small number of samples was analysed, genes were identified that could potentially be specific, though insensitive, markers of higher ccRCC tumour grades.


Introduction
Kidney cancer is the 14th most common cancer worldwide, with clear cell renal cell carcinoma (ccRCC) accounting for 70% of total cases in adult patients (Bray et al., 2018;Moch et al., 2016). The global incidence of ccRCC was over 403,262 in 2018 with a mortality rate of 1.8 per population of 100,000 (Bray et al., 2018). Conventionally, tumour grading is used to guide clinical decisions on ccRCC patient management and is based on the highest dedifferentiated morphological features in the tumour (Dagher et al., 2017;Delahunt et al., 2014). Due to the rapid developments in molecular diagnostics, mutational profiling is being introduced clinically for a range of tumour types. Although tumour morphology has been shown to correlate with both tumour progression and prognosis (Dagher et al., 2017), the prognostic value of mutational profiling in ccRCC is less clear.
In the next generation sequencing era, the most frequently mutated genes in ccRCC have been described in several large studies such as The Cancer Genome Atlas project (Network, 2013). Subsequently, the mutation rank system was introduced, in which mutations were scored and ranked by their functional impact based on in silico analysis using gene variant predictor tools (Gonzalez-Perez et al., 2013). Using this ranking, some genes were found to harbour more frequent mutations with higher predicted functional impact compared to other genes, and were referred to as cancer driver genes. For ccRCC, several cancer driver genes were identified by these studies, including VHL, PBMR1, SETD2, BAP1, and MTOR. However, knowledge on the associations between https://doi.org/10.1016/j.yexmp.2020.104431 Received 8 January 2020; Received in revised form 20 February 2020; Accepted 28 March 2020 these driver genes and the prognosis for ccRCC patients is still very incomplete. So far, in a few studies, inactivating mutations in TP53, and the ccRCC driver genes BAP1, SETD2 have been shown to be correlated with the prognosis of ccRCC patients (Hakimi et al., 2013;Manley et al., 2017).
One of the challenges in studying the relationship between tumour characteristics, including the somatic driver gene variants, and clinical outcome, is intratumour heterogeneity. This heterogeneity exists both at the level of histomorphological tumour grading and at the level of the mutational profile (Singh et al., 2015). In the last decade, several studies have shown intratumour heterogeneity in ccRCC at the genomic level (Gerlinger et al., 2014;Gerlinger et al., 2012;Mitchell et al., 2018;Turajlic et al., 2018). However, these studies did not examine this mutational heterogeneity in the context of regional tumour grading. A limited number of studies have reported on the association of the mutation profile with the specific morphological features; e.g. rhabdoid and sarcomatoid differentiation (Bi et al., 2016;Singh et al., 2015). As it is of interest to know if regional tumour grading, which plays an important role in current tumour classification and ccRCC prognosis, matches certain regional mutational profiles, we set out to study this in a series of six ccRCC patients using targeted sequencing and in one patient using whole exome sequencing (WES). Our objective was to investigate whether mutational profiles would be distinct between regions with distinct grades within individual tumours and if, in the total series of tumours, specific genes would uniquely be mutated in higher tumour grades.

Patient materials and tumour grading
Formalin-fixed paraffin-embedded (FFPE) tissues were collected from nephrectomy specimens of seven patients with ccRCC. Hematoxylin and Eosin (H&E) stained tissue sections were examined by two pathologists (P.F. and G.K-U.). Each section within those tumours with a homogenous histomorphology was assigned a regional tumour grade using the 4-grade system of the World Health Organization/ International Society of Urological Pathology (WHO/ISUP grading) (Dagher et al., 2017;Moch et al., 2016). For each grade in a patient, homogeneous tissue sections were collected using macrodissection from all available blocks without degraded tissue. In total, DNA was isolated from 31 tumour regions and seven matched normal kidney FFPE samples from seven patients (Table 1, Fig. 1). From each patient, the healthy part of the kidney, as indicated by the absence of tumour cells in the selected FFPE block, was taken as a source of normal DNA. The use of human material and clinical data from patients in this study has been approved by the Board Medical and Health Research Ethics Committee at the Faculty of Medicine, Universitas Gadjah Mada (UGM), Indonesia, as recognized by the Forum for Ethical Review Committee in Asia and the Western Pacific (FERCAP). Approval was given by the committee on September 4, 2013 (Reference: KE/FK/795/ EC). The study was carried out in accordance with the University Medical Centre Groningen Medical Ethical Review board guidelines (project number 20190251, filed January 4, 2016) and Dutch ethical guidelines and laws, and complied with the regulations stated in the Declaration of Helsinki.

Targeted sequencing and whole exome sequencing
The targeted sequencing protocol is based on the Single Primed Enrichment Technique as implemented in the Ovation™ Target Enrichment System (NuGEN, San Carlos, CA, USA). Our in-house designed set of landing probe areas covers the entire consensus coding region of the 42 gene-panel which include the most frequently mutated genes in ccRCC (Forbes et al., 2011;Network, 2013) (https://cancer. sanger.ac.uk/cosmic, database accessed 23 June 2014) and genes associated with the VHL/HIF pathway and the (PI3K)-AKT-MTOR pathway in ccRCC with a frequency ≥ 1% (Sato et al., 2013) (see  Supplementary materials and methods, Supplementary Table S1). This set includes 12 ccRCC cancer driver genes in a total of 31 acknowledged cancer driver genes (https://www.intogen.org/search, release 2014.12) (Gonzalez-Perez et al., 2013). Library preparation was done according to the manufacturer's protocol, starting from 500 ng of DNA. Enriched libraries were sequenced with the Illumina HISEQ 2500™ (Illumina, San Diego, CA, USA) using single-end sequencing with 100 bp reads.

Sequence data analysis and somatic mutation identification
Data analysis was done using a pipeline based on the Genome Analysis Toolkit (GATK) best practice recommendation (McKenna et al., 2010). For both targeted sequencing and WES, two different variant calling algorithms have been used: HaplotypeCaller from GATK and FreeBayes (Garrison and Gabor, 2012;McKenna et al., 2010). Called variants were annotated and filtered to identify true somatic mutations (see Supplementary materials and methods).
For each patient, the mutations present in each tumour region were classified into major clonal, minor clonal, absent, or inconclusive using the following approach. For each sample, we first determined the somatic variant with the highest mutant read frequency (MRF). Variants with a total number of mutant reads ≥ 5 and an MRF ≥ 50% of the highest MRF seen for that sample, were considered to be major clonal variants likely to be present in the majority of the tumour cells. Mutations with a total number of mutant reads equal to 3 or 4, or with a total number of mutant reads ≥ 5 and an MRF < 50% of the highest MRF, were defined as minor clonal variants likely to be present in a minority of the tumour cells. If none of the above criteria were met and the total read count was ≥ 10, the mutation was defined as absent. If none of the above criteria were met and the total read count was < 10, the mutation was considered inconclusive. Only variants with a major clone in at least one sample from a patient were included in the final analysis. Matched normal kidney samples were included to remove personal variants in the targeted sequencing and WES data. The Integrative Genomic Viewer was used to confirm the authenticity of the identified somatic mutations (Robinson et al., 2011). Sequencing data are available in the European Nucleotide Archive repository (accession number PRJEB37556).

Copy number variation analysis
ArrayCGH was carried out using 500 ng genomic DNA from four tumour samples from Patient 7 using the Complete Genomic SureTag DNA Enzymatic Labelling Kit protocol and an OligoaCGH/ChIP-on-Chip Hybridization kit (Agilent, Santa Clara, CA, USA). Details of the protocol have been reported previously (Ferronika et al., 2019).

Data interpretation
First, we analysed the distribution of mutations between different tumour regions per patient. We then grouped all 31 samples based on their regional tumour grade and analysed the distribution of mutations between the four different tumour grade groups. A Venn diagram was made to visualize the association of the mutated genes with specific tumour grades. P. Ferronika, et al. Experimental and Molecular Pathology 115 (2020) 104431 3. Results

Mutations identified by targeted sequencing
The 42 genes of our targeted panel were selected based on their high mutational frequency in ccRCC, as described in several studies and in The Cancer Genome Atlas Network database (Dalgliesh et al., 2010;Duns et al., 2012;Forbes et al., 2011;Network, 2013;Pena-Llopis et al., 2012;Sato et al., 2013), and from genes associated with the VHL/HIF pathway and the (PI3K)-AKT-MTOR pathway in ccRCC (Sato et al., 2013). The landing probes designed for the 42 genes are given in Supplementary Table S1. The average distance between different landing probes was optimized for FFPE based on preliminary data obtained with a catalogue assay from NuGEN (data not shown).
With the targeted sequencing assay, we analysed 27 tumour regions and six matched normal kidney tissues from six patients. The average tumour cell percentage within the samples was 70% in Patient 1, 85% (range 80-90%) in Patient 2, 85% (range 80-85%) in Patient 3, 78% (range 75-80%) in Patient 4, 83% (range 70-90%) in Patient 5, 70% in Patient 6, and 75% (range 70-80%) in Patient 7. Sequencing resulted in a mean target coverage of 45× (Supplementary Table S2). An overview of the detected mutations is given in Fig. 2 and Supplementary Table  S3. For each patient, we focused on the variants that were identified as a major mutation in at least one of the sections from that patient.
We analysed the distribution of mutations among different tumour regions in each patient. For Patient 1, four grade 2 regions and one grade 3 region were analysed from the primary tumour (Table 1). Although the samples were from five different locations within the primary tumour, we observed little variation in their mutational spectrum for our selection of 42 genes ( Fig. 2A, Supplementary Table S3). All five regions shared identical mutations in two cancer driver genes: VHL and ERBB3. The VHL mutation was present as a major clone in all regions, whereas the ERBB3 mutation was present as a major clone in one grade 2 region and as a minor clone in all other regions.
In Patient 2, we assessed two grade 2 regions, one grade 3 region, and four grade 4 regions from the available FFPE blocks. We identified rhabdoid and sarcomatoid dedifferentiation in the grade 4 regions. Mutations in VHL, BAP1, and ROS1 were shared by all regions, and were major clonal in the majority of them (Fig. 2B, Supplementary Table S3). The VHL and ROS1 mutations were each identified as a minor clone in a single region. Mutations in TSC1 and KDM5C were present in some of the regions, but were never unique for one grade.
For Patient 3, all three FFPE blocks were of tumour grade 3. Three cancer driver genes, ARID1A, VHL, and PBRM1, were mutated in all tumour regions, with ARID1A and VHL being major clonal. The PBRM1 mutation was major clonal in two samples and appeared to be a minor clone in the third sample. Six additional genes, including the cancer driver genes MTOR, KMT2C, and ZFHX3, were mutated in only a single region (Fig. 2C, Supplementary Table S3).
From Patient 4, two grade 3 regions and three grade 4 regions were analysed. Surprisingly no single mutation was shared between all regions, even though the sequence depth appeared sufficient in all instances where we did not detect a mutation. Interestingly, a VHL mutation was shared between only one of the two grade 3 regions and all grade 4 regions (Fig. 2D, Supplementary Table S3). A PBRM1 mutation was shared by one grade 3 and two grade 4 regions.
For Patient 5, we analysed two grade 2 regions and two grade 3 regions. Mutations were detected for two ccRCC driver genes from our panel: VHL and KDM5C. The mutations were shared by all samples as a major or minor clone (Fig. 2E, Supplementary Table S3).
For Patient 6, all available FFPE blocks were of grade 2, and we analysed three tumour regions. A VHL mutation was detected in all samples. A PBRM1 mutation was detected in two samples but was inconclusive in the third due to low coverage. Three additional genes were mutated in only a subset of the samples (Fig. 2F, Supplementary Table S3). The VHL gene was mutated in all the patients' tumours as a major clone and was thus clearly a trunk mutation in all cases. For patients 3 and 6, PBRM1 appeared to be a second trunk mutation. BAP1, ARID1A, and KDM5C were second trunk mutations in one case each.

Mutational profiles identified by whole exome sequencing
To look beyond the 42 panel genes, we carried out WES on tumour material of one additional patient. For this patient, we analysed four tumour regions and matched normal kidney cortex of this patient as control. We analysed one grade 1 region, one grade 2 region, and two grade 4 regions. Sequencing resulted in a mean target coverage of 41× in which 93% of the target region had > 10× coverage (Supplementary Table S2).
We identified 12 nonsynonymous somatic mutations that were shared by all regions, including two as a major mutation. In addition to that, 27 nonsynonymous somatic mutations located in 26 genes were variously distributed among tumour regions and responsible for intratumour heterogeneity (Fig. 3 and Supplementary Table S4). For the four regions, we observed 14 (grade 1), 11 (grade 2), 15 (grade 4) and 19 (grade 4) major clonal mutations, respectively. Thus the mutational load of the grade 4 samples is slightly higher compared to the low-grade samples. This also holds true when the minor clonal mutations are taken into consideration; 19 (grade 1), 23 (grade 2), 23 (grade 4) and 28 (grade 4). Eight of the 38 mutated genes are known cancer driver genes, of which five are included in our custom panel for targeted sequencing. Three of these cancer driver genes (VHL, PBRM1, and ARID1A) are related to ccRCC development. Whereas all regions harboured a VHL mutation, a PBRM1 mutation was only observed in one grade 2 region. The two grade 4 tumour regions shared the highest Fig. 1. Sample selection and DNA isolation. Multiple FFPE blocks containing tumour and normal tissue were collected from each patient (A). Each FFPE block was serially cut in 10 μm sections (B). The first and the last section (3 μm in thickness) were stained with H&E and used as reference in identifying tumour nests with different WHO/ISUP tumour grade (C). DNA was isolated from each tumour region with different WHO/ISUP grade (G1-4) in multiple block sections. P. Ferronika, et al. Experimental and Molecular Pathology 115 (2020) 104431 number of mutations including major clonal mutations in MYO7B, ARID1A, and FMN1. The two grade 4 regions also uniquely share mutations in the cancer driver genes TP53, and PIK3CA. Intratumour heterogeneity was also observed for the pattern of copy number variations (CNVs) of patient 7 (Fig. 4). All samples shared loss of the short arm of chromosome 3 (3p), which is characteristic for ccRCC. The grade 1 and grade 2 tumour samples, P7T3 and P7T1, are near diploid and differ from each other at five genomic segments. P7T3 has gain of 9q, and appears to be mosaic for del18. For P7T1, we observed loss of 9p, 10q, and 16q. A small fraction of the cells of P7T1 appear to have a gain of 9q and a loss of chr20. The two grade 4 tumour samples, P7T2 and P7T4, had virtually identical and complex copy number profiles that showed much more CNVs than the lower grade tumour samples. Chromosomal alterations which clearly distinguished these two grade-4 tumour samples from the lower grade include the amplification of 3q, chr5, and chr8, and loss of chromosomes 4, 13, and The samples are arranged in order from left to right based on tumour grade (G1-G4). Well-known cancer driver genes in ccRCC and other cancer(s) are highlighted in colors. The classification of the mutations is indicated by colors as shown in the legend at the bottom of the Figure. Mutations are included only if they represent a major mutation in at least one of the subregions of a patient's tumour, and if the read depth at their location was sufficient in at least all but one samples of the patient. Read counts and minor allele frequencies for the mutations are listed in Supplementary Table S3. Ferronika, et al. Experimental and Molecular Pathology 115 (2020) 104431 14. The copy number profile in these two samples suggests an aneuploidic nature and includes a characteristic chromotripsis-like structure for chromosome 5.

Tumour-grade specificity of mutated genes
To obtain more insight into possible grade-related mutation profiles, we pooled the observed mutated genes of all six patients per tumour grade (Fig. 5A). We included the mutations observed for patient 7, but limited ourselves to the genes present in our targeted sequencing panel.
We also restricted ourselves to the major clones of each tumour sample.
To represent the overlap between grades and genes uniquely mutated per grade, we constructed a Venn diagram (Fig. 5B). VHL mutations were present in all tumour grades. The mutated genes that were present  Supplementary Table 4. Ferronika, et al. Experimental and Molecular Pathology 115 (2020) 104431 in multiple tumour grades, except tumour grade 1, were BAP1, ROS1, KDM5C, and PBRM1 (Fig. 5B). PBRM1 was present in multiple tumour grade regions from four patients, whereas BAP1 and ROS1 were present in multiple tumour regions of only one patient (Fig. 2). A number of genes were mutated in regions with only one specific tumour grade (Fig. 5B), namely ERBB3, LRRK and DROSHA in grade 2; six genes including MTOR in grade 3; and TP53 and PIK3CA in grade 4. It is important to recognise that mutations in these genes were not identified in all samples with a particular tumour grade, so, for example, not all samples with grade 2 showed mutations in ERBB3.

Discussion
For the six patients subjected to targeted sequencing, most mutations were indeed observed in the two most frequently mutated genes in ccRCC, e.g. VHL and PBRM1 (Forbes et al., 2017;Network, 2013). Six . The X-axis of the plots represents the genomic position starting from 1pter until Xqter, while the Y-axis indicates log2 intensity ratio between tumour and reference. Samples P7T2 and P7T4 which were classified as tumour grade 4, show a virtual identical, and highly aberrant CNV pattern, including a chromotripsis-like pattern for chromosome 5. P. Ferronika, et al. Experimental and Molecular Pathology 115 (2020)  If tumour phenotypes in ccRCC were consistently driven by mutations in specific genes, these genes would very likely have ended up high in the ranks of the mutation databases and thus be captured by our panel, which is based on published gene mutation frequencies and cancer driver-gene rankings. Indeed, we observed genomic  Ferronika, et al. Experimental and Molecular Pathology 115 (2020) 104431 heterogeneity testing those genes. Eighteen genes of our 42 panel genes, including five ccRCC driver genes, contributed to the intratumour heterogeneity in our group of seven patients. Perhaps surprisingly, given its known role as an early driver of ccRCC development, this included VHL. VHL was qualified as a minor clone in one grade 4 sample from patient 2 and one grade 3 sample from patient 5, and found not to be mutated in one grade 3 sample from patient 4. In the last case, functional loss of VHL cannot be ruled out because it could have been inactivated by promoter methylation. This mechanism is known to occur frequently for VHL in ccRCC (Sato et al., 2013), and would not have been detected by our sequencing protocol. Using a targeted panel of the 42 most frequently mutated genes may not be sufficient to determine the full spectrum of subclonal mutations. WES did reveal a number of unique mutations for each sample from patient 7. Moreover, we observed a small increase in the number of major clone mutations for the samples with a high tumour grade. A previous study also reported an increase of the mutational load with the tumour grade of ccRCC tumours (Turajlic et al., 2018), although they did not make an intratumour comparison. An increase in the mutational load became more apparent in our WES patient when looking at largescale genomic changes, i.e. CNVs (Fig. 4). We identified some CNVs which are present in grade 4 samples but not in the lower grade, including the loss of 4p, chr13, and chr14. Loss of 4p, chr13, and chr14 has also been observed previously in higher ccRCC tumour grades (Arai et al., 2008;Moore et al., 2012). In line with our findings, previous studies found an increase in genomic instability to be correlated with higher ccRCC tumour grades (Ball et al., 2017;Kluzek et al., 2017;Moore et al., 2012). In contrast to a previous report that observed mainly copy number gains in the higher tumour grades, we saw a balanced increase of both gains and losses (Thiesen, 2017).
We saw a marked increase in chromosomal instability from the lowgrade to high-grade samples in this patient, including a characteristic chromotripsis-like event for the grade 4 samples that carry a TP53 mutation. Chromotripsis events have been shown to have a strong association with TP53 mutations in medulloblastoma and acute myeloid leukemia (Rausch et al., 2012). Chromotripsis for chromosome 5 has previously been found in renal cancer cell line TK10 (Stephens et al., 2011). These observations are also in line with the intratumour heterogeneity that we observed in another ccRCC patient (Ferronika et al., 2019). If confirmed in a larger series, the extend of chromosomal instability might be a phenomenon that could be used to distinguish individual tumour grades.
Our study suggested 11 mutated genes as potential markers for tumour grade, but the value of these findings may be limited by the small size of our series. Moreover, these 11 genes were found to be wild type in several samples with the same tumour grade, so they cannot fully reflect the genomic changes that drive ccRCC evolution and heterogeneity at tumour-grade level. Apparently, variations in ccRCC tumour grade are not accomplished by single mutations in one of the 42 mostfrequently-mutated genes in ccRCC, or even in groups of genes representing specific pathways. Instead, the underlying mechanism of tumour grade heterogeneity may be an expressional shift in pathways caused by an interplay between rare mutations and epigenetic changes. In addition to DNA sequencing, gene expression (and gene methylation) data might help in understanding the biology behind the tumour grades. The ultimate goal would be to help improve ccRCC patient prognostics, and this requires larger series of patients and, ideally, RNA as well as DNA analysis of their tumours.
Our study confirms the phenomenon of intratumour heterogeneity of ccRCC. Through multi-region sampling we were able to detect a higher number of mutations, as a subset of the mutations were present in local tumour regions only. The identification of these local mutations, including CNVs, may be crucial in the future when these alterations become predictive for prognosis as markers for metastatic potential and indicators for targeted therapy. Using data from multiregional samples would in fact help to distinguish trunk (driver) mutations from local mutations, which could help in selecting the optimal targeted therapy. Thus, either way, multiregion sampling may guide selection of the treatment that can most effectively target all, or the most threatening, tumour populations.