Ethics

This study was approved by the Institutional Review Board of Yokohama City University (A200100004) and RIKEN (H20-11), and each patient has provided written informed consent for publication.

RNA sequencing (RNA-seq)

In this study, RNA-seq of twenty-one sporadic ChRCCs and seven oncocytomas was performed (Supplementary Dataset 1). Bulk RNA-seq of sixteen BHD-associated renal tumours and five adjacent normal kidneys from ten BHD patients has been carried out in our previous study.17 RNAs were extracted from frozen tissues using the RNA extraction reagent ISOGEN (Nippon-gene) following the manufacturer's protocol, qualified and quantified by Bioanalyzer (Agilent Technologies). Library preparation was done using TruSeq Stranded mRNA Library Prep Kit and TruSeq RNA CD Index Kit according to the manufacturer protocol (Illumina). Sequencing was performed with HiSeq 2500 system with 126-bp reads (Illumina). Using STAR (v2.7.9a), the FASTQ files were mapped to a reference sequence generated from the GRCh37 fasta and GENCODE Release 37 comprehensive gene annotation for reference chromosomes (gencode.v19.annotation.gtf). Aligned SAM files were sorted and duplicates were marked with picard MarkDuplicates. Mapped reads were counted by featureCounts. Based on the read count data, gene expression levels were calculated as fragments per kilobase of exon per million mapped reads (fpkm). For the series of analyses from FASTQ file data to gene expression level normalization, the riboduct pipeline was used (https://github.com/msfuji/riboduct), and the details of the commands and parameter settings are available in the GitHub page. Then, these read count data and normalized data were used for further expression analysis. Gene Set Enrichment Analysis (GSEA) was done using GSEA software (v.4.2.3) (https://www.gsea-msigdb.org/gsea/index.jsp); default parameters were applied.19 Plots of principal component analysis (PCA) and L1CAM and FOXI1 expressions were created using iDEP (integrated Differential Expression and Pathway analysis) (v.0.96) (http://bioinformatics.sdstate.edu/idep96/), a tool for RNA-seq analysis on a web browser.20 In the iDEP analyses, raw count data created from featureCounts was used, and default settings were used for each parameter.

Analysis of single-cell transcriptomes

Using the single-cell RNA-seq data from our previous work, we performed gene expression analysis in BHD-associated renal tumour cells.17 The R package Seurat (v3.1.2) was used for the analysis.21 Parameter settings for clustering were resolution 2.0 for the FindClusters function, dims 1:4 for the RunUMAP function, and the package default values for the other parameters. BHD-associated HOCT cells were clustered into L1CAM expressing cell cluster (cL1CAM) or FOXI1 expressing cell cluster (cFOXI1), and gene expressions in each cell cluster were investigated.

Whole genome sequencing (WGS)

WGS was performed on sixteen BHD-associated renal tumours and nine adjacent normal kidneys from nine BHD patients (Table 1 and Supplementary Dataset 1). Genomic DNA was extracted from frozen tissues using the QIAGEN DNA Mini Kit (QIAGEN) following the manufacturer's protocol, qualified and quantified by Bioanalyzer (Agilent Technologies). The library construction was done using the TruSeq Nano DNA Library Prep Kit (Illumina) following the manufacturer's protocol. Paired-end sequencing of 150-bp reads was performed using NovaSeq6000. Reads mapping and further analyses were done using Genomon2 pipeline (v.2.6.3) (https://github.com/Genomon/genomon). Briefly, sequence reads were mapped to the reference human genome GRCh37 using BWA-mem. GATK (The Genome Analysis Toolkit) was used to realign bam files. For the Genomon2 pipeline analysis, the default settings were used for each parameter, and a control panel of normal samples was used to eliminate germline mutations and errors. We obtained the BAM file data of 49 samples of sporadic ChRCC from the data resource of The Pan-Cancer Analysis of Whole Genomes (PCAWG) study (KICH).12 In tumour mutation burden (TMB) analysis and mutational signature analysis, for the comparison of BHD-associated renal tumours and KICH, the output data (VCF files) using the same method described above were used. In TMB analysis, SP123975 was excluded because of the exceptionally higher number of mutations (54,708 mutations) due to DNA mismatch repair deficiency.

Germline and somatic variant calling

Single-nucleotide variants (SNV), multi-nucleotide variants (MNV) and small insertion-deletions (Indels) were called using SAGE v2.8 (Somatic Alterations in Genome) from the Hartwig Medical Foundation (HMF: https://github.com/hartwigmedical/hmftools/releases/tag/sage-v2.8). Germline calling parameters were based on the "germline mode" in the SAGE documentation. The panel-of-normals (PON) VCFs used in the pipeline for annotation and filtering of somatic variants called by SAGE were created from 149 Japanese germline samples.22 Genome-wide somatic variant calls were performed with a quality cutoff of 70 and then filtered using bcftools with a quality cutoff of 70 for hotspots and 100 for other genomic locations. The maximum acceptable variant allele frequency (VAF) for the reference sample was set at 4% and the minimum VAF for the tumour was set at 2.5%. The output VCFs were filtered for PON variant sites, sites that were variant or likely fixed (MAF >0.01 and AF < 0.99; AF ≥ 0.99) in dbSNP or ALFA population samples.22 The resulting VCFs were annotated using ANNOVAR and SnpEff.

Analyses of gene and structural variants (SV), and copy number alterations (CNA)

Structural Variation (SV) was called using GRIDSS 2.12.0 (the Genomic Rearrangement IDentification Software Suite)23 with the BAM file for each sample as the input. Scripts were used from GRIDSS- PURPLE-LINX (GPL) pipeline (v.1.3.2) (https://github.com/hartwigmedical/gridss-purple-linx). Analysis of copy number alteration (CNA) was done using CNApp (https://tools.idibaps.org/CNApp/), a tool for the quantification of copy number alterations.24 TUMOUR.purple.cnv.somatic.tsv, output files from the PURPLE, were used as input files. In order to minimize the copy number noise, the optimal solution for ploidy and purity in the PURPLE was selected and the noise was removed. After this process, F43-T1, F59-T11 and F59-T12 were excluded because of the high copy number noise. In the CNApp analysis, the default settings were used for each parameter.

Mutational signature analysis

Analysis of somatic mutation signatures was performed using SigProfiler (v.3.3) (https://cancer.sanger.ac.uk/signatures/tools/) with the VCFs created from the Genomon pipeline. Somatic mutation signature profiles were decomposed to COSMIC 3.2 signatures.25 The default settings were used for each parameter in all of analyses.

Putative cell-of-origin (COO) estimation

COOBoostR is a xgboost-based machine learning algorithm which aims to predict putative tissue-of-origin (TOO) or cell-of-origin (COO) (https://github.com/SWJ9385/COOBoostR).26 In our study, COOBoostR was applied to predict putative COO of BHD-associated renal tumours and ChRCCs. Briefly, 1-megabase window-based somatic point mutation frequencies from WGS data were calculated for individual tumour samples along with processing single-cell ATAC-seq data from normal kidney tissue27 by calculating 1-megabase window sequencing count density per each cell type (prespecified inside the metadata). For the parameters, default values were used except for the eta and depth (set to 0.3 and 10).

Putative variant acquisition timing analysis

R package MutationTimeR (v.1.0.2) (https://github.com/gerstung-lab/MutationTimeR) was used to examine the timing of a second hit occurrence of the remaining FLCN allele in BHD-associated renal tumours.28 MutationTimeR is an R package used by the PCAWG consortium to calculate mutation times for 2778 whole genome sequence samples.28 VCFs of BHD-associated renal tumour were used and each parameter was set as default setting. A small “clusters” data.frame, which includes the estimated number of subclonal clusters, SNVs in each cluster and clonal proportions of each cluster, created from Clip pipeline (v.1.2.1) (CliP: Clonal structure identification through penalizing pairwise differences) (https://github.com/wwylab/CliP) (BioRxiv: https://doi.org/10.1101/2021.03.31.437383) was used as an input for MutationTimeR.

Mitochondrial DNA analysis

Mitochondrial genome (mtDNA) copy number per cell in tumour and normal tissues was calculated as follows. Germline SAGE variants were extracted, the mean depth of coverage for variant sites on autosomes (mean.DPauto) or on chrM (mean.DPchrM) was calculated, and the mtDNA copy number in each tissue was estimated using the formula.29

This formula scales as:

for purity (f: 0–1 for tumours, 1 for normal tissues) and ploidy (2 in normal tissues).

Also, variants in the mtDNA of BHD-associated renal tumours were analyzed using VarScan and Mutect2, from which variant unions were extracted and further filtered as described below. Filtering criteria were as follows. Position was not in (302–315, 513–525, 568–573, 956–965, 3105–3109, 5895–5899, 8270–8289, 16180–16195 (indel)). Variant reads were≥10, variant base quality score was≥20. Forward read/Total variant read ratio: >0.1 and <0.9. Mean read position: >0.15 and <0.85 (from 5′ and 3’ both). Reference read length - Variant read length was <25. Reference read mapping quality - Variant read mapping quality was <10. Reads with normal context were >80% of variant reads (except variant next to germline mutation). Variant NM - Ref NM was≤3. Remove Germline Polymorphism (germline list of SAMPLES and PCAWG (VAF >95%), dbSNP (Frequency >1%)). Calculate Noise Level in normal sample in SAMPLE and PCAWG - > remove VAF lower than noise level. Variants that passed the filter were labeled as somatic or tumour-specific dysplastic variants based on VAF differences observed in normal tissue. We analyzed mtDNA variants of sporadic ChRCCs in TCGA cohort in the same method as we stated above .12

Immunohistochemistry

Immunohistochemistry was done on 4-mm-thick formalin-fixed paraffin-embedded sections of five each of L1CAM expressing sporadic ChRCCs (gL1CAM) and FOXI expressing sporadic ChRCCs (gFOXI1). In immunohistochemistry of FOXI1 protein, the tissue slides were autoclaved in Tris–EDTA buffer (pH 6.0) for antigen retrieval, then incubated with goat polyclonal FOXI1 antibody at 1:400 dilution (#ab20454; Abcam, Cambridge, United Kingdom) for 1 h at room temperature and visualized with N-Histofine Simple Stain MAX PO(G) (NICHIREI BIOSCIENCE INC., Tokyo, Japan). The FOXI1 protein was considered to be positive when nucleic and/or cytoplasmic FOXI1 staining was observed. In immunohistochemistry of L1CAM protein, the tissue slides were autoclaved in Tris–EDTA buffer (pH 9.0) for antigen retrieval, then incubated with mouse monoclonal L1CAM antibody at 1:200 dilution (clone 2C2, #ab24345; Abcam) for 1 h at room temperature and visualized with the Real EnVision Detection System, Peroxidase/DAB+, Mouse (Agilent Technologies). The L1CAM protein was considered to be positive when membranous L1CAM staining was observed.

Statistics

Welch's two sample t-test was applied to determine whether the means of two populations were different, and differences were considered to be statistically significant at a value of P < 0.05. All statistical tests were 2-sided. For the analysis of the correlation between the two gene expressions, Pearson's correlation coefficient and p-value were calculated.

BHD-associated renal tumours most commonly harbour biallelic FLCN alterations but fewer mutations and copy number alterations

We conducted WGS analysis of sixteen BHD-associated renal tumours and nine adjacent normal kidneys from nine BHD patients (Table 1 and Supplementary Dataset 1). WGS analysis in this study, together with our previous WES analysis, confirmed somatic second hit FLCN alterations in all of the BHD-associated renal tumours (Fig. 3a) (Supplementary Dataset 4).16 In addition, WGS analysis of BHD-associated renal tumours demonstrated that very few genes were commonly mutated except for FLCN, as was observed in our previous WES study, suggesting that FLCN alteration should be the main driver for BHD-associated renal tumourigenesis.16 In F107 case, interestingly, amplification of VEGFA and CCND3 on chromosome 6 was observed, which may lead to the formation of clear cell renal cell carcinoma. F133-T1, whose histology was unclassified renal cell carcinoma, harboured a variant in succinate dehydrogenase D (SDHD), a gene responsible for SDH-associated tumour predisposition syndrome.32^,33

Open in a separate window

Fig. 3

BHD-associated renal tumours harbour less numbers of gene variants compared to sporadic ChRCCs. (a) Mutation matrices show histology, clinical characteristics (age, size and stage) and somatic variants revealed by WGS of sixteen BHD-associated renal tumours. (b) Bar plots (left panel) and box plots (right panel) show the number of mutations in sixteen BHD-associated renal tumours and forty-eight sporadic ChRCCs in TCGA cohort (KICH). Welch's two sample t-test was used.

The average number of mutations per tumour was 1107, 1090 and 1650 in BHD-associated ChRCCs, BHD-associated HOCTs and sporadic ChRCCs in the TCGA cohort (KICH), respectively, supporting that BHD-associated renal tumours have less inter-tumour heterogeneity than sporadic ChRCCs (Fig. 3b, and Supplementary Dataset 5). In fact, BHD-associated renal tumours did not show increased number of mutations in an age-dependent manner, whereas sporadic ChRCCs in the TCGA cohort did exhibit an age-dependent increase in number of mutations, further supporting the notion that BHD-associated renal tumours have less inter-tumour heterogeneity compared to ChRCCs (Supplementary Fig. S2a).

Of note, we observed no common structural variants in BHD-associated renal tumours, suggesting that structural variants may not be key drivers for BHD-associated renal tumourigenesis (Supplementary Fig. S2b). High-resolution copy number alteration (CNA) analysis using WGS data exhibited that CNAs of BHD-associated renal tumours were considerably different from those already reported in sporadic cases; fewer CNAs were observed in BHD-associated HOCT and BHD-associated ChRCC compared to those in BHD-associated clear cell renal cell carcinoma, which might reflect the indolent nature of BHD-associated HOCT and BHD-associated ChRCC. Loss of chromosomes 1, 2, 6, 10, 13 and 17, which is characteristic of sporadic ChRCC, was not observed in BHD-associated renal tumours (Supplementary Fig. S2c).1^,12^,16^,31

BHD-associated renal tumours demonstrate unique mutational signatures

Mutational signatures are imprints of DNA damage and repair processes that have been operative during tumourigenesis and can provide insights into environmental and endogenous causes of cancer and biological implications of somatic mutations.25^,34 In this study, mutational signature analysis was performed using WGS data from 16 BHD-associated renal tumours and 49 sporadic ChRCCs in TCGA cohort (Supplementary Fig. S3a). BHD-associated renal tumours and sporadic ChRCCs in the TCGA cohort generally demonstrated similar patterns of mutational signatures. However, some of the signatures were different between these tumours, including SBS1 (Spontaneous deamination of 5-methylcytosine, clock-like signature), SBS2 (APOBEC activity), SBS5 (Unknown, clock-like signature), SBS13 (APOBEC activity), insertions and deletions (ID)2 (Slippage during DNA replication of the template DNA strand) and ID8 (DBS repair by non-homologous end-joining), highlighting molecular differences of these tumours.

BHD-associated renal tumours harbour mtDNA with higher copy numbers and fewer variants

A notable proportion (22.2%) of sporadic ChRCCs harbor non-silent mtDNA variants without any variants in nuclear-encoded kidney cancer-associated genes, whereas BHD-associated renal tumours harbor no disruptive alterations in mtDNA, suggesting that non-silent mtDNA variants in sporadic ChRCCs may be drivers of tumourigenesis of sporadic ChRCCs whereas mitochondrial function is most likely retained in BHD-associated renal tumours.12^,29^,35 Consistent with these results, heteroplasmic level or variant allele frequencies (VAFs) of total and pathogenic variants in mtDNA of BHD-associated renal tumours were significantly lower than those of sporadic ChRCCs (p = 0.0368 and 0.0117, respectively) (Fig. 4a) (Supplementary Dataset 6). mtDNA copy numbers were higher in BHD-associated renal tumours compared to sporadic ChRCCs and variable depending on each histology (Fig. 4b and c) (Supplementary Fig. S4).

Open in a separate window

Fig. 4

BHD-associated renal tumours harbour higher mtDNA copy number with fewer variants compared to sporadic ChRCCs. (a) Bar plots of number of mutations (left upper panel) and scatter plots of variant allele frequencies (VAF) (left lower panel) in mtDNA from sixteen BHD-associated renal tumours and forty-nine sporadic ChRCCs in TCGA cohort (KICH) are shown. Bar graphs show average number of variants shown in left upper panel (right upper panel). Box plots show average VAF shown in left lower panel (right lower panel). Welch's two sample t-test was used. (b) Box plot shows the mtDNA copy number in BHD-associated renal tumours and various sporadic renal tumours. BHD-Clear cell: BHD-associated clear cell renal cell carcinoma. BHD-Unclassified: BHD-associated unclassified renal cell carcinoma. KICH: sporadic ChRCCs in TCGA cohort. KIRC: sporadic clear cell renal cell carcinomas in TCGA cohort. KIRP: sporadic papillary renal cell carcinomas in TCGA cohort. Welch's two sample t-test was used. (c) Box plots show mitochondrial DNA copy number in each histological type of BHD-associated renal tumours and their adjacent normal kidneys (Paired normal). Welch's two sample t-test was used.

BHD-associated renal tumours and sporadic ChRCC may arise from different normal nephron cells

A recent machine learning algorithm integrating bulk or single cell-level ATAC-seq (Assay for Transposase-Accessible Chromatin using sequencing) data of normal cell types with somatic point mutation frequency derived from the WGS data of tumours enables us to predict putative cancer cell-of-origin (COO), based on a hypothesis that the somatic variants might preferentially accumulate outside of the open chromatin regions of the cancer cell of origin.36 In this study, we applied the COOBoostR algorithm to predict a putative cancer COO of BHD-associated renal tumours and sporadic ChRCCs.26 Interestingly, the COOBoostR predicted that BHD-associated renal tumours and sporadic ChRCCs may arise from a variety of normal nephron cells. The most frequently predicted cell-of-origin was proximal tubules for BHD-associated renal tumours and intercalated cells for sporadic ChRCCs, implying the possibility that these two tumour types might arise from different cells of origin (Fig. 5a) (Supplementary Dataset 7). Of note, none of L1CAM expressing tumours in BHD-associated renal tumours and sporadic ChRCCs might arise from intercalated cells, suggesting that cell-of-origin may determine tumour characteristics of these tumours.

Open in a separate window

Fig. 5

Bioinformatic algorithms depict natural history of BHD-associated renal tumourigenesis. (a) The bar plots show putative cell of origin of BHD-associated renal tumours (left panel) and sporadic ChRCCs in TCGA cohort (KICH) (right panel) predicted by COOBoost algorithm. Groups of gFOXI1, gNeutral and gL1CAM correspond to those in Fig. 2d and e. (b) The plots show the results of MutationTimeR algorithm. The dot plots show variant allele frequencies (VAF) of point mutations; purple dot indicates clonal mutation occurred at late stage, green dot indicates clonal mutation occurred at early stage, blue dot indicates clonal mutation occurred at unknown stage, red dot indicates subclonal mutation (left upper panels in each case). The stacked barplots show copy number of each allele; dark and light grey bars indicate major and minor allele, respectively (left middle panels in each case). The box plots with median values and 95% confidence intervals indicate predicted mutation timings of primary and secondary gains; blue bar indicates mono-allelic gains (N:1), pink bar indicates copy neutral (CN)-LOH/loss + gain (N:0) and green bar indicates bi-allelic gains (N:2) (left lower panels in each case). The histograms show magnified timelines of left lower panels in each case (right panels in each case).

We thank Dr. Noboru Nakaigawa, Dr. Keiichi Kondo, Dr. Narihiko Hayashi, Dr. Kimito Osaka and Dr. Yumiko Yokomizo in Department of Urology, Yokohama City University for surgical resection of renal tumours and Dr. Akira Nishiyama and Dr. Koichi Murakami in Department of Immunology, Yokohama City University for their excellent works in scRNA-seq analysis.

Authors were supported by JSPS KAKENHI Grant Number as following; H.Hasumi. by 19K09694, I.K. by 19K16563, Y.Ito. by 19K18591, H.I. by 20K18121 and 19K23781, M.K. by 21H03068, M.Furuya. by 20K07395, M.Y. by 19K09717, S.F. by 19K07769, H.N. by 18H04049 and M.B. by 18H02938. This work was also supported by the Intramural Research Program of the National Institutes of Health (NIH), National Cancer Institute (NCI), Center for Cancer Research. This project has been funded in whole or in part with Federal funds from the National Cancer Institute, National Institutes of Health, under Contract No. HHSN261201500003I. The content of this publication does not necessarily reflect the views or policies of the Department of Health and Human Services, nor does mention of trade names, commercial products or organizations imply endorsement by the US Government.