The study was approved by the Institutional Review Board of Institute of Post Graduate Medical Education & Research (IPGME&R), Kolkata, India. Written informed consent was obtained from all participants and all methods in this study were performed in accordance with the ethical principles founded in the Declaration of Helsinki.

Affected individuals were recruited at Dept. of Urology, IPGME&R, Kolkata. A standard TRUS-guided 18-core prostate biopsy with subsequent histopathological analysis was carried out in those with abnormal serum PSA level (>4.0 ng/ml). All cases included were newly diagnosed patients without any prior treatment before surgery. Individuals with clinical evidence of prostatic lesions receiving radiotherapy and/or undergoing androgen deprivation therapy (ADT) were excluded. Digital images of biopsy slides were reviewed by genitourinary oncology pathologist. Samples were collected into two cohorts. Cohort-1 (Discovery Cohort) containing 13 BPH and 33 PCa tissue biopsy samples were subjected to 16S rRNA amplicon sequencing. Cohort-2 (Validation Cohort) containing 15 BPH and 16 PCa tissue biopsy samples together with Cohort-1 were used for validation using real-time qPCR analyses. Clinical characteristics and other information related to both sample cohorts are described in Table S1 .

A total of 77 formalin-fixed paraffin-embedded (FFPE) tissue samples in both cohorts including 28 BPH and 49 PCa tissue biopsy specimens were received as 10 µm sections on non-charged glass slides. Genomic DNA was isolated according to the standard procedure as previously described (Mitra et al., 2010). Isolated genomic DNA was stored at -20°C for future use. The quality and quantity of extracted DNA was determined by agarose gel electrophoresis and the A260/280 ratio using Synergy H1 Multimode Microplate Reader (BioTek Instruments, Inc., Winooski, VT, USA). Approximately 50 ng and 1 ng of DNA from each sample were used for 16S rRNA amplicon sequencing and real-time qPCR analyses, respectively.

For characterization of bacterial populations, taxonomical analysis, and species identification, different hypervariable regions of 16S rRNA genes were amplified and sequenced on an Ion GeneStudio S5 System (Thermo Fisher Scientific Inc., Waltham, MA, USA) using Ion 550 Chip. Ion AmpliSeq Libraries were prepared using an automated Ion Chef system (Thermo Fisher Scientific Inc., Waltham, MA, USA) and subsequently placed in the Ion Chef System for emulsion PCR. The concentration and quality of the amplicons were measured using Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).

To assign taxonomy, each unique sequence was mapped using two comprehensive 16S rRNA reference databases – Greengenes v13.5 and the Thermo Fisher Scientific in-house MicroSeq 500 v2013.1. Reads were aligned against the databases using MegaBLAST. The expectation value (E-value) for the searches was set to 0.01, and the max target hits value was set to 100. To assign taxonomy, the minimum alignment percentage read to a subject sequence was set to a threshold of 90. A read was assigned to a genus only when the identity score of the sequence alignment was at 97% or higher. For species assignment, the minimum percentage identity of the alignment was set to 99%. The taxonomy distribution counts, or abundance derived from the clustered reads, was subsequently transformed into the relative abundance of the individual species.

For statistical analysis, MicrobiomeAnalyst, a web-based tool, was utilized. QIIME (v1.9.0) was used to evaluate alpha diversity including Observed Species, Chao1, and Shannon indexes. While Observed Species and Chao1 indexes are the indicators of species richness, Shannon index evaluates species diversity. The difference of alpha diversity between groups was evaluated by Wilcoxon Rank-Sum Test using SPSS (version 22).

To compare microbial compositions between different groups, beta diversity was evaluated by calculating weighted UniFrac distances using Bray-Curtis method from the OTU abundance and utilized in Principal Component Analysis (PCoA). PERMANOVA algorithm on weighted UniFrac distance matrices for statistical significance between groups using 999 permutations in QIIME was applied to generate PCoA plots. The Differential abundance analysis between groups was performed using Metastats and P-values were adjusted for multiple hypotheses testing using the False Discovery Rate (FDR) based on the Benjamini-Hochberg procedure. The differential abundances of OTUs and specific OTU enrichment between different groups were determined using LEfSe based on Kruskal–Wallis H test. P-value and FDR were adjusted to 0.05. The unique bacterial composition among sample groups was identified using Random Forest classification algorithm within MicrobiomeAnalyst. Galaxy, an online tool for metagenomics analysis, was used to build the Cladogram.

Functional compositions of the bacterial communities among different groups were predicted using Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) according to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database using STAMP v2.1.3. P-value and FDR cutoff were adjusted to 0.05 level of significance.

Gene-specific primers were designed using Primer-BLAST tool in NCBI database for real-time qPCR analyses and are listed in Table S3 . qPCR primers were obtained from IDT, Inc. (Coralville, IA, USA). The optimum primer melting temperature (Tm) was set at 60°C and the maximum GC content was kept at 55%. qPCR analysis was performed using iTaq Universal SYBR Green Supermix (BIO-RAD, Hercules, CA, USA) in CFX Connect Real-Time PCR detection System (BIO-RAD, CA, USA) with the following thermal profile – one cycle: 95°C for 10 min; 40 cycles: 95°C for 10 s followed by 60°C for 10 s; and finally the dissociation curve at –95°C for 1 min, 55°C 10 s, and 95°C for 10 s. Unless otherwise stated, each sample was performed in duplicate and calculation was made using a -ΔCT method to quantify relative abundance compared with human genomic GAPDH control. The -ΔCt values of each sample were plotted using GraphPad Prism 8.0.1 for data output.

Correlation among bacteria and oncogenic viruses were analyzed using R software. Spearman’s rank test was performed using the ΔCt value of both bacteria and viruses. P-value cutoff was adjusted to 0.05 significance level.

To understand the changes of microbial compositions associated with prostate cancer (PCa) development, we prospectively collected a total of 46 tissue biopsy samples from 13 benign prostatic hyperplasia (BPH) and 33 PCa patients (‘Discovery Cohort,’ Table S1 ). The total DNA of each specimen from the ‘Discovery Cohort’ was extracted and subjected to 16S rRNA amplicon-based sequencing on an Ion GeneStudio S5 System targeting the hypervariable regions of 16S rRNA gene using two distinct primer sets. While set-I primer was used for amplification of hypervariable regions including V2, V4, and V8, set-II primer was used for hypervariable regions including V3, V6-7, and V9. Additionally, as described later, to validate the sequencing data, we also collected another 31 tissue biopsy samples from 16 BPH and 15 PCa patients (‘Validation Cohort,’ Table S1 ).

A total of 35,320,332 raw reads were generated after sequencing of the 46 samples from the ‘Discovery Cohort’. The number of raw sequence reads varied by approximately 10 fold across samples. After quality trimming and chimera checking, 30,641,026 high-quality reads were mapped, while 6,301,956 reads were excluded due to low copy number reads (<10 copy numbers) and 41,312 reads were found to be un-mapped ( Table S2 ). Subsequently, 24,297,758 mapped reads identified Operational Taxonomic Units (OTUs) at the level of family (>95% similarity index), genus (>97% similarity index), and species (>99% similarity index). Overall, a total of six phylum, 127 genera, and 291 bacterial species were identified from all 46 samples in the Discovery Cohort.

Rarefaction curves of species richness against sequences per sample were plotted for BPH and PCa samples to determine the efficiency of the sequencing process ( Figures 1A, B ). Most of the samples, though not completely, reached a saturated plateau phase, indicating that the depth of sequencing was sufficient for the diversity analysis ( Figures 1A, B ). Both Observed Species (p = 0.015) and Chao1 index (p = 0.034) showed that species richness was significantly decreased in PCa samples as compared to BPH samples ( Figures 1C, D ). In addition, the diversity estimators for both Shannon index (p = 0.158) and Simpson index (p = 0.411) indicated a trend of depletion in relative diversity of species composition in PCa samples as compared to BPH samples, although the data were not statistically significant ( Figures 1E, F ).

To assess the diversity among two groups, we evaluated weighted UniFrac distance matrix from the OTU abundance through utilizing PERMANOVA algorithm and subsequently applied in Principal Component Analysis (PCoA) ( Figure 1G ). Given that PCoA analyses revealed no significant difference (p = 0.082) in the bacterial compositions between BPH and PCa samples ( Figure 1G ), a ‘Random Forest’ algorithm was applied to further confirm the difference in bacterial community among the BPH and PCa biopsy samples ( Figure 1H ). The decision trees extracted from the random forest classification identified distinct bacterial composition in two BPH samples and 11 BPH samples exhibited overlapping species with PCa samples (class error: 0.850; Figure 1H ). In contrast, all 33 PCa samples demonstrated unique bacterial compositions ( Figure 1H ).

The bacterial communities associated with the diseased prostate lesions were further analyzed at different taxonomic levels. Six phyla including Firmicutes, Proteobacteria, Bacteroidetes, Actinobacteria, Fusobacteria, and Deinococcus-Thermus collectively comprised the entire sequences in both groups ( Figures S1A, B ). In general, Proteobacteria was the most abundant phylum in the prostate microbial ecology, contributing ~40.6% in diseased prostate lesions ( Figures S1A, B ). Overall, the results indicated that only Actinobacteria phylum was significantly depleted in PCa samples as compared to the BPH category ( Figures S1A, B ). At the genus level, Prevotella, Cupriavidus, Propionibacterium, Acinetobacter, and Corynebacterium represented the top five genera in both samples ( Figures S1C, D and S2A ). Of all genera detected, both groups shared approximately half of the total genera identified, i.e., 56/107 in BPH category and 56/118 in PCa lesions ( Figure S2B ).

To identify the differentially enriched genera within BPH and PCa samples, we employed Linear Discrimination Analysis (LDA) Effect Size (LEfSe) method ( Figure S1E ). LEfSe analysis of top 20 bacterial genera identified Kocuria, Staphylococcus, Corynebacterium, Cellvibrio, Pseudomonas, Paracoccus, Brachybacterium, Pseudoxanthomonas, Anaerococcus, Stenotrophomonas, Microvirga, Empedobacter, Lysobacter, Brevibacterium, Comamonas, Serinicoccus, Rhodobacter, Chryseobacterium, and Aeromicrobium as enriched genera in BPH samples, whereas only Bradyrhizobium genus was found to be significantly elevated in PCa samples ( Figure S1E ).

Detailed analyses of bacterial compositions at species level demonstrated Prevotella copri, Cupriavidus campinensis, Propionibacterium acnes, and Paracoccus sp. covered almost 50% of species variety in both BPH and PCa samples ( Figures 2A, B ). Next, cladogram and LEfSe analyses were conducted to further uncover the differential species compositions between BPH and PCa tissue biopsy samples ( Figures 2C, D ). The LDA scores demonstrated that among the top 20 most significantly enriched species Kocuria palustris, Cellvibrio mixtus, Pseudomonas stutzeri, Paracoccus sp, Staphylococcus hominis, Corynebacterium tuberculostearicum, Brachybacterium paraconglomeratum, Staphylococcus arlettae, Staphylococcus cohnii, and Anaerococcus octavius were notably elevated in the BPH samples, while Cupriavidus taiwanensis, Methylobacterium organophilum, Brevundimonas vancanneytii, Neisseria flavescens, Acinetobacter junii, Bradyrhizobium cytisi, Cupriavidus basilensis, Caulobacter segnis, Leclercia adecarboxylata, and Neisseria elongata were significantly increased in the PCa samples ( Figure 2D ).

To confirm the 16S rRNA sequencing data, quantitative real-time PCR (qPCR) analyses were performed of top two bacteria in each category – BPH (K. palustris and C. mixtus) and PCa (C. taiwanensis and Methylobacterium organophilum), along with three most abundant bacteria (P. copri, C. campinensis, and P. acnes) identified in diseased prostate lesions, using two distinct species-specific primers ( Figure 3 ). Moreover, in order to further corroborate the results, we used two sample cohorts – ‘Cohort-1’ (Discovery Cohort) and ‘Cohort-2’ (Validation Cohort) as described in Table S1 . As similar to LEfSe analyses between BPH and PCa biopsy samples, qPCR data further confirmed Kocuria palustris and Cellvibrio mixtus as BPH specific and Cupriavidus taiwanensis, and Methylobacterium organophilum as PCa specific bacterial species in both sample cohorts ( Figures 3A, D ). In contrast to 16S rRNA sequencing results, qPCR analyses demonstrated Prevotella copri, Cupriavidus campinensis, and Propionibacterium acnes were significantly enriched in PCa samples as compared to the BPH group in both sample cohorts ( Figures 3E–G ).

In order to visualize the functional effects resulting from the altered microbial community associated with disease progression in prostate, we employed Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) software. PICRUSt analyses can predict the functional Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways related to the composition of a metagenome and have been demonstrated to provide a good representation of metagenomic prediction. The LEfSe outputs of KEGG pathways among BPH and PCa tissue biopsy samples identified functions related to starch and sucrose metabolism, galactose metabolism, carbohydrate metabolism, primary immunodeficiency, ubiquitin system, Ion channels, proteasome, phenylpropanoid biosynthesis, electron transfer carriers, glycan degradation, and N-Glycan biosynthesis were significantly associated with BPH condition, while pathways such as nitrotoluene degradation, steroid hormone biosynthesis, non-homologous end-joining, and primary bile acid biosynthesis were directly linked to PCa development ( Figure S3A ). PCoA also demonstrated that the predicted functions of bacterial compositions among BPH and PCa were significantly clustered (p < 0.05) ( Figure S3B ).

It is known that prostatic disease commonly affects middle-aged and elderly people (Perdana et al., 2016). In order to determine whether there is any link between age and microbial dysbiosis among BPH and PCa samples, correlation studies among the top 10 bacterial species identified in LEfSe analyses in each category with the patient’s age were further conducted ( Figures S4A, B ). While no correlation was established among BPH-specific bacteria with patient’s age, a number of bacterial species specific to PCa samples were moderately positively correlated with increasing patient’s age ( Figure S4B ). These include Neisseria elongate (r = 0.262), Caulobacter segnis (r = 0.277), and Bradyrhizobium cytisi (r = 0.346) ( Figure S4B ). Among the top two bacteria in LEfSe analyses along with the three most abundant bacteria in diseased prostatic lesion, C. taiwanensis (r = 0.176), C. campinensis (r = 0.162), and P. acnes (r = 0.191) were weakly positively correlated with PCa pateint’s age ( Figure S4C ). In contrast, while C. taiwanensis (r = -0.239) and M. organophilum (r = -0.131) were somewhat negatively correlated with PCa grade based on Gleason score, P. copri was found to positively correlated (r = 0.196) with PCa grade ( Figure S4D ).

Studies suggest that a number of human oncogenic viruses including high-risk HPVs and EBV are associated with PCa development (Pascale et al., 2013; Whitaker et al., 2013; Smelov et al., 2016). To evaluate the potential involvement of viral etiology in our samples we designed qPCR primers for seven human tumor viruses including EBV, two high risk HPVs – HPV-16 and HPV-18 - HBV, HTLV-1, HCV, KSHV, and MCPyV along with two more human polyomaviruses – JCV and BKV ( Table S3 ). qPCR analyses of Cohort-1 with primer set-I against all virus-specific antigens demonstrated that only EBV, HPV-16, HPV-18, and HBV are significantly associated with PCa samples as compared to BPH controls ( Figures 4A–D and S5A–F ). The housekeeping gene, human GAPDH gene, was utilized as control assuming the genomic segment bearing GAPDH gene remained unaffected in both BPH and PCa samples. A higher negative -ΔCt (average GAPDH Ct value – average target primer Ct value) indicated elevated presence of the virus in the sample as detected by specific primer set targeting specific viral gene. We further validated the association of these four selected tumor viruses in Cohort-2 using a second set of primers ( Figures 4A–D and Table S3 ). The results clearly demonstrated that all four tumor viruses were significantly associated with the PCa lesions in comparison to BPH samples ( Figures 4A–D ).

To further corroborate the connection of these tumor viruses with the identified microbial signature associated with BPH and PCa lesions, the co-occurrence and co-exclusion patterns of EBV, HPV-16, HPV-18, and HBV with the most abundant bacterial species identified in LEfSe and qPCR analyses in each category were further investigated ( Figure 4E ). As depicted in both LEfSe and qPCR analyses, BPH-specific bacteria K. palustris and C. mixtus and PCa-specific bacteria C. taiwanensis and M. organophylum were moderately correlated with each other (r = 0.362 and 0.394, respectively), indicating that they fell into two distinct groups ( Figure 4E ). Among BPH-specific bacteria, interestingly only K. palustris was positively correlated with EBV (r = 0.446) and as expected negatively correlated with HPV-16 (r = -0.461) and to a lesser extent with HBV (r = -0.263) ( Figure 4E ). Among PCa-specific bacteria, while C. taiwanensis was positively correlated with HPV-16 (r = 0.411) and to lesser extents with HPV-18 (r = 0.328) and HBV (r = 0.246), M. organophylum were positively correlated with EBV (r = 0.456) and to lesser extents with HPV-18 (r = 0.325) and HPV-16 (r = 0.261) ( Figure 4E ). Among the most abundant bacteria in both BPH and PCa tissue samples, all three species P. copri, C. campinensis, and P. acnes were in general found to be somewhat positively correlated with the PCa-specific bacteria but not with BPH-specific bacteria ( Figure 4E ). Among these three species P. copri and P. acnes were particularly strongly correlated (r = 0.543) ( Figure 4E ). In addition, while P. copri was positively correlated with both C. taiwanensis (r = 0.369) and to a lesser extent M. organophylum (r = 0.246), C. campinensis and P. acnes were correlated with only M. organophylum (r = 0.318, and 0.282, respectively) ( Figure 4E ). While P. copri and P. acnes were robustly correlated with two tumor viruses – EBV (r = 0.551 and 0.599, respectively) and HPV-18 (r = 0.509 and 0.444, respectively) - C. campinensis was weakly correlated with only HPV-18 (r = 0.279) ( Figure 4E ). Among the four tumor viruses, EBV and HPV-18 were positively correlated with each other (r = 0.551), whereas HPV-16 and HBV were grouped together (r = 0.370) ( Figure 4E ).