Prostate Cancer-Induced Remodeling of the Topologies and Interplay of the TLR (Toll-Like Receptor) and Chemokine Signaling Pathways

1. Introduction

Toll-like receptors are a key family of pattern recognition receptors in the innate immune system and play a fundamental role in host defenses [1]. The activation of TLR signaling pathways prevents not only pathogen infection by recognizing pathogen-associated molecular patterns (PAMPs) [2] but also cancer by recognizing dangerous endogenous damage-associated molecular patterns (DAMPs) [3]. Once TLRs recognize their ligands, an intracellular signaling cascade is triggered to activate master transcriptional factors including NF-κB and IRFs to induce pro-inflammatory cytokines, upregulate costimulatory molecules, and secrete interferons [4]. These activations control tumor immune microenvironments and affect T cell anti-tumor activity.

Prostate cancer (PCa) is the most common cancer in men worldwide [5]. Althoughmost PCa in the early stages can be managed by hormone therapies, some of them inevitably will become hormone deprivation-resistant, metastatic, and deadly, a condition known as metastatic castration-resistant PCa (mCRPC). Patients with mCRPC often fail to respond to current treatments, including cancer immunotherapy. We recently found that prostate cancer cells express a repertoire of TLRs and activation of TLR5 by the bacterial protein flagellin, a known TLR5 agonist, induces robust anti-tumor immune responses along with NF-κB activation in prostate cancer [6].

This study complements two previous publications [7,8] with the analysis of the cancer-induced remodeling of the Toll-like receptor signaling pathway (hereafter denoted TLR) and its interplay with the chemokine signaling pathway (CHS) in the three cancer nodules harbored by a surgically removed prostate tumor.

Such analysis provides insight into the transcriptomic behavior of the prostate cancer cells and innate and adaptive immunities. Understanding the remodeling of the TLR might provide a new strategy for developing novel personalized gene therapies for this malignancy [9], including the use of TLR agonists [10].

2. Materials and Methods

2.1. PROSTATE tissues

We re-analyzed the gene expression profiles of three cancer nodules and the normal surrounding tissue of a surgically removed prostate tumor from a 65y old black man deposited (together with the experimental protocols) on the publicly accessible NCBI Gene Expression Omnibus [11,12,13]. The three cancer nodules were denoted as “A”, “B”, and “C”; “A” is the primary cancer nodule (Gleason Score GS = 4 + 5 = 9), while nodules “B” and “C” (each with GS = 4 + 4 = 8). Point (~ 2mm) biopsies from each cancer nodule and surrounding normal tissue (denoted as “N”) were split into four biological replicas and the total RNA extracted and hybridized with Agilent-026652 Whole Human Genome Microarray 4x44K v2 [14].

2.2. Transcriptomic Characterization of Individual Genes (Appendix A)

The expression levels of adequately quantified genes (by spots with no corrupt pixels and foreground fluorescence larger than twice the background) were normalized to that of the median gene of that microarray. Each gene “i” was characterized in every prostate region “c” by the primary independent variables: AVE_i^(c), REV _i^(c) and COR _ij^(c), that are mathematically defined in Appendix A for each condition/region “c”.

AVE_i^(c) (formula A1) is the average normalized expression level over biological replicas normalized to the median gene expression.

REV _i^(c) (formula A2) is the mid-interval of the chi-square estimate of the expression coefficient of variation of the normalized to the median expression levels of a gene among the biological replicas of a given region.

COR _ij^(c) (formula A3) is the correlation coefficient of the (log₂) normalized to the median expression levels of genes “i” and “j” in the region “c”. As proved recently [15], given the technical noise of the microarray technology, the configuration function of the transcriptome is suitable for approximating the set of all possible pair-wise correlations among gene expressions in biological replicas.

There are three statistically significant types of expression correlation:

i. 1 ≥ COR _ij^(c) ≥ 0.95, synergistic expression (i.e. expressions of both genes fluctuate synchronously across biological replicas);

ii. -0.95 ≥ COR _ij^(c) ≥ -1, antagonistic expression (i.e. expressions of the two genes fluctuate in antiphase (when is up the other is down) across biological replicas);

iii. |COR _ij^(c)| <0.05, independent expressions (i.e. expression change of one gene has no consequences for the other one).

We used also the derived characteristics Relative Control Strength “RCS” (A4), Relative Expression Control “REC” (A5), Coordination degree “COORD”, and Gene Commanding Height “GCH” (A6).

RCS indicates how stable is that gene expressed among the biological replicas of the profiled region with respect to the median gene, the log₂ ratio in the definition of REC dividing the genes into equally two subsets with symmetric positive and negative values. High positive RECs point out genes whose expression control is very strict, most likely because their right expression levels are critical for the cell survival and phenotypic expression. By contrast, the most negative RECs pinpoint genes whose expressions may loosely fluctuate, as they might be vectors of adaptation to the changing environment.

COORD = percentage of synergistically expressed gene pairs + percentage of antagonistically expressed pairs - percentage of independently expressed genes within the analyzed pathways. Together, the percentage of synergistically expressed gene pairs and the percentage of antagonistically expressed pairs indicate how much coupled (i.e. influencing each-other the expression) are the genes in the functional pathway

GCH combines the expression control and the expression coordination with all other genes and serves to identify the Gene Master Regulators (GMRs, top GCH) of the entire transcriptome or just of selected functional pathways in the profiled region.

2.3. Transcriptomic Characterization of Functional Pathways and Their Interplay

In this study, we resume the pathway analyses from previous papers ([7,8]) with the KEGG [16]-derived signaling pathways CHS ([17]) and TLR ([18]), and their transcriptomic interplay through the shared genes. All transcriptomic characteristics of individual genes were after that ranked and averaged over the pathway(s) they belong.

2.4. Quantification of Transcriptomic Regulation (Appendix A)

A gene is significantly regulated in the cancer nodule (A, B or C) with respect to the normal tissue (N) when it satisfies the composite criterion (A7). It requires that |x| (the absolute fold-change of the expression ratio x (A8), negative for down-regulation) exceeds the corresponding cut-off value (“CUT”) and the p-value of the heteroscedastic t-test for the means’ equality is less than 0.05. CUT (A9), computed for each gene in the two-region comparison, accounts for both biological variability across the replicas and technical noise of the probing spots. Fold-change formulas (A10 & A11) quantify the changes of REC and GCH.

3. Results

135 out of 191 CHS and 76 out of 109 TLR KEGG-identified genes were profiled and analyzed in this report. Transcripts of missing genes were either not probed by the used Agilent microarray platform, or amount below the detectable limit (i.e. yielding a foreground fluorescence signal less than twice the background in the hybridized spot), or were hybridized to microarray spots with either corrupted or saturated pixels in at least one sample. The analyzed pathways share 33 genes, out of which we quantified 28: AKT1/2/3, CCL3/3L3, CCL4L2, CCL5. CHUCK, CXCL9/10/11, IKBKB/G, MAP2K1, MAPK1/10, NFKB1, NFKBIA, PIK3CA/B/D, PIK3R1/2/3, RAC1, RELA, STAT1/2.

3.1. The three CHS and TLR genes with the highest expression level (largest AVE)

We began by identifying the top three genes with the highest average expression level (AVE) in three prostate cancer nodules (A, B, C) and the surrounding normal tissue (N) (Table 1, expression ratios in Table S1 from the Supplementary Materials).

In the surrounding normal tissue “N”, the top three CHS genes were GNAI2, RAC3, and STAT2 (common to both pathways). Interestingly, GNAI2 and STAT2 were also expressed at levels close to the normal one (no significant regulation) in the cancer nodules, suggesting their important roles in the immune response in the entire prostate tissue. However, RAC3 expression was lower in all three cancer nodules compared to the normal tissue (expression ratios: x^{(A vs N)} = -2.11; x^{(B vs N)} = -1.51, x^{(C vs N)} = -1.35), indicating that RAC3 down-regulation may contribute to decreased cell proliferation, migration, and invasion in prostate cancer [19].

Among the CHS genes, CXCL14 (positive fold-changes: x^{(A vs N)} = 6.69; x^{(B vs N)} = 7.17, x^{(C vs N)} = 4.58) and NFKBIA (x^{(A vs N)} = 2.23; x^{(B vs N)} = 2.90, x^{(C vs N)} = 1.60) had high expression levels in all three cancer nodules but not in the surrounding normal tissue. The up-regulation of these genes underscores their important roles in prostate cancer development. In the analysis of TLR genes, IRF7 was also highly expressed in all three cancer nodules (although twice more in B than in C) but not in the surrounding normal tissue (expression ratios: x^{(A vs N)} = 3.02; x^{(B vs N)} = 4.15, x^{(C vs N)} = 2.05).

3.2. The Three Most Controlled CHS and TLR Genes in Each Profiled Region

The REC analysis identifies the most critical CHS and TLR genes for cell survival and phenotypic expression in the cancer nodules “A”, “B”, “C”, and the normal tissue “N”. Table 2 presents the REC values (formula A5) of the top three genes in each profiled region. The fold-changes (negative for down-regulation, formula A10) of the Relative Control Strength (formula A4) in the cancer nodules with respect to “N” are presented in Table S2 from the Supplementary Materials.

Within the CHS pathway, the most controlled genes (highest positive REC values) in the cancer nodules were: GNB1 (“A”), MAPK1 (“B”) and CCL16 (“C”). According to the FC scores, each of these genes had an impressive increase of the control strengh with respect to the normal tissue in one nodule but not in the other two. Thus, GNB1 had:

FC^{(A vs N)} = 23.99, FC^{(B vs N)} = 1.38, FC^{(C vs N)} = 2.04, MAPK1 had: FC^{(A vs N)} = 2.14, FC^{(B vs N)} = 9.92, FC^{(C vs N)} = 2.07, while CCL16 had: FC^{(A vs N)} = -1.11, FC^{(B vs N)} = 2.08, FC^{(C vs N)} = 9.53.

The situation is similar within the TLR pathway where the most controlled gene in one cancer nodule had lower fold-change increases in the other nodules. Thus, TYK2 had FC^{(A vs N)} = 27.98, FC^{(B vs N)} = 6.05, FC^{(C vs N)} = 1.11, TOLLIP had: FC^{(A vs N)} = 5.07, FC^{(B vs N)} = 13.19, FC^{(C vs N)} = -1.10, more control in “B” but even less control in “C” than in “N”), and IRF5 had: FC^{(A vs N)} = 2.19, FC^{(B vs N)} = 1.45, FC^{(C vs N)} = 6.82.

On the opposite side, the control strength of the top three CHS genes in N exhibited a substantial reduction in each of the three cancer nodules. Thus, the control strength of the top CHS gene, CCL15, was reduced by: FC^{(A vs N)} = -25.22, FC^{(B vs N)} = -5.21, FC^{(C vs N)} = -6.29. The control stength of VAV2 was reduced by: FC^{(A vs N)} = -8.23, FC^{(B vs N)} = -7.91, that of KRAS by FC^{(A vs N)} = -3.29, FC^{(B vs N)} = -6.66 and FC^{(C vs N)} = -2.60. Although less impressive, the control strength of the top TLR genes was also reduced in the cancer nodules. Thus, PIK3R1 had FC^{(A vs N)} = -3.49, TBK1 had FC^{(B vs N)} = -3.12, FC^{(C vs N)} = -2.45, PIK3CA had: FC^{(A vs N)} = -1.60, FC^{(B vs N)} = -2.29, FC^{(C vs N)} = -3.04.

The most important conclusion of this analysis is that the homeostatic expression control mechanisms are different among the cancer nodules, even equally pathologicaly graded and harbored by the same tumor.

3.3. The Three Least Controlled CHS and TLR Genes in Each Profiled Region

The REC analysis identified also the CHS and TLR genes whose loose expression control (largest negative REC values) make them suitable as adaptation vectors to the environmental change. Table 3 presents the least controlled three genes in each profiled region, while Table S3 from the Supplementary Materials presents the fold-changes of the control strength with respect to “N”.

Interestingly, genes such as CCL15 jumped from being very strictly controlled in the normal tissue (REC = 2.75) to being very loosely controlled in nodule C (REC = -1.91), meaning an over 25x increase in the Relative Expression Variation (REV). CCL19 was the least controlled in the equally graded nodules “B” and “C”, but not in “A” and even less controlled in “N”. Within the TLR pathway, the least controlled genes in nodules “A” (TICAM2) and “B” (CXCL11) had also relaxed control in the normal tissue.

It is important to note that in the primary cancer nodule A, the expression control of TLR4 is substantially reduced (by 4.14x, i.e. by 314%). However, the control in the equally graded nodules is increased by 1.40x (i.e. 40%) in B and by 2.28x (i.e. 128%) in C (Table S3). These results indicate major switches in the cell’s priorities in controlling the expression levels of critical genes.

Since the least controlled genes might be used by the cell as adaptation vectors to accommodate various local environmental conditions, these results indicate that the transcriptomic molecular mashinery might be adjusted in various ways to respond the environmental challenges.

3.4. Regulation of the TLR-Signaling Pathway In The Cancer Nodules Compared To The Surrounding Normal Tissue

Figure 1 presents the statistically significant regulation of the genes included in the KEGG-determined Toll-like receptor signaling pathway [16] in the cancer nodules “A”, “B” and “C” with respect to the surrounding normal tissue “N”.

Overall, the common features among the three nodules are the upregulation of the cell surface TLRs (TLR1/2/6 and TLR4) and of the associated CD14. Furthermore, these cell surface TLRs activate via intracellular signaling cascades the transcriptional factors NF-kB and IFNs that trigger the response of the innate and then adaptive immune system. Of note is that even the equally graded cancer nodules “B” and “C” from the same tumor exhibit differences in the gene expression regulation, proving the tumor transcriptomic heterogeneity. Thus, out of the 76 quantified genes, 17 (22.37%) were up- and 10 (13.16%) were down-regulated in “B”, while 13 (17.11%) were up- and 4 (5.26%) were down-reguled in “C”. Nodule “A” prensented 20 (26.32%) up- and 9 (11.84%) down-regulated genes. Although the significat reglation profiles were different among the three cancer nodules, no TLR gene was oppositelly regulated in two nodes.

3.5. Regulation of the Chemokine Signaling (CHS) Pathway

Figure 2 presents the statistically significant regulation of the genes included in the KEGG-determined Chemokine signaling pathway [17] in the cancer nodules “A”, “B” and “C” with respect to the surrounding normal tissue “N”. Of note is that even the equally graded cancer nodules “B” and “C” from the same tumor exhibit differences in the gene expression regulation, proving the tumor transcriptomic heterogeneity. Thus, out of the 135 quantified genes in this pathway, 24 (17.78%) were up- and 21 (15.56%) were down-regulated in “B”, while 25 (18.52%) were up- and 8 (5.93%) were down-regulated in “C”. Nodule “A” presented 32 (23.70%) up- and 20 (14.81%) down-regulated genes. Although the significant regulation profiles were different among the three cancer nodules, no CHS gene was oppositely regulated in two nodes. A number of 28 genes were included in both TLR and CHS signaling pathways: AKT1/2/3, CCL3, CCL3L3, CCL4L2, CCL5, CHUK, CXCL9/10/11, IBKB/G, MAP2K1, MAPK1/10, NFKB1, NFKBIA, PIK3CA/B/D, PIK3R1/2/3, RAC1, RELA, STAT1/2. We considered these common genes as connectors between the two pathways.

Yellow background of the gene symbols in Figure 1 and Figure 2 does not mean that expression level did not change in the cancer nodules with respect to the normal tissue but that the change was not statistically significant from the perspective of our composite criterion (A7).

3.6. Remodeling of the TLR-CHS Interplay

3.6.1. Significant Expression Coordination

Figure 3 presents the (p < 0.05) statistically significant positive (i.e. synergistic) and negative (antagonistic) expression coordination of the TLR-signaling pathway genes with the CHS signaling pathway genes through the 28 shared genes in the surrounding normal prostate tissue “N” (a) and the three cancer nodules. Of note are the significant increases of the coordination scores (COORD) from 1.09% and 2.50% in “N” to 8.97% and 13.51% in “A”, 4.52% and 6.79% in “B”, and 3.43% and 6.73% in “C”. It is also important to observe the substantial differences among the cancer nodules collected from the same tumor even when they were equally graded. The coordination differences indicate that cancer-triggered remodeling of the TLR-CHS transcriptomic interplay is dependent on the local, never exactly repeating conditions.

3.6.2. Significant Expression Independence

Figure 4 presents the (p < 0.05) statistically significant expression independence of the TLR and CHS genes with their 28 common hubs in the normal tissue “N” (a) and the three cancer nodules. Significant expression independence means that there is no functional pathway (including, but not limited to TRL and CHS signaling ones) where the encoded proteins by the two genes might directly interact.

We found substantial decreases of the percentages of independently expressed genes, from 5.90% and 4,94% in “N” to 2.84% and 2.02% in “A”, 3.64% and 2.80% in “B”, and 3.43% and 2.80% in “C”. The percentage of the independently expressed gene pairs has negative contribution to the overall coordination score “COORD” provided in Figure 3 for the transcriptomic interplay of the CHS and TLR signaling pathways through the shared genes.

3.6.3. Prostate Cancer Couples Several Normally Independently Expressed TLR and CHS genes

Figure 5 presents the (p < 0.05) statistically significant independently expressed TLR-CHS gene pairs in the normal tissue that were turned into (p < 0.05) statistically significant synergistically or antagonistically expressed genes. The 28 genes shared by the two pathways were excluded from this analysis and only the pairs with significant changes in their expression coordination are presented. Of note is that cancer made genes that are normally independently expressed to coordinate their expression fluctuations, indicating deep remodeling of the molecular mechanisms. It is important to observe the inter-nodule differences in the significant gene expression coupling, even between equally graded nodules from the same tumor. Thus, in nodule “A”, 11 independently expressed gene-pairs were turned into synergistically and 11 into antagonistically expressed, in “B”, 6 were turned synergistically and 4 antagonistically, while in “C”, 10 were turned synergistically and 7 antagonistically expressed.

3.6.4. Prostate Cancer Decouples Several Coordinately Expressed TLR and CHS Gene-Pairs in the Normal Tissue

We found that cancer not only couples genes that are normally independenly expressed but also that it decouples genes whose expression is normally significantly synergistically or antagonistically coordinated. Figure 6 presents part of the significantly coordinatelly expressed TLR-CHS gene-pairs (shared genes excluded) in the normal tissue that were decoupled in the cancer nodules. Of note are again the substantial differences among the three cancer nodules with 6 coordinatelly expressed gene-pairs turned into independently expressed in “A”, 8 in “B” and 5 in “C”.

3.6.5. Prostate Cancer Switches to the Opposite the Normal Coupling of Several Gene Pairs

We found that in addition to coupling some normally independently expressed genes or decoupling genes whose expressions are significantly correlated, cancer also switches the normal type of expression coordination. Figure 7 presents part of the TLR-CHS gene pairs whose normal significant synergistic or antagonistic expression was switched to the opposite by cancer. Of note are again the differences among the three cancer nodules. Thus, “A” turned 6 synergistic pairs into antagonistic and 3 antagonistic into synergistic, “B” turned 3 synergistic into antagonistic and 3 antagonistic into synergistic, and “C” turned 5 synergistic into antagonistic and only one antagonistic into synergistic.

3.7. Key Influential Genes in the CHS and TLR Signaling Pathways

In addition to being strctly controlled because its right level is critical for the cell physiology, expression of a key gene is also a major regulator of the transcription of many other genes. Analysis from the GCH score perspective (definition formula A6) revealed distinct patterns of influential genes across the three prostate cancer nodules and the surrounding normal tissue “N”, highlighting the heterogeneous nature of the disease. In nodule “A”, the top three influential genes were GNB1, JAK2, and TYK2, while nodule “B” showed a predominance of CCR6, CCR1, and BCAR1. Nodule “C” exhibited a different pattern with CCL16, PLCB1, and PREX1 as the most influential genes. The surrounding normal tissue “N” displayed yet another distinct profile with CCL15, PRKCD, and VAV2 emerging as the top influential genes (Table 4). Table S4 in the Supplementary material gives the fold-changes (negative for down-regulation) of the GCH scores in the cancer nodules with respect to “N”.

4. Discussion

Although limited to three cancer nodules and the surrounding normal tissue from a surgically removed prostate tumor, our all-inclusive transcriptomic analysis provides new insights into the dynamic regulation of the toll-like receptors and chemokine signaling pathways. The investigation was carried out from the perspective of the Genomic Fabric Perspective that provides the most theoretically possible comprehensive characterization of the transcriptome [20].

The observed heterogeneity of the gene expression profile, control and networking across the cancer nodules of the tumor highlights the complexity of prostate cancer biology. These results may inform future personalized therapeutic strategies aimed at targeting specific signaling pathways that are dysregulated in prostate cancer, with the potential to develop more personalized and effective treatments.

The present report confirms our prior observations regarding the unrepeatability of the cancer-induced remodeling of the transcriptomic topology not only from person to person but also among the cancer nodules (even those equally graded) from the same tumor. The transcriptomic uniqueness of the investigated prostate cancer nodules was before disclosed for the P53-signaling, apoptosis, block of differentiation, evading apoptosis, immortality, insensitivity to anti-growth signals, proliferation, resistance to chemotherapy, and sustained angiogenesis [7]. A similar conclusion resulted from the analysis of: mTOR signaling, metabolic genes and expression coordination of AKT2 with its partners central to the KEGG-constructed prostate cancer pathway [21] in another prostate tumor [8]. The standard immortalized prostate cancer cell lines LNCaP [22,23] and DU145 [24,25] presented also substantial differences in the transcriptomic organization. The singleness of the cancer nodules seems to be not limited to the prostate tumors but extended to all poly-clonal malignancies. Meaningful inter-nodule differences were found through analysis of: chemokine and VEGF signaling, apoptosis, basal transcription factors, cell cycle, oxidative phosphorylation, and RNA polymerase pathways in a case of metastatic clear cell renal cell carcinoma [15,26].

The observation that GNAI2 maintained the highest expression levels across tumor (normal and cancerous) regions (Table 1) suggests that the malignancy does not alter the fundamental role of the encoded protein (Giα2) in the chemokine signaling. Although not statistically significant, GNAI2 slight over-expression in all nodules: x^{(A vs N)} = 1.19, x^{(B vs N)} = 1.65, x^{(C vs N)} = 1.18 (Table S1) justifies why Giα2 is essential also for the prostate cancer cell migration [27]. An interesting paper [28] has shown that the GNAI2 expression decreased in the early stage of ovarian cancer but increased compared to normal in advanced cancers.

The differential expression of several key genes between normal and cancerous tissues provides insight into potential cancer-specific mechanisms. For instance, the downregulation of RAC3 in cancer nodules, contrary to its high expression in normal tissue, indicates a possible tumor-suppressive role that contradicts some previous findings in other cancer types [19]. This discrepancy might reflect tissue-specific functions or temporal dynamics of cancer progression [29].

The upregulation of CXCL14 in all cancer nodules aligns with previous findings [30], while the elevated expression of NFKBIA (Table 1) is somehow surprising, suggesting a more complex regulation of CHS and TLR signaling pathways in prostate cancer. NFKBIA encodes IκBα, which binds and inhibits the NF-κB complex [30]. This could indicate an intricate feedback mechanism attempting to maintain cellular homeostasis in the face of persistent inflammatory signaling. The up-regulation of IRF7 (Table S1), a critical transcription factor in type I interferon responses, is still not enough explored in the context of potential gene therapy of prostate cancer [31].

The results from Table 2 listing the most strictly controlled genes indicate that the phenotypic expression of histo-pathologically distinct regions in the same tumor depends on most likely different sets of critical genes, even for equally pathologically graded nodules. GNB1, shown to promote malignancy in the colon [32], was the most controlled gene in “A” (but not in “B” and “C”) while MAPK1, a valuable target in anti-cancer therapy [33], was the most controlled in “B” (but not in “A” and “C”). Effectiveness of intra-tumoral injection of CCL16 [34], the most controlled in “C” (but not in “A” and “B”) might be explained by its significant downregulation in prostate cancer. The control of GRK5, known for its regulating power in prostate cancer metastases [35], increased by 6.62x in “B” with respect to “N” but stayed practically the same in “C” (FC^{(C vs N)} = 1.12).

Table 3 shows that the most likely adaptation genes are also different among the cancer nodules. Our analysis revealed dramatic shifts in gene expression control between normal and cancerous tissues. The observation that CCL15 transformed from being strictly controlled in normal tissue to loosely controlled in cancer nodules suggests a significant disruption of normal regulatory immune response mechanisms [36,37]. Conversely, we found a spectacular increase of the expression control of TYK2 (FC^{(A vs N)} = 27.98) and GNB1 (FC^{(A vs N)} = 23.99) in nodule “A” (but to a much lesser extent in “B” and “C”). These results indicate that each cancer clone selectively enforce strict regulation of certain signaling components [35,38] while allowing others to vary widely.

The identification of 28 shared genes between TLR and CHS signaling pathways provides insight into the interconnected nature of these inflammatory response networks [39,40,41,42]. The consistent upregulation of cell surface TLRs (TLR1/2/6 and TLR4) and CD14 across all three nodules suggests a common mechanism for enhanced immune surveillance or inflammation in prostate cancer. This is further supported by the activation of downstream NF-κB and interferon signaling pathways [43,44,45,46,47].

Given the prominent role of TLR4 in immune system signaling [48] and its potential use as a prostate cancer biomarker [49], it is puzzling to note its loosed control in nodule “C” while its control is even stricter in the other two nodules with respect to “N”. This is another indication of the dependency of cellular homeostasis on the local environmental conditions.

One of the most striking observations is the evidence for substantial tumor heterogeneity, even between equally graded nodules from the same tumor, as illustrated in Figure 1bc and Figure 2bc. The differential significant regulation patterns observed in nodules B and C not only as affected genes but also as percentages (22.37% vs. 17.11% up-regulated and 13.16% vs. 5.26% down-regulated TLR genes, respectively) emphasize the importance of considering intra-tumoral heterogeneity in therapeutic strategies [50,51,52,53,54].

For instance, while three TLR signaling mitogen-activated kinases (MAP2K2, MAP2K3, MAPK13) were upregulated and four (MAP2K4, MAP2K6, MAPK9, MAPK10, MAPK11) were down-regulated in ”B”, only MAP2K3 was upregulated and MAP2K6 was down-regulated in “C” (Figure 1). The differential regulation of the mitogen-activated kinases in the two equally graded nodules indicates that participation of the MAPK signaling [55] to the “development, invasion, metastasis, and drug resistance” [56] of cancer cells is not uniform even across the same prostate. Also, while four CHS signaling guanine nucleotide binding proteins (GNB1, GNG10, GNG4, GNG5) were upregulated and three (GNB4, GNG12, GNG7) were down-regulated in “A”, only two genes (GNG7, GNG8) were down-regulated in “B” and no gene was up- or down-regulated in “C” (Figure 2).

The coordination analysis is based on the “Principle of Transcriptomic Stoichiometry” (PTS) [57] requiring expression correlation of the genes whose encoding products are linked in a functional pathway. This principle, theoretically ensuring the efficiency of the functional pathways, extends to the gene networking the classical Dalton’s Law of Multiple Proportion from chemistry [58]. Like the chemical elements that can combine in multiple proportions to produce various molecules, the genes can also be networked in several ways to respond the stringent necessities of the cells in their actual environment.

Figure 3 presents a remarkable increase in coordination scores from normal tissue (1.09% and 2.50%) to cancer nodules (up to 8.97% and 13.51% in nodule A), suggesting that cancer development involves extensive remodeling of pathway interactions [59,60]. This increased coordination degree (higher in the primary nodule “A” than in the secondary “B” and “C”) might represent a more rigid signaling state that could potentially be exploited on novel therapeutic avenues. Of note are also the different ways by which the genes are networked in the four investigated regions, especially between nodules “B” and “C”, a strong argument that the personalization of the cancer gene therapy should go beyond the individual to the most aggressive cancer clones from the tumor.

Identification of the gene pairs whose expression coordination is statistically significant was complemented in each profiled prostate region with the independently expressed gene-pairs. The independent expression analysis selects genes that are not linked within or between any functional pathway(s) in the profiled region, although they can be coordinated in other tissues. Figure 4 shows substantial decreases of the percentages of independently expressed genes, from 5.90% and 4,94% in “N” to 2.84% and 2.02% in “A”, 3.64% and 2.80% in “B”, and 3.43% and 2.80% in “C”. Like the spectrum of the expression coordination in Figure 3, the spectrum of the expression independence in Figure 4 is not repeatable even among the cancer nodules from the same tumor, making less likely to be repeated from person to person.

The transformation of independently expressed gene-pairs in normal tissue into coordinately expressed pairs in cancer, and vice versa, reveals the extensive rewiring of cellular signaling networks during cancer development [61,62,63].

Of particular interest is the switching of gene pair relationships from synergistic to antagonistic expression (or vice versa) in cancer tissue. For instance, the antagonistic pairs MAPK11-NFKBIA and SOS2-NFKBIA in “N” are turned into synergic ones in “A”, while the synergist pair CCL17-CXCL10 in “N” is reversed into an antagonistic one in “A”. In a synergistic pair, the expressions of the two genes are stimulating each-other either directly or through a common transcription factor, while in an antagonistic pair the two genes are inhibiting each-other. These changes suggest that cancerization not only alter individual gene expression but fundamentally reorganize the relationships between signaling components [64,65]. The findings of this analysis question the “universality” of the functional pathways as designed by various specialized software, KEGG [16] and Ingenuity PathAssist [66] included.

The prostate cancer intratumor heterogeneity was revealed in previous studies through coordination analysis of gene networks involved in metabolism [7], enzymatic activity [8], as well in expression coordination of AKT2, MTOR [8], and PTEN with their partners [20]. The unrepeatable remodeling of the gene network among cancer nodules of was also reported in a case of clear cell renal cell carcinoma (ccRCC) for all KEGG-designed functional pathways of the excretory system [15,67], chemokine signaling, basal transcription factors, RNA polymerase, and cell cycle [26]. In thyroid cancer investigations, we analyzed the remodeling of the functional pathways involved in the genetic information processing [68], oxidative phosphorylation [69], or those by which the ncRNAs regulate apoptosis [70].

A significant result of the present study is the hierarchization of the CHS and TLR signaling genes in the four regions with respect to their gene commanding height (GCH) score (Table 4, fold-changes in Table S4). GCH combines the strength of the expression control and expression coordination. High expression control stipulates the importance of that gene for the cell physiology, while high coordination degree indicates how influential it is on the expressions of other genes [71].

As expected, we found different hierarchies in the four regions (Table 4), as well as substantial GCH score changes for the same gene from one region to another (Table S4). For instance, the GCH score of CCL15, a chemokine protective factor for both liver and lung neoplasms [72], was reduced by 122.87x in “A”, 12.42x in “B” and 10.65x in “C”. This substantial reduction of the GCH score means that while experimental manipulation of this gene might have strong (even lethal) repercussions on the healthy tissue, it will have practically no effect on the cancer cells.

The downgrade of VAV2, known to activate RhoA that mediates prostate cancer cell–cell repulsion required for metastasis [73], mainly by substantially reducing its expression control (Table 2), increases the flexibility of the cancer cells to adapt to different local environments. We found also that the EPH receptor A4 (EPHA4), a mediator of the contact inhibition locomotion [74] through regulating Vav2-RhoA was significantly down-regulated in all three cancer nodules by: x^{(A vs N)} = -3.35, x^{(B vs N)} = -7.20, x^{(C vs N)} = -6.34.

The relevance of FADD, an important receptor protein involved in transmitting apoptosis initiation signals [75,76] was significantly reduced in nodule “C” (FC^{(C vs N)} = -17.41), to a lesser extent in nodule “B” (FC^{(B vs N)} = -3.25) and much lesser in “A” (FC^{(A vs N)} = -1.45). In contrast, the GCH score of RIPK1, an essential regulator for necrosis induced by tumor necrosis factor α (TNF-α) [77], increased more in nodule “A” (FC^{(A vs N)} = 4.56) than in the other two nodules (FC^{(B vs N)} = 1.16; FC^{(C vs N)} = 2.04). Our results indicate that the molecular mechanisms for cancer cell survival and proliferation are not uniform across the malignant part of the surgically removed prostate. Furthermore, in a previous work (Figure S1 in Ref. [7]), we reported the inter-nodule heterogeneity of the regulation of genes associated with: block of differentiation, evading apoptosis, immortality, insensitivity to anti-growth signals, proliferation, resistance to chemotherapy, and sustained angiogenesis.

Although not regulated in any cancer nodule with respect to “N”, TOLLIP, encoding an autophagy product protein produced during prostate cancer development [78], had a substantial increase of the expression control strength (and implicit relevance) in nodule B. Interestingly, the same gene exhibited a decrease of the expression control strength in nodule C (Table S2). In contrast, the significantly upregulated IRF5 in all three cancer regions (x^{(A vs N)} = 1.57, x^{(B vs N)} = 2.20, x^{(C vs N)} = 1.94), part of the so-called “seven gene signature of the neuroendocrine prostate cancer” [79], had a higher relevance increase in nodule “C” (FC^{(C vs N)} = 6.20) than in “B” (FC^{(B vs N)} = 1.63) or “A” (FC^{(A vs N)} = 2.40) (Table S4). –

5. Conclusions

The substantial differences in the expression level, control and inter-coordination even between two equally graded cancer nodules from the same tumor invalidates the cancer gene biomarker paradigm [80], regardless of its wide usage [81,82]. Our Genomic Fabric approach [83] extends to novel characteristics (expression control and gene networking) of the recognized transcriptomic heterogeneity of prostate cancer [51,52,53]. Moreover, the reported results confirm the “Stem Cell Theory” [54,84,85] claiming that “the same genetic abnormalities and microenvironmental aberrations lead to different biological effects and clinical outcomes”. On top of this, the GCH analysis identifies also the Gene Master Regulators (GMR) [26] whose manipulation is expected to be the most consequential for the persistence of that phenotype [86]. Given that the set of major cancer nodules GMRs are singled out for each person, they are the most legitimate targets of a personalized gene therapy [7,8].

Comparing the top CHS and TLR genes with the GMRs of these regions (determined in [7]) indicates the importance of the two analyzed immune response pathways for the entire transcriptomes. Thus, in “N”, with GCH = 16.84 (CCL15) and GCH = 7.28 (FADD) both CHS and TLR are far below GCH = 74 of TOR1A (torsin family 1, member A), an important regulator in lipid metabolism [87]. With GCH = 34.15 (GNB1) and GCH = 11.86 (RIPK1) they are far below GCH = 149 of ENTPD2 (ectonucleoside triphosphate diphosphohydrolase 2), in “A”, part of 7-gene prognosis signature of gastric cancer [88]. In “B”, with GCH = 5.46 (CCR6) and GCH = 7.21 (TOLLIP) are below AP5M1 (adaptor-related protein complex 5, mu 1 subunit, a.k.a. MUDENG), recognized for its pro-apoptotic function [89]. Finally, in “C”, with GCH = 14.72 (CCL16) and GCH = 11.34 (RIPK1) both CHS and TLR are below GCH = 49 of BAIAP2L1 (BAI1-associated protein 2-like 1, a.k.a. IRTKS), reported as up-regulated in many tumors [90]. The comparison of the GCHs shows that, at least for this patient, manipulation of the CHS and TLR would not be lethal in any region of the prostate.

The idea of gene master regulators of prostate cancer was used by several other groups, but they are limited to only the expression levels of certain transcription factors [91,92,93,94] or tumor suppressors involved in major processes like adipogenesis [95] or cell-cycle [96]. Nevertheless, our GMRs offer a distinct perspective by considering how critical is the right transcription level for the cell survival (revealed by the expression control across biological replicas) and how influential that gene (through significant expression correlations with other genes). However, because of the unrepeatable nature of the transcriptomic topology, the real value of the GMRs of one person can be tested only on that person. In addition, given the dynamics of the local conditions, the GMRs change in time, so their validation test should be performed not too late from the identification moment.