This study was approved by the Board of Ethical Committee of Clinical Hospital Center Split (2181-147-01/06/M.S.-20-4). All patients have signed a written informed consent. A total of 12 male patients who underwent TURBT and were treated with adjuvant BCG immunotherapy for high-risk NMIBC were included in this study. Inclusion criteria for patient selection were male gender, no previous history of urothelial cancer, no other site malignancies, no active infection, and no antibiotic usage for 20 days prior to the first BCG dose application and during BCG cycles. Immunomodulating conditions such as autoimmune diseases and diabetes were considered as exclusion criteria since they can modulate urobiota composition (Shaheen et al., 2022). Also, all patients were preoperatively treated with the same prophylaxis (norfloxacin). The median follow-up for the whole cohort is 38.5 months. The recurrence was diagnosed with a mean follow-up of 14 months. One patient’s follow-up was limited to 10 months due to a death unrelated to bladder cancer and was excluded from the analysis of response-to-therapy groups, but remained included in all other analyses (labeled Limited FU). The diagram of patient selection with the exclusion process is presented in Figure 1A , whereas Table 1 shows the patient’s demographics and relevant clinical data.

Urine samples were obtained with aseptic bladder catheterization before each of the six induction BCG doses. Therefore, the first sample (timepoint T0) represents the patient’s urinary microbiota prior to BCG instillation, whereas all other samples (timepoints T1-T5) represent the urobiota composition during the course of the BCG induction therapy (shown in Figure 1B ). Urine samples were stored in a sterile container and were transported to the Laboratory for Cancer Research at the University of Split, School of Medicine, within 4 h for further processing.

Approximately 30 to 50 mL of collected urine was centrifuged for 10 min at 7,500 g. The supernatant was discarded, and pellets were resuspended in 100 µL of molecular-grade water and stored at −20°C. Alternatively, whole urine samples were frozen at −80°C, and centrifugation was performed after thawing.

DNA isolation was performed using urine pellets as described previously (Bučević Popović et al., 2018). We used a DNA PowerSoil Kit (Qiagen, Germany) according to the manufacturer’s protocol. Extracted DNA was eluted in 20 µL of water and stored at −20°C. DNA concentrations were determined by NanoDrop (Thermo Fisher, USA) and Qubit fluorometric assay (Thermo Fisher, USA) and ranged from too low for detection to 500 ng/µL. To check for DNA contamination during DNA extraction and sample collection, we performed negative controls using DNA-free molecular-grade water that was passed through the sterilized catheter, followed by the same DNA extraction procedure. Negative controls did not give DNA bands in the PCR reactions.

The PCR reaction was performed using the universal primers 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) that amplify the hypervariable region 4 (V4) of the bacterial 16S ribosomal gene (Caporaso et al., 2011; Caporaso et al., 2012). The second PCR step was performed to add Illumina adapters and sample indices. The same amount of PCR products from each sample was pooled; the mix was purified and then size selected using Agencourt AMPure XP-PCR magnetic beads (Beckman Coulter, USA).

Sequencing was performed using a paired-end 2 bp × 250 bp sequencing reagent cartridge, according to manufacturer’s protocol instructions (Illumina, USA). The library was quantified with Qubit and real-time PCR and checked with a bioanalyzer for size distribution detection. PCR and sequencing were performed at Novogene Company (Novogene Co., China).

Paired-end reads were assigned to samples according to their unique indices and truncated by cutting off the barcode and primer sequences. Demultiplexed reads were imported into the Qiime2 platform (McKinney, 2010; McDonald et al., 2012; Bolyen et al., 2019). Trimmed reads were denoised using the DADA2 (2023.5.0) package using standard parameters (Callahan et al., 2016). This included filtering out of low-quality, PhiX, and chimera reads and finally generation of amplicon sequence variants (ASVs). Next, we used a classifier based on the SILVA database (v.138) to assign taxonomy to the sequences in the ASV table (Pruesse et al., 2007; Quast et al., 2013; Yilmaz et al., 2014). Naïve Bayes pretrained classifier Silva 138 99% OTUs using the 515F/806R region of sequences from Qiime2 were used for this purpose (Bokulich et al., 2018; Robeson et al., 2021).

Microbiome composition was investigated using species relative abundance and alpha and beta diversity indices (Pielou, 1966; Vázquez-Baeza et al., 2013; Weiss et al., 2017). Phylogenetic analysis was performed using MAFFT multiple-sequence alignment and FastTree phylogenetic tree building (Price et al., 2010; Katoh and Standley, 2013).

Data analysis in R using RStudio was performed using packages mia and vegan for microbiome analysis (Oksanen et al., 2012; Ernst et al., 2024). To make graphical representations, package ggplot2 was used (Wickham, 2016). Differential abundance testing was performed in R version 4.2.2 using the ANCOMBC-2 and Linda packages on phylum, family, and genus taxonomic levels (Lin and Peddada, 2020; Zhou et al., 2022). Paired-sample differential abundance was performed using the ALDEx2 package in R (Gloor et al., 2016; Wickham et al., 2019). Alpha diversity (Shannon index, species richness) and beta diversity (Bray–Curtis) indices were calculated using a rarefied table with the same number of reads for all samples (Gloor et al., 2016,241). This depth was chosen according to lowest read number per sample in order to keep all samples in the analysis. From rarefaction curves of alpha diversity measures, this depth was estimated as satisfactory ( Supplementary Figure 1 ).

A total of 5,596,063 reads were obtained from 72 urine samples (mean frequency of reads per sample: 77,723, range: 30,241-101,073). Sequences were assigned to 42,212 amplicon sequence variants or ASVs (mean frequency of reads per ASV: 132.57, range: 1-870,084).

We performed data analysis using MedCalc statistical software (version 23.0.6 MedCalc Software Ltd., Ostend, Belgium).

The distribution of Shannon’s Biodiversity Index (SI) as well as ASV at different time points was summarized with the median and range due to the limited sample size, and the associated 95% confidence intervals for these changes were also reported.

To test for significant differences in both Shannon’s Index and ASV across different time points, we used the Friedman test, with the Wilcoxon test for paired samples applied as post hoc analysis. The signed rank sum test was additionally used to determine whether the median of changes relative to baseline differed from zero. The association between weekly change in SI and the initial biodiversity level of an individual was investigated using simple linear regression. The two-sided significance level was set at 0.05.

The relative abundance of the most abundant genera and phyla of bacteria recorded in all patients at all timepoints is shown in Figure 2A . Enterococcus sp. (17.6%), Serratia (7.1%), Pseudomonas (4.9%), and Escherichia–Shigella (2.8%) were overall the most prevalent genera ( Figures 2A–C ). Firmicutes (49.1%) and Proteobacteria (28.7%) were overall the most prevalent phyla, followed by Actinobacteriota (10.4%) and Bacteriodota (6.3%) ( Figure 2D ).

When looking at the relative composition in urinary microbiota before the BCG therapy (timepoint T0), the most prevalent genera for all patients were Enterococcus sp. (11.35%), Lactobacillus sp. (7.83%), Serratia sp. (5.64%), and Escherichia–Shigella sp. (4.14%) ( Figure 2B ).

There were differences in most abundant taxa in the two response-to-therapy groups when looking at all timepoints. Enterococcus (14.7%), Serratia (9.8%), and Lactobacillus (6.3%) were the most prevalent genera in the responder group ( Figure 2C ). In the non-responder group, most prevalent genera were Lactobacillus (10.3%), Peptoniphilus (7.0%), and Pseudomonas (5.0%). The most prevalent phyla ( Figure 2E ) were similar in both response-to-therapy groups, the most abundant being Firmicutes in the responder group (46.7%) and in the non-responder group (42.2%), followed by Proteobacteria (responder—30.9%, non-responder—29.8%) and Actinobacteriota (responder—10.9%, non-responder—11.9%).

Before the first BCG administration (timepoint T0), the most prevalent genera in responders were Enterococcus (13.9%), Serratia (7.8%), and Lactobacillus (7.6%), whereas in non-responders, Lactobacillus (8.7%), Escherichia–Shigella (6.2%), and Clostridium (5.2%) ( Figure 2D ). Most abundant phyla before the first BCG administration were similar in both response groups: Firmicutes (R—49,2%, NR—42.2%), Proteobacteria (R—29,9%, NR—28.3%), Actinobacteriota (R—9.4%, NR—8.8%), and Bacteroidota (R—5.7%, NR—12.1%) ( Figure 2E ).

In the last timepoint (T5—before the last induction BCG dose), the three most abundant genera were the same as in T0 for BCG responders: Serratia (15.8%), Enterococcus (14.8%), and Lactobacillus (5.4%), whereas for non-responders, the three most abundant genera changed to Enhydrobacter (6.9%), Lactobacillus (6.5%), and Pseudomonas (5.4%) ( Figure 2D ). The most abundant phyla in the last timepoint were Firmicutes (R—43.7%, NR—32.6%), Proteobacteria (R—36.9%, NR—40.5%), Actinobacteriota (R—9.2%, NR—9.8%), and Bacteroidota (R—4.8%, NR—6.1%) ( Figure 2E ).

The baseline Shannon’s Biodiversity Index (SI) for most respondents (50%) was between 7.5 and 8.5 (median SI 8.0, range 6.2-9.6) ( Figure 3A ). We observed a significant change in biodiversity, as measured by the SI, across the six therapy time points (Friedman test, p=0.037). Specifically, there was a notable decrease in the SI compared with the baseline level after the first and third weeks of therapy (p ≤ 0.041), whereas during the other therapy weeks (weeks 2, 4, and 5), the SI showed no significant deviation from the baseline ( Figures 3A, B ). In the initial week of therapy, following the application of BCG, a decrease in biodiversity was noted in 10 out of 12 patients (signed rank sum test, p=0.021), with a median SI decrease of 1.3 (95% CI 0.2, 2.9).

Regarding the species, more precisely ASVs, richness of the bacterial community, we also observed significant changes in the number of observed features relative to the baseline value at the 0.1 significance level (Friedman test, p=0.081). The number of bacterial community features was notably lower only during the first week of therapy (p=0.005), returning to levels comparable with baseline thereafter ( Figures 3C, D ).

When assessing weekly biodiversity changes (between successive therapy doses), we made an interesting observation. In the second week of treatment (T2), nearly three-quarters of patients (8 out of 12 subjects) exhibited an increase in SI compared with their SI levels in the first week of therapy (T1) ( Figure 3E ). However, the median increase did not reach statistical significance (signed rank sum test, p=0.176). In the subsequent weeks of therapy, microbial biodiversity in response to therapy was not so pronounced as we did not observe significance in the direction of weekly changes in terms of growth or decline (signed rank test, p≥0.204), nor did we find a significant correlation between weekly changes and the initial state of biodiversity of each individual (R2 ≤ 13%, p-value for regression coefficient ≥0.244).

When considering alpha diversity variance attributed to response-to-therapy status, responders had non-significantly lower mean Shannon diversity before the first BCG instillation, (median(R) = 7.40, median (NR) = 8.55, p-value = 0.29, Wilcoxon test with Holm adjustment, Supplementary Figure 2A ). Species richness was similar in responders and non-responders in the first time point (median (R)= 1582.5, median (NR) = 1674.5, p-value = 1, Wilcoxon test with Holm adjustment, Supplementary Figure 2B ). Overall microbial community dissimilarity, calculated as Bray–Curtis dissimilarity, did not show any differentiation based on response status (p-value 1, Permanova test, stratified per patient, number of permutations 9999, Supplementary Figure 3A ); however, when comparing only the urobiome community before the start of the therapy (T0), the differentiation was significant at level 0.1 (p-value 0.0539, Permanova test, number of permutations 9999) ( Supplementary Figure 3B ), whereas in the last timepoint (T5), no significant clustering was detected (p-value 0.1856, Permanova test, number of permutations 9999, Supplementary Figure 3C ).

Differential abundance analysis revealed that before the onset of the BCG treatment, non-responders had 12 times higher abundance of genus Aureispira compared with ones that did respond to treatment (p < 0,001) ( Figure 4A ). On a species level, Negativicoccus succinicivorans showed 27 times lower abundance among non-responders (p < 0,001) ( Figure 4A ). At the end of the BCG treatment ( Figure 4B ), genera and species related to the phylum Pirellulaceae had three (p < 0,001) times higher abundance in the non-responder’s group, along with species the Parabacteroides johnsonii, showing the highest increase in abundance (p < 0,001), almost three times.

In patients that responded to the BCG treatment, the only difference during the course of treatment was observed at the final stage of therapy, showing an increase in abundance of almost two times in several species (p < 0,001, ANCOMBC2) compared with a starting point, but all of them were uncultured ( Figure 5 ).

In a group of patients that did not respond to therapy, Linda analysis revealed a significant decline in the abundance of phylum Nanoarchaeota, up to 60 times (p = 0,025) after the first administration of the BCG treatment and an additional decline to 135 times by the end of the therapy (p = 0.005) ( Supplementary Figure 4 ). Figure 6 shows the dynamics of abundance changes on a genus and species levels throughout the course of the therapy in the non-responder group (ANCOMBC2 analysis). Genus Nitrosospira had 10 times higher abundance at the end of the BCG treatment compared with the starting point (p < 0.001). Furthermore, genus Aureispira, except for being more abundant in non-responders before the onset of the therapy, also showed a decreasing trend in abundance throughout the treatment period among non-responders. At the end of the BCG therapy, genus Aureispira was eight times less abundant (p < 0.001) compared with the values before the onset of the treatment.

Furthermore, we did pair sample testing on samples from the same patient to detect differences in taxon abundance between two therapy timepoints (T0 vs. T1, T1 vs. T2, and T0 vs. T5) at levels of species, genus, family, and phylum, but no taxons were detected to be significantly changed.

Regarding BCG bacterial detection, we were not able to determine up to species-level sequences belonging to the genus Mycobacterium. This genus represented 0.18% sequences in both responders and non-responders before the first instillation (T0). In the last timepoint, its relative abundance in responders was higher, 0.33%, whereas in non-responders it was somewhat lower at 0.14% ( Supplementary Figure 5 ).