Genome-Wide Profiling of TRACK Kidneys Shows Similarity to the Human ccRCC Transcriptome

Renal cell carcinoma (RCC) is the most common primary cancer arising from the kidney in adults, with clear cell RCC (ccRCC) representing approximately 75% of all RCCs (2, 3). Histologically, ccRCC cells are characterized by a transparent cytoplasm caused by deposition of glycogen, phospholipids, and neutral lipids, including cholesterol esters (4, 5). This phenotype suggests that there are metabolic changes in ccRCC cells, resulting in abnormal deposition of glycogen and lipids. Seven genes commonly mutated in human kidney cancer, including VHL (NCBI Gene ID: 7428), MET (4233), FLCN (201163), TSC1 (7248), TSC2 (7249), FH (2271), and SDH (6390, 6391, and 6392), have been identified to date (6). Interestingly, all seven genes are involved in the regulation of metabolic pathways (6). These data support the theory that kidney cancer is a metabolic disease (6, 7).

Loss of expression or mutation of the von Hippel–Lindau (VHL) tumor suppressor gene is found in hereditary and most sporadic ccRCCs (2, 8). This suggests an etiological role for VHL gene loss in renal carcinogenesis. However, the exact pathway by which loss of VHL leads to ccRCC has not been definitively elucidated. The best studied and likely most important effect of VHL loss is the resulting increase in expression of the alpha subunits of hypoxia-inducible factors 1 (HIF1α) and 2 (HIF2α; refs. 9–11). Increased expression of these two transcription factors has been proposed as a key step in ccRCC carcinogenesis (refs. 9, 10; for review, see ref. 12). We previously generated the murine TRAnsgenic model of Cancer of the Kidney (TRACK) that expresses a triple-mutant (P402A, P564A, and N803A) human HIF1α construct in murine proximal tubule cells (PTCs) and showed that this model mimics the histologic alterations found in early stage human ccRCC (1). The cellular histologies displayed in TRACK mice are also similar to those observed in the kidneys of individuals with VHL disease, including areas of distorted tubular structure, cells with clear cytoplasm and increased glycogen and lipid deposition, multiple renal cysts, and areas of multifocal ccRCC (1).

To delineate the gene expression pattern in TRACK transgenic positive (TG⁺) kidneys, we sequenced the entire transcriptome of the TRACK TG⁺ kidney cortex by high-throughput sequencing of cDNA libraries (RNAseq) and compared it with both the gene expression profile in WT/transgenic negative (TG⁻) kidneys and the gene expression profile observed in sporadic human ccRCC using the Oncomine database (Compendia Bioscience) and The Cancer Genome Atlas (TCGA) database. We report that the pattern of gene expression in TRACK TG⁺ kidneys is similar to that of human ccRCC, including expression of transcripts associated with glycolysis/gluconeogenesis, the pentose phosphate pathway, and the lipid metabolism pathway. These data provide evidence that constitutive activation of HIF1α in kidney PTCs transcriptionally reprograms the regulation of metabolic pathways in a manner similar to that observed in human ccRCC.

TRACK TG⁺ and TG⁻ mice were housed five per cage in a 12-hour light/dark cycle in the Research Animal Resource Center of Weill Cornell Medical College (WCMC). The care and use of animals in this study were approved by the Institutional Animal Care and Use Committee of WCMC.

Total RNA from one thin, outer slice of kidney cortex per kidney was used for whole transcriptome sequencing. Total RNA was extracted using mini-RNAeasy columns (Qiagen). The complete transcriptomes of kidney cortices from three γ-HIF1αM3 TRACK (#43 line) TG⁺ and three γ-HIF1αM3 TG⁻ male C57BL/6 mice (about 13 months old) were sequenced (51-bp single-end reads) on an Illumina HiSeq2000 Sequencer following standard protocols in the Genomics Resources Core Facility at WCMC. Three lanes were used to sequence all 12 samples (4 samples/lane, 2 kidney samples/mouse). The TRACK RNAseq data have been deposited in the GEO database (accession no. GSE54390).

Data analysis was mainly performed with the Tuxedo tools software (13). In brief, the RNAseq reads were first aligned to the Mus musculus genome (UCSC version mm10) using Tophat version 2.0.6 (13, 14). The aligned reads were assembled into transcripts, their abundance was estimated, and they were tested for differential expression using Cufflinks version 2.1.1 (13, 15). CummeRbund version 2.0.0 (13) was used to analyze the differential expression analysis results. To identify pathways changed in the TRACK TG⁺ versus TG⁻ kidneys, functional enrichment analysis was performed using the goseq package in R software (16). A stringent threshold in selecting differentially expressed (DE) genes (fold change > 3, q < 0.01) was used to reduce the false positive ratio. Heatmaps of log₂ transformed RPKM (Reads Per Kilobase per Million mapped reads) values were created in R using the heatmap.2 command of the gplots package.

Human ccRCC gene expression changes were retrieved from Oncomine (Compendia Bioscience) by combining five different datasets of human ccRCC patient samples (17–20). The same five Oncomine datasets of Cancer versus Normal Analysis of ccRCC that we used in the γ-HIF2αM3 TG⁺ RNAseq analysis (21) were used in this analysis.

Human ccRCC mRNA data were downloaded from TCGA (http://cancergenome.nih.gov/). Only tumor patient data with matched normal and normal patient data with matched tumor were downloaded. All data satisfying this requirement were downloaded, including a total of 470 tumor samples and 68 normal samples. Differential expression between ccRCC and normal kidneys was calculated using the downloaded RPKM values. Statistical analyses were performed by the Student t test followed by FDR adjustment (q value). Statistical significance was defined as q < 0.05.

Expression of a mutated, constitutively active HIF1α in the PTCs of the γ-HIF1αM3 TRACK mice results in early stage tumors morphologically similar to human ccRCC (1). To identify changes in gene expression associated with ccRCC carcinogenesis, we examined the whole transcriptome from cells in TRACK kidney cortex slices compared with TG⁻ controls. At the time of sacrifice (∼13 months old), about 50% of the proximal tubules in TRACK mice show clear cell abnormalities. However, further abnormalities, e.g., carcinoma in situ, are not ubiquitously seen in the kidney cortex. The average number of reads per sample was approximately 45.6 million, and approximately 96% of reads mapped to the Mus musculus genome (Supplementary Table S1). A scatter plot of the RPKM values of TRACK TG⁺ versus TG⁻ kidneys shows that the majority of transcripts evaluated display no change between these two samples (Fig. 1A), but there are transcripts that show increased or decreased levels in the TRACK TG⁺ versus TG⁻ kidneys (Fig. 1A). Changes in some of these transcripts have been confirmed by semiquantitative RT-PCR (ref. 1 and Supplementary Fig. S1). Principal component analysis shows that there is a clear distinction between TRACK TG⁺ and TG⁻ kidneys (Fig. 1B).

We have shown that the high expression of CA-IX, Glut1, and VEGF proteins in the TRACK kidneys is mainly localized in the clear cell proximal tubules (1). Here, we also used immunohistochemistry to examine the protein levels of NDUFA4L2, SLC16A3, and HK2, three of the top genes overexpressed in TRACK kidneys by RNAseq. We detected high expression of NDUFA4L2, SLC16A3, and HK2 primarily in the clear cell proximal tubules (Supplementary Fig. S2). These transcripts are also highly expressed in human ccRCC (see next section).

We first examined the overrepresentation of 274 DE genes (259 overexpressed and 15 underexpressed) in KEGG (Kyoto Encyclopedia of Genes and Genomes). KEGG is a database resource for understanding high-level functions and utilities of biologic systems (22, 23). The KEGG PATHWAY is a collection of manually drawn pathway maps of molecular interactions and reaction networks. The 274 DE genes are overrepresented in 5 KEGG pathways (q < 0.05, top 10 pathways are shown in Table 1). Interestingly, seven of the 10 KEGG pathways are related to metabolism, and five of these seven KEGG pathways have the lowest P values in the list (Table 1), indicating that metabolic changes are the most prominent changes in TRACK TG⁺ kidneys. The most significant KEGG pathway is mmu00010: glycolysis/gluconeogenesis (P = 3.87E−06). The pentose phosphate and PPAR signaling pathways are also emphasized.

We also examined the overrepresentation of DE genes in the Gene Ontology (GO) consortium (24). The GO covers three domains: cellular component, molecular function, and biologic process, and we focused on the biologic process domain. The 274 DE genes are overrepresented in 98 GO biologic process ontologies (q < 0.05; Supplementary Table S2). Similar to the KEGG pathway analysis, seven of the top ten overrepresented GO ontologies are metabolism related, e.g., glucose metabolic process (P = 3.57E−10), glucose catabolic process (P = 2.36E−08).

The most significant metabolic change that occurs in ccRCC, and probably in most cancers, is the shift from oxidative phosphorylation of glucose to glycolysis under normoxic conditions, known as the Warburg effect (25). We examined the transcript levels of key genes involved in glycolysis and the tricarboxylic acid (TCA) cycle in the TRACK TG⁺ kidneys versus TG⁻ kidneys (Fig. 2). Most glycolysis gene transcripts are increased in the TRACK TG⁺ kidneys compared with the TG⁻ kidneys, e.g., hexokinase, phosphofructokinase, and pyruvate kinase. TCA cycle genes do not show large decreases in transcript levels in TRACK TG⁺ kidneys compared with TG⁻ kidneys (fold change between 0.8 and 0.9), presumably as a result of the presence of nontransformed cells in the cortex samples. Transcripts encoding pyruvate dehydrogenase kinase, the enzyme that inactivates pyruvate dehydrogenase through phosphorylation (26), are increased in the TRACK TG⁺ kidneys compared with TG⁻ kidneys (fold change = 7.8). Inactivation of pyruvate dehydrogenase by pyruvate dehydrogenase kinase prevents the conversion of pyruvate to acetyl-CoA (26), which is oxidized in mitochondria in the TCA cycle. The lower level of acetyl-CoA thereby inhibits the TCA cycle. Consistent with this, the transcript for lactate dehydrogenase, which converts pyruvate to lactate, is increased in the TRACK TG⁺ relative to TG⁻ kidneys. All of these changes indicate an important role for HIF1α in changing the metabolism of kidney cells to glycolysis under normoxic conditions. Heatmaps created from genes involved in glycolysis and in the TCA cycle also confirmed these changes (Fig. 3).

Using a fold-expression change of >2 or < 0.5 and q < 0.05, 655 genes showed increased mRNA levels and 55 genes showed decreased mRNA levels in TRACK TG⁺ kidney samples compared with the TG⁻ control kidney cortex samples. As expected, several HIF1α target genes were highly expressed at the transcript level in TRACK kidney samples, including hexokinase 2 (HK2, 70X), carbonic anhydrase IX (CA-IX, 9.3X), glucose transporter 1 (Glut-1 or Slc2a1, 4.8X) etc., which were confirmed by semiquantitative RT-PCR (ref. 1; Supplementary Fig. S1).

We first compared the expression profiles of human ccRCC, as reported in the Oncomine database (27), with our Track TG⁺ expression profile. We identified the 20 genes most highly overexpressed and 20 genes most highly underexpressed at the RNA level in human ccRCC by combining five different datasets of human ccRCC patient samples (17–20). The total number of ccRCC patients in all five data sets is 175. Five of the 20 genes most highly overexpressed in human ccRCC were similarly expressed in TRACK TG⁺ kidneys (fold change > 2; Table 2): NDUFA4L2, C7ORF68, EGLN3, SLC16A3, and CA-IX. None of the 20 genes highly underexpressed in human ccRCC shows significant downregulation in our TRACK TG⁺ kidneys (fold change < 0.5; Table 2). However, the number of genes that shows reduced expression in TRACK TG⁺ versus TG⁻ kidneys is most likely an underestimate because only 30% to 50% of PTs in TRACK TG⁺ kidneys demonstrate morphologic changes (1).

We also compared the expression profiles of human ccRCC mRNA downloaded from the TCGA database with our Track TG⁺ expression profile. We identified the 20 genes most highly overexpressed and the 20 genes most highly underexpressed at the RNA level in human ccRCC from the TCGA data. Five of the 20 genes most highly overexpressed in human ccRCC were similarly expressed in TRACK TG⁺ kidneys (fold change > 2; Table 3). None of the 20 genes highly underexpressed in human ccRCC shows significant downregulation in our TRACK TG⁺ kidneys (fold change < 0.5; Table 3).

We then performed the reverse comparison by identifying the top genes overexpressed or underexpressed at the mRNA level in the TRACK TG⁺ versus TG⁻ WT kidneys and comparing these transcripts to those in the TCGA database and the combined Oncomine human ccRCC datasets. Eleven of the 20 genes highly overexpressed in TRACK TG⁺ kidneys show overexpression (fold change > 2) in the TCGA data (Table 4). Ten of the 20 genes highly underexpressed in TRACK TG⁺ kidneys show underexpression (fold change < 0.5) in the TCGA dataset (Table 4). Four of the 20 genes highly overexpressed in TRACK TG⁺ kidneys show overexpression (fold change > 2) in the combined Oncomine datasets (Table 4). Eight of the 20 genes highly underexpressed in TRACK TG⁺ kidneys show underexpression (fold change < −2) in the combined Oncomine datasets (Table 4). We conclude from analysis of these data that expression of a mutant, constitutively active HIF1α protein in kidney PTs results in a transcriptome that partially resembles those of human ccRCCs.

We previously established the TRACK model, which mimics early stage human ccRCC through expression of a mutated, constitutively active human HIF1α construct in the PTCs (1). Here, we present genome-wide transcriptome analysis of the TRACK TG⁺ kidney cortex cells, which are mainly PTCs. We identified 655 upregulated genes and 55 downregulated genes that are differentially expressed in the TRACK TG⁺ kidneys compared with the TG⁻ kidneys (fold change > 2). Some of these genes also show increased or decreased transcript levels in human ccRCC specimens compared with normal kidneys, respectively (Table 4). For example, NDUFA4L2 is the gene overexpressed to the greatest extent in the human ccRCC datasets in Oncomine. The median fold change in human ccRCC compared with normal kidney from all 5 datasets is 53.9 (Table 2). Similarly, NDUFA4L2 transcript levels are about 68-fold higher in TRACK TG⁺ kidneys compared with TG⁻ kidneys (Table 2). Similarities can also be seen in genes that are downregulated in TRACK TG⁺ versus TG⁻ and human ccRCC versus normal kidneys, e.g., 4-hydroxyphenylpyruvic acid dioxygenase (HPD; Table 4). These results suggest that the transcriptome of the TRACK TG⁺ kidneys partially resembles that of human ccRCC cells. Furthermore, we have examined the expression patterns of some of these top genes overexpressed in the TRACK kidneys by immunohistochemistry (ref. 1; Supplementary Fig. S2). The high expression of these proteins occurs primarily in the clear cell proximal tubules. These data support our contention that the changes in the transcriptome of TRACK kidneys that we have seen by RNAseq are mainly, if not all, caused by changes in these clear cell proximal tubules.

New evidence suggests that ccRCC is a metabolic disease (6, 7). All known kidney cancer susceptible genes are involved in the regulation of metabolic pathways (6). ccRCC cells contain increased amounts of glycogen and lipid in their cytoplasm (4, 5). Increased glycogen and lipid are also seen in the TRACK TG⁺ kidney cells (1). Our analysis of the entire transcriptome of the TRACK TG⁺ kidney PTCs identified some altered metabolic pathways, including increased transcript levels of genes involved in glycolysis/gluconeogenesis (Fig. 2A) and decreased transcript levels of genes involved in the TCA cycle (Fig. 2B). Furthermore, increased mRNA levels of the protein pyruvate dehydrogenase kinase can inactivate pyruvate dehydrogenase, which converts pyruvate to acetyl-CoA (26). Decreased levels of acetyl-CoA, together with decreased transcript levels of genes involved in the TCA cycle, should result in decreased activity of the TCA cycle. The products of the TCA cycle, NADH and succinate, are used in the oxidative phosphorylation pathway. Decreased TCA cycle activity and decreased levels of NADH and succinate should result in lower levels of substrates, and as a result, decreased oxidative phosphorylation. This phenotype recapitulates much of what is described for the Warburg effect in tumor cells (25). HIF1α, rather than HIF2α, is the main regulator of these pathways (28, 29), further emphasizing the importance of HIF1α in modulating the metabolic alterations in ccRCC.

The histone deacetylase, Sirtuin-6 (SIRT6), can suppress the transcription of the HIF1α gene (30). Increased glucose uptake and increased HIF1α activity were shown in SIRT6-deficient cells and mice, and this increase in HIF1α activity was sufficient to cause tumorigenesis (30). SIRT6 has been reported to be a potential tumor suppressor gene (31, 32). Thus, either the absence of VHL or the lack of SIRT6 can result in an increase in HIF1α activity and tumorigenesis. In summary, our results suggest that constitutive activation of HIF1α in kidney PTCs in our TRACK model induces a phenotype similar to the tumor phenotype associated with the Warburg effect (25).

In addition to inactivating mutations in the VHL gene, inactivating mutations in other genes, e.g., PBRM1 (33), SETD2 (34), etc., are commonly seen in patients with advanced ccRCC (35–38). Mutations and/or changes in the expression of these genes probably play important roles in the development of advanced ccRCC. This might explain why TRACK mice only mimic early stage ccRCC. This hypothesis is being tested in the TRACK mouse model by mutation/knockout of additional genes, such as PBRM1 and SETD2.

Importantly, the transgenic mice we generated that overexpress a constitutively active HIF2α do not show major transcript changes that reflect those in human ccRCC (21). Thus, our data reported here and those of many other researchers (reviewed in ref. 12) indicate that the activity of HIF1α is more closely linked to ccRCC tumorigenesis than the activity of HIF2α.

No potential conflicts of interest were disclosed.

Conception and design: L. Fu, D.R. Minton, D.M. Nanus, L.J. Gudas

Development of methodology: L. Fu, L.J. Gudas

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): L. Fu, D.R. Minton, L.J. Gudas

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): L. Fu, T. Zhang, L.J. Gudas

Writing, review, and/or revision of the manuscript: L. Fu, D.R. Minton, T. Zhang, D.M. Nanus, L.J. Gudas

Study supervision: D.M. Nanus, L.J. Gudas

The authors thank Daniel Stummer for editorial assistance and the Gudas and Nanus labs for thoughtful discussions of data.

This research was supported by WCMC, the Genitourinary Oncology Research Fund, and the Turobiner Kidney Cancer Research Fund. L. Fu held the Robert H. McCooey Genitourinary Oncology Research Fellowship and is supported by NCI grant NIH T32 CA062948. D.R. Minton is supported by an NIH T32 Training Grant (5T32CA062948). The results published here are in part based upon data generated by the TCGA Research Network: http://cancergenome.nih.gov/.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.