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

**Ivana Samaržija 1,2**

<sup>1</sup> Laboratory for Epigenomics, Division of Molecular Medicine, Ruder Boškovi´c Institute, Bijeniˇcka 54, ¯ 10000 Zagreb, Croatia; Ivana.Samarzija@irb.hr

<sup>2</sup> Laboratory for Cell Biology and Signalling, Division of Molecular Biology, Ruder Boškovi´c Institute, ¯ Bijeniˇcka 54, 10000 Zagreb, Croatia

**Abstract:** 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.

**Keywords:** prostate cancer; bone metastasis; lymph node metastasis; liver metastasis; gene expression; meta-analysis; focal adhesion; protein filament; androgen receptor signaling

## **1. 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].

**Citation:** Samaržija, I. Site-Specific and Common Prostate Cancer Metastasis Genes as Suggested by Meta-Analysis of Gene Expression Data. *Life* **2021**, *11*, 636. https:// doi.org/10.3390/life11070636

Academic Editors: Ana Faustino and Paula A. Oliveira

Received: 15 May 2021 Accepted: 28 June 2021 Published: 30 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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.

#### **2. Materials and Methods**

#### *2.1. 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].

**Table 1.** Description of datasets used in this study.


#### *2.2. 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].

#### *2.3. Visualization of Metastasis Genes Networks and Identification of the Hub Genes*

The networks representing up-regulated and down-regulated genes among the 434 shared metastasis genes were retrieved by STRING [21] with default settings (default medium confidence 0.4) and visualized in Cytoscape [22]. The genes that showed interactions are depicted as a network, while the genes with no connections were omitted

from the Figures. Twenty hub genes from a set of 434 genes that are shared among all metastatic sites were detected with the use of cytoHubba application [23].

#### *2.4. Survival and Expression Analysis for the Hub Genes*

The prognostic significance of each hub gene was performed by gene expression profiling interactive analysis (GEPIA [24]) taking into account disease-free survival. *p* < 0.05 was considered to indicate a statistically significant difference.

The analysis of mRNA expression for the first 10 hub genes was performed using GEPIA software, including TCGA (The Cancer Genome Atlas) and GTEx (Genotype-Tissue Expression) data.

#### **3. Results**

*3.1. The Top Ten Most Changed Genes Shared by Metastasis from All Analyzed Sites Belong to the Class of Filaments and Proteins Involved in Bone and Prostate Biology*

As shown in Table 2, a substantial proportion (40 of 260) of genes that are among the 50 most up-regulated or 50 most down-regulated in prostate cancer bone, lymph node, and liver metastasis groups are shared at least between two groups (30 genes) or between all three groups (10 genes). It is interesting to note that all 40 genes that are found overlapping in two or three groups change in the same direction (up- or down-regulated), adding to the hypothesis that the change in their expression is functional, "driving" change rather than "passenger" change. The ten genes found in all groups are SPP1 (up-regulated) and MYH11, MSMB, ACTG2, CNN1, PCP4, KRT15, NEFH, DES, and CHRDL1 (down-regulated).

**Table 2.** 50 most up- (first row) and down-regulated (second row) genes in prostate cancer metastasis originating from bones, lymph nodes, and liver. The genes are ordered by decreasing |fold change|. Genes shared among three lists are depicted in red; genes shared among bone and LN metastasis are depicted in blue; genes shared among bone and liver metastasis are depicted in orange; genes shared among LN and liver metastasis are depicted in green.


According to the KEGG enrichment analysis of the nine shared down-regulated genes MYH11 and ACTG2 belong to the enriched process of vascular smooth muscle contraction (*p*-value = 0.0015). Additionally, five of those nine genes are also either filaments or they are involved in their processes, as listed below. CNN1 is a thin filament-associated protein

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].

#### *3.2. 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.

**Figure 1.** Venn diagram of differentially expressed genes in prostate cancer bone, lymph node, and liver metastasis in comparison to primary prostate tumors. The number in parenthesis indicates the total number of differentially expressed genes and the number in the intersection indicates the number of overlapping genes.

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.


**Table 3.** Gene ontology and pathway enrichment analysis of the genes common to all metastatic sites analyzed.


**Figure 2.** Network of up-regulated genes in prostate cancer metastasis from an intersection of bone, lymph node, and liver metastasis.


**Figure 3.** Network of down-regulated genes in prostate cancer metastasis from an intersection of bone, lymph node, and liver metastasis.

#### *3.3. 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).

#### *3.4. 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.


**Table 4.** Gene ontology and pathway enrichment analysis of the genes differentially expressed in bone metastasis.

**Table 5.** Gene ontology and pathway enrichment analysis of the genes differentially expressed in lymph node metastasis.


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".

**Table 6.** Gene ontology and pathway enrichment analysis of the genes differentially expressed in liver metastasis.


## *3.5. 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.

**Figure 4.** Hub genes that are found among 434 overlapping genes that are changed in metastases from all sites analyzed. The highest-ranked node is in red, and the lowest in yellow.

**Figure 5.** *Cont.*

1

**Figure 5.** Disease-free survival for hub genes as retrieved by GEPIA software using TCGA data.



2


**Table 8.** Gene ontology and pathway enrichment analysis of the hub genes.

**Figure 6.** The expression of the first 10 hub genes by rank in TCGA and GTEx datasets. Statistically significant differences are marked with an \*.

#### **4. 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–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.

#### **5. 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.

**Funding:** This research was funded by MY ZABA START 2019 donation from Zagrebaˇcka banka.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** The author would like to thank Robert Bakari´c and Andreja Ambriovi´c Ristov for their advice, support, and critical reading of the manuscript.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**

