Site-Specific and Common Prostate Cancer Metastasis Genes as Suggested by Meta-Analysis of Gene Expression Data

Anticancer therapies mainly target primary tumor growth and little attention is given to the events driving metastasis formation. Metastatic prostate cancer, in comparison to localized disease, has a much worse prognosis. In the work presented here, groups of genes that are common to prostate cancer metastatic cells from bones, lymph nodes, and liver and those that are site-specific were delineated. The purpose of the study was to dissect potential markers and targets of anticancer therapies considering the common characteristics and differences in transcriptional programs of metastatic cells from different secondary sites. To that end, a meta-analysis of gene expression data of prostate cancer datasets from the GEO database was conducted. Genes with differential expression in all metastatic sites analyzed belong to the class of filaments, focal adhesion, and androgen receptor signaling. Bone metastases undergo the largest transcriptional changes that are highly enriched for the term of the chemokine signaling pathway, while lymph node metastasis show perturbation in signaling cascades. Liver metastases change the expression of genes in a way that is reminiscent of processes that take place in the target organ. Survival analysis for the common hub genes revealed involvements in prostate cancer prognosis and suggested potential biomarkers.


Introduction
The vast majority (90%) of cancer-related deaths is not caused by primary tumors but metastasis on distal organs [1]. Still, cancer chemotherapies are mainly designed in a way that they consider events that drive primary tumor growth only, and little attention is given to pathways governing metastatic outgrowth [2]. While the literature and research on gene expression differences between primary tumors and metastasis is abundant [3], little is known on which signaling pathways are shared among metastatic cells colonizing distinct secondary sites, and which strategies are site-specific. Several publications exist on such topics, including the recent reports by Hartung et al. [4,5]. This type of knowledge is essential to suggest signaling pathways that could be therapeutic targets in metastatic cells in cancer types that spread to multiple organs.
Prostate cancer (PCa) is the second most commonly occurring cancer in men and the fifth leading cause of death worldwide [6]. The 5-year survival rate of non-metastatic prostate cancer is 98.9%, but, the rate in patients who are initially diagnosed with metastatic prostate cancer is less than 30% and does not show improvements [7]. The most usual sites colonized by prostate cancer cells are bones, lymph nodes, lungs, liver, and brain while rare locations include adrenal glands, breasts, eyes, kidneys, muscles, pancreas, salivary glands, and spleen. Recently, it was shown that patients with liver-only metastasis have worse cancer-specific and overall survival than patients with bone-only and lung-only metastasis [8]. The driving events in prostate cancer dissemination include entangled actions of several signaling pathways that are potentiated by changes in gene expression, genetic alterations [9], and post-translational modifications [10]. To better understand signaling events that are site-specific and common to prostate cancer metastasis, herein a metaanalysis of publicly available GEO gene expression datasets that consist of primary prostate cancer samples as well as samples of metastasis from bones, lymph nodes, and liver was conducted. Differentially expressed genes that are shared among all sites analyzed are presented. A substantial number of differentially expressed genes are secondary site-specific which emphasizes the need to study metastasis separately according to the secondary site. Survival analysis for hub genes found among the genes commonly changed in all metastatic sites was conducted and has revealed potential biomarkers.

Gene Expression Datasets
The list and description of gene expression datasets downloaded from the GEO database [11] and used in the meta-analysis are provided in Table 1. All the datasets are from microarray chips. The samples within datasets belong to four categories: primary tumors, bone metastasis (51 samples versus 175 primary tumor samples, from three datasets), lymph node metastasis (103 samples versus 232 primary tumor samples, from four datasets), and liver metastasis (26 samples versus 83 primary tumor samples, from two datasets). Depending on the platform, the annotation was either downloaded from the file deposited on the GEO database or obtained using gProfiler [12].

Meta-Analysis of Gene Expression Data and Enrichment Analysis
Meta-analysis of the gene expression data was performed using ImaGEO software that displayed good performance in the comprehensive gene expression meta-analysis [18]. The p-value method (minP) and default settings were used. This meta-analysis method was chosen as combining p-values provides an advantage over effect size combination for standardization of the associations from genomic studies to a common scale allowing to compare very heterogeneous datasets, for example, datasets from different tissues [18]. The criteria for differential gene expression were false discovery rate (FDR) p-value < 0.05 and |log2fold change| > 1. The intersection and the list of genes differentially expressed among bone, lymph node, and liver metastasis were obtained using GeneVenn [19]. Overview of enrichment analysis was obtained by Enrichr [20] using default settings and GO Biological Process (BP), GO Cellular Component (CC), GO Molecular Function (MF), KEGG (Kyoto Encyclopedia of Genes and Genomes), and WikiPathways data are presented. The top five processes are listed. The background used is set by default by Enrichr software, as Enrichr cannot upload a background list [20].
that is also implicated in the regulation and modulation of smooth muscle contraction. DES is a muscle-specific type III intermediate filament essential for proper muscular structure and function. NEFH is an intermediate filament protein, part of neurofilaments. KRT15 is a keratin that belongs to intermediate filament proteins responsible for the structural integrity of epithelial cells. PCP4 is a calmodulin regulator protein that functions as a modulator of calcium-binding by calmodulin. Among other roles, it was shown that Ca2+ and calmodulin regulate the binding of FLNA to actin filaments [25]. In summary, MYH11 and ACTG2 genes belong to genes involved in the assembly of actin fibers, while CNN1 and PCP4 are proteins capable of binding actin or influence its association with partner proteins. KRT15, NEFH, and DES belong to the three (keratins, neurofilaments, and desmin, respectively) of five classes of intermediate filaments.
SPP1 gene codes for osteopontin, the protein that is involved in the attachment of osteoclasts to the mineralized bone matrix. CHRDL1 protein has recently been shown to improve osteogenesis of bone marrow mesenchymal stem cells [26]. CNN1 gene also plays a role in osteoblast and osteoclast function and formation [27].
MSMB is one of the three major proteins secreted by the epithelial cells of the prostate. It is also secreted by epithelial cells in many other organs. The protein inhibits the growth of cancer cells in an experimental model of prostate cancer, but this property was shown to be cell line-specific [28].

Reorganization of Focal Adhesions and Changes in Androgen Receptor Signaling Are Common Characteristics of Prostate Cancer Metastasis Regardless of the Target Organ
As shown in Figure 1, the meta-analysis revealed that the intersection of differentially expressed genes among prostate cancer bone, lymph node, and liver metastasis consists of 434 genes. This gene list was analyzed to establish prostate cancer metastasis "core transcriptional program". Table 3 is listing the results of enrichment analysis for all 434 shared genes. It is clear from Table 3 that "Focal adhesion" is the most changed enrichment term in the intersection of metastasis from all sites as it is among the top enriched terms from the three lists-KEGG, WikiPathways, and GO CC. The up-regulated and down-regulated gene networks are depicted in Figures 2 and 3, respectively. The enrichment term "Prostate cancer" with 10/97 genes, "miRNA regulation of prostate cancer signaling pathways" with 8/33 genes, and "Androgen receptor" (12/90) are also among the most changed enrichment terms (Table 3) confirming the specificity of the results.

Results Suggest That Transcriptional Landscape of Bone Metastasis Is Profoundly Rewired in Comparison to Lymph Node and Liver Metastasis
In comparison to the number of genes differentially expressed in lymph nodes (2509) and liver (1269) metastasis, the changes in bone (7871 genes) metastasis are several times higher, suggesting the profound change in the transcriptional repertoire. Enrichment analysis was done for all the genes that are changed in bone metastasis (Table 4). On the top of the list of genes whose expression changes in bones is the term "VEGFA-VEGFR2 signaling pathway" followed by genes from focal adhesions and other signaling pathways. The BP category showed enrichment in genes involved in neutrophil biology. The top of the enrichment list of differentially expressed genes found in bone metastasis only (data not shown) is dominated by genes involved in immune system function, especially "Chemokine signaling pathway", "T-Cell antigen receptor signaling pathway" and "B-cell receptor signaling pathway".
The enrichment analysis of 100 most up-or down-regulated genes ( Table 2) in bone metastasis revealed enriched term of "Regulators of bone mineralization" (IBSP, COL4A1, and SPP1) being on the top of the list (BioCarta 2016).

Lymph Node Metastasis Are Characterized by Changes in Signaling Networks While the Liver Metastasis Transcriptional Program Is Reminiscent of Processes That Take Place in the Target Organ
On the top of the list of genes that are changed in lymph node metastasis (Table 5) is the term "VEGFA-VEGFR2 signaling pathway", followed by "Focal adhesion". This was found to be similar to the result of the genes differentially expressed in bone metastasis.
Other signaling pathways whose components are found to have significant differential expression in lymph node metastasis include PI3K-Akt, EGF/EGFR, MAPK, and TGF-beta signaling pathways.  On the top of the list of genes that are changed in liver metastasis (Table 6) is the term "VEGFA-VEGFR2 signaling pathway". When the enrichment analysis for the 100 most upor down-regulated genes (Table 2) in liver metastasis was done, the term "Complement and coagulation cascades" was found to be highly enriched with 13/79 genes being present on the list (data not shown). Genes specifically changed in liver metastasis only (data not shown) are reminiscent of processes taking place in the target organ according to the terms that are enriched and that include "Folate metabolism", "Selenium micronutrient network", "Fat digestion and absorption", "Cholesterol metabolism", "Phenylalanine metabolism" and fore-mentioned "Human complement system" and "Complement and coagulation cascades".

Prostate Cancer Metastasis Hub Genes and Their Involvement in Patient Disease-Free Survival
To detect the potential driving network of prostate cancer metastasis, hub genes among the 434 common genes of all metastatic sites were singled out and are shown in Figure 4. The disease-free survival analysis revealed the statistically significant involvement of the following genes ( Figure 5): AURKA, BUB1, CCNB2, CDC20, CDKN3, CENPF, CHEK1, FOXM1, HMMR, MELK, PTTG1, TOP2A, TPX2, TRIP13, TYMS, UBE2C. Their biological roles are listed in Table 7 and among the enriched terms listed in Table 8 are "Cell-cycle" and "Microtubule cytoskeleton". The analysis of mRNA expression for the first ten hub genes is shown in Figure 6, confirming their up-regulation in prostate cancer.

Discussion
Metastasis formation is a complex process driven by a variety of genes and molecular events [3]. Meta-analysis on differential gene expression from primary tumors and metastasis are common for different cancer types. However, the analyses that stratify metastatic samples according to the secondary sites are rare with more interest shown in very recent years [29][30][31][32][33]. To the best of my knowledge, the work that is presented here is the first such study for prostate cancer and its most common metastasis sites. These types of studies are important as they accumulate information that could be of interest when designing anti-metastatic therapies.
The main finding of the presented study is that a group of differentially expressed genes encoding for filaments or associated proteins is the most differentially expressed by fold change and the processes of focal adhesion and androgen receptor signaling are among the most changed in metastasis from all sites analyzed. Moreover, there is a substantial difference in expression programs from metastasis from different sites. In the following chapters, based on the results, several questions are considered: what are the site-specific transcriptional programs that predispose or characterize the metastasis from bone, lymph node, and liver and could potentially be used as targets to treat metastasis from those particular sites; what are the common genes that could be used as targets to treat metastasis from all sites; what are the strategies that are used by cancer cells to colonize different organs?
Although the lymph nodes are the first sites encountered by prostate cancer cells that enter the circulation system, the bones are the most common sites that are homing them [34]. This is very intriguing as this study shows that, among the three most common distal sites, metastasis found in bones underwent a much more profound change of transcriptional program (7871 genes with changed expression) in comparison to metastasis from lymph nodes (2509) and liver (1269). The question arises as to what facilitates the colonization of the bones when, despite the complete reorganization of the transcriptional program that is needed, they are still the first choice for prostate cancer cell homing? The suggested concept of pre-metastatic niche offers the explanation that the primary cancer cells, possibly through exosomes, prime the distal organs which increases chances and enables their homing [35]. However, the question is whether there is a predisposition already within primary prostate cancer cells and their transcriptional program that is kept and subsequently upgraded in metastatic cells that makes them prone to bone colonization ("seed pre-selection" [36]). From a list of differentially expressed genes, the driving events in bone-specific metastases that could be directing them to the target organ are extracted. Upregulation of SPP1 and downregulation of CHRDL1 and CNN1 (genes involved in bone biology) in metastasis from all sites analyzed are found which suggests they are changed early on and could belong to the genes that make prostate cancer metastatic cells bone-gravitating. The expression level of SPP1 is elevated in cancers, particularly those that spread preferentially to the skeleton. "Osteomimicry" of malignant cells is partially conferred by transmembrane receptors bound by SPP1. Binding of integrins on malignant cells by SPP1 results in activation of signaling cascades within the cell that promotes metastasis [37]. CHRDL1 gene encodes an antagonist of bone morphogenetic protein 4 and may play a role in embryonic bone formation, while overexpression of CNN1 in osteoblast lineage cells was shown to regulate bone mass [27]. Because of their prominent role in bone biology as listed above, these genes could contribute to the site-preference during the metastatic process. Also, an extensive change in transcription of immune cell-related genes (chemokines, T-and B-cells, and neutrophil-related genes) was recorded in genes differentially expressed in bone metastasis. This transformation supports the role of prostate cancer cell and immune system crosstalk which is crucial for the formation of metastasis in this organ [38]. The contribution of chemokines to the metastatic process has been well documented [39,40]. For example, the prostate cancer metastasis-promoting role of CXCL5 has been recently shown [41]. Some other chemokines from the list of differentially expressed genes (e.g CCL5, CCR2, reviewed in [40]) are also implicated in the metastatic process from different cancer types. On the list of differentially expressed chemokines in bone metastasis (data not shown), 13 out of 14 chemokines or chemokine receptors with changed expression in bone metastasis only, are up-regulated, indicating a strong increase of activity in the network of chemokines with known (minor part) and yet to be investigated (major part) roles in prostate cancer metastasis. This finding suggests that targeting the chemokine signaling pathway could alleviate the bone metastasis burden in prostate cancer patients.
As noted above, to colonize lymph nodes and liver, prostate cancer cells undergo changes in transcription that are several times less extensive in the number of affected genes than in bone metastasis. From this data, it could be suggested that the brute force that drives prostate cancer cells through circulation is probably more involved in regional lymph node and liver colonization than in bones. However, enrichment analysis of genes whose expression is changed in lymph nodes only revealed extensive changes in the expression of kinases, suggesting rewiring of signaling cascades. While the term "Focal adhesion" was on the top of the enrichment lists of genes shared among all metastasis, this list was extended to 60 genes, including 7 integrins, and was on the top of the list shared among bone and lymph node but not liver metastasis (data not shown). This indicates that bone and lymph node metastasis highly rewire the expression of adhesion molecules and related pathways, while for liver metastasis this change is much less prominent. Taken together, these results suggest that there are parts of focal adhesion (integrin) network that are commonly changed in prostate cancer metastasis from all analyzed sites, but also a part that is specific to metastasis from each site, giving them an integrin code that is, apparently, important during every step of the metastatic cascade [42] and a subject to a frequent change [43].
The observation from this work which indicates that the expression of genes encoding for microfilaments and intermediate filaments or associated proteins (MYH11, ACTG2, KRT15, NEFH, DES, CNN1, and PCP4) are the most extensively down-regulated in prostate cancer metastasis from all sites suggests fundamental reorganization of their cytoskeleton. KRT15 is downregulated in the progression of normal prostate tissue to prostate cancer and further to lymph node metastasis [44]. Other mentioned genes influence tumors either by inhibiting [45,46] or promoting their pathogenicity [47,48] or displaying a dual role [49,50]. Those studies investigated the roles of individual proteins, however, in this system, it is very likely that they act in concert which could be envisioned from their involvement in the same process of cytoskeleton assembly.
It is interesting to note that liver metastases are highly enriched for genes that are involved in the processes of a target organ. This phenomenon of metastatic cells adapting to the host site has been recently described [5], but the potential contamination with the host tissue should also be taken into account.
The role of androgen receptor signaling in prostate cancer progression is multilayered and has been extensively studied [51,52]. The change in androgen receptor signaling pathway (up-regulation of androgen receptor and changes in expression of related genes) that are documented here are in agreement with recent findings that elevation of androgen receptor promotes prostate cancer metastasis by induction of epithelial-mesenchymal transition [53]. Also, this result is in agreement with a study by Guo et al. [54] who suggest androgen receptor as one of the hub genes in metastatic prostate cancer. However, the study presented here encompasses a larger sample size and presents a more extended approach. In addition, herein differences in differential gene expression between different secondary sites are in focus which is, to the best of my knowledge, a unique approach for prostate cancer metastasis and not so common for other cancer types.
In this study, 20 genes were singled out as hub genes among those that change expression in all metastatic sites analyzed. For most of these genes involvement in diseasefree survival was shown suggesting that these genes might be considered in potential targeted therapies.
Finally, heterogeneity of prostate cancer calls for studies that include an even larger number of samples than this study did. The experimental validation of these results would bring the confirmation of the importance of the genes that are suggested here to play a role in the prostate cancer metastatic process. Also, this analysis involved only the data for protein-coding genes and their subsequent mRNA expression. However, the roles of non-coding RNAs are also known to be highly important in prostate cancer progression, both as regulatory elements and biomarkers [55,56]. It would be interesting to reveal their roles in prostate cancer metastasis from different secondary sites.

Conclusions
Although sharing changes in the expression of basic groups of genes belonging to the class of filaments and focal adhesions and androgen signaling pathway, metastasis from different sites differ profoundly in their transcriptional program. Based on this finding, it can be concluded that it is important to study separately cancer cells originating from different secondary sites because the results of gene expression data are expected to be skewed when metastasis samples are not stratified. In addition, data on the differentially expressed genes that are site-specific and common to all metastases provide potentially useful information for targeting metastatic cells. It would be interesting to see whether metastatic cells from different primary organs use similar strategies when colonizing the same secondary site. AURKA, BUB1, CCNB2, CDC20, CDKN3, CENPF, CHEK1, FOXM1, HMMR, MELK, PTTG1, TOP2A, TPX2, TRIP13, TYMS, UBE2C are the hub genes identified in this study that show involvement in the disease-free survival of prostate cancer patients.