Genomic Deregulation during Metastasis of Renal Cell Carcinoma Implements a Myofibroblast-Like Program of Gene Expression

Cancer involves multistep changes in the genome that lead to deregulation of gene profiles and disruption of molecular networks. A variety of genomic approaches have been used to identify the molecular profiles that contribute to and reflect cancer progression. Moreover, gene expression profiling has been useful in the prediction of clinical outcome, disease progression, and metastatic recurrence (1). During carcinogenesis, genetic and epigenetic alterations drive tumor evolution toward higher grades of malignancy; however, the extent to which these changes alter gene expression, and how these alterations influence each other, remains incompletely understood. Association between gene dosage and gene expression levels has been reported in a variety of tumors (2), and approximately 12% of gene expression variation can be explained by differences in DNA copy number (3). In addition, integrated analyses of genome copy number and gene expression have provided insights into the etiology and mechanisms by which DNA copy number changes contribute to tumor evolution (2). Genomic methylation patterns are clearly deranged in cancer cells (4); however, the impact of promoter methylation on genome-wide gene expression changes during the pathogenesis of individual cancers still is poorly defined. Recently, the impact of global cancer-related changes in DNA methylation, genomic imbalance, and gene expression has been analyzed using whole genomic profiling approaches (5, 6). Alternative approaches using recent advances in computational biology have successfully deciphered the functional impact of transcription factor networks on tumor progression (7).

Kidney cancer is the seventh most common cancer in the United States. Clear cell renal cell carcinoma (RCC) arises from the renal epithelium and accounts for approximately 85% adult kidney cancers (8). Hereditary type RCC is associated with von Hippel-Lindau disease (VHL), which arises as a result of mutations in the VHL gene; the latter is also inactivated by deletion, mutation, or gene silencing in most sporadic cases (9). In approximately 30% of patients, metastases are detected at the time of primary diagnosis, and an additional 30% to 50% of initially localized tumors progress to distant metastases (10). The genetic mechanisms of RCC progression, especially those which underlie metastatic conversion, are poorly understood.

The prevailing multistep progression model of tumorigenesis assumes that random mutations followed by selection eventually give rise to the metastatic phenotype in a rare subpopulation of tumor cells (11, 12). Early genomic studies suggested that the primary tumor may or may not carry a molecular signature indicative of high metastatic potential (13). More recently, Massague and colleagues extensively characterized molecular signatures that allow breast tumor cells to colonize specific secondary sites (14–16). It has also been proposed that metastasis could arise from early disseminated tumor cells based on studies that showed that disseminated tumor cells can be identified in the secondary site very early on, even before the primary tumor becomes clinically apparent (17, 18). Following up on Paget's “seed and soil” hypothesis (19), proposed more than a hundred years ago, it has been suggested that metastatic growth requires the existence of a “metastatic niche” (20).

We used an integrated genomic approach and performed a comprehensive analysis of clinical and xenograft models of RCC metastasis. For this, we used the RCC-derived and nonmetastatic SN12C cell line and a collection of derived cell variants capable of efficient spontaneous metastasis to the lung when injected orthotopically. Gene expression profiling of nonmetastatic and metastatic tumor cell lines revealed that metastatic activity was associated with the acquisition of a mesenchymal signature. Many of the genes in this signature are specifically enriched in the prototypical fibrotic cell. Notably, some of these genes are also know to control metastatic activity. By associating gene expression profiles and genetic and epigenetic alterations in the metastatic cells, we were able to gain new insights into the mechanisms by which the genome of RCC tumor cells is rearranged during metastasis, and to identify functionally relevant genes in this model.

Full methods are described in supplemental methods.

The SN12C cell line was established in culture from a primary RCC from a 43-year-old male patient and was described previously (21). Derivation of the lung-metastatic LM1 cells was also previously described (22). LM2 cells were isolated from lung metastasis of NOD/SCID mice injected with 10⁶ LM1 cells in the kidney subcapsule and sacrificed 2 months later. Lung nodules from each mouse were dispersed into separate culture wells to generate individual cell lines. The cells were cultured in minimum essential medium supplemented with 10% fetal bovine serum. All animal studies were performed in accordance with a protocol approved by the Institutional Animal Care and Use Committee; 6- to 8-week-old NOD/SCID mice (NCI) were used for all xenograft studies.

Mouse lungs were harvested at necropsy. For hematoxylin and eosin staining, tissues were fixed overnight in 10% neutral buffered formalin, washed with PBS, and dehydrated in 70% ethanol prior to paraffin embedding (Histoserv Inc.).

Fresh-frozen paired primary and metastatic tumor biopsies were obtained from 16 patients (32 samples) with RCC evaluated at the Memorial Sloan-Kettering Cancer Center following an Institutional Review Board–approved protocol.

RNA was isolated from cell lines and tumor tissues using an RNeasy kit (Qiagen). For expression profiling, the HG-U133 Plus 2.0 array (Affymetrix) was hybridized with cRNA from the samples using standard protocols. For microRNA (miRNA) analysis, total RNA was isolated using the Qiagen miRNeasy kit and hybridized to Agilent Human miRNA Microarray kit V2-G4470B (Agilent Technologies) following labeling and hybridization protocols recommended by the manufacturer.

Cells were lysed using sample buffer and proteins were separated by SDS-PAGE and blotted onto nitrocellulose membranes and hybridized with the anti-S100A4 antibody ab27957 (Abcom).

Genomic DNA was isolated using the Qiagen DNA extraction kit. DNA digestion, labeling, and hybridization were performed following Agilent's protocol version 4.0 for Agilent Human Genome CGH 244A oligo microarrays.

Cells were seeded (1×10⁶) in culture medium and maintained for 24 hours before treatment with 5 μmol/L 5 aza 2′ deoxycytidine (5′-AZA; Sigma-Aldrich) for 3 days (23). Medium containing 5′-AZA was replaced every 24 hours during the treatment. Control cells were treated in the same way, but without the addition of 5′-AZA.

For overexpression, S100A4 cDNA was subcloned into the pQCXIN retroviral-based vector (Clontech). For knockdown, the pLKO.1 plasmid (clone TRCN0000053608) encoding short hairpin RNA (shRNA)-targeting S100A4 was obtained from Open Biosystems. Viral supernatants were generated by transfecting 293-FT cells with the shRNA constructs in combination with the packaging vectors pVSVG and pDR2.

SN12C cells, when implanted in the renal subcapsule of NOD/SCID mice, showed limited lung metastatic activity (3–10 nodules). Lung-derived nodules, when expanded in culture and reinoculated into the kidney capsule showed higher lung metastatic activity (80 nodules), designated LM1 (22). Another round of selection in vivo yielded a second generation of highly metastatic cells (400 nodules), designated LM2 (Fig. 1B). Interestingly, the LM cells also showed an increased metastatic spread to other organs (ref. 22; our results). Of note, the growth rate of the highly metastatic cells at the primary site was similar to that of the parental SN12C cells (Fig. 1A). Histologic review of lung metastasis showed that the tumor cells, in addition to colonizing the lung parenchyma, associated in thrombus-like structures inside lymphatic and blood vessels (Fig. 1C).

To evaluate genome wide changes in gene expression during metastatic conversion, we subjected the 3 cell types, SN12C, LM1, and LM2, to gene expression profiling using the Affymetrix U133 plus 2.0 arrays. Class comparison between SN12C and LM2 cells identified 349 unique genes that showed differential expression: 191 were upregulated and 158 were downregulated in the most aggressive metastatic cells. Analysis of the gene expression patterns revealed a mesenchymal signature as the dominant feature of the gene list that characterized the metastatic cells (Table 1; Supplementary Fig. S1). Epithelial–mesenchymal transition (EMT) has been proposed to be one of the initial steps during metastatic conversion that endows tumor cells with the ability to invade surrounding tissues (24). We therefore investigated whether the observed mesenchymal phenotype of LM cells correlated with an increase in invasiveness. Surprisingly, the ability of LM cells to invade through both an artificial basement membrane (Boyden chamber assay) and an artificial endothelium (human umbilical vein endothelial cell) was not significantly altered when compared with parental SN12C cells (Supplementary Fig. S2). This result suggested that the SN12C cells may already have undergone the initial stages of EMT. Consistent with this assumption, we found that epithelial markers such as E-cadherin were downregulated whereas mesenchymal markers such as N-cadherin and vimentin were already upregulated in the SN12C poorly metastatic cells (Supplementary Fig. S2C). Emerging evidence suggests that epithelial cells are also an important source of myofibroblasts in fibrosis and cancer (25). Therefore, we reasoned that a significant function of tumor EMT would be to provide tumor cells with stromal myofibroblastic properties, which would bypass the requirement of an activated stroma during metastatic colonization. Detailed analysis of mesenchymal genes upregulated in the LM cells revealed a high enrichment for expression of genes related to fibrosis and desmoplastic response (S100A4, COL11A1, COL1A1, COL3A1, CALD1, and TNNC1; refs. 26, 27). Conversely, the expression of some of the genes characteristic of the epithelial lineage (KRT34, DSC3, and KRT33A) was lost in the metastatic cells (Supplementary Fig. S3).

To validate the observations from SN12C and its metastatic derivatives in the clinical disease described earlier, we subjected a cohort of paired primary and metastatic biopsies from 16 RCC patients to gene expression profiling using the Affymetrix U133 plus 2.0 arrays. The metastatic components in these pairs were also highly enriched for genes typically associated with the myofibroblast phenotype (SFRP2, COL11A1, OGN, COL14A1, COL1A2, DCN, ASPN, and CALD1; refs. 26, 27; Table 2). Collectively, these results indicate that during RCC metastatic conversion, tumor cells acquire a myofibroblast-like phenotype.

To formally demonstrate that the myofibroblast-like phenotype of the highly metastatic cells affects the metastatic behavior, we selected, from among a subset of genes biologically relevant to the mesenchymal phenotype, S100A4, a fibroblast marker for functional validation. We were particularly interested in S100A4 because this molecule could serve as a functional link between fibrosis and metastasis during kidney cancer progression, as strong evidence suggests S100A4 to be a crucial myofibroblast-expressed factor that regulates metastasis (28). To determine whether this gene plays a causal role in lung metastasis, we stably decreased its expression in LM2 cells using shRNA and found that S100A4 levels had a profound effect on metastatic activity (Fig. 2A). To investigate whether S100A4 was sufficient to induce metastatic spread, we engineered SN12C cells to overexpress S100A4 at levels comparable to those seen in LM2 cells. Injection of SN12C-S100A4 cells into the kidney subcapsule showed a very modest lung colonization advantage in comparison with control SN12C cells (Fig. 2B) suggesting that S100A4 is required but not sufficient for LM2 cells to metastasize, and that cooperation of other genes may be required for the attainment of full metastatic competence. Globally, these results support a model in which RCC cells co-opt a myofibroblastic gene expression program to colonize secondary sites during tumor progression.

We performed high-resolution whole genome profiling of SN12C and derived metastatic LM1 and LM2 cells using the Agilent High-Density a-CGH Human 244A array to cross reference gene dosage and gene expression changes. The resultant data were analyzed with the Partek Genomic Suite. This analysis identified approximately 50 significant regions of chromosome copy number differences that distinguished the nonmetastatic and metastatic cells, spanning from 2 Kbp to 80 Mbp. These regional changes in chromosome copy number were scattered over 14 chromosomes with a certain degree of enrichment on chromosomes 1, 2, and 4. Of the 50 aberrations, 10 were gains and 40 were losses (Supplementary Fig. S4). Comparison of probe heat maps from the 3 different genomes showed complex differences (Fig. 3A). To identify genes belonging to regions of significant copy number change, the results of the segmentation algorithm of array comparative genomic hybridization (a-CGH) profiles were annotated with the Agilent human genome annotation library and further aligned with the list of genes differentially expressed between the 3 groups using the Partek Genomic Suite. This analysis yielded a list of 57 genes that mapped to chromosomal regions of DNA copy number changes (Fig. 3B). Only 1 gene mapped to a significant region of gain on chromosome 16, whereas 56 genes mapped to regions of significant loss on chromosomes 1, 2, 4, 10, 13, and 16 (Supplementary Fig. S5). Interestingly, 50 of the 57 genes mapped to chromosome 4. None of the genes defining the myofibroblast-like phenotype mapped to regions of DNA copy number variation, arguing against a role for genomic rearrangements during the acquisition of the myofibroblastic phenotype.

Because acquired changes in gene expression may be influenced by both genetic and epigenetic factors, we set out to study the contribution of promoter methylation in the establishment of the metastatic phenotype in our model of tumor progression. To identify epigenetically silenced genes, we used an epigenetic reactivation strategy that combined treatment of cancer cells in vitro with the DNA methyltransferase inhibitor 5′-AZA, followed by global gene expression analysis using microarrays. Treatment of SN12C cells resulted in the upregulation of more than 600 unique genes. To uncover genes upregulated by promoter demethylation during the SN12C-LM transition, we performed a Venn analysis using a dataset of genes upregulated in SN12C cells after 5′-AZA treatment and a dataset of genes spontaneously upregulated in LM cells (Fig. 4A). By this method we identified 21 genes, including S100A4, CALD1, COL9A3, CRIP1, CRIP2, COL1A1, and TNNC1, as hypomethylated in the LM cells, supporting a metastasis promotion role for those genes. Conversely, treatment of LM cells with the same inhibitor followed by analysis of reexpressed genes and Venn analysis of overlapping genes with the list of genes downregulated during SN12C-LM transition uncovered 31 hypermethylated candidate genes (Fig. 4B). Remarkably, some of the hypermethylated and silenced genes in the LM cells (KRT34 and KRT33A) have been known to be associated with an epithelial phenotype (29). Overall, 52 of 300 genes (∼20%) deregulated during metastatic progression were upregulated after demethylation agent treatment, and therefore had their expression likely regulated by changes in promoter methylation. A total of 48 genes harbored CpG islands in their promoters. Moreover, 70 and 10 genes were downregulated by drug treatment of SN12C and LM cells, respectively, probably reflecting the reexpression of a repressor (data not shown). To validate the results from the 5′-AZA reexpression experiment, we selected 4 genes predicted to be demethylated during metastatic conversion (S100A4, COL1A1, CRIP1, and CRIP2) and subjected them to quantitation of CpG methylation using the EpiTYPER (Sequenom) mass-spectrometry analysis of bisulfite-converted DNA (Fig. 4D). Of these 4 genes, 3 showed decreased CpG methylation in the LM cells, as predicted. SN12C was also demethylated in a cohort of clinical samples as the tumor progressed to a metastatic stage (Supplementary Fig. S6).

Taken together, these results suggest that a number of genes undergo methylation-associated silencing or reexpression in the course of RCC tumorigenesis and metastasis. Even more interestingly, many of the methylated genes are associated with a myofibroblastic phenotype, suggesting that promoter methylation may be one of the pivotal events controlling metastatic conversion.

Tumor-acquired promoter methylation often coincides with allele loss. To determine whether any of the 52 hypo- and hypermethylated candidate genes were also subject to DNA copy number changes, we cross-referenced this list of genes with the dataset of genes mapped to regions of significant DNA dosage changes, using the Venn diagram tool. We determined that only 3 of the 52 genes whose expression was regulated by methylation belonged to regions of significant copy number losses (Fig. 4C). However, we did not identify any candidate-methylated genes mapping to any of the regions of gene dosage gain. These results suggest a complex pattern in the genetic control of gene expression in which a subset of genes is synergistically impacted by both promoter methylation and DNA copy variation. Parenthetically, the approach utilized here provides a new way to sort and prioritize functionally relevant genes.

miRNAs are attractive candidates as upstream regulators of metastatic progression because they can regulate entire sets of genes. Therefore, we reasoned that miRNA deregulation may account for at least some of the gene expression changes that distinguish the highly metastatic LM cells from the poorly metastatic SN12C cells. We performed array-based miRNA profiling of SN12C, LM1, and LM2 cells. Hierarchical clustering based on the expression of the differentially expressed miRNAs correctly classified the cell populations into 2 groups (Fig. 5). Of 700 human miRNAs assayed, 38 were upregulated and 34 were downregulated in the highly metastatic cells, using a 1.5-fold threshold. On the assumption that the expression profiles of miRNA genes and their targets are inversely correlated, we identified their putative targets using TargetScan 4.2 and integrated the target genes with expression data from the same experiment. Using this methodology, we identified a network of putative functional miRNA-target regulatory relations involving 100 genes (Supplementary Fig. S7). We also validated the changes in expression of the 2 top miRNA, LM2 versus SN12C (Supplementary Fig. S8). A large proportion of genes included in the miRNAdown-mRNAup-predicted network have reported mesenchymal functions (CALD1, COL1A1, COL1A2, COL3A1, COL9A3, GREM1, and MYLK). Next, we sought to identify the underlying mechanisms that can explain the miRNA expression changes associated with the transition to a highly metastatic state in our model. miRNAs are frequently located in cancer-associated regions of the human genome (30). Thus, to determine whether genomic gains and losses might have impacted on miRNA expression, we correlated gains and losses of genomic regions, as determined by a-CGH on each cell line, with relative miRNA expression values. Notably, we found that only 2 miRNAs (hsa-miR-574-3p and hsa-miR-95) mapped to a region of gene-dosage change (Supplementary Fig. S9). In addition to genetic alterations, the expression of specific miRNAs could be dysregulated by epigenetic aberrations. To identify miRNAs with putative DNA methylation-related inactivation involved in metastasis, we treated SN12C and LM2 cell lines with 5′-AZA followed by hybridization to an expression microarray. By comparing the list of genes reexpressed by 5′-AZA and the list of genes differentially expressed in SN12C and LM2, we found that hsa-miR-224 and hsa-miR-34c-5p miRNAs undergo specific hypermethylation-associated silencing (Supplementary Fig. S10).

This study describes a comprehensive analysis of gene expression and genomic alterations in in vitro and xenograft models of RCC metastasis to elucidate the genetic mechanisms and associated biological functions that empower tumor cells with the ability to colonize distant sites. Molecular profiling of tumor cells in the xenograft model and RCC tumor biopsies uncovered a mesenchymal signature as the most significant feature linked to high metastatic activity. Unexpectedly, tumor cell invasion was not significantly increased in the highly metastatic populations. These results prompted us to search for alternative mechanisms by which a mesenchymal phenotype may be of advantage at the metastatic site. Detailed analysis of the gene functions associated with the mesenchymal signature revealed a substantial enrichment for profibrotic genes, suggesting that metastatic spread imposes a selective pressure for cells with a myofibroblastic makeup. We therefore suggest that an important function of tumor EMT is to generate cells with properties of stromal fibroblasts. It is logical to presume that a tumor cell endowed with such a phenotype would be highly effective at colonizing distant organs because it could bypass the requirement of an activated “compatible” stroma in the initial stages of metastatic establishment. Myofibroblasts produce and modify the extracellular matrix (ECM), secrete angiogenic and proinflammatory factors, and stimulate epithelial cell proliferation and invasion. Myofibroblasts were originally characterized by their role in wound healing (31) and its induction has been associated with diverse types of organ fibrosis (32). Myofibroblasts are also abundant in the reactive tumor stroma and are an established source of tumor promoting factors such as cell surface proteins, secreted growth factors, and ECM proteins (33, 34). A direct prometastatic effect of tumor-associated myofibroblasts was also recently implied (35). Consistent with our hypothesis, it has also been suggested that myofibroblasts are derived from malignant or normal epithelial cells undergoing EMT (36). Additional support for this hypothesis comes from the fact that increased levels of collagens and laminins have been associated with an increased likelihood of clinical metastasis of multiple human solid tumors (13). To substantiate the relevance of the mesenchymal signature, we selected S100A4 for functional validation. Enforced repression of this factor profoundly affected metastatic activity in agreement with previous studies that demonstrated that S100A4 is a crucial myofibroblast-expressed factor regulating metastasis (37). Furthermore, the expression of S100A4 in renal tubular epithelium undergoing EMT during fibrosis (38), together with the fact that proximal tubules arise from the differentiation of mesenchymal cells, suggest a model wherein RCC epithelial cells reverse this developmental process back to the original mesenchymal state, in order to metastasize. The plasticity of the renal epithelial cell has been further supported by studies of the role of the developmental gene, GREM1, which is reactivated in adult renal fibrotic disease (39). Remarkably, GREM1 is also a component of the myofibroblastic signature identified in our in vitro model of RCC progression.

Although deregulation of gene expression during tumor progression is relatively well characterized, the associated genetic and molecular mechanisms are largely unknown. In addition, although the sequence of genetic events that drive the establishment of primary tumors has begun to emerge, the subsequent events that lead to metastasis have remained generally obscure. As a consequence, few mutations, genomic alterations, or allelic imbalances are currently known to distinctively endow tumor cells with metastatic functions (40). We reasoned that by studying the mechanisms by which RCC genomes are deranged during metastatic conversion, we would be able to improve the ability to pinpoint critical genes involved in this process. By analyzing mRNA expression levels and DNA copy number changes in parallel, we found that changes in gene dosage have a modest effect on gene expression. These results suggest that during the evolution toward a metastatic phenotype, these cells acquire a set of chromosomal losses and gains that are perhaps associated with retention or otherwise of specific cell clones within the cell population. Although, none of the genes mapped to regions of genomic imbalance related to mesenchymal functions, we cannot exclude that these changes do not indirectly impact the mesenchymal signature.

In addition to changes in DNA sequence, gene expression can be modulated by aberrations in the patterns of DNA methylation. However, the pattern of DNA methylation abnormalities in cancer cells seems paradoxical. Compared with normal cells, cancer cells are concomitantly hypermethylated at specific CpG island sequences but hypomethylated at CpGs found in most other sites, resulting in a net loss of genomic 5meC content (41). Understanding the relative timing of DNA hypermethylation and hypomethylation alterations in cancer is crucial to understanding the importance of these changes during tumor progression. We found that a number of mesenchymal-associated genes were reactivated on 5′-AZA treatment suggesting that they were hypermethylated in the poorly metastatic cells. S100A4 showed complete methylation of the intronic CpG sites analyzed in the poorly metastatic cells, consistent with its epigenetic transcriptional silencing. Some of those genes, including S100A4, COL1A1, and GREM1, have previously been shown to be methylated in other tumor models and in RCC (42–44). The proposed impact of promoter methylation on EMT is not unprecedented in the literature, as altered methylation of some gene promoters have been reported to be one of the principal causes of EMT during tumor progression (45) or during epithelial to myofibroblast transition (46).

Several bioinformatic algorithms have been constructed to predict miRNA gene targets. These algorithms predict hundreds of potential gene targets, which cannot all be experimentally validated. Integration of putative targets with mRNA expression date provides a rational method to prioritize functionally relevant targets. Interestingly, the most significant correlation in this study was found between Let-7 and COL1A1, COL1A2, COL3A1, COL9A3. In consonance with our results, Let-7 is widely viewed as a tumor suppressor miRNA (47, 48).

Our findings, revealing the lack of major genetic and epigenetic alterations accounting for the direct deregulation of miRNA expression, suggest that miRNAs alterations are likely to be induced indirectly as a consequence of the dysregulation of specific transcription factors. Moreover, this observation argues that miRNA might be downstream targets of pathways that are commonly dysregulated in cancer and not initiating events during tumor progression.

In summary, in this study, we have proposed and provided proof of principle for a new mechanism of tumor progression in RCC based on acquisition of a myofibroblastic trait by tumor cells. We also showed that this approach that cross-references multiple whole-genome datasets can identify targets and genetic mechanisms important for tumor progression.

No potential conflicts of interest were disclosed.

We thank I. Fidler for providing the SN12C cell lines.

Grant Support This study was supported by grants from the V-Foundation, a Syms kidney cancer award, Pfizer, Inc., and the NIH (CA-121327).