Recurrent exon-deleting activating mutations in AHR act as drivers of urinary tract cancer

Bladder cancer has a high recurrence rate and low survival of advanced stage patients. Few genetic drivers of bladder cancer have thus far been identified. We performed in-depth structural variant analysis on whole-genome sequencing data of 206 metastasized urinary tract cancers. In ~ 10% of the patients, we identified recurrent in-frame deletions of exons 8 and 9 in the aryl hydrocarbon receptor gene (AHRΔe8-9), which codes for a ligand-activated transcription factor. Pan-cancer analyses show that AHRΔe8-9 is highly specific to urinary tract cancer and mutually exclusive with other bladder cancer drivers. The ligand-binding domain of the AHRΔe8-9 protein is disrupted and we show that this results in ligand-independent AHR-pathway activation. In bladder organoids, AHRΔe8-9 induces a transformed phenotype that is characterized by upregulation of AHR target genes, downregulation of differentiation markers and upregulation of genes associated with stemness and urothelial cancer. Furthermore, AHRΔe8-9 expression results in anchorage independent growth of bladder organoids, indicating tumorigenic potential. DNA-binding deficient AHRΔe8-9 fails to induce transformation, suggesting a role for AHR target genes in the acquisition of the oncogenic phenotype. In conclusion, we show that AHRΔe8-9 is a novel driver of urinary tract cancer and that the AHR pathway could be an interesting therapeutic target.

www.nature.com/scientificreports/ However, the relevance and functional consequences of cancer-associated AHR mutations in bladder cancer are still poorly understood.
Here, we describe the identification of novel recurrent exon-deleting AHR alterations in the pan-cancer whole-genome sequencing database of metastatic cancers from Hartwig Medical Foundation. Using state-ofthe-art data analysis tools with improved functionality for detecting structural variants we identified a previously unnoticed recurrent in-frame deletion of exons 8 and 9 in AHR. Together with the Q383H point mutation and AHR gene amplifications, AHR variants comprise ~ 19% of patients with urinary tract cancer. We demonstrate that the recurrent alterations lead to constitutively activated AHR signaling and induce an oncogenic phenotype in bladder cells. Our results suggest that aberrant AHR signaling is an important driver of urothelial tumorigenesis.

AHR is frequently mutated in urinary tract cancers. To identify genetic factors involved in urinary
tract cancers, we analyzed the Hartwig database, which represents whole-genome sequencing data of solid metastatic tumors and normal tissue of 4500 patients 14 . Utilizing improved structural variant detection algorithms 15,16 , we detected a novel deletion spanning exons 8 and 9 of the AHR gene (AHR Δe8-9 ) in 22 (10.7%) out of the 206 urinary tract samples (Fig. 1a,b). The AHR Δe8-9 is identified in several urinary tract cancer subtypes, most of which are in bladder cancer (12/22) and pyelum (7/22) (Supplementary Table 1). The AHR Δe8-9 mutant is highly specific to urinary tract cancer, because in the complete pan-cancer data set of 4500 patients it was identified in only three non-urinary tract samples (2 × lung, 1 × skin; Supplementary Table 2). The identified exonic deletion was validated in RNA-sequencing, which was available for 129 out of 206 urinary tract samples. Read support for the AHR exon 7-10 splice junction was identified in all urinary tract samples harboring AHR Δe8-9 (n = 16; Supplementary Table 1) and confirmed the in-frame loss of exons 8 and 9 for six samples in which the genomic start position of the deletion is located within the coding sequence of exon 8.
The relatively small size of the AHR deletion (~ 3 kb) may explain why it has not been identified in previous studies 9,17 . Most of the 3' breakpoint junctions are positioned in a narrow window of 30 bp between 2 Alu elements (AluYa5 and AluY) and overlap with the 3' site of the AluYa5 element (Fig. 1a). Such inverted Alu repeats have been identified as strong initiators of genetic instability 18 and we hypothesize that the inverted Alu repeats in the AHR gene facilitate the deletion of exons 8 and 9. The high prevalence in urinary tract cancer of such a rare event strongly suggests positive selection. Analysis of an independent RNA-sequencing dataset of urinary tract cancers from Weill Cornell Medicine confirmed the presence of the 7-10 splice junction in 8% (2/24) of samples. One of the bladder cancer samples with the detected splice junction was of primary origin (out of 8 primary tumors in this dataset). Together, these results demonstrate that the AHR Δe8-9 is present in primary and metastatic urinary tract cancer (Supplementary Table 3).
To identify additional events in the AHR gene we performed a targeted analysis of the Hartwig database. We identified recurrent gene amplifications and a recurrent point mutation, in 16 (7.8%) and 2 (1.0%) of the urinary tract cancer patients, respectively (Fig. 1a,b). The c.1149G>C (AHR Q383H ) point mutation was previously reported as an APOBEC-associated hotspot mutation in bladder cancer based on the TCGA PanCancer Atlas data 11 . The AHR Q383H point mutation was found in 11 patients out of ~ 11,000 patients in the TCGA PanCancer Atlas, which are mainly of primary tumor origin 19,20 . The majority (n = 8) occurred in bladder urothelial carcinomas with the other three in hepatocellular carcinoma, papillary renal cell carcinoma, and lung adenocarcinoma. Thus, like the AHR Δe8-9 variant, the AHR Q383H mutation is highly specific for urinary tract cancer.
AHR is a ligand-activated transcription factor that contains a basic Helix-Loop-Helix/PER-ARNT-SIM (bHLH/PAS) motif. AHR contains two PAS domains of which the PAS-B domain contains the ligand binding domain (LBD) (Fig. 1c) 12,21 . AHR is part of a cytoplasmic protein complex containing HSP90, p23, and XAP2. Upon ligand binding, AHR dissociates from the complex and translocates to the nucleus where it forms a heterodimer with ARNT 22 . The AHR/ARNT complex subsequently binds to Xenobiotic Response Elements (XREs) in the genome to activate the transcription of target genes such as those coding for the phase I and II drug metabolizing enzymes CYP1A1 and CYP1B1 12 . At the protein level, the AHR Δe8-9 deletion results in an in-frame protein coding sequence with the loss of 84 amino acids (p. 303-387), disrupting the PAS-B domain  www.nature.com/scientificreports/ and the C-terminal part of the ligand-binding domain of the protein. The HSP90 and XAP2 protein binding domains are predicted to be affected, while the dimerization region of the protein likely remains intact (Fig. 1c). Among the urinary tract cancer patients in the Hartwig database, AHR alterations are mutually exclusive with FGFR3 alterations and PIK3CA mutations, although the significance is impacted by low sample numbers (FGFR3: Odds ratio = 0.24 p = 0.21, PIK3CA: Odds ratio = 0.30 p = 0.32). The tendency of mutual exclusivity with genes from the RTK-Ras-PI3K pathway strengthens the idea that AHR is an independent driver of urinary tract cancer, with potential convergence on the RTK-Ras-PI3K pathway. AHR alterations do not show mutual exclusivity with genes involved in chromatin remodeling or the tumor suppressors TP53 and Rb1 (Fig. 1d, Supplementary Table 4).
For tumors with AHR amplifications, AHR expression is increased compared to the non-amplified urinary tract cancer samples (Fig. 1e). Urinary tract cancers with AHR Δe8-9 showed increased expression of AHR target genes, which was not observed in the AHR amplified or AHR non-affected samples (Fig. 1f). No RNA sequencing data was available for the 2 patients with the AHR Q383H mutation.
The hAHR Q383H mutation affects ligand binding affinity and specificity. The mouse ortholog of AHR Q383H (Ahr Q377 ) can form hydrogen bonds with Ahr ligands in the ligand-binding domain and mutations that change this residue affect ligand binding affinities [23][24][25] . To examine if ligand binding affinities are also affected for the human AHR (hAHR), we measured the transcriptional activity of hAHR Q383H and hAHR WT in an XRE-luciferase reporter assay. Both the hAHR WT and hAHR Q383H expressing cells showed strong transcriptional activity with AHR ligands TCDD and B[a]P (Fig. 2a). Stimulation with B[a]P resulted in a higher induction of luciferase transcription for the hAHR Q383H mutant than for the hAHR WT , suggesting a difference in ligand affinity between both variants. The AHR antagonist CH-223191 reduced transcriptional activation in hAHR WT cells stimulated with TCDD or B[a]P ( Fig. 2a) 26,27 . Surprisingly, incubation of hAHR Q383H with CH-223191 resulted in transcriptional activation and no antagonism was observed when CH-223191 was combined with TCDD or B[a]P (Fig. 2a). To examine this altered ligand-binding specificity, we analyzed protein localization in RPE1 cells, which is a genetically stable near diploid cell line well-suited for microscopy. In the absence of exogenous ligands, both hAHR Q383H and hAHR WT proteins are mainly localized to the cytoplasm and both proteins show nuclear translocation upon incubation with TCDD or B[a]P (Fig. 2b, Supplementary Fig. 1). Incubation with the antagonist CH-223191 also resulted in nuclear translocation of hAHR Q383H , while the hAHR WT remains localized to the cytoplasm (Fig. 2b). The difference in responses between the hAHR Q383H and hAHR WT to AHR agonists and antagonists demonstrate differential sensitivity to AHR pathway modulation.
The hAHR Δe8-9 mutant leads to constitutive AHR pathway activation. To study the functional consequences of AHR Δe8-9 , we also created transgenic hAHR Δe8-9 RPE1 cells. In contrast with hAHR WT , the hAHR Δe8-9 protein is localized in the nucleus regardless of the presence or absence of a ligand (Fig. 2b, Supplementary  Fig. 1). In XRE-luciferase reporter assays, the degree of transcriptional activation by hAHR Δe8-9 was similar for all the conditions, irrespective of the presence of AHR ligands. Moreover, the AHR antagonist CH-223191 did not affect the transcriptional activation activity of hAHR Δe8-9 (Fig. 2a). Increased concentrations of hAHR WT and hAHR Δe8-9 constructs in the transfections resulting in increased activation of the luciferase reporter and revealed that in untreated conditions hAHR Δe8-9 induced higher transcriptional activation of the Luciferase reporter than hAHR WT (Fig. 2c).
To further characterize the consequences of hAHR Δe8-9 on a disease-relevant cell type, we created transgenic mouse bladder organoids expressing hAHR WT , hAHR Δe8-9 and a DNA binding deficient variant of hAHR Δe8-9 (hAHR Δe8-9 DBD ). We expected similar protein expression levels between the cell lines, because these were produced in an identical manner. However, we observed weak protein expression of hAHR Δe8-9 , but robust expression for the wildtype and the DNA binding deficient variant (Fig. 3a). These results suggest increased degradation of the DNA binding proficient hAHR e8-9 protein, possibly the result of negative feedback loops that depend on the DNA binding of AHR.
The hAHR WT and hAHR Δe8-9 bladder organoid lines were exposed to TCDD and together with untreated lines subjected to bulk RNA sequencing. Untreated hAHR Δe8-9 organoids differentially expressed several genes compared to the hAHR WT organoids, including the AHR target genes Cyp1a1, Cyp1b1, Ahrr, Gsta1, and Tiparp demonstrating constitutive activation of the AHR pathway (Fig. 3b,c). The expression pattern of up and down regulated genes is similar for hAHR Δe8-9 in untreated and TCDD treated conditions, confirming that loss of exons 8 and 9 leads to ligand-independent activation of the AHR pathway, which is in line with the upregulation of these genes in AHR Δe8-9 positive urinary tract cancers (Fig. 1f). Thus, AHR Δe8-9 leads to constitutive activation of the AHR pathway.
The hAHR WT organoids treated with TCDD also show upregulation of the canonical AHR target genes, but do not show the same downregulated genes (Fig. 3b). This observation indicates a different effect on the transcriptome between the constitutively active hAHR Δe8-9 and the 24 h TCDD stimulated hAHR WT condition. Moreover, overexpression of hAHR WT and hAHR Q383H without the addition of exogenous ligands already results in modest pathway activation when compared with control organoids that are not transgenic for hAHR (Fig. 3c). These observations are in line with the higher induction of luciferase transcription in the untreated condition for the hAHR WT and hAHR Q383H constructs compared to the control in the luciferase assay (Fig. 2a).
The hAHR Δe8-9 mutant induces cellular transformation. The bladder is a stratified epithelium, with stem cells that reside in the basal cell layer that support organ regeneration and renewal. Upon differentiation, the stem cells give rise to intermediate cells and luminal umbrella cells 28   www.nature.com/scientificreports/ when compared to organoids composed of differentiated cells 29 . We observed a mixture of cystic and compact organoids in the mouse bladder organoids expressing hAHR WT , indicating a heterogeneous population of differentiated and undifferentiated cells. No cystic organoids were observed in the organoids expressing hAHR Δe8-9 indicating a more basal stem-cell like phenotype (Fig. 3d).
We performed GO enrichment analysis on all differentially expressed genes in RNA-sequencing data of the mouse organoids to understand which processes are affected by hAHR Δe8-9 . Most outstanding is the downregulation of genes related to the extracellular matrix and cell periphery organization (Fig. 3e, Supplementary Table 5). Moreover, hAHR Δe8-9 expressing mouse organoids show downregulation of the differentiation markers Upk3a, Upk1a, and Krt18 and upregulation of the stem cell markers Wnt5a, and Krt1 when compared to hAHR WT organoids (Fig. 3f) [29][30][31] . Together, these observations indicate that hAHR Δe8-9 promotes a basal stem-cell like phenotype in bladder cells.
Because RNA-seq was performed on bulk cultures, the transcriptional changes induced by hAHR Δe8-9 may reflect a shift in the composition of the cell types towards a more homogeneous population, or the transformation of cells towards a novel phenotype. To discriminate between these scenarios, we performed scRNA-seq on the mouse bladder organoid lines. Dimensional reduction and unsupervised clustering revealed the presence of 7 clusters (Fig. 4a). Based on the genes enriched in the different clusters, cluster 0 represents a basal phenotype (characterized by Krt14, Trp63, Bcam and Agrn), cluster 1 represents an intermediate/luminal phenotype (characterized by Krt19, Krt18, Upk1b, Cldn4, Cldn7, Ceacam1 and Alcam), and clusters 3 and 5 represent cells that are in the S-phase and M-phase of the cell cycle, respectively (Fig. 4a,b, Supplementary Table 5). Together, these clusters represent a classical stem cell system where stem cells divide to give rise to new stem cells or to cells that differentiate. The vast majority of all cells of the control lines fall into these clusters and in cluster 4, which we were not able to link to a particular cellular phenotype. Strikingly, two clusters (clusters 2 and 6) were almost exclusively occupied by mouse bladder cells expressing hAHR Δe8-9 (Fig. 4a,c). In addition to the canonical AHR target genes such as Cyp1b1 and Tiparp, cluster 2 expresses basal stem cell markers such as Krt17, Wnt5a, Itga6, Wnt4 29 . In addition, cluster 2 is characterized by genes that are associated with urothelial cancer, such as: Htra1 32 , Cyb5r1 33 , Steap1 33,34 , Ptgs2 35,36 , and Trib3 37 (Fig. 4b). Cells in cluster 6 show upregulation of Dsp, Pkp1, Ppl, Jup and Krt13, which are involved in desmosome and intermediate filament cytoskeleton organization. This cluster has some overlap with markers expressed in the luminal-intermediate cluster (cluster 1) of the control cells, but with less apparent expression of the umbrella cell markers. This indicates that the cells in cluster 6 represent an intermediate cell type that however fails to differentiate towards luminal umbrella cells. Thus, the constitutive activation of the AHR pathway by hAHR Δe8-9 leads to a transformation of cells towards a less differentiated phenotype and the activation of genes linked to urothelial cancer.
Cells of a DNA binding deficient variant of the oncogenic hAHR Δe8-9 (hAHR Δe8-9DBD ) clustered among cells of the other control lines and apart from hAHR Δe8-9 cells (Fig. 4a). This shows that the cellular transformation of the hAHR Δe8-9 organoids depends on the transcriptional activation caused by AHR and is not the result of other transcription-independent effects caused by the deletion. To functionally explore if the transformed phenotype is accompanied by oncogenic properties, a soft agar growth assay was performed, a technique to measure cellular transformation in vitro. A higher number of colonies were counted for hAHR Δe8-9 organoids compared to hAHR WT organoids or control organoids demonstrating that hAHR Δe8-9 confers anchorage-independent growth to bladder cells (Fig. 4d).
Together, these observations demonstrate that transcriptional changes driven by the constitutively active hAHR Δe8-9 lead to a transformation of cells towards a cancerous phenotype that fails to differentiate and is able to grow independent of anchorage to the extracellular matrix.

Discussion
In this study, we show high prevalence of AHR alterations in urinary tract cancers and provide functional validations to support that these aberrations are oncogenic and drive urinary tract cancer. These findings are in line with mouse studies that have shown that AHR overexpression or pathway activation can lead to malignant transformation of epithelial cells [38][39][40] . A variety of molecular signaling pathways have been previously linked to AHR-mediated tumorigenesis whether or not in a cancer type-specific background and/or driven by AHR agonists 13,41 . This includes the biotransformation of hydrocarbons by CYP enzymes to mutagenic intermediates which can induce DNA adducts, cross-talk of AHR with other signaling pathways, and interaction of AHR with other binding partners than ARNT to promote transcription of non-canonical genes [42][43][44][45][46][47] . However, the underlying molecular signaling pathways driving the tumorigenesis of urinary tract cancers in patients with activating mutations in AHR remains unclear. The tendency of mutual exclusivity between AHR, FGFR3, and PIK3CA alterations may indicate convergence on the same pathway, although a fully independent parallel oncogenic pathway in bladder cancer driven by AHR cannot be excluded. The hAHR Q383H and hAHR Δe8-9 variants are almost uniquely detected in urothelial cancers, which suggests this tissue is particularly sensitive to deregulated AHR signaling. However, our results, nor information in literature, provides clues why AHR pathway activation is so specific for bladder cancer, so this will require further investigations.
We demonstrate that the hAHR Q383H mutation leads to increased ligand activation sensitivity which may lead to overactivity of the AHR pathway in the urinary tract, thereby driving tumorigenesis. This is mechanistically different from hAHR Δe8-9 which is constitutively active in a ligand-independent way. We presume that for hAHR Δe8-9 loss of the ligand binding domain leads to a conformational change that exposes a nuclear import signal that in the native protein is only exposed upon ligand binding dependent dissociation from the cytoplasmic complex. As for the hAHR Q383H mutation, increased sensitivity of the AHR pathway likely also holds true for AHR amplifications since overexpression of hAHR wt in different models shows a modest increase of AHR pathway activation compared to the controls. It is not known if the ligands that induce overactive AHR signaling   www.nature.com/scientificreports/ in the tumors with AHR amplified and AHR Q383H backgrounds are of environmental origin (like components in tobacco smoke) or have an endogenous source (like metabolites) 48 . Surprisingly, we observed hAHR Q383H activation upon treatment with the AHR antagonist CH-223191. This contrasts with a recent study that demonstrated reduced viability upon treatment with CH-223191 of the bladder cancer cell line KMBC2 harboring the AHR Q383H mutation 11 . A possible explanation for the apparent discrepancy in results may lie in the different types of experiments that were conducted. The reduced viability of the KMBC2 cell line after CH-223191 incubation may also be independent of alterations in AHR pathway activity. Patients that harbor AHR activating mutations could potentially benefit from AHR targeted therapies 49 . Here we show that the mode of action of AHR activation differs between the different mutations, which implies that tailored therapies depending on the underlying mutational event are required. Most classic AHR antagonists function by interference with the ligand-binding domain, but this domain is not targetable for the AHR Δe8-9 mutant as this domain is deleted. Therefore, functional screens to identify specific AHR Δe8-9 targeting compounds or the identification of essential downstream activated processes could be a next step towards the identification of novel treatment strategies for selected urinary tract cancer patients in the context of precision medicine.

Material and methods
Driver gene status urinary tract cohort. Variant detection and driver likelihood status are based on the Hartwig database 14 . Structural variation was detected and classified with the GRIDSS, PURPLE, LINX algorithms, using pairs of matched tumour/normal BAM files as input 15,16 . Driver likelihood is introduced to select for a sample specific likelihood of each gene based on the type of variant and taking mutation load per sample into account. To select for affected genes in the urinary tract cancer samples, the driver likelihood score is set to > 0.8 and detected gene fusions with a high impact are included. Visualized are the genes affected in more than 12.5% of the samples (top 12 genes). Mutually exclusivity is pairwise calculated with Fisher Exact Test based on Odd Ratio cut-offs as described in Gao et al. 19 .

Gene constructs.
A plasmid containing the hAHR WT sequence was purchased from Origene (RC209832).
The Q383H point mutation was introduced with site-directed mutagenesis with forward primer cattgtaactcacagaccactaacagatg and reverse primer gttagtggtctgtgagttacaatgatataatc. The hAHR WT , hAHR Q383H , and hAHR Δe8-9 sequences were cloned in a lentiviral plasmid (Addgene #52961) and subsequently to pcDNA3.1 vector with the In-fusion HD Eco-dry cloning (Takarabio). Primers for lentiviral plasmid are: forward N-flag caggaccggttctaggatatcgccaccatggattacaaagacgatgacgataagaacagcagcagcgccaac, forward caggaccggttctaggatatcgccaccatgaacagcagcagcgcc, reverse ttgttgcgccggatcgcaggaatccactggatgtcaaatcag, reverse C-flag ttgttgcgccggatcgcttatcgtcatcgtctttgtaatccaggaatccactggatgtcaaatcag, forward deletion tggttgtgatgccaaagatgaggaaggaacagagca and reverse deletion gttccttcctcatctttggcatcacaaccaatagg. Primers cloning of AHR sequences to pcDNA3.1 vector are forward taccgagctcggatcatatcgccaccatgaacag and reverse gatatctgcagaattttacaggaatccactggatgtcaaat. For the DNA binding deficient variant of AHR Δ8-9 , mutations were introduced to substitute the amino acids Histidine 39 and Arginine 40 with Alanines 50 .
Lentivirus production and transduction. Lentivirus particles containing hAHR constructs were produced by transient calcium phosphate transfection of lentiviral transfer, packaging, and envelope plasmids into HEK293T cells. Virus particles were harvested 48 h posttransfection and concentrated with Lenti-X Concentrator (Takarabio) according to the manufacturer's directions. Concentrated virus was resuspended in 250-500 µl Advanced DMEM/F-12 supplemented with HEPES, Glutamax, and 1% Pen-Strep. The day before transduction, RPE-1 cells were seeded in 6-well plates and mouse bladder organoids were seeded as single cells onto a layer of solidified matrigel in 96 well plates. Organoid culture medium was supplemented with 10 μM Y-27632 to enhance survival of single cells. The day after, RPE-1 cells and mouse bladder organoids were transduced with a dilution series ranging from 10-150 µl concentrated lentivirus. To enhance transduction efficiencies, the medium was supplemented with 8 µg/ml polybrene. After transduction, RPE1 cells and the mouse bladder organoids were placed on 10-15 µg/ml and 1 µg/ml puromycin selection, respectively. We continued with the transgenic lines that showed approximately 50% survival upon selection. RNA-sequencing. The RNA from the transgenic mouse bladder organoids was isolated with Trizol (Ther-moFisher) according to the manufacturer's instructions. RNA-seq libraries were prepared using TruSeq Stranded Total RNA Library Prep Kit (Illumina) according to the manufacturer's protocol. RNA-seq libraries were pooled and sequenced on a NextSeq2000 (Illumina) in 1 × 50 bp single end mode. RNA sequencing reads were aligned against mouse reference genome GRCm38 using STAR and the number of reads mapping to genes was counted using featureCounts all by using a custom in-house pipeline (https:// github. com/ UMCUG eneti cs/ RNASeq-NF). The Bioconductor package DESeq2 was used to normalize raw read counts and to perform differential gene expression analysis with apeglm shrinkage 52,53 . The analyses were performed with significant (P.adjust < 0.05) and differentially expressed (Log2FoldChange > 2.5) genes with exception of luminal and basal marker analyses where smaller differences were included (Log2FoldChange > 1.5). The Bioconductor package clusterProfiler and Revigo were used for GO enrichment analysis 54,55 . scRNA-seq. The transgenic mouse bladder organoids were cultured in 2D and dissociated to single cells using a 10 min incubation with TryLE. Cells were sorted into 384-well capture plates (1 plate per condition) and the scRNA library preparation and sequencing were performed according to the SORT-seq protocol by Single Cell Discoveries B.V. 56 . For all single cells, reads were aligned to the mouse reference genome GRCm38 and Sort-seq read counts were filtered to exclude reads with identical library-, cell-and molecules 56 . With the Seurat R package, low quality cells were removed by a cut-off of 10,000 transcripts per cell and the data was normalized and scaled 57 . The top 2000 most variable genes in the dataset were identified and used for principal component analysis to determine dimensionality and clustering of the dataset. Cluster gene markers were detected using a Wilcoxon rank sum test between each cluster and the rest of the cells in the dataset with a bonferroni correction for multiple testing.
Soft agar assay. 3% agarose (REF11388991001) was dissolved and autoclaved in 100 ml EBSS. One volume of melted 3% agarose was mixed with four volumes of Advanced DMEM/F12 to obtain a 0.6% solution. This mixture was added to 6-well plates (1.5 ml/well) in which the gels were allowed to solidify. Subsequently, Tryple was used to prepare single-cell suspensions of the mouse bladder organoids. The cells were counted and for each condition, 2.5 ml of cell suspension was prepared at a concentration of 1 × 10 e4 cells/ml in complete medium. The cell suspension was mixed 1:1 with a warm 0.6% agarose solution to get 0.5 × 10 e4 cells/ml in 0.3% agarose solution. Per well 1.5 ml agarose/cell mixture was plated. The next day 300 µl of the medium was added and the cells were refed every 2-4 days. After 3 weeks the cells were stained with nitroblue tetrazolium chloride solution and pictures were made. Colonies were counted using ImageJ.

Data availability
The bulk RNA-sequencing and the scRNA-seq data of the mouse bladder organoids have been deposited in ENA with the accession code PRJEB49233.