Multilayer-omics analysis of renal cell carcinoma, including the whole exome, methylome and transcriptome.

The aim of this study was to identify pathways that have a significant impact during renal carcinogenesis. Sixty-seven paired samples of both noncancerous renal cortex tissue and cancerous tissue from patients with clear cell renal cell carcinomas (RCCs) were subjected to whole-exome, methylome and transcriptome analyses using Agilent SureSelect All Exon capture followed by sequencing on an Illumina HiSeq 2000 platform, Illumina Infinium HumanMethylation27 BeadArray and Agilent SurePrint Human Gene Expression microarray, respectively. Sanger sequencing and quantitative reverse transcription-PCR were performed for technical verification. MetaCore software was used for pathway analysis. Somatic nonsynonymous single-nucleotide mutations, insertions/deletions and intragenic breaks of 2,153, 359 and 8 genes were detected, respectively. Mutations of GCN1L1, MED12 and CCNC, which are members of CDK8 mediator complex directly regulating β-catenin-driven transcription, were identified in 16% of the RCCs. Mutations of MACF1, which functions in the Wnt/β-catenin signaling pathway, were identified in 4% of the RCCs. A combination of methylome and transcriptome analyses further highlighted the significant role of the Wnt/β-catenin signaling pathway in renal carcinogenesis. Genetic aberrations and reduced expression of ERC2 and ABCA13 were frequent in RCCs, and MTOR mutations were identified as one of the major disrupters of cell signaling during renal carcinogenesis. Our results confirm that multilayer-omics analysis can be a powerful tool for revealing pathways that play a significant role in carcinogenesis.

What’s new?

Large-scale systems biology approaches are currently reshaping biomedical research identifying new pathways or reinforcing significance of previously discovered pathways in cancer biology. Here the authors performed multilayer -omics analyses in clear renal carcinoma or healthy control samples. They found frequent tumor-associated genetic aberrations of GCN1L1, MED12, and CCNC, all members of the CDK8 Mediator complex involved in regulating β-catenin-driven transcription, as well as alterations in MACF1, also a member of the Wnt/β-catenin signaling pathway. These findings underscore the significance of the Wnt/β-catenin signaling pathway during renal carcinogenesis and confirm the power of large-scale sequencing efforts in revealing pathways that may become therapeutic targets in specific cancers.

Clear cell renal cell carcinoma (RCC) is the most common histological subtype of adult kidney cancer and frequently affects working-age adults in midlife.1 Recently, large-scale PCR-based exon resequencing and whole-exome analysis by exon capturing have revealed that renal carcinogenesis involves inactivation of histone-modifying genes such as SETD2,2 a histone H3 lysine 36 methyltransferase, KDM5C,2 a histone H3 lysine 4 demethylase and UTX,3 a histone H3 lysine 27 demethylase, as well as the SWI/SNF chromatin remodeling complex gene, PBRM1.4 Moreover, it is well known that clear cell RCCs are characterized by inactivation of the VHL tumor-suppressor gene encoding a component of the protein complex that possesses ubiquitin ligase E3 activity.5 Another exome analysis study has revealed frequent mutation of a further component of the ubiquitin-mediated proteolysis pathway, BAP1.6 Non-synonymous mutations of the NF2 gene and truncating mutations of the MLL2 gene have also been reported.2

Not only genetic, but also epigenetic events appear to accumulate during carcinogenesis, and DNA methylation alterations are one of the most consistent epigenetic changes in human cancers.7,8 In fact, we have shown that noncancerous renal tissue obtained from patients with RCCs is already at the precancerous stage associated with DNA methylation alterations, even though no remarkable histological changes are evident and there is no association with chronic inflammation or persistent infection with viruses or other pathogens.9,10 Furthermore, using single-CpG resolution methylome analysis with the Infinium array, we have demonstrated that DNA methylation alterations at precancerous stages may determine tumor aggressiveness and patient outcome.11

It is well known that DNA methylation alterations around promoter regions affect the expression levels of tumor-related genes.7 Once the DNA methylation status has been altered, such alterations are stably preserved on the DNA double strands by covalent bonds through maintenance-methylation mechanisms by DNMT1 during carcinogenesis.7 Therefore, tumor-related genes showing alterations of both expression level and DNA methylation may have a larger impact on carcinogenesis than those showing only alterations of expression. Therefore, subjecting tissue specimens to a combination of both methylome and transcriptome analyses may be a powerful approach for revealing genes that are involved in carcinogenetic pathways.

Although one article reporting the use of an integrated multilayer-omics approach including exome analysis to examine human clear cell RCCs was published while this manuscript was in preparation,12 the entire pathway of carcinogenesis in the kidney may not yet be fully explained. In this study, to identify pathways having a significant impact during renal carcinogenesis, we subjected paired samples of both noncancerous renal cortex tissue (N) and cancerous tissue (T) from patients with clear cell RCCs to whole-exome, methylome and transcriptome analyses.

Material and Methods

Patients and tissue samples

Sixty-seven paired T and N samples were obtained from materials that had been surgically resected from 67 patients with primary clear cell RCCs. N mainly consists of proximal tubules, which are the origin of clear cell RCCs. These patients had not received any preoperative treatment and had undergone nephrectomy at the National Cancer Center Hospital, Tokyo. Tissue specimens were provided by the National Cancer Center Biobank, Tokyo. Histological diagnosis was made in accordance with the World Health Organization classification.13 All the tumors were graded on the basis of previously described criteria14 and classified according to the pathological Tumor-Node-Metastasis classification.15 The clinicopathological parameters of these RCCs are summarized in Supporting Information Table S1.

All patients included in this study provided written informed consent. This study was approved by the Ethics Committee of the National Cancer Center, Tokyo and was performed in accordance with the Declaration of Helsinki.

Exome analysis

High-molecular-weight DNA was extracted using phenol-chloroform, followed by dialysis. Three-microgram aliquots of genomic DNA from the 67 paired samples were fragmented by a Covaris-S2 instrument (Covaris, Woburn, MA) to provide DNA fragments with a base pair peak at 150–200 bp. The DNA fragments were end-repaired and ligated with paired-end adaptors (NEBNext DNA sample prep, New England Biolabs, Ipswich, MA). The resulting DNA library was purified using Agencourt AMPure XP Reagent (Beckman Coulter Genomics, Danvers, MA) and amplified by PCR (4 cycles). Five-hundred-nanogram aliquots of the adaptor-ligated libraries were hybridized for 24 hr at 65°C with biotinylated oligo RNA bait, SureSelect Human All Exon 50 Mb (Agilent Technologies, Santa Clara, CA). The hybridized genomic DNA was subjected to 10 cycles of PCR reamplification. Following the manufacturer’s standard protocols, the whole-exome DNA library was sequenced on an Illumina HiSeq 2000 (Illumina, San Diego, CA) using 75-bp paired-end reads.

After completion of the entire run, image analyses, error estimation and base calling were performed using the Illumina Pipeline (version 1.3.4) to generate primary data. First, the reads were aligned against the reference human genome from UCSC human genome 19 (Hg19) using the Burrows Wheeler Aligner Multi-Vision software package.16 Because duplicated reads had been generated during the PCR amplification process, paired-end reads that were aligned to the same genomic positions were removed using SAMtools. Second, the following loci were removed: (i) read depth <6 and (ii) base quality score <3 in the T sample. Third, we used the following Bayesian data analysis pipeline developed in our laboratory: (i) single nucleotide polymorphism (SNP) array analysis was performed on each paired cancerous and noncancerous tissue samples using Illumina HumanOmni1-Quad BeadChip (see “SNP microarray analysis”) and the genomic region, which is considered to be 1 copy in the pure cancerous genome was identified by the visual inspection of the log R ratio and B allele frequency plots on the Illumina Genome Viewer in the GenomeStudio software. (ii) Heterozygous SNP loci were selected from the above 1-copy region using GATK UnifiedGenotyper (Broad Institute, MA). (iii) At the SNP loci, which were 1 copy in the pure cancerous genome but heterozygous in the noncancerous genome, the ratio of the contaminating non-cancerous cells in the cancerous tissue was estimated from the allele frequencies of the cancerous genome by fitting to a binominal mixture model. (iv) Considering the estimated ratio of the contaminating noncancerous cells, the posterior probability of the genotypes of the cancer cells was calculated. Mutation was called if the posterior probability of being homozygous for the allele recorded in the reference human genome sequence was 0.001 or lower, and the ratio of the nonreference allele was 0.02 or lower in the noncancerous tissue sample, which had a read depth of at least 15. Fourth, Annovar extracted candidates that were nonsynonymous and did not correspond to the refSNP number. Fifth, candidates were discarded if the frequency of the nonreference allele was >2% in the N sample. Somatic mutations were also removed from the candidates if the root mean square mapping quality score of the reads covering the somatic mutation was <20. Finally, if the Blast search did not detect homologous regions for which the edit distance was 7 or <7 within the neighboring 151-bp stretch (75 bp both up- and downstream), the candidate was considered as a somatic mutation. Somatic insertions/deletions (indels) were called using both SAMtools and Pindel17 as described previously.18 Effects of amino acid substitutions on protein function due to single nucleotide nonsynonymous mutations have been estimated using the Sorting Intolerant from Tolerant (SIFT) (http://sift.jcvi.org)19 and polymorphism phenotyping (PolyPhen)-2 (http://genetics.bwh.harvard.edu/pph2/),20 and those due to indels have been estimated using SIFT.21 All data from exome analysis will be submitted to the Genome Medicine Database of Japan (GeMDBJ, https://gemdbj.nibio.go.jp/dgdb/).

Sanger sequencing

To verify the nonsynonymous single-nucleotide mutations and indels detected by the exome analysis and described in Table 1, the target sites and the flanking sequences of each patient’s DNA template were amplified individually with specific primers designed using Primer6.0. The PCR products were then sequenced with an ABI 3730 DNA Analyzer using the BigDye Terminator v1.1 Cycle Sequencing kit (Life Technologies, Carlsbad, CA).

Table 1

Genes showing 3 or more genetic aberration scores in clear cell RCCs

							Predicted protein function
			Genetic aberration score				Nonsynonymous single-nucleotide mutation2		Indel	Copy number aberration (%)3
Genes	Chr1	Entrez Gene ID	Non-synonymous single-nucleotide mutation	Indel	Intragenic break	Total	SIFT	PolyPhen-2	SIFT	Loss	Gain
VHL	3	7,428	22	14	0	36	0	1	Damaging	77.61	11.94
PBRM1	3	55,193	11	10	1	22	0	1	Damaging	73.13	10.45
TTN	2	7,273	9	3	0	12	0.75	0.387878	Neutral	0.00	38.81
KDM5C	X	8,242	4	4	0	8	0	0.998	Damaging	53.73	26.87
MUC16	19	94,025	6	0	0	6	0	NA	–	2.99	29.85
CUBN	10	8,029	5	1	0	6	0.32	0.987	Damaging	0.00	26.87
SETD2	3	29,072	3	3	0	6	0	0.99	Damaging	76.12	7.46
ABCA13	7	154,664	5	0	0	5	0	NA	–	0.00	44.78
BIRC6	2	57,448	4	1	0	5	0.02	NA	Damaging	4.48	35.82
GCN1L1	12	10,985	3	2	0	5	0	0.735079	Damaging	0.00	37.31
HERC2	15	8,924	5	0	0	5	0.01	0.902	–	1.49	25.37
BAP1	3	8,314	4	0	0	4	0	1	–	74.63	10.45
KIAA0100	17	9,703	4	0	0	4	0.05	0.999	–	0.00	29.85
MTOR	1	2,475	4	0	0	4	0	0.999	–	7.46	25.37
SPTBN1	2	6,711	3	1	0	4	0	0.993	NA	0.00	35.82
SPTA1	1	6,708	2	2	0	4	0.09	0.513	Damaging	0.00	34.33
CADM2	3	253,559	1	0	3	4	0.09	0.012	–	29.85	25.37
ERC2	3	26,059	1	0	3	4	0.01	NA	–	71.64	10.45
ADAM23	2	8,745	3	0	0	3	0	0.998	–	2.99	37.31
AKAP9	7	10,142	3	0	0	3	0	0.986	–	0.00	46.27
ANKRD26	10	22,852	3	0	0	3	0	0.995	–	2.99	28.36
ARHGEF33	2	100,271,715	3	0	0	3	0	NA	–	2.99	35.82
BRD4	19	23,476	3	0	0	3	0	0.997	–	0.00	29.85
C1orf112	1	55,732	3	0	0	3	0	0.952	–	0.00	34.33
CCNC	6	892	3	0	0	3	0	0.876	–	2.99	22.39
CPAMD8	19	27,151	3	0	0	3	0	0.439286	–	0.00	29.85
CSMD3	8	114,788	3	0	0	3	0	0.999	–	1.49	31.34
DNAH5	5	1,767	3	0	0	3	0.1	0.169	–	0.00	46.27
FAT1	4	2,195	3	0	0	3	0	NA	–	1.49	22.39
FAT2	5	2,196	3	0	0	3	0	0.999	–	0.00	71.64
FMN2	1	56,776	3	0	0	3	0	0.957	–	0.00	34.33
FNIP1	5	96,459	3	0	0	3	0.1	0.45171	–	0.00	65.67
KIF26B	1	55,083	3	0	0	3	0	NA	–	0.00	34.33
LIMCH1	4	22,998	3	0	0	3	0	0.992	–	2.99	20.90
LRBA	4	987	3	0	0	3	0.01	0.939	–	1.49	23.88
MACF1	1	23,499	3	0	0	3	0	0.791225	–	4.48	25.37
MADD	11	8,567	3	0	0	3	0	0.999	–	0.00	29.85
MED12	X	9,968	3	0	0	3	0.01	0.576	–	55.22	25.37
MGAM	7	8,972	3	0	0	3	0	NA	–	0.00	46.27
OBSCN	1	84,033	2	1	0	3	0	NA	Neutral	1.49	34.33
PLCE1	10	51,196	3	0	0	3	0	0.999	–	10.45	25.37
PREX2	8	80,243	3	0	0	3	0	1	–	5.97	31.34
PTPN4	2	5,775	3	0	0	3	0	0.999	–	0.00	35.82
ROR2	9	4,920	3	0	0	3	0	1	–	7.46	20.90
RP1	8	6,101	3	0	0	3	0.01	0.992	–	5.97	31.34
RYR2	1	6,262	3	0	0	3	0	NA	–	0.00	34.33
SYNE1	6	23,345	3	0	0	3	0.04	0.918	–	4.48	20.90
TTI1	20	9,675	3	0	0	3	0	0.999	–	0.00	29.85
VWDE	7	221,806	3	0	0	3	0.04	NA	–	0.00	44.78
ATM	11	472	2	1	0	3	0	1	NA	7.46	28.36
DNAH2	17	146,754	2	1	0	3	0.14	0.048	Damaging	0.00	29.85
FOXN2	2	3,344	2	1	0	3	0.08	0.255	Neutral	1.49	35.82
PTEN	10	5,728	2	1	0	3	0.01	0.988	Damaging	8.96	25.37
SAMD9L	7	219,285	2	1	0	3	0	0.968	Damaging	0.00	46.27
SI	3	6,476	2	1	0	3	0.01	0.992	Damaging	10.45	32.84
TCHH	1	7,062	2	1	0	3	NA	0.998	Damaging	0.00	34.33
TUBGCP6	22	610,053	2	1	0	3	0	0.993	NA	1.49	29.85
UGGT2	13	55,757	2	1	0	3	0.01	0.726	Neutral	0.00	25.37
CCDC178	18	374,864	1	2	0	3	0	0.235	Damaging	2.99	22.39
HGSNAT	8	138,050	1	2	0	3	0	NA	Damaging	16.42	25.37
NIPBL	5	25,836	1	2	0	3	0.05	0.98	Damaging	0.00	46.27

Chromosome.

Minimum SIFT score and maximum PolyPhen-2 score among all detected mutations of each gene (A SIFT score of <0.05 means “damaging.”19 PolyPhen-2 scores of >0.85 and 0.15–0.85 mean “probably damaging” and “possibly damaging,” respectively).20 NA: not available using SIFT or PolyPhen-2; –: indels of the gene were not detected.

The incidence of loss (1 or less copy number) or gain (3 or more copy number) detected using ASCAT or GPHMM in all 67 tumors. SIFT and PolyPhen-2 scores and copy numbers of each gene in each RCC were described in Supporting Information Table S3.

SNP microarray analysis

Two-hundred-nanogram aliquots of DNA from the 67 paired samples were genotyped with the HumanOmni1-Quad BeadChip (Illumina) in accordance with the manufacturer’s protocols. The data were assembled using GenomeStudio software (Illumina). For the single-nucleotide mutation detection, we developed the Bayesian data analysis pipeline using SNP microarray data (see “Exome analysis”). Localization of intragenic breakpoints, in which the end point of a deletion or duplication lies within a gene, in each of the T samples was clearly identified by the visual inspection of the B allele frequency plots on the Illumina Genome Viewer in the GenomeStudio software (Supporting Information Fig. S1). Copy number data has been obtained using Allele-Specific Copy Number Analysis of Tumors (ASCAT; http://heim.ifi.uio.no/bioinf/Projects/ASCAT/)22 and Global Parameter Hidden Markov Model (GPHMM; http://bioinformatics.ustc.edu.cn/gphmm/)23 software.

Infinium analysis

Five-hundred-nanogram aliquots of DNA from the 67 paired samples were subjected to bisulfite conversion using an EZ DNA Methylation-Gold™ Kit (Zymo Research, Irvine, CA). Subsequently the DNA methylation status at 27,578 CpG loci was examined at single-CpG resolution using the Infinium HumanMethylation27 Bead Array (Illumina). The data were assembled using GenomeStudio methylation software (Illumina). At each CpG site, the ratio of the fluorescent signal was measured using a methylated probe relative to the sum of the methylated and unmethylated probes, that is, the so-called β-value, which ranges from 0.00 to 1.00, reflecting the methylation level of an individual CpG site. All data of Infinium analysis will be submitted to GeMDBJ.

Pyrosequencing

DNA methylation levels of Infinium probe sites of the RAB25, GGT6, C3 and CHI3L2 genes and the 5′-region of the ABCA13 gene were measured by pyrosequencing. The PCR and sequencing primers were designed using Pyrosequencing Assay Design Software ver.1.0 (QIAGEN, Hilden, Germany). To overcome any PCR bias, we optimized the annealing temperature as described previously.24 Each of the primer sequences and PCR conditions are given in Supporting Information Figure S2. The PCR product was generated from bisulfite-treated DNA and subsequently captured on streptavidin-coated beads. Quantitative sequencing was performed on a PyroMark Q24 (QIAGEN) using the Pyro Gold Reagents (QIAGEN) in accordance with the manufacturer’s protocol.

Expression microarray analysis

Total RNA was isolated using TRIzol reagent (Life Technologies). From the 67 paired samples, 29 pairs, from which a sufficient amount of total RNA for both N and T samples was available, were subjected to expression microarray analysis. Two-hundred-nanogram aliquots of total RNA from the 29 paired samples were used for the production of fluorescent complementary RNA, and all samples were hybridized to the SurePrint G3 Human Gene Expression 8 × 60 K microarray (Agilent Technologies). The signal values were extracted using the Feature Extraction software (Agilent Technologies). All data of Expression microarray analysis will be submitted to GeMDBJ.

Quantitative RT-PCR analysis

cDNA was reverse-transcribed from total RNA using random primers and Superscript III RNase H− Reverse Transcriptase (Life Technologies). From the 67 paired samples, 66 pairs, from which a sufficient amount of cDNA for both N and T samples was available, were subjected to quantitative RT-PCR analysis. mRNA expression was analyzed using custom TaqMan Expression Assays (probe and PCR primer sets, Supporting Information Table S2) on the 7500 Fast Real-Time PCR System employing the relative standard curve method. All CT values were normalized to that of GAPDH in the same sample.

Multilayer-omics scoring

If any of the somatic nonsynonymous single-nucleotide mutations, indels or intragenic breaks was observed in one of the T samples, a genetic aberration score of one was assigned for the gene. If the Δβ (βΤ − βΝ) was 0.2 or more, the gene was considered to be hypermethylated in the T sample relative to the corresponding N sample. If the Δβ (βΤ − βΝ) was −0.2 or less, the gene was considered to be hypomethylated in the T sample relative to the corresponding N sample.

The expression level (E value) of each gene was expressed as the log2-signal intensity normalized by the median for all probes in the sample. If the ΔΕ (EΤ − EΝ) was 4 or more, the expression of the gene was considered to be elevated in the T sample relative to the corresponding N sample. If the ΔΕ (EΤ − EΝ) was −4 or less, the expression of the gene was considered to be reduced in the T sample relative to the corresponding N sample.

All probes of the Infinium HumanMethylation27 Bead Array and SurePrint G3 Human Gene Expression 8 × 60 K microarray were aligned against the reference human genome from Hg19. Infinium array probe and expression microarray probe pairs were annotated to each individual gene. If the probe of the Infinium array was designed for the upstream region including the promoter region, exon 1 or intron 1 of the gene, if Δβ (βΤ − βΝ) of the gene was 0.2 or more (DNA hypermethylation), and if ΔΕ (EΤ − EΝ) based on the expression microarray was −4 or less (reduced expression) in one paired sample of T and N, then a gene downregulation score of one was assigned. If the probe of the Infinium array was designed for the upstream region including the promoter region, exon 1 or intron 1 of the gene, if Δβ of the gene was −0.2 or less (DNA hypomethylation), and if ΔΕ (EΤ − EΝ) based on the expression microarray was 4 or more (overexpression) in one paired sample of T and N, then a gene upregulation score of one was assigned.

Pathway analysis

MetaCore software (http://www.genego.com) is a pathway analysis tool based on a proprietary manually curated database of human protein–protein, protein–DNA and protein compound interactions. The MetaCore pathway analysis by GeneGo was performed among genes showing genetic scores of 3 or more or showing downregulation or upregulation scores of 5 or more. Pathways for which the p value was <0.05 were considered to play a significant role in renal carcinogenesis.

Results

Genetic aberrations

Exome analysis detected somatic non-synonymous single-nucleotide mutations and indels of 2,153 and 359 genes among the 67 clear cell RCCs, respectively. SNP array analysis revealed intragenic breaks in 8 genes among the 67 RCC samples. In total, 2,440 genes showed non-synonymous single-nucleotide mutations, indels and/or intragenic breaks in RCCs and were assigned genetic aberration scores (described in “Multilayer-omics scoring” in the Material and Methods section) of 1 or more. Genetic alterations in each RCC are summarized in Supporting Information Table S3. The 2,131 and 248 genes that were assigned a genetic aberration score of 1 and 2 are listed in Supporting Information Table S4, and the 61 genes that were assigned genetic aberration scores of 3 or more are listed in Table 1. All 256 mutations (209 somatic nonsynonymous single-nucleotide mutations and 57 indels with 10 exceptions, for which Sanger sequencing failed due to difficulties with PCR primer design) listed in Table 1 were verified by Sanger sequencing. In addition, mutations of 54 (89%) of the 61 genes included in Table 1 were also found in the clear cell RCC database in The Cancer Genome Atlas (http://cancergenome.nih.gov/; Supporting Information Table S5), indicating the reliability of our whole-exome analysis results.

Effects of amino acid substitutions due to genetic aberrations on protein function estimated using SIFT19,21 and PolyPhen-220 software are shown in Table 1 and Supporting Information Table S3. In 60 of 61 genes listed in Table 1, SIFT and PolyPhen-2 analyses (less than 0.05 SIFT score19 or more than 0.15 PolyPhen-2 score20 for nonsynonymous single-nucleotide mutations and “damaging” SIFT score21 for indels) indicated that amino acid substitutions due to genetic aberrations impair the functions of proteins.

The incidence of copy number loss (1 or less) and gain (3 or more), detected using ASCAT22 and GPHMM23 software, of the genes that were assigned genetic aberration scores of 3 or more is described in Table 1. The copy number of each gene showing genetic aberrations in each RCC is described in Supporting Information Table S3. Nonsynonymous single-nucleotide mutations and indels were frequently concordant with copy number alterations (Table 1), suggesting that such genetic aberrations may actually result in dysfunction of proteins in RCCs.

In addition to recurrent genetic aberrations, expression microarray analysis revealed reduced mRNA expression [ΔΕ (EΤ − EΝ) was -4 or less as described in “Expression microarray analysis” in the Material and Methods section] of the ERC2 and ABCA13 genes in 21 and 31% of RCCs, respectively. These mRNA expression alterations were verified quantitatively by real-time RT-PCR analysis [mean ERC2 expression levels in T samples (n = 66): 8.91 ± 29.72; those in N samples (n = 66): 110.02 ± 75.31 (p < 1.00 × 10−12, Mann-Whitney U-test) and mean ABCA13 expression levels in T samples (n = 66): 8.43 ± 45.12; those in N samples (n = 66): 47.82 ± 89.51 (p < 1.00 × 10−12, Mann-Whitney U-test)].

Probes for the ERC2 gene were designed for the Infinium array, and DNA hypermethylation around the 5′-region of the ERC2 gene was detected in only 6% of RCCs, indicating that reduced expression of the ERC2 gene may not be attributable to DNA methylation alterations during renal carcinogenesis. Since the probes for the ABCA13 gene were not designed for the Infinium array, we examined DNA methylation levels in the 5′-region of the ABCA13 gene by pyrosequencing. No significant differences in the DNA methylation levels of the ABCA13 gene between T samples (0.528 ± 0.060, n = 67) and N samples (0.510 ± 0.149, n = 67) were observed (Supporting Information Fig. S2a). Our data for RCCs were consistent with the data in the public database Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/): no significant differences in DNA methylation levels of the ABCA13 gene were evident between bile duct cancer and normal bile duct tissue (Accession number: GSE49656) and between breast cancer and normal breast tissue (GSE37754), indicating that reduced expression of the ABCA13 gene may not be attributable to DNA methylation alterations during renal carcinogenesis.

Alterations of expression associated with DNA hypermethylation or hypomethylation

All genes showing DNA methylation alterations [0.2 or more Δβ (βΤ − βΝ) or −0.2 or less Δβ (βΤ − βΝ)] or mRNA expression alterations [4 or more ΔΕ (EΤ − EΝ) or −4 or less ΔΕ (EΤ − EΝ)] in each RCC are summarized in Supporting Information Table S6 along with genes showing genetic aberration scores of 1 or more. The DNA methylation status of the 5′-region can regulate the mRNA expression level of each gene. DNA methylation status is stably preserved on DNA double strands by covalent bonds and inherited through cell division by maintenance-methylation mechanisms by DNMT1. Therefore, altered mRNA expression due to DNA methylation alterations may be more stably fixed during multistage human carcinogenesis in comparison to mRNA expression alterations without DNA methylation alterations. Therefore, we have calculated upregulation and downregulation scores based on both DNA methylation status and expression levels described in the Material and Methods section: 86 genes showed reduced expression [−4 or less ΔΕ (EΤ − EΝ)] associated with DNA hypermethylation [0.2 or more Δβ (βΤ − βΝ)] in 5 or more patients (downregulation scores of 5 or more; Table 2) and 28 genes showed overexpression [4 or more ΔΕ (EΤ − EΝ)] associated with DNA hypomethylation [−0.2 or less Δβ (βΤ − βΝ)] in 5 or more patients (upregulation scores of 5 or more; Table 2).

Table 2

Genes showing downregulation or upregulation scores of 5 or more in clear cell RCCs

Gene	Chromo-some	Entrez GeneID	Downregulation score1
(a) Genes showing reduced mRNA expression associated with DNA hypemethylation in their 5′-regions
CLCNKB	1	1,188	24
SCNN1A	12	6,337	24
RAB25	1	57,111	22
TMEM213	7	155,006	22
ATP6V0A4	7	50,617	22
NR0B2	1	8,431	21
KCNJ1	11	3,758	21
GGT6	17	124,975	21
CLDN8	21	9,073	20
CLDN19	1	149,461	19
MUC15	11	143,662	16
RANBP3L	5	202,151	15
HRG	3	3,273	14
TSPAN8	12	7,103	14
RGS7	1	6,000	11
PTH1R	3	5,745	11
CWH43	4	80,157	11
F11	4	2,160	11
IRX2	5	153,572	11
EHF	11	26,298	11
CBLC	19	23,624	11
ATP6V1B1	2	525	10
LRRC2	3	79,442	10
CLDN16	3	10,686	10
EGF	4	1,950	10
WISP3	6	8,838	10
PHYHD1	9	254,295	10
FLJ45983	10	399,717	10
WIT-AS	11	51,352	10
ACSF2	17	80,221	10
ALDOB	9	229	9
ANKRD2	10	26,287	9
WT1	11	7,490	9
SOST	17	50,964	9
CYP4F3	19	4,051	9
COL18A1-AS1	21	378,832	9
BSND	1	7,809	8
TACSTD2	1	4,070	8
SLC44A4	6	80,736	8
KHDRBS2	6	202,559	8
VWC2	7	375,567	8
CHRM1	11	1,128	8
COL4A6	X	1,288	8
XPNPEP2	X	7,512	8
PROM2	2	150,696	7
ACPP	3	55	7
CKMT2	5	1,160	7
NEFM	8	4,741	7
KCNA4	11	3,739	7
FLRT1	11	23,769	7
OLFM4	13	10,562	7
SERPINA4	14	5,267	7
STRA6	15	64,220	7
CRABP1	15	1,381	7
SLC7A10	19	56,301	7
CSDC2	22	27,254	7
VWA5B1	1	127,731	6
LAD1	1	3,898	6
SYN2	3	6,854	6
SLC22A13	3	9,390	6
ABHD14A	3	25,864	6
UPK1B	3	7,348	6
KCTD8	4	386,617	6
SFRP1	8	6,422	6
GATA3	10	2,625	6
DAO	12	1,610	6
TMPRSS3	21	64,699	6
CHD5	1	26,038	5
PRELP	1	5,549	5
PLD5	1	200,150	5
MAL	2	4,118	5
ENTPD3	3	956	5
TNNC1	3	7,134	5
ANK2	4	287	5
PART1	5	25,859	5
SVOPL	7	136,306	5
DMRT2	9	10,655	5
AMBP	9	259	5
RBP4	10	5,950	5
SLC22A12	11	116,085	5
PDZRN4	12	29,951	5
PROZ	13	8,858	5
RHCG	15	51,458	5
KLK6	19	5,653	5
BEX1	X	55,859	5
ZCCHC16	X	340,595	5

If the probe of the Infinium array was designed in the 5′-region of the gene, if Δβ (βΤ − βΝ) was 0.2 or more (DNA hypermethylation) and if ΔΕ (EΤ − EΝ) based on the expression microarray was −4 or less (reduced expression) in one paired sample (T and N), then a gene downregulation score of 1 was assigned.

If the probe of the Infinium array was designed in the 5′-region of the gene, if Δβ (βΤ − βΝ) was −0.2 or less (DNA hypomethylation) and if ΔΕ (EΤ − EΝ) based on the expression microarray was 4 or more (overexpression) in one paired sample (T and N), then a gene upregulation score of 1 was assigned.

Expression alterations of genes included in Table 2 were validated using the clear cell RCC database in the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/; Supporting Information Table S7): reduced or increased mRNA expression of 97 (89%) of the109 genes, which are included in Table 2 and for which probes were designed in the expression microarrays described in the database, were found, indicating the reliability of our expression analysis. Since genome-wide DNA methylation data for RCCs obtained using array-based analysis with appropriate resolution were not available in the public database, Infinium assay data for other human malignant tumors deposited in the Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo/) were used instead for validation (Supporting Information Table S8). In addition, DNA methylation levels of the representative genes, RAB25, GGT6, C3 and CHI3L2, included in Table 2 based on the Infinium assay were successfully verified using pyrosequencing (Supporting Information Figs. S2b–S2e), indicating the reliability of our Infinium assay.

MetaCore pathway analysis by GeneGo was performed for 61 genes assigned genetic aberration scores of 3 or more, 86 genes assigned downregulation scores of 5 or more (frequent reduction of expression associated with DNA hypermethylation) and 28 genes assigned upregulation scores of 5 or more (frequent overexpression associated with DNA hypomethylation; total 174 genes). Twenty potentially significant GeneGo pathways (p < 0.05) and the affected genes are listed in Table 3. Mutations of 5 (100%) of the 5 genes included in Table 3 were found in the clear cell RCC database of The Cancer Genome Atlas (Supporting Information Table S5). Reduced or increased mRNA expression of 11 (92%) of the 12 genes, which are included in Table 3 and for which probes had been designed in expression microarrays described in the clear cell RCC database of the Gene Expression Omnibus, were found (Supporting Information Table S7), supporting the participation of these genes in renal carcinogenesis.

Table 3

Statistically significant GeneGo pathway maps revealed by MetaCore pathway analysis

		Involved genes
Pathway	P-value	Genes	Entrez Gene ID	Multilayer-omics scoring (exome, methylome and transcriptome)
Cell adhesion_tight junctions	9.98 × 10−4	CLDN8	9073	Downregulation score 20
		CLDN16	10686	Downregulation score 10
		CLDN19	149461	Downregulation score 19
Blood coagulation	1.26 × 10−3	VWF	7450	Upregulation score 6
		F11	2160	Downregulation score 11
		FGG	2266	Upregulation score 6
Translation_non-genomic (rapid) action of androgen receptor	1.36 × 10−3	MTOR	2475	Genetic score 4
		PTEN	5728	Genetic score 3
		EGF	1950	Downregulation score 10
Signal transduction_PTEN pathway	2.04 × 10−3	MTOR	2475	Genetic score 4
		PTEN	5728	Genetic score 3
		EGF	1950	Downregulation score 10
Development_EGFR signaling via PIP3	7.04 × 10−3	PTEN	5728	Genetic score 3
		EGF	1950	Downregulation score 10
Protein folding and maturation_Bradykinin/ Kallidin maturation	1.34 × 10−2	KLK6	5653	Downregulation score 5
		XPNPEP2	7512	Downregulation score 8
Transcription_receptor-mediated HIF regulation	1.95 × 10−2	MTOR	2475	Genetic score 4
		PTEN	5728	Genetic score 3
Serotonin modulation of dopamine release in nicotine addiction	2.24 × 10−2	PTEN	5728	Genetic score 3
		CHRM1	1128	Downregulation score 8
Signal transduction_AKT signaling	2.34 × 10−2	MTOR	2475	Genetic score 4
		PTEN	5728	Genetic score 3
cAMP/ Ca(2+)-dependent Insulin secretion	2.34 × 10−2	PLCE1	51196	Genetic score 3
		RYR2	6262	Genetic score 3
Immune response_interleukin-4 signaling pathway	2.45 × 10−2	MTOR	2475	Genetic score 4
		GATA3	2625	Downregulation score 6
Role of alpha-6/beta-4 integrins in carcinoma progression	2.55 × 10−2	MTOR	2475	Genetic score 4
		EGF	1950	Downregulation score 10
G-protein signaling_regulation of cAMP levels by muscarinic acetylcholine receptor	2.55 × 10−2	PLCE1	51196	Genetic score 3
		CHRM1	1128	Downregulation score 8
Development_PIP3 signaling in cardiac myocytes	2.77 × 10−2	MTOR	2475	Genetic score 4
		PTEN	5728	Genetic score 3
Some pathways of EMT in cancer cells	3.22 × 10−2	MTOR	2475	Genetic score 4
		EGF	1950	Downregulation score 10
Development_beta-adrenergic receptors signaling via cAMP	3.34 × 10−2	RYR2	6262	Genetic score 3
		TNNC1	7134	Downregulation score 5
Development_IGF-1 receptor signaling	3.34 × 10−2	MTOR	2475	Genetic score 4
		PTEN	5728	Genetic score 3
Translation _regulation of EIF4F activity	3.45 × 10−2	MTOR	2475	Genetic score 4
		EGF	1950	Downregulation score 10
G-protein signaling_RAP2B regulation pathway	3.81 × 10−2	PLCE1	51196	Genetic score 3
DNA damage_DNA-damage-induced responses	4.87 × 10−2	ATM	472	Genetic score 3

Genes for which correlation with Wnt/β-catenin signaling was indicated by MetaCore pathway analysis, together with their genetic aberration, DNA methylation alterations and mRNA expression alterations, are illustrated schematically in Figure 1. Mutations, mRNA expression alterations or DNA methylation alterations of 32 (89%) of the 36 genes included in Figure 1 were found in Supporting Information Tables S5, S7 or S8, supporting the participation of the Wnt/β-catenin signaling pathway in renal carcinogenesis. In addition, MetaCore pathway analysis was separately performed for RCCs with and without genetic aberrations and/or DNA hypermethylation [Δβ (βT − βN) >0.2] of the VHL gene (Supporting Information Table S9 and Fig. S3).

Figure 1

Genes for which a correlation with Wnt/β-catenin signaling was indicated by MetaCore pathway analysis. The numbers of genetic aberrations, DNA hyper- or hypo-methylation and/or increased or reduced mRNA expression (shown in Supporting Information Table S6) detected among the 67 examined RCCs are indicated schematically: legends are shown at the left of the panel. The 36 marked genes that showed genetic aberration, DNA methylation alterations and/or mRNA expression alterations in one or more RCCs were correlated with Wnt/β-catenin signaling.

Discussion

High frequencies of genetic aberrations of the VHL (53%), PBRM1 (33%), KDM5C (12%) and SETD2 (9%) genes, which have been highlighted in previous resequencing2 and exome analyses,4,6 supported the reliability of our approach. In addition to PBRM1, somatic mutation of another member of the SWI/SNF complex, SMARCA4, was detected. In addition to SETD2 and KDM5C, somatic mutation of another histone modification protein, JARID2, was also detected. The significance of aberrations of chromatin remodeling and histone modification-related proteins in RCCs was confirmed.

Among genes showing frequent genetic aberrations (genetic aberration score of 4 or more in Table 1), GCN1L1 has recently been reported to be associated with the CDK8 mediator complex, which includes CDK8, cyclin C (also known as CCNC), MED12 and MED13.25 CDK8 directly regulates β-catenin-driven transcription25 and human CDK8 is known to be an oncogene that is amplified in a subset of colon cancers.26 In addition, our quantitative RT-PCR analysis revealed a tendency for down regulation of β-catenin after knockdown of CDK8 by siRNA in RCC cell lines A-498 and ACHN (Supporting Information Fig. S4). These results are consistent with those of previous studies showing that knockdown of CDK8 in the human colon cancer cell line HCT11627 and the human gastric cancer cell line SNU-63828 resulted in significant reduction of β-catenin, indicating correlations between CDK8 and the Wnt/β-catenin pathway.

The fly MED12 and MED 13 homologs, kohtalo and skuld, respectively activate Wnt/β-catenin target genes through direct interaction with the Wnt pathway component Pygopus.29 However, let-19 and doy-22, homologs of human MED12 and MED13, respectively, in Caenorhabditis elegans, suppress the transcription of Wnt/β-catenin target genes.30 Frequent mutation of human MED12 has been reported in human uterine leiomyomas.31 Deletion of the CCNC gene is frequently detected in human lymphoid malignancies32 and sarcomas.33 Wnt/β-catenin signaling is constitutively active in RCCs and activates their cell growth and metastasis.34 However, unlike other human carcinomas, the incidence of mutation of exon 3 of the β-catenin gene is not so high in RCCs.34 Analogously with other members of the CDK8 mediator complex, mutations of GCN1L1 may participate in renal carcinogenesis via Wnt/β-catenin signaling.

All 5 amino acid substitutions of the GCN1L1 occurred within or near to Huntingtin protein, eEF3, protein phosphatase 2A and TOR (HEAT) repeats, which are crucial for protein-protein interaction35 (Supporting Information Fig. S5). In addition, SIFT and PolyPhen-2 software predicted that amino acid substitutions due to mutations of the GCN1L1 gene result in dysfunction of GCN1L1 protein (Table 1). The present study demonstrated not only a genetic aberration score of 5 for GCN1L1, but also a genetic aberration score of 3 for MED12 and CCNC (Table 1). SIFT and PolyPhen-2 analyses have predicted that amino acid substitutions due to mutations of the MED12 and CCNC genes also result in dysfunction of the proteins (Table 1). Taken together, the present data indicate that the function of the CDK8 mediator complex may have been disturbed in 16% of the examined 67 RCCs. Genetic aberrations in members of the CDK8 mediator complex may thus participate in the Wnt/β-catenin-related carcinogenetic pathway in clear cell RCCs.

MACF 1, a member of the plakin family of cytoskeletal linker proteins, regulates dynamic interactions between actin and microtubules to sustain directional cell movement.36 MACF1 is known to function in the Wnt signaling pathway through association with a complex containing axin, β-catenin, GSK3β and APC during mouse embryogenesis.36 Somatic mutation of MACF1 (Table 1) may also participate in the Wnt/β-catenin-related carcinogenetic pathway in clear cell RCCs. With respect to 29 RCCs for which transcriptome analysis was performed, mRNA expression levels of the targets genes of the Wnt/β-catenin signaling, such as MYC,37 MYCN,37 IGF2,38 POU5F1,39 SOX9,40 CYR61,41 ENPP242 and MITF,43 tended to be higher in the 8 RCCs with mutations of any of the GCN1L1, MED12, CCNC and MACF1 genes than in 21 RCCs without them (Supporting Information Table S10), indicating that such mutations may result in activation of Wnt/β-catenin signaling.

The downregulation score for the SFRP1 gene was 6: reduced expression associated with DNA hypermethylation of SFRP1 was frequent in clear cell RCCs. Members of the secreted frizzled-related protein (SFRP) family contain an N-terminal domain homologous to the cysteine-rich domain of the Wnt receptor Frizzled and lack a transmembrane region and the cytoplasmic domain required for transduction of signals into the cells.44 This enables SFRPs to downregulate Wnt/β-catenin signaling by competing with Frizzled for Wnt binding via their cysteine-rich domain. Silencing of SFRP1 due to DNA hypermethylation is known to result in activation of Wnt/β-catenin signaling.44

Since this study indicated possible alternative activation mechanisms (mutations of the GCN1L1, MED12, CCNC and MACF1 genes and reduced expression of SFRP1 due to DNA hypermethylation), we extensively examined Wnt/β-catenin signaling. MetaCore pathway analysis revealed that the 36 genes (marked in Fig. 1 and included in Supporting Information Table S6), which showed genetic aberration, DNA hypermethylation or hypomethylation and/or increased or reduced mRNA expression in one or more RCCs, are included in the Wnt/β-catenin signaling pathway. The present multilayer-omics analysis revealed that the Wnt/β-catenin signaling pathway may be of greater significance in renal carcinogenesis than was realized previously.

ERC2, which had a genetic aberration score of 4, is localized in presynaptic active zones and plays a critical role in neurotransmitter release.45 Interaction between ERC2 and the tandem PDZ protein syntenin-1, which is known to associate with many synaptic proteins, together with multimerization of ERC2 both promote the localization of syntenin-1 at presynaptic ERC2 clusters and contribute to the molecular organization of active zones.45 Although the significance of ERC2 in human cancers has remained unclear, frequent intragenic breaks in the ERC2 gene indicated disruption of ERC2 function in RCCs. In addition to recurrent genetic aberration, the present quantitative RT-PCR revealed frequent reduction of ERC2 expression in clear cell RCCs relative to the corresponding N samples. Although frequent genetic and transcriptional inactivation of ERC2 may be involved in renal carcinogenesis, further functional analysis of ERC2 in RCCs is needed.

ABCA13 is a member of ATP-binding cassette sub-family A (ABC1) and a transmembrane transporter.46 Xenobiotics, including anticancer drugs, are extensively metabolized by activation enzymes such as cytochromes P450 and conjugation enzymes such as glutathione S-transferases or glucuronide transferases. Biotransformation represented by ABC transporters represents another important component of xenobiotic metabolism. In addition, ABC transporters play a crucial role in the development of resistance through efflux of anticancer agents from cancer cells.46 The disease-free interval of patients with colorectal cancers treated by adjuvant chemotherapy is significantly shorter in patients with low ABCA13 transcript levels.47 In addition to recurrent genetic aberration (Table 1), the present quantitative RT-PCR revealed frequently reduced expression of ABCA13 in RCCs relative to the corresponding N samples. Our findings suggest that it may be necessary to pay more attention to aberrations of ABCA13 at both the genetic and expressional levels when deciding the indications for chemotherapy in patients with clear cell RCCs.

In Table 3 based on MetaCore pathway analysis, it is feasible that expression of CLDNs required for generating cation-selective paracellular channels48 was reduced in clear cell RCCs, which lack the original absorptive function of the renal tubule. Moreover, MTOR mutations were highlighted as one of the major disrupters of multiple cell signaling during renal carcinogenesis: the MTOR gene participated in 10 (50%) of the 20 significant pathways in Table 3. The mammalian target of rapamycin (mTOR) encoded by the MTOR gene is a serine/threonine kinase that regulates cell growth, proliferation and autophagy.49 mTOR inhibitors, such as rapamycin and its derivatives, are being introduced for targeted therapy of clear cell RCCs. Overactivation of mTOR is generally considered to be due to homozygous deletion of the PTEN tumor suppressor gene.50 However, all 4 mutations of the MTOR gene detected in this cohort were located close to the kinase domain (data not shown) and may be activating mutations, as a previous in vitro study has suggested that mutations located close to the kinase domain activate the mutant form of mTOR.50 In addition, all detected mutations of the MTOR gene showed a SIFT score of 0 and PolyPhen-2 scores of 0.998 or 0.999, strongly suggesting that all MTOR mutations affect protein function (Table 1 and Supporting Information Table S3). MTOR mutation may be a marker for predicting the sensitivity of clear cell RCCs to rapamycin therapy.

In summary, the present exome analysis has revealed frequent genetic aberrations of GCN1L1, MED12, CCNC, MACF1, ERC2, ABCA13 and MTOR in clear cell RCCs. In addition to confirming the significance of aberrations of chromatin remodeling and histone modification-related proteins, the present multilayer-omics analysis has highlighted the significance of dysregulation of the Wnt/β-catenin signaling pathway including CDK8 mediator function, as well as the need to pay closer attention to MTOR mutations, causing major disruption of cell signaling during renal carcinogenesis, in relation to chemosensitivity. Multilayer-omics analysis can be considered a powerful tool for revealing significant carcinogenetic pathways in human cancers.