Single‐cell transcriptomes of mouse bladder urothelium uncover novel cell type markers and urothelial differentiation characteristics

Abstract Objectives Much of the information to date in terms of subtypes and function of bladder urothelial cells were derived from anatomical location or by the expression of a small number of marker genes. To have a comprehensive map of the cellular anatomy of bladder urothelial cells, we performed single‐cell RNA sequencing to thoroughly characterize mouse bladder urothelium. Materials and methods A total of 18,917 single cells from mouse bladder urothelium were analysed by unbiased single‐cell RNA sequencing. The expression of the novel cell marker was confirmed by immunofluorescence using urinary tract infection models. Results Unsupervised clustering analysis identified 8 transcriptionally distinct cell subpopulations from mouse bladder urothelial cells. We discovered a novel type of bladder urothelial cells marked by Plxna4 that may be involved with host response and wound healing. We also found a group of basal‐like cells labelled by ASPM that could be the progenitor cells of adult bladder urothelium. ASPM+ urothelial cells are significantly increased after injury by UPEC. In addition, specific transcription factors were found to be associated with urothelial cell differentiation. At the last, a number of interstitial cystitis/bladder pain syndrome–regulating genes were found differentially expressed among different urothelial cell subpopulations. Conclusions Our study provides a comprehensive characterization of bladder urothelial cells, which is fundamental to understanding the biology of bladder urothelium and associated bladder disease.


| INTRODUC TI ON
The bladder urothelium is a stratified epithelium that serves as a crucial barrier between the blood and urine. Beside barrier function, the bladder urothelium also has a range of other important functions such as sensory function and immune function. 1 Investigating the bladder urothelial cell types, specific cell markers, signalling receptors and genes will help us to learn more about the relationship between bladder urothelial cells and bladder diseases including BPS/IC, bladder cancer and urinary tract infections (UTIs).
Since the pathophysiology and treatment strategies of bladder diseases are highly relevant to the bladder urothelium, it is vital to understand the biology and physiology of bladder urothelium. Much of the information to date in terms of subtypes and function of bladder urothelial cells were derived from anatomical location or by the expression of a small number of marker genes. Common technique to study the bladder urothelium is using cell culture preparations and enzyme digested or surgically isolated tissue. 5,6 However, it seems difficult, or even impossible, to study cell type-specific characteristics and functions of bladder urothelial cells by these traditional methods. The heterogeneity and hierarchy of bladder urothelial cells have not been well revealed.
Single-cell RNA sequencing is an advanced technology that can analyse mRNA transcripts from thousands of individual cells simultaneously, and so greatly facilitates the studies on cell function under normal physiological activities as well as during disease processes. 7 By utilizing this powerful tool, the cellular hierarchy and developmental trajectory of multiple organs have been discovered, such as lung, kidney and intestine. [8][9][10] Although previous studies have mapped all major mouse organs, including mouse bladder, 11,12 high-throughput single-cell RNA sequencing of bladder urothelial cells have not been reported so far. The present study aims to investigate the unsupervised classification and function of mouse bladder urothelial cells precisely by single-cell RNA sequencing. The transcriptomic map of bladder urothelial cells will provide a resource for studying bladder urothelial cell types, specific cell markers, signalling receptors and genes that can help us to learn more about the relationship between bladder urothelial cells and bladder diseases.

| Bladder urothelial cell isolation
Ten-week-old female C57BL/6 mice were ordered from the animal experimental centre of Shandong University. All mice were maintained in the specific pathogen-free animal facility of our institution. Five mice were sacrificed by cervical vertebra dislocation and bladder was removed. The bladder was cut open through the bladder neck with sharp scissor and was nailed to the Sylgard-coated plate (Supplement Figure S1). The lumen of the bladder was exposed using four pins to form a parallelogram. The bladder was incubated in dispase II solution (04 942 078 001, Sigma) at 37°C for 1.5 hours. The bladder was then enlarged to achieve a maximum stretch (Supplement Figure S1E). After dispase II solution was discarded, 10-15ml PBS-EDTA (10mmol/L, CZ0026, Leagene) was added to Sylgard-coated plate and shaking for 1.5 hours. The PBS-EDTA was collected and centrifuged at 300g for 5 minutes. The

| Single-cell RNA sequencing
Single-cell RNA-seq libraries were prepared with Chromium Single cell 3' Reagent (v3) Kits according to the manufacturer's protocol. Singlecell suspensions were loaded on the Chromium Single Cell Controller Instrument (10 × Genomics) to generate single-cell gel beads in emulsions (GEMs). Briefly, 10 5 -10 6 single cells were suspended in calcium-and magnesium-free PBS containing 0.04% weight/volume BSA. After generation of GEMs, reverse transcription reactions were engaged barcoded full-length cDNA followed by the disruption of emulsions using the recovery agent and cDNA clean up with DynaBeads Myone Silane Beads (Thermo Fisher Scientific). cDNA was then amplified by PCR with appropriate cycles which depend on the recovery cells. Subsequently, the amplified cDNA was fragmented, end-repaired, A-tailed, index adaptor ligated and library amplification.
Then, these libraries were sequenced on the Illumina sequencing platform (HiSeq X Ten) and 150 bp paired-end reads were generated. The single-cell RNA-seq data are available at GEO: GSE163029.

| Quality control
The Cell Ranger software pipeline (version 3.0.0) provided by 10 × Genomics was used to demultiplex cellular barcodes, map reads to the genome and transcriptome using the STAR aligner and downsample reads as required to generate normalized aggregate data across samples, producing a matrix of gene counts versus cells. We processed the unique molecular identifier (UMI) count matrix using the Seurat R package (version 2.3.4; https://satij alab.org/seura t/). To remove lowquality cells and likely multiplet captures, which is a major concern in microdroplet-based experiments, we apply a criterion to filter out cells with UMI/gene numbers out of the limit of mean value ± 2-fold of standard deviations assuming a Guassian distribution of each cells' UMI/gene numbers. Following visual inspection of the distribution of cells by the fraction of mitochondrial genes expressed, we further discarded low-quality cells where > 10% of the counts belonged to mitochondrial genes. After applying these QC criteria, 18,917 single cells and 31,053 genes in total remained and were included in downstream analyses. Library size normalization was performed in Seurat on the filtered matrix to obtain the normalized count.

| Dimensional reduction, clustering and cell type identification
Top variable genes across single cells were identified using the method as previously described. 7 Briefly, the average expression and dispersion were calculated for each gene, and genes were subsequently placed into 8 bins based on expression. Principal component analysis (PCA) was performed to reduce the dimensionality on the log transformed gene-barcode matrices of top variable genes.
Cells were clustered based on a graph-based clustering approach and were visualized in 2-dimension using t-SNE. Likelihood ratio test that simultaneously test for changes in mean expression and in the percentage of expressed cells was used to identify significantly differentially expressed genes between clusters. Here, we use the R package SingleR, a novel computational method for unbiased cell type recognition of scRNA-seq, with reference transcriptomic data sets 'Immgen' to infer the cell of origin of each of the single cells independently and identify cell types. 14

| RNA velocity analysis
We performed RNA velocity analysis using the R package velocyto.R v0.6. 15 The RNA velocity was calculated on the basis of spliced and unspliced transcript reads and estimated using gene-relative model.
The resulting velocity estimates were projected onto the t-SNE embedding obtained in Seurat.

| SCENIC analysis
We used SCENIC (Single Cell rEgulatory Network Inference and Clustering) for the simultaneous reconstruction of gene regulatory networks (GRNs) and the identification of stable cell states. 16 SCENIC contains three main steps, including co-expression analysis, To evaluate the specific regulon for each cell type, we calculated the regulon specificity score which is based on the Jensen-Shannon divergence, a measure of the similarity between two probability distributions. To systematically characterize the combinatorial patterns, we compared the atlas-wide similarity of regulon activity scores of every regulon pair based on the Connection Specificity Index (CSI). [17][18][19] The CSI for all regulons was calculated with scFunctions package.

| Pseudotime analysis
We determined the developmental pseudotime with the Monocle2 package. 20 The data, previously scaled and clustered by the Seurat tool, were loaded into a monocle object with default parameters. We obtained variable genes with Monocle2 and ordered the cells onto a pseudotime trajectory based on the union of highly variable genes obtained from all cells.

| HE stain and immunofluorescence
The mice bladder tissues were kept in 10% formalin and embedded in paraffin. Paraffin sections were stained for haematoxylin and eosin.

| Single-cell profiling and unbiased clustering of mouse bladder urothelial cells
Single bladder urothelial cell suspension was obtained from the mouse bladder urothelium using trypsin digestion (Supplement Figure S1). Trypan blue staining and Calcein-AM/PI staining showed that the per cent of alive cell is more than 90% (Supplement Figure   S2). Mouse bladder urothelial cell types were catalogued in an unbiased manner using droplet-based single-cell RNA sequencing ( Figure 1A). Using stringent quality controls (Supplement Figure S3), 18 917 bladder urothelial cells were further analysed.
Unsupervised clustering analysis identified 8 distinct cell clusters consisting of as few as 127 cells to as many as 6117 cells per cluster (the clusters were censored for a minimum of 100 cells) ( Figure 1B and 1C).

| Classification of urothelial cells based on cell type-specific marker genes
To further investigate the identity of each cell cluster, we performed differential gene expression analysis ( Figure 1D Figure 2B). In addition, cluster 4 has higher expression of genes including Foxa1, Gata3 and Grhl3, which are typically associated with differentiated urothelial cell phenotypes. 21,23 To further investigate the differences between these two clusters, we identified their differential genes. Five genes (Ivl, c-Jun, Klf6, Psca and Rab11fip1) increased ≥ 2-fold change in cluster 4 and two genes increased ≥ 2-fold change in cluster 3 (Supplement Figure S4). Cluster 5 and cluster 6 were classified into basal cell category because of expression of keratin-5, ITGb4, Trp63 and Shh. In addition, cluster 5 was enriched for DNA helicase gene markers ( Figure 2C) and cluster 6 was enriched for cell cycling regulatory gene ( Figure 2D). Further analysis also showed that was enriched for tricarboxylic acid (TCA) cycle regulating genes, suggesting those cell subtypes possess active metabolism function ( Figure 2E). Cluster 5 and cluster 6 were named B2 and B3, respectively. Cluster 7 was not enriched with any of the known subtype-specific markers except for Upk IIIb. Cluster 7 was named as Plxna4 + urothelial cell according to the most specific gene-Plxna4 ( Figure 2E). Altogether, the present single-cell transcriptome atlas provides a molecular definition of 8 bladder urothelial cell types.

| Reconstructing the developmental trajectory of mouse bladder urothelial cells
It is usually believed basal cells can divide and produce superficial cell daughters. To further understand the relationship among different types of bladder urothelial cells, we utilized reverse graph embedding to generate a trajectory plot ( Figure 3A) and perform cell trajectory analysis by using pseudotime analysis ( Figure 3B). Cluster 6 is the starting point of the trajectory, and followed by cluster 5 and cluster 0. BTS cells were in the middle position which suggested its transitional state. Superficial cell (cluster 3 and cluster 4) and Plxna4 + urothelial cells appeared as the two ends of the trajectory branch 1. It suggested cluster 6 might be the progenitor cell and Plxna4 + urothelial cell might play as a superficial cell. To further analyse the genes in terms of the changes in pseudotime analysis, we clustered genes via a pseudotemporal expression pattern. The top 50 genes that vary as a function of pseudotime were clustered, as shown by the heat map ( Figure 3C). We found that Malat1, Ftl1 and Tpt1 were highly expressed at the early stage of developmental trajectory, while Ly6d, S100a6, Fau, Gstm1 and Tmsb4x were highly expressed at the end stage of developmental trajectory. These genes were suggested to regulate cellular pluripotency or identified as growth inducible genes and might be linked to the process of urothelial cell proliferation and differentiation.

| RNA velocity analysis revealed the developmental lineages
RNA velocity was analysed to study the developmental lineage of    Figure 5A). So we named this group as the ASPM + urothelial cells. Genome-wide analysis has suggested ASPM as a possible stem/progenitor cell marker. 26,27 Immunofluorescence results showed that few ASPM + cell was found in healthy mouse urothelium ( Figure 5B). However, the numbers of ASPM + basal cells were significantly increases in basal layer after UPEC injury ( Figure 5C).

| The expression patterns of transcription factors in bladder urothelial cells
The maintenance of cell identity involves the coordinated action of many regulators, among which transcription factors have been long recognized to play a central role. 18 SCENIC was used to identify transcription factors in different bladder urothelial cell subpopulations and to predict essential regulators for urothelial cell type. Our network analysis identified top 3 specific regulons for each cluster ( Figure 6A). Mybl1, Rad21 and pole4 may cooperate in regulating ASPM + cell differentiation. Nfib, sp5 and Gli1 may cooperate in regulating Plxna4 + cell differentiation. t-SNE plot provides additional support that the activities of these regulons are highly specific to ASPM + cell subpopulations ( Figure 6B) and Plxna4 + cell subpopulations ( Figure 6C). This analysis provides an opportunity to systematically identify critical regulators for urothelial cell identity.
Transcription factors often work in combination to coordinate gene expression levels. To systematically characterize the combinatorial patterns, we compared the atlas-wide similarity of regulon activity scores of every regulon pair based on the CSI. We identified co-expression modules associated with 165 transcription factors, ranging in size from 10 to 2172 target genes, with a median size of 50 genes. The regulons which have similar functions were categorized into a module. Strikingly, regulons were organized into 4 major modules. The difference of CSI activity in each cluster was not substantial. Module M2 contains 12 regulators, in which Foxa1, Grhl3, Foxq1, Grhl1, Htatip2, Ikzf2, Irf8 and Creb3l2 were considered to play important role in cell differentiation ( Figure 7A). The activity of those regulons was higher in basal-like cell subpopulations (clusters 0, 5 and 6) but lower in S2 cell subpopulations (cluster 4) ( Figure 7B).
On the contrary, the regulon specificity score (RSS) were higher in S2 cell (cluster 4) and BTS cell (cluster 1) ( Figure 7C), which was similar to gene expression level ( Figure 7D). It indicated that combinations of these regulons play important roles in urothelial cell differentiation.
The RSS and gene expression of transcription factors, rather than activity of regulons, may more accurately reflect its identity. In order to investigate the function of these 8 regulons, their target genes were searched by KEGG pathways database. It was found that the target genes of these 8 regulons are associated with epidermal barrier formation, such as regulation of actin cytoskeleton, bacterial invasion of epithelial cells, endocytosis and lysosome (Supplement Figure S5).

| Cell type-specific expression of IC/BPSassociated genes
Since the pathogeneses of IC/BPS is theorized to be associated with impairment of the urothelium, we asked whether specific bladder urothelial cell subtypes are enriched for IC/BPS-associated genes.
We investigated the expression profiles of genes known to contribute to pathogeneses of IC/BPS. [28][29][30] We found a number of IC/BPSassociated genes differentially expressed among different subtypes  be analysed to define cell subtypes. In addition, the discrepancies between the present study and the previous studies as to the subtype of bladder urothelial cells could be due to the fact that the mouse bladder urothelium used in our study was in static state, rather than post-injury. Microbial, chemical or processing/harvesting urothelial stimuli/injury can shift the bladder urothelium from near-quiescence to a highly proliferative state, which can favour the lineage tracing study for the progenitor cells or stem cells.

| D ISCUSS I ON
But it also broke the homeostasis of bladder urothelium and may impact the identities of bladder urothelial cells, which will provide the misleading information about the cell types. Harvesting and processing urothelial cells may also impact the transcriptomic profiles of the urothelial cells which could lead to spurious expression of genes.
As a crucial barrier between the blood and urine, the bladder urothelium has to recover quickly after injury. Conflicting evidence exists as to the location of progenitor cells responsible for urothelial restoration after injury. Previously, it is generally accepted further analyse the specific gene expressed in the B3 subpopulation and found that ASPM was the most highly expressed genes in F I G U R E 8 Heat map showing basal cell was enriched for IC/BPS regulating genes B3 compared to other subpopulations. So we named this group as the ASPM + urothelial cells. Genome-wide analysis has suggested ASPM as a possible gastric stem/progenitor cell marker. 26,27 In the UTI model we established, we found that ASPM + urothelial cells were significantly increased after acute injury, suggesting that ASPM + urothelial cells are involved in the bladder urothelial regeneration. These evidence indicated that ASPM may also be a marker of progenitor cells in bladder urothelium. BTS subpopulation was regarded as the transitional state. We found that both BTS cells and I cells are responsible for the development of superficial cells (S1 and S2). These results are in line with previous lineage studies that both intermediate cells and basal cells can divide and produce superficial cells. 35,38 Bladder urothelial cells are closely related to the host response to UTIs and they are also a major source of proinflammatory cytokines and chemokines following bacterial infection. [39][40][41] By using the scRNA-seq of bladder urothelium, we discovered a novel bladder urothelial cell type characterized by the specific expression of Plxna4. IF results indicated that Plxna4 + urothelial cells were located on the apical layer of bladder urothelium in mouse, rat and human beings. Plexins are cell surface receptors that play an important role in axon guidance, cell migration, wound repair, mechanosensation, immune-cell regulation and cancer progression. 42-4 4 Until now, the expression and function of Plxna4 in bladder urothelium has not been reported. Wen et al demonstrate that Plxna4 is required for optimal cytokine production upon Tolllike receptor (TLR) stimulation and bacterial challenge, suggesting a critical role of Plxna4 in mediating the host response to infection. 45,46 In addition, we found that Plxna4 + urothelial cells were highly enriched for WFDC1, which is a key modulator of inflammatory and wound repair responses. 47 Hence, it is hypothesized that these Plxna4 + urothelial cells may play a role in the physiological process of host response to UTIs. It is hypothesized that these  In the present study, the bladder urothelium of human beings was not used for the reason that it is impossible to get normal bladder urothelium from healthy volunteers. Since fresh cell suspensions are needed for the single-cell RNA sequencing, bladder samples from the body donated take some time and this will greatly impact the cell activity and cell transcriptome. Another possible way to get the bladder urothelium is to use the bladder samples from cystectomy. However, patients undergoing cystectomies are usually diagnosed with muscle-invasive bladder cancer. The bladder urothelium from these patients cannot be defined as normal bladder urothelium, especially if patients have been previously treated intravesical bacillus Calmette-Guérin (BCG) or cytotoxic chemotherapy. Cystoscopic bladder urothelial biopsies could also be considered as another option to obtain normal human bladder urothelium. However, these biopsies are typically small and the number of urothelial cells obtained is low. These biopsies may lose some urothelial cell subpopulations, especially the basal cells.
Hence, after careful deliberation, we decided to only perform the single-cell sequencing from mouse bladder urothelial cells. Another limitation of this study is that there is no way to prove that these cluster profiles were not somehow contributed or affected by the way the cells were harvested/processed.
In conclusion, our study provides the most comprehensive information on the cell types of mouse bladders urothelium and opens up a range of possible avenues for future research. This work provides an important reference providing unprecedented and granular insights into the transcriptional landscape of bladder urothelial cells.
This information is fundamental to understanding bladder biology and bladder diseases.

ACK N OWLED G EM ENTS
This work was technically supported by OE Biotech Co., Ltd.
(Shanghai, China). This work was supported by grants from the National Natural Science Foundation of China (81970661), the National Natural Science Foundation of China (81900637) and the Tai Shan Scholar Foundation to Benkang Shi.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are avail-