Microbiome Dysbiosis Associates to Castration Resistance and Cancer Stemness in Metastatic Prostate Cancer

1. Introduction

Prostate cancer is the second leading cause of death in males in America [1,2]. It is estimated that roughly one in eight men will be diagnosed with prostate cancer in his lifetime [1,2]. These cancers are classified as either localized, regional, or distant based on the degree by which they’ve metastasized to other bodily sites. Localized and regional prostate cancers are rarely fatal, though distant prostate cancers exhibit a 5-year survival rate of only 32% [1]. Understanding the factors that promote the acquisition of the more aggressive stage is crucial to improving these patients’ survival.

The growth and proliferation of prostate cancers are known to be heavily mediated by testosterone signaling [3]. In a healthy individual, testosterone is capable of binding to and activating the androgen receptor (AR) protein, which then stimulates the production of secretory proteins in the prostate [4,5]. Through mechanisms less understood, AR signaling is also known to cause specific genomic deletions, amplifications, and translocations that promote the growth and proliferation of prostate cancers [4,5]. As such, androgen deprivation therapies are commonly implemented as treatments for high-risk prostate cancers [6]. These largely include medical and surgical castrations [6]. Nonetheless, these therapies yield fairly poor survival rates, as the patients who undergo castration likely have more advanced diseases [6,7]. Many patients develop castration-resistance (CR), which is defined by sustained growth of a cancer despite serum testosterone levels being at or below the level expected with castration. The mechanisms by which prostate cancers acquire this resistance are well characterized, with AR signaling being integral to many [8]. These involve mutations in AR, mutations in AR coactivators and corepressors, androgen-independent activation of AR, and alternate means of androgen biosynthesis [8]. In these cases, AR signaling inhibitors are often prescribed as an adjuvant treatment to androgen deprivation therapy [9,10,11,12,13,14]. However, many patients remain insensitive to these as well [15], with few treatment options available thereafter. The median survival length for non-metastatic castration-resistant prostate cancer (CRPC) cases is estimated to be 30.3 months, and that for metastatic CRPC cases is only 13.3 months. Understanding the factors that influence the acquisition of CR in prostate cancers may significantly improve our ability to treat these diseases. Understanding the influences of metastasis in CRPC may prove doubly useful.

Numerous genetic factors have been identified for their implications in CRPC, primarily comprised of mutations in the genes of the AR signaling pathway [8]. Metastatic contributors have been identified, too, and include the loss of PTEN, aberrations in the PI3K-AKT signaling pathway, and the acquisition of DNA repair defects [16]. The degree of stemness observed of a cancer is also known to influence CR and metastasis [17,18,19]. Cancer stem cells (CSC) are thought to compose only about 1% of tumor’s mass, though they are crucial toward the tumor’s growth and proliferation [20]. CSCs are believed to originate from epithelial cells through a process known as epithelial-mesenchymal transition (EMT) [21]. In this process, malignant epithelial cells gain mesenchymal-like traits, and become highly invasive in doing so [22]. The formation of CSCs through the process of EMT provides a tumor with a high capacity for colonization, ultimately promoting the cancer’s metastasis [17,18,19].

The influences of epigenetic factors in CRPC are less explored, though the human microbiome may be highly relevant. The human microbiome is a collection of microorganisms that populate the gastrointestinal system [23]. The microbiome has become increasingly implicated in human diseases, including inflammatory bowel disease, psoriasis, and diabetes among others [24,25]. The microbiome is thought to influence an array of biological pathways, largely through metabolite-mediated immune modulation [26,27]. As such, studies have also characterized the microbiome for its implication in various cancers, particularly colorectal cancers [28,29,30,31]. Less is known regarding the microbiome’s influence beyond the gastrointestinal system, though numerous studies have demonstrated the importance of the microbiome in the development and progression of prostate cancer [32,33]. Specific dysbioses of the gut microbiome have been identified between castration-sensitive (CS) and CR prostate cancers [33]. Moreover, antibiotic therapies and fecal transplants in mouse models have demonstrated the gut microbiome’s ability to modulate the effectiveness of androgen deprivation therapy [34]. The mechanisms of these relations are less understood.

Hence, this study attempts to characterize microbiome dysbioses for correlations to both CR and cancer stemness. RNA sequencing data was downloaded for bone and soft tissue biopsies of patients with metastatic CRPC across two studies: phs000915 (n = 147) and phs001141 (n = 143). These sequences were mapped to bacterial sequences to yield species-level abundance approximations in each sample. We identified exact correlations of these species to known transcriptional markers of CR and cancer stemness. Further, we observed enrichment of the AR, PI3K-AKT, and endocrine resistance signaling pathways with respect to these species’ abundances. Specific enrichment of EMT and pluripotency signaling was also observed. We propose that the human microbiome is heavily associated with CR and metastasis in prostate cancer. Through this investigation, we may better understand the pathology of metastatic CRPC, creating new avenues of research as to the treatment of this disease.

2. Results

2.1. Cross-Study Normalization and Contamination Correction

In order to account for variation in the sequencing procedures of the chosen datasets, cumulative sum scaling was performed as a means of normalization (see Materials and Methods). PCoA was conducted to demonstrate the effectiveness of this technique. In this analysis, the abundance values of all species in a sample are reduced into several arbitrary dimensions. The proximity of one sample to another is reflective of wholistic similarities in the samples’ abundance profiles. The samples were analyzed both before (Figure 1A) and after (Figure 1B) the normalization procedures were performed. After normalization, the samples were observed to be of closer proximity, with fewer outliers present. This served to confirm the compatibility of the chosen samples for further analyses. Cumulative sum scaling was similarly performed on the samples’ gene expression counts.

In tissue extraction and sequencing procedures, there remains the possibility for contaminant species to be introduced into a sample [35]. Contamination correction was performed to identify and exclude these species (see Materials and Methods). To visualize the phylogenic division of these species, they were grouped by class of phylum and plotted in a phylogenic tree (Figure 1C). The bone metastases were found to contain the greatest number of contaminant species, followed by the liver and lymph node metastases.

2.2. Species Abundance Correlates to Castration-Resistance Biomarker Expression

A list of markers known to be implicated in CR at a transcriptional level was collected from literature [16]. Among others, this includes AR, PTEN, several PI3K family genes, and several AKT family genes. Spearman’s correlations were used to assess the relationship between each marker and each species. Of the bone and lymph node cohorts, a remarkable number of species were significantly correlated to these markers (Figure 2A). Among others, greater abundance of Staphylococcus xylosus corresponded to an increase in the expression of all but few of the markers. Lesser abundance of Streptococcus pneumoniae and Veillonella parvula also corresponded to an increase in the markers’ expression. Significance in the liver cohort was less pronounced, though numerous correlations were still observed. The overlap in species-marker correlations between each cohort was further plotted (Figure 2B). 35 significant correlations were common to all sites of metastasis. The abundance values of each species were further simplified to binary “high” or “low” classifications, based on each sample’s relation to the median abundance of that species. The expression counts of several of markers were plotted with respect to the abundance of Klebsiella pneumoniae and Pseudomonas savastanoi (Figure 2C). Lesser abundance of these species corresponded to significantly greater expression of the AKT and PI3K family genes.

2.3. Castration-Resistance Pathway Enrichment

Select KEGG gene sets were chosen to model the cellular pathways of CR. The Prostate Cancer pathway contains many genes of which aberrations are indicative of CR, including AR, CREBBP, and PTEN [16]. PI3K-AKT signaling is known to be dysregulated in CRPC [16]. The Endocrine Resistance pathway was chosen to model resistance to androgen deprivation therapy.

Gene set enrichment analysis was conducted to assess the extent by which each of these pathways was enriched with respect to each species’ abundance. Of the lymph node cohort, enrichment plots were created to visualize the running enrichment score of these pathways with each species (Figure 3A). Only the twenty species with the greatest number of significant correlations to the above CR markers are shown. The peak of each curve indicates the total enrichment score of the pathway. These pathways were found to be negatively enriched with respect to most species, meaning that lesser abundance of these species yields an increase in AR signaling, PIK3-AKT signaling, and endocrine resistance. The individual genes of each pathway were further assessed for correlation to these species’ abundances. This served to determine the significance by which each of the pathways’ components were enriched. Spearman’s correlations were computed between each gene’s expression and each species’ abundance. The resultant correlation coefficients and test statistics were combined for all cohorts and plotted (Figure 3B). Only the twenty species with the greatest number of significant correlations to the above CR markers are shown. When viewed individually, most genes can be noted for their consistently positive or consistently negative correlations to these species. AR and PTEN, for instance, are downregulated with greater abundance of all but few species. This may suggest that the microbiome is generally less diverse in high-risk CRPC patients. Whether these observations are mechanistically linked is unclear, as is the direction of their potential regulation.

2.4. Species Abundance Correlates to Cancer Stem Cell Marker Expression

Cancer stemness is known to be heavily implicated in CR and metastasis [17,18,19]. To assess the microbiome’s relation to these factors, a list of transcriptional CSC markers was first collected from literature [36]. Among others, this includes CD44, CXCR4, SOX2, and NANOG. Spearman’s correlations were used to assess the relationship between each marker and each species. Of all cohorts, select markers appeared to be broadly upregulated or broadly downregulated with respect to these species (Figure 4A). EGFR expression was negatively correlated to the abundance of all but a few species, as was SLC3A2. BMI1, CD24, and TDGF1 were widely overexpressed with greater abundance of these species. It is unclear whether these genes are linked to the microbiome mechanistically, though the consistency in the correlations’ directions should be noted. The overlap in species-marker correlations between each cohort was further plotted (Figure 4B). Only 2 significant correlations were common to all sites of metastasis. The abundance values of each species were further simplified to binary “high” or “low” classifications, based on each sample’s relation to the median abundance of that species. The expression counts of EGFR and SLC3A2 were plotted with respect to the abundance of Brevundimonas subvibrioides and Geobacillus thermodenitrificans (Figure 4C). Lesser abundance of these species corresponded to significantly greater expression of these markers.

2.5. Epithelial Mesenchymal Transition and Pluripotency Regulation Pathway Enrichment

CSC formation is believed to be dependent on the invasion of epithelial cells through the process of EMT [21]. The microbiome’s relevance to these processes has been demonstrated [37,38], though not with regards to CR or metastasis in prostate cancer. Select KEGG gene sets were first chosen to model EMT and stemness. The Adherens Junction pathway contains many of the genes involved in EMT, and was chosen to model this process [39]. The Signaling Pathways Regulating Pluripotency of Stem Cells pathway was chosen to model stemness.

Gene set enrichment analysis was again conducted to assess these pathways for enrichment with respect to each species’ abundance. Of the lymph node cohort, enrichment plots were created to visualize the running enrichment score of these pathways with each species (Figure 5A). Only the twenty species with the greatest number of significant correlations to the above CSC markers are shown. Similar to the above pathways of AR signaling, PI3K-AKT signaling, and endocrine resistance, these pathways were negatively enriched with respect to all but few species. This may indicate that EMT and stemness are greater in patients with a lack of microbial diversity. The individual components of each pathway were further assessed for correlation to these species’ abundances. This served to determine the significance by which each of the pathways’ components were enriched. Spearman’s correlations were computed between each gene’s expression and each species’ abundance. The resultant correlation coefficients and test statistics were plotted for all cohorts (Figure 5B). Only the twenty species with the greatest number of significant correlations to the above CSC markers are shown. Of the bone and lymph node cohorts, the consistency in the genes’ correlational directions can again be noticed. Many of the Wnt family genes of the Adherens Junction pathway are negatively correlated to a majority of the species. This again may suggest the relation between decreased microbial diversity and an increased capacity for invasion and metastasis.

3. Discussion

Our results suggest that the tumoral microbiome is strongly connected to CR in prostate cancer. Of all cohorts, numerous species were observed to correlate significantly to the expression of the chosen CR markers. AR expression was significantly greater in samples with lesser abundance of Veillonella parvula and Streptococcus pneumoniae. The reverse was true of the genus Staphylococcus, which was previously reported to be of greater abundance in prostate cancers [40]. Species of the genus Shewanella were also observed to correlate positively to AR expression, and have been shown to be enriched in malignant prostate cancers [41]. The microbiome is known to be implicated in an array of human diseases, and is thought to exert its effects through the release of metabolites [24,25,26,27]. It is unknown whether these metabolites interact directly with AR or related proteins, though metabolomic analyses may speak to this influence.

Several PI3K and AKT family genes, too, were observed to correlate negatively to many of the species studied, namely Klebsiella pneumoniae and Pseudomonas savastanoi. The PI3K-AKT pathway as a whole was negatively enriched with these species, as well. The effects of probiotics on this pathway have been demonstrated [42]. It is thought that that metabolites of select species are capable of suppressing aberrant activation of this pathway, ultimately suppressing a cancer’s growth [42]. Further metabolomic analysis might confirm the extent by which this is true in CRPC. Nonetheless, the tumoral microbiome appears to closely follow the AR, PI3K-AKT, and endocrine resistance signaling pathways. Specific dysbioses of the microbiome may ultimately be involved in the acquisition of CR in prostate cancer.

We also observed similar microbial relations to cancer stemness and pluripotency. Numerous CSC markers were found to correlate significantly to the abundance of these species, including Brevundimonas subvibrioides and Geobacillus thermodenitrificans. Among these markers, EGFR and SLC3A2 were significantly downregulated with respect to greater abundance of most species. Hence, lesser abundance of these species may correlate to an increase in cancer stemness. Interestingly, androgen deprivation therapy is known to decrease the diversity of the gut microbiome, ultimately yielding a broad decrease in abundance [43,44]. Given the above species-marker correlations, this decrease may correlate to increased EGFR and SLC3A2 expression, increased EGFR and SLC3A2 signaling, and ultimately increased pluripotency [45,46]. Overexpression of EGFR has been shown to be implicated in the metastasis of prostate cancers to bone [47]. SLC3A2 has similarly been shown to regulate proliferation, migration, and therapy resistance in cancer cells [48]. Our analyses demonstrate the microbiome’s correlation to these factors, suggesting it may be involved in the rapid progression of CRPC.

We observed the microbiome’s correlation to the EMT and pluripotency regulation pathways studied, as well. These pathways were negatively enriched with respect to a majority of the species studied. Moreover, many of the pathways’ component genes were consistently downregulated with greater abundance of these species. This was especially pronounced for the Wnt family genes of the EMT pathway. Androgen deprivation therapy is known to decrease diversity of the microbiome [43,44]. Our analyses suggest that this decrease may correlate to enrichment of the EMT pathway, ultimately promoting the migration of epithelial cells into the mesenchyme. In this way, the microbiome may encourage the formation of CSCs, ultimately yielding an increase in a tumor’s colonizing capacity [17,18,19,21]. With the assumption that cancer stem cells only originate after EMT, the observed correlations of these species to pluripotency may only be coincidental and mediated by EMT. The microbiome and its metabolites are known regulators of EMT [37,38], though less is known of the microbiome’s relation to pluripotency. Investigation the interaction of microbial metabolites with the above CSC markers may be used to further test this hypothesis.

Our results are limited due to the correlational nature of this study. We are unable to claim a causal relationship between microbiome dysbiosis with CR. Further metabolomic analysis may serve to confirm these relations. The sample collection and sequencing procedures used by these studies also inherently differed. We attempted to mitigate these differences through the above normalization procedures, though the potentially confounding effects should be noted. Additionally, species-level profiling performed with the use of a reference sequence database will only capture culturable species and may omit species otherwise present. This is common to most microbiome studies that employ direct sequence alignment.

4. Materials and Methods

4.1. Data Acquisition

RNA sequencing data was downloaded from the dbGaP Data Browser (https://www.ncbi.nlm.nih.gov/gap/) for bone and soft tissue biopsies of patients with metastatic CRPC. These samples spanned two studies: phs000915 (n = 147) and phs001141 (n = 143). Only bone (n = 159), lymph node (n = 92), and liver (n = 39) metastases were analyzed, as the remaining sites contained insufficient sample sizes. Samples of each site were considered independently throughout the remainder of this study.

4.2. Bacterial Read Mapping

Sequences were mapped to bacterial species using the Pathoscope 2.0 software [49], with reference sequences sourced from the NCBI Nucleotide Database (https://www.ncbi.nlm.nih.gov/nucleotide/). Bacterial reads were targeted for quantification, and human reads were filtered. Standard parameters were chosen.

4.3. Gene Read Mapping

Sequences were mapped to the hg38 reference genome using the STAR software [50], with reference sequences sourced from the NCBI Nucleotide Database (https://www.ncbi.nlm.nih.gov/nucleotide/). Standard parameters were chosen.

4.4. Cross-Study Normalization

Cumulative sum scaling was performed to increase the validity of cross-study comparisons. In this technique, the expression value of a gene is divided by the sum of all genes’ expressions in that sample. This technique was similarly performed on species abundance values. Principle coordinate analysis (PCoA) was conducted to confirm the effectiveness of this technique. Sample dissimilarities were calculated by a Euclidean distance.

4.5. Microbial Contamination Correction

In tissue extraction and sequencing procedures, there remains the possibility for contaminant species to be introduced into a sample [35]. These species are not reflective of a patient’s tumoral microbiome, and are likely introduced in a fixed amount. As such, these species are expected to be of similar abundances in all samples, regardless of the total abundance of taxa in a sample. To identify and exclude these species, the total abundance of all species was summed in each sample. Spearman’s correlations were used to assess the relationship of each individual species to the total abundance of all species in each sample. Species that did not exhibit a significant relation (p < 0.05) were deemed contaminants and were excluded from further analyses.

4.6. Expression Correlation Analyses

A list of genes known to be implicated in CR at a transcriptional level was collected from literature [16]. A list of known CSC markers was similarly collected [36,47,48,51]. Spearman’s correlations were used to assess the relation of each species’ abundance to each gene’s expression.

4.7. Gene Set Enrichment Analyses

The clusterProfileR R package was used to assess pathway enrichment with respect to each species’ abundance [52]. Pathways were sourced from the KEGG PATHWAY Database (https://www.genome.jp/kegg/pathway.html). The Prostate Cancer (hsa05215), PI3K-AKT (hsa04151), and Endocrine Resistance (hsa01522) pathways were used to model CR. The Adherens Junction (hsa04520) and Signaling Pathways Regulating Pluripotency of Stem Cells (hsa04550) pathways were used to model EMT and stemness, respectively.

5. Conclusions

Our observations suggest that the microbiome is heavily associated with CR and stemness in prostate cancer. Numerous species were observed to correlate strongly to the expression of the CR and CSC markers studied. These species further correlated to enrichment of select signaling pathways involved in CR and cancer stemness. The potential regulation of these factors by the microbiome should be further investigated to explore the microbiome’s relevance in CRPC. Ultimately, this may provide great nuance as to the aggressive nature of this disease. By exploring the microbiome’s implications in CR, stemness, and metastasis, we may better understand the pathology of high-risk prostate cancers, allowing for needed improvement in the treatment of patients with these diseases.