Single-cell RNAseq and longitudinal proteomic analysis of a novel semi-spontaneous urothelial cancer model reveals tumor cell heterogeneity and pretumoral urine protein alterations

Bladder cancer, one of the most prevalent malignancies worldwide, remains hard to classify due to a staggering molecular complexity. Despite a plethora of diagnostic tools and therapies, it is hard to outline the key steps leading up to the transition from high-risk non–muscle-invasive bladder cancer (NMIBC) to muscle-invasive bladder cancer (MIBC). Carcinogen-induced murine models can recapitulate urothelial carcinogenesis and natural anti-tumor immunity. Herein, we have developed and profiled a novel model of progressive NMIBC based on 10 weeks of OH-BBN exposure in hepatocyte growth factor/cyclin dependent kinase 4 (R24C) (Hgf-Cdk4R24C) mice. The profiling of the model was performed by histology grading, single cell transcriptomic and proteomic analysis, while the derivation of a tumorigenic cell line was validated and used to assess in vivo anti-tumor effects in response to immunotherapy. Established NMIBC was present in females at 10 weeks post OH-BBN exposure while neoplasia was not as advanced in male mice, however all mice progressed to MIBC. Single cell RNA sequencing analysis revealed an intratumoral heterogeneity also described in the human disease trajectory. Moreover, although immune activation biomarkers were elevated in urine during carcinogen exposure, anti-programmed cell death protein 1 (anti-PD1) monotherapy did not prevent tumor progression. Furthermore, anti-PD1 immunotherapy did not control the growth of subcutaneous tumors formed by the newly derived urothelial cancer cell line. However, treatment with CpG-oligodeoxynucleotides (ODN) significantly decreased tumor volume, but only in females. In conclusion, the molecular map of this novel preclinical model of bladder cancer provides an opportunity to further investigate pharmacological therapies ahead with regards to both targeted drugs and immunotherapies to improve the strategies of how we should tackle the heterogeneous tumor microenvironment in urothelial bladder cancer to improve responses rates in the clinic.

Introduction combines two pro-oncogenic factors, overexpression on Hgf and an activating mutation in Cdk4 (R24C). According to the Pathology Atlas of the Human Cancer Transcriptome high mRNA expression of HGF or CDK4 in urothelial bladder cancer correlates with a worse overall survival compared to patients with low expression of these genes (P-scores: HGF 0.021 and CDK4 0.016) [29]. Moreover, moderate protein expression for both genes is observed in urothelial cells as well as other cell types of the urinary bladder [29].
In urothelial bladder cancer, the alterations of the Hgf/c-met axis are commonly found in high-risk tumors [30][31][32][33][34] as well as in animal models [35,36]. The activating mutation R24C of Cdk4 allows cells to overcome senescence and to acquire a proliferation phenotype by promoting the aberrant phosphorylation of Rb which inactivates its tumor suppressor function. The expression of the mutated Cdk4 (R24C) results in spontaneous development of epithelial tumors or sarcomas in mice at an older age (>8 months) or melanoma upon skin exposure to carcinogens [37][38][39]. Cell-cycle dysregulation of urothelial cancer in humans is commonly identified as either loss of function of RB [40][41][42] along with CDK4 upregulation [43][44][45][46] or genetic alterations of CDKN2A that usually hamper the function of its transcript p16 which is the inhibitor of CDK4 [47,48], herein modelled in mice by overactivation of Cdk4 due to a point mutation [37,49,50]. Additionally, scattered mutations were introduced to the developing urothelial tumors by administration of the carcinogen OH-BBN [22].
The Hgf-Cdk4 R24C strain has been previously used to study melanoma-genesis upon exposure to carcinogenic stimuli [49,51], but is susceptible to spontaneous melanoma and to a lesser extent hepatocellular carcinoma development at an average age older than 8 months [52]. Moreover, this strain has previously been validated as potential model of spontaneous conjunctival or uveal melanoma, however there was no evidence of such spontaneously occurring pathology [53].
Herein, we subjected the animals to experiments of carcinogen induction or as healthy controls during the first 8 months of their life without observing high rates of spontaneous malignancy development, while using fully developed animals. Despite the presence of the genetic modifications in all cells, the carcinogen exposure in our study yielded bladder tumors that were histologically identified as urothelial carcinoma with squamous differentiation. The simultaneous overexpression of Hgf and increased activity of Cdk4 (R24C) led to accelerated progression of OH-BBN induced urothelial tumors compared to Cdk4 (R24C) alone. Using this novel model, we characterized it in comparison to commonly used models (OH-BBN exposed C57BL/6 mice and the syngeneic MB49 model) using transcriptomics, proteomics and assessed response to immunotherapy in both, male and female mice.

Carcinogen induction
Mice were housed in groups of five per cage and N-Butyl-N-(4-hydroxybutyl) nitrosamine (OH-BBN) (Sigma-Aldrich cat. B8061) was added into drinking water at the final concentration of 0.05%. The drinking solution was freshly made and replaced twice weekly for 10 weeks. During this period the animals were weighed and examined twice weekly. The monitoring was performed every 24 hours after establishment of advanced palpable tumors and hematuria development. The humane endpoint was determined as 10 days of visible light hematuria or 5 days of complete hematuria or 15% body weight loss. The two experimental endpoints that were used to sacrifice the mice and perform several studies on the bladder tissue were 1) 10 weeks of carcinogens exposure and 2) palpable bladder tumor and humane (MIBC) endpoint.

Bladder measurement and histopathology analysis
At 10-weeks carcinogen endpoint OH-BBN (female n = 5, male n = 6) and at MIBC endpoint (females n = 6, males n = 3) mice were sacrificed, bladders were harvested and emptied of urine, measured with a caliper and the ellipsoid formula "(4/3) x π x r1 x r2 x r3" was used to determine their volumes. Bladders as well as lungs and kidneys were fixed in Formaldehyde 4% buffered (Q Path1 cat. 11699455) and embedded in paraffin. Sections of 4 μm thick were cut and stained with HE. The sections were examined by a pathologist for determining the presence and characteristics of malignant tissue, using healthy mouse tissue as reference in a genotype-, age-and sex-matched manner. The same criteria as for staging human bladder cancer were used. Digital microscopy images were acquired with a 3.3 M pixel CMOS camera (Olympus SC30, Olympus corp. Tokyo, Japan) and the image acquisition software CellSense v1.18 at 10x (200 μm scale bar), 20x (100 μm scale bar) and 40x (50 μm scale bar) magnifications.

Generation of murine Hgf-Cdk4 R24C urothelial cancer cell lines
After 10 weeks of OH-BBN exposure mice were maintained while given only normal drinking water for one (n = 2 mice), five (n = 2 mice) or 10 (n = 1 mouse) weeks. At the designated time points mice were sacrificed, bladders were resected and digested under sterile conditions with Liberase TL (Roche, 2 ml PS/2.5 mg cat. 291959) for 15 min. Single cell suspensions were acquired via filtering in 70-μm cell strainers (Sigma cat. 431751) and seeded in T75 tissue culture treated flasks (Sarstedt) in RPMI 1640 medium supplemented with Glutamax (Thermo Fisher cat. 61870-010), 10% FBS (Gibco cat.11550356), Penicillin/ Streptomycin (Thermo Fisher cat. 15140-122), and cultured at 37˚C, 5% CO2. Adherent cell fractions were separated from floating cells and cell cultures maintained for additional 16 weeks. Single cell colonies were isolated by limiting dilution in 96 well plates, expanded and mycoplasma tested. For determination of the doubling time cell lines were seeded in triplicates and cell counts acquired with an automated cell counter (Bio-Rad TC20) during the first four days were applied to the function y = a2Tx-b (a, b = parameters, T = doubling time). Additionally, CDK4 (R24C) mutation was confirmed by restriction enzyme digestion as previously described [26] with HindIII-HF endonuclease (New England Biolabs, US-MA, cat. R3104S).

qRT-PCR of cell lines and urothelium
The expression of selected genes was analyzed in six new Hgf-Cdk4 R24C cell lines, MB49 cells and urothelium scrape from healthy Hgf-Cdk4 R24C bladder. Briefly, RNA was isolated using the Total RNA Purification Kit (Norgen Biotek cat. 17200) and DNase column treatment (Norgen Biotek cat. 25710) was performed according to the manufacturer's instructions. The RNA concentration was quantified using NanoDrop Spectrophotometer (Thermo Fisher Scientific) and 1.5 μg of RNA was reverse-transcribed into cDNA in a 30 μl reaction using the iScript RT kit (Bio-Rad cat. 170-8890) according to the manufacturer's instructions. Gene expression was determined using qRT-PCR with Sybr Green detection and PowerUp Master mix (Thermo Fisher cat. A25742) on a CFX384 Touch detection system (Bio-Rad). The primer pairs (S1 Table) were designed by us and ordered from Integrated DNA Technologies. The qPCR conditions consisted of 2 minutes initial activation at 50˚C, 2 minutes denaturation at 95˚C, followed by 40 cycles of 20 seconds denaturation at 95˚C, 20 sec annealing at 55˚C and 60 sec of extension at 72˚C. After the completion of all cycles a final extension step for 10 minutes at 72˚C was applied and melting curves were generated.

Flow cytometry analysis of cell lines
Cells were harvested, washed in buffer (PBS containing 1.5% BSA) and antibody dilution 1:100 in 50 ul buffer were added. Incubation of 30 minutes in 4˚C in the dark, was followed by 3 washes and 30 000 events were acquired in a CytoFLEX flow cytometer (Beckman Coulter) with a three-laser configuration (408, 488 and 640 nm). For the analysis of IFNγ induced stimulation, the 6 new cell lines as well as MB49 cells were cultured with or without IFNγ (10ng/ml) for 24h and expression of selected immune recognition cell surface markers (S2 Table) was analyzed by flow cytometry.

The orthotopic MB49 model
The tumorigenic cell line MB49 (Mouse Bladder-49) was kindly provided by Dr K. Esuvaranathan, National University of Singapore, Singapore. The original stock vial was expanded in culture to a low passage clone and cryopreserved in liquid nitrogen as a primary cell line stock for rederivation purposes. Fig 2B, 2C and 2E shows the expression pattern of known cell adherence, urothelial and immune associated markers on the MB49 cell line. Cells used in the experiments were thawed and cultured for not more than 7 passages. The cell culture and the orthotopic MB49 model instillation were performed as previously described [54], the animals were sacrificed and tumors were harvested after 5 days. The anesthesia used was 75 mg/kg ketamine and 1 mg/kg medetomidine administered intraperitoneally and reversal was induced by 1 mg/kg atipamezole administered subcutaneously, both at 0.1 ml per 10 g of body weight. During the duration of the experiment the animals were weighed and examined twice weekly. Urine and serum were collected 1 day prior to and 5 days post intravesical instillation of MB49 cells. The humane endpoint for the orthotopic MB49 model was the same as for carcinogen induced bladder cancer: 10 days of visible low hematuria or 5 days of complete hematuria or 15% body weight loss. In the experiments performed no MB49 tumor bearing animals reached the humane endpoint.

Single cell RNA sequencing (single-cell RNAseq) of whole bladders
Single cell transcriptomic analysis was performed to study the transcriptomic heterogeneity of tumor cells. Single cell analysis was performed for whole mouse bladders after a) 10 weeks of OH-BBN exposure and after b) progression to MIBC, as well as of c) healthy Hgf-Cdk4 R24C bladders, d) orthotopic MB49 tumor bearing (C57BL/6) bladders and e) healthy C57BL/6 bladders. On the 10-week endpoint Hgf-Cdk4 R24C tumors required pooling of samples due to low cell numbers per bladder (one pool of n = 4 females and one of n = 5 males). For the MIBC endpoint, one whole tumor bearing bladder per male and female were analyzed. Healthy Hgf-Cdk4 R24C bladders were pooled together (n = 4 males and n = 4 females). A pool of day 5 orthotopic MB49 tumors (pool of n = 5) and a pool of 5 healthy C57BL/6 bladders as control were also analyzed. Bladders were dissociated with the mouse tumor dissociation kit (Miltenyi cat. 130-096-730) according to the manufacturer's instructions to create single cell suspensions. Samples were submitted to downstream procedures if viability was 75% or higher. Samples were processed using Chromium Single Cell 3' GEM, Library & Gel Bead Kit v3 (10x Genomics cat. 1000092) according to the manufacturer's instructions and sequencing performed with a NovaSeq6000 sequencer. Appropriate QC and integrity analysis were performed. Read alignment, barcode matrices and cell clustering were performed with the Cell Ranger 3.0 software. Gene expression analysis and visualization of the data was performed in R studio [55] using the Seurat 3.0 package [56] for R. One large urinary epithelium cell cluster (14386 cells) was identified in the dataset on the basis of single cell clustering at low resolution and use of the cell identity prediction tool of the PanglaoDB database [57] using differential gene expression. Subsequently, the urinary epithelium cells were extracted from the dataset as a subset and clustered anew. Then, cell clusters found exclusively in tumor bearing samples were extracted as subsets and clustering was performed in order to investigate whether heterogeneous tumor cells could be identified. The expression of Mki67 was used as a cancer cell proliferation marker gene and its expression on the chosen cell clusters was assessed to determine highly proliferating cells. The expression of Ceacam1 was used as a marker of cancer cell stemness. MB49 cells were identified by expression of Mki67, Vimentin and Y-linked genes as these cells are originally derived from a male mouse and are implanted in female C57BL/6 bladders. The expression levels of a set of molecular subtype relevant genes were used to determine the molecular subtype profile tendency of each urinary epithelium cell cluster.

Microhematuria detection and Proximity Extension Assay (PEA)
Collection of all urine samples was performed at the same time (between 9-11AM) to facilitate alignment of the data. Microhematuria was assessed at timepoints day 24, day 48 and day 70 of carcinogen exposure when visible hematuria was absent by using a semi-quantitative test (Combur Test UX, Roche) and the remaining urine was cryopreserved along with urine from the MIBC endpoint at -80˚C for PEA analysis. Blood was collected at the 10-week endpoint and at the MIBC endpoint via tail vein incision, centrifuged at 10000 x g for 5 minutes, serum was collected and cryopreserved at -80˚C. Urine and serum samples were randomized on 96-well plates and blinded PEA analysis was performed (Mouse Exploratory panel, Olink Proteomics, Uppsala, Sweden). Healthy control samples from age-, genotype-, sex-matched animals were used as reference for NPX values (Normalized Protein eXpression, Olink's arbitrary unit). In the text, protein names are denoted in non-italic with a first capital letter and the rest lower cases for murine proteins, while human proteins are denoted in non-italic capital letters.

Anti-PD1 treatment in animals with autochthonous bladder tumors
To assess the therapeutic effect of anti-PD1 therapy on OH-BBN exposed Hgf-Cdk4 R24C male and female animals, one week after the 10-week exposure period, animals received six intraperitoneal administrations of 100 μg InVivoMAb anti-mouse PD-1 (CD279), Clone RMP1-14 (BioXCell cat. BE0146) in 50 μl PBS every 4.5 days for 3 weeks. Animals were monitored twice weekly or every 24h when signs of sickness appeared. Mice were sacrificed at humane endpoints (maximum 10 days of low hematuria or 5 days of complete hematuria or 15% body weight loss). Six out of 10 animals (3 out of 4 males, 3 out of 6 females) that received anti-PD1 therapy suffered sudden death prior reaching the humane endpoint. The survival of the treated animals was compared to historically recorded OH-BBN exposed, but untreated (6 female and 7 male) control Hgf-Cdk4 R24C mice which were also monitored daily upon appearance of palpable tumors and sacrificed at the humane endpoint.

Tumorigenicity of new tumor-derived cell lines and in vivo anti-PD-1 or CpG ODN 1668 treatment of subcutaneous tumors
The tumorigenicity of the 6 new Hgf-Cdk4 R24C urothelial cancer cell lines was assessed in 2 or 4 mice by subcutaneous inoculation of 5x10 6 cells to C57BL/6 female mice or 2x10 6 cells to Hgf-Cdk4 R24C female mice. Animal weight and tumor growth were assessed twice a week by caliper measurement for two months (experimental endpoint set to assess tumorigenicity of cell lines). To assess the therapeutic potential of anti-PD-1 or CpG ODN 1668, the only tumorigenic cell line, 74/7, was inoculated to the flanks of Hgf-Cdk4 R24C mice (8 females and 7 males per group, 3-5 months old). On day 0 the mice received subcutaneous inoculum of 2x10 6 cells. Animals were monitored, weighed and tumors were measured with a caliper twice a week and the ellipsoid formula "(4/3) x π x r1 x r2 x r3" was used to determine their volumes. From day 10 the mice were treated with 5 doses of 50 μg anti PD-1 (BioXCell cat. BE0146) peritumorally every third day or with 6 doses of 50 μg CpG ODN 1668 (tlrl-1668-5, Invivogen) peritumorally every third day. The control group received an equal number of PBS injections. For survival analysis, the animals were monitored until the tumor size reached the volume of 1 cm 3 which was the determined humane endpoint measure for subcutaneous tumors.

Statistical analysis and illustrations
Statistical significance of survival curves was determined by Log-rank (Mantel-Cox) testing. For direct comparisons of groups non-normal distribution of the data was assumed (due to small numbers of groups) and therefore non-parametric comparisons were performed (Mann-Whitney). Error bars represent standard error of the mean. Power analysis to determine the number of animals required for experimental assessment of anti-tumor therapeutic effect (increased survival/reduced tumor size) (alpha 0.05, power 80%), indicated that the sample size should be 5 mice per group. Hence, the study was designed to be sufficiently powered for the primary endpoint of determining if there is a difference among groups in anti-tumor response. In occasions where transgenic mice were used this number was sometimes lower due to the low ratio (below Mendelian ratio) of birth of Hgf-Cdk4 R24C mice. For statistical analysis of proteomic data in the form of NPX values, Mann-Whitney tests were used for comparisons between groups. For analysis of longitudinal proteomic alterations across time the Generalized Least Square model was applied. Volcano plots were generated with R studio [55]. KEGG pathway analysis was performed with the function 'Proteins with Values/Ranks' in the online database STRING [58]. Graphical workflow illustrations were created with BioRender. com. Graphs were generated with GraphPad Prism 7.0 and R studio 1.1.463 [55].

Establishment of OH-BBN induced urothelial tumorigenesis in Hgf-Cdk4 R24C mice
Administration of the pro-carcinogen OH-BBN (metabolized to active carcinogen in the liver and in the bladder) is used to selectively induce urothelial cancer in animal models [17].
Herein, we used tumor induction by 10 weeks of exposure to 0.05% OH-BBN in drinking water in Hgf-Cdk4 R24C or Cdk4 R24C mice. For some parameters we also compared the tumorigenesis in C57BL/6 female mice using either OH-BBN induced tumors or orthotopically implanted syngeneic MB49 tumors. In a pilot experiment, we compared tumor induction in Hgf-Cdk4 R24C weaned (4-6 weeks old; before full immune system development) and adult (10-14 weeks old, with fully developed immune system) female mice. Histopathology examination of bladders showed that carcinogen exposure of older mice was associated with increased immune cell infiltration, in contrast to mice exposed at weaning age that presented with low or no inflammation (S1A Fig). Since inflammation is an important disease component in humans, we used adult mice for further studies. Bladders were collected for histology analyses at two timepoints: 1) 10 weeks of carcinogen exposure and 2) at development of symptoms of advanced urinary bladder tumor (humane endpoint due to tumor burden) ( Fig 1A). Pathology assessment revealed that non-carcinogen exposed female Hgf-Cdk4 R24C mouse bladder histology was normal (Fig 1B), whereas carcinogen treated Hgf-Cdk4 R24C females displayed a urothelial cancer of Tis-T1 stage with squamous differentiation (Fig 1C) which was more advanced than for transgenic Hgf-Cdk4 R24C males ( Fig 1D) or Cdk4 R24C littermates ( Fig 1E).
Next, progression and survival of the tumor-bearing Hgf-Cdk4 R24C and Cdk4 R24C mice were investigated under active surveillance until reaching the tumor burden related humane endpoint. Histopathology analysis upon active surveillance revealed that in transgenic mice, tumors progressed to MIBC with squamous differentiation (Fig 1F and 1G). For wildtype mice, MIBC rarely develops after 10 weeks of OH-BBN exposure [17], thus C57BL/6 mice were not investigated towards progression to MIBC. However, for the first endpoint of 10 weeks exposure as a comparison we found that female C57BL/6 mice developed urothelial dysplasia with or without focal carcinoma in situ (CIS) (Fig 1I), which was not as advanced as in transgenic female or male mice (Fig 1C and 1D). In contrast, orthotopic MB49 tumors presented with more consistent and homogenous histology than OH-BBN induced tumors, with low degree of inflammation ( Fig 1J). The transgenic mice reached a humane endpoint within 5-11 weeks of active surveillance after ending carcinogen exposure, while the survival difference among sexes did not reach statistical significance (Median survival: Female: 56 days; Male: 43 days). The survival of Cdk4 R24C mice upon progression to MIBC was significantly longer (p = 0.0377) compared to their Hgf-Cdk4 R24C littermates ( Fig 1K) (Median survival: Hgf-Cdk4 R24C : 56 days; Cdk4 R24C : 83 days). Based on this we pursued our studies with the Hgf-Cdk4 R24C mice which seem to progress faster than mice carrying the Cdk4 R24C mutation alone without overexpression of Hgf. Interestingly, after only 10 weeks of carcinogen exposure, female Hgf-Cdk4 R24C bladders presented with a noted enlarged volume than males, whereas at the MIBC timepoint male bladders tended to be larger than for females without reaching statistical significance (Fig 1L). The summary of detailed histopathology analysis for Hgf-Cdk4 R24C mice is presented in Fig 1M. In half of the female Hgf-Cdk4 R24C bladders with MIBC, lymphovascular invasion by tumor cells was observed, but none of the mice developed distant metastasis in lungs or kidney tumors. Despite the advanced tumor stage noted in a few mice (T4), no hydronephrosis was observed.

Development of Hgf-Cdk4 R24C urothelial cancer cell lines
To develop a novel in vitro model of bladder cancer, we generated cell lines originating from Hgf-Cdk4 R24C urothelial tumors using whole bladder cell pools (Fig 2A). Thereafter, to exclude the potential presence of cells of non-malignant origin in the cell cultures, we analyzed the expression of epithelial-mesenchymal transition (Icam1 and E-cadherin) and fibroblast  (Pdgfra and Cd90.2) markers by flow cytometry (Fig 2B). To further determine the cell origin and the molecular profile of the new cell lines, we examined the expression of genes that characterize muscle (Des, Myh11, Tagln), stroma (Vim, Pdgfra, Col1a2), epithelial mesenchymal transition (EMT) (Cdh1, Cdh2), urothelial (Krt7), basal (Krt5, Krt14, Trp63) and luminal (Krt18, Krt20) cells and compared the novel cell lines with established cell lines and healthy Hgf-Cdk4 R24C urothelium (Fig 2C). All of the novel cell lines displayed a lack of or very low expression of genes related to muscle and stroma when compared to healthy bladder, the urothelial cell line MB49 or the non-urothelial cell lines TRAMP-C2, NIH/3T3 and CT26. All cell lines showed that they are undergoing or have completed EMT through expression of lower levels of Cdh1 and higher levels of Cdh2 compared to healthy bladder tissue, with the exception of 73/2 cell line. Expression of two general urothelial markers (Upk2, Krt7) was stronger in healthy bladder followed by the new cell lines 74/7, 55ad/1, 51/4 and 73/2, while a tendency towards basal molecular profile was displayed by these cell lines as assessed by expression of Krt5 or Krt14. Trp63 was strongly expressed by non-urothelial cell lines followed by the new 43/1 and 55fl/1 cell line. All new cell lines and MB49 cells expressed moderate or lower levels of the luminal gene markers Krt18 and Krt20 in comparison to healthy bladder. Hierarchical clustering based on the expression of these genes revealed that the novel cell lines and MB49 cells were more similar to each other than to healthy bladder (Fig 2C).  (Fig 2D). To explain the lack of tumorigenicity in immunocompetent mice of all novel cell lines except for 74/7, the expression of immune recognition receptors was examined before and after IFNγ stimulation in vitro using flow cytometry ( Fig 2E). As previously shown [59], MB49 cells that were used for comparison in this experiment were MHC class I positive, and upon IFNγ stimulation upregulated MHC class I, II and PD-L1 (Fig 2E). Without IFNγ stimulation the five non-tumorigenic new cell lines were MHC class I and II deficient, PD-L1 low/neg , CD80 deficient (except for 73/2) and FasL low . However, they displayed increased MHC class I and/or PD-L1 expression after exposure to IFNγ for 24h in vitro ( Fig 2E). Conversely, the 74/7 tumorigenic cell line showed no response to in vitro IFNγ stimulation (Fig 2E) providing evidence that the other cell lines did not form tumors in vivo due to immune rejection.

Single cell transcriptomic analysis identified distinct transformed tumor cell clusters
To our knowledge there is very limited number of studies that have characterized healthy epithelial cells of the human and mouse bladders [60], as well as tumor, stroma and immune cells from human bladder tumors [61] using single cell transcriptomics. Molecular subtyping for human bladder cancer is performed by bulk DNA or RNA sequencing of tumor biopsies. To better visualize the intratumoral transcriptomic heterogeneity in bladder cancer, we performed single-cell RNAseq analysis on murine OH-BBN induced Hgf-Cdk4 R24C urothelial tumors as well as on orthotopic MB49 tumors (Fig 3A). The urinary epithelium cells were extracted from the dataset and analyzed as described in Materials and Methods. Initially eight cell clusters of urinary epithelium origin with transcriptomic heterogeneity were derived, cl00-cl07 number of cells from each cluster in each bladder analyzed for Hgf-Cdk4 R24C or C57BL/6 mice. Out of all cell clusters in Hgf-Cdk4 R24C and MB49 tumor bearing bladders (B, C, D), cl01 which is present only in OH-BBN induced Hgf-Cdk4 R24C bladders and cl04 which is present only in MB49 tumor bearing bladders, were further subsetted/ subclustered and the expression of (D) tumor or cancer stem cell marker genes, as well as molecular subtype relevant genes is visualized for each model. https://doi.org/10.1371/journal.pone.0253178.g003 (Fig 3B and 3C, S1 File). Cluster cl01 was enriched only in OH-BBN induced Hgf-Cdk4 R24C bladders and cl04 only in MB49 tumor-bearing bladders (Fig 3D). Cells belonging to cl01 and cl04 were subclustered again and then characterized exclusively in terms of gene expression of proliferation markers, cancer cell stemness and molecular subtype profile markers. Upon reclustering, four clusters emerged from cl01 ( Fig 3E): 1) a highly proliferating cluster (based on Mki67 expression) with basal/squamous and cell cycle deregulated profile gene expression, while some cells of this cluster in NMIBC also displayed Fgfr3 expression, 2) a cancer stem cell cluster (based on Ceacam1 expression) with mixed expression of subtype-related genes and higher Pparg expression than the other clusters, 3) a cancer stem cell cluster enriched in female MIBC that displayed a heterogeneous profile of subtype-related genes expression, and 4) a cluster with a heterogeneous Ceacam1 expression that was present mainly in male MIBC including cells with high and low Egfr expression, but overall low expression of subtype-related genes. The MB49 cell cluster (Fig 3E, right) displayed a strong cell cycle dysregulation gene expression profile, similar but higher in the expression scale compared to the highly proliferating Hgf-Cdk4 R24C tumor cell cluster, and with visibly higher expression of Tp53 compared to all other Hgf-Cdk4 R24C clusters. These findings suggest that urothelial cells that are present only in carcinogen induced tumors, but not in healthy bladders, may or may not present high proliferation or cancer stemness, while still contributing to the tumor burden and possibly to different subtype profiles.

Distinct protein profiles of serum and urine among bladder cancer models
Histopathology analysis depicted strain specific differences in OH-BBN induced bladder cancer development in wild type transgenic mice, while single-cell RNAseq showed differences in molecular profile between the MB49 and the OH-BBN induced Hgf-Cdk4 R24C model. To further investigate differences of the systemic protein profile and local tumor microenvironment among the different models, we performed a targeted proteomic profiling on liquid biopsies.
Serum. When comparing serum signatures from the two OH-BBN induced strains (10-week endpoint) and the orthotopic MB49 model (day 5) (Fig 4) a mutual decrease of Ghrl was identified among the three models. Decreased levels of Plin1 were seen in wild-type OH-BBN exposed animals and MB49 day 5 orthotopic tumor bearing animals when compared to healthy controls of respective model (Fig 4A-4C). Moreover, the orthotopic MB49 model displayed decreased Ccl5, Cxcl9 and Cxcl1 levels, suggesting a less immunogenic systemic impact of the tumors. In this model, we further observed increased levels of the anti-proliferative cell cycle regulator Map2k6 and the apoptosis mediator Casp3 (Fig 4A). Additionally, an increased concentration of the apoptosis marker Parp1 was observed in wild type mice ( Fig  4B) and increase of the tumor suppressor transcription factor Eda2r was seen in transgenic animals when compared to healthy controls of respective model (Fig 4C). Interestingly, Il6 levels were decreased in C57BL/6, but increased in the transgenic animals. This inverse correlation could suggest a polarized early systemic immune response during the time of tumor establishment in the different strains. In line with this observation, the pro-inflammatory , carcinogen induced wild type female C57BL/6 (10-week endpoint) and carcinogen induced male and female Hgf-Cdk4 R24C (10-week and MIBC endpoint). Significantly altered serum proteins in (A) MB49 orthotopic model (n = 6), (B) Carcinogen induced female C57BL/6 mice (n = 6), (C) 10-week endpoint Hgf-Cdk4 R24C mice (n = 13), (D) MIBC endpoint Hgf-Cdk4 R24C mice (n = 8) compared to respective healthy controls (n = 6-10). Levels of significantly changed protein levels in the urine samples of (E) MB49 orthotopic model (n = 4), (F) Carcinogen induced female C57BL/6 mice (n = 6), (G) 10-week endpoint Hgf-Cdk4 R24C mice (n = 18), (H) MIBC endpoint Hgf-Cdk4 R24C mice (n = 8) compared to respective healthy controls (n = 4-13). Fold change for each protein was calculated as 2 (mean NPX cancer-mean NPX healthy control) . Colored circles indicate significance level of p<0.05; red for increased proteins and blue for decreased proteins for cancer model compared to healthy controls, grey for proteins below the significance level. cytokines Csf2, Il1β and Il5 were downregulated in the serum of C57BL/6 mice. At the MIBC endpoint serum of transgenic mice displayed increased Epo, nitric oxide and tumor growth regulator Ddah1, and Fslt3.
Urine. No significantly altered protein levels were detected in the inoculated MB49 mice at day 5 in urine when compared to urine before instillation (Fig 4E). However, the urine profile of OH-BBN induced transgenic and wild type mice were distinct at the 10-week endpoint (Fig 4F and 4G). Contrary to the serum profile (Fig 4B), urine from OH-BBN induced C57BL/ 6 mice displayed increased levels of the pro-inflammatory cytokines Ccl3 and Il1β, as well as an upregulation of tumor promoting Ppp1r2, Pak4, tumor growth marker Nadk, and Dctn2 ( Fig 4F). Also, increased level of Ahr was detected in wild type OH-BBN induced mice. The transgenic mice displayed increased pro-inflammatory response markers linked to the fatty acid metabolism/prostaglandin pathway at 10-week endpoint, such as S100a4 and Il1β ( Fig  4G). In MIBC the fatty acid/prostaglandin pathways remained as an upregulated signature along with other pro-tumorigenic factors (Fig 4H), such as Pdgfb and Tgfb1 as well as signals of an active DNA repair machinery, such as Parp1 and Apbb1ip.

Proteomic analysis of urine and serum from early carcinogenesis to progression in male and female mice
An obstacle in studying human clinical material is performing repeated sample analysis due to logistic, ethical and feasibility constraints. Studies capturing early tumorigenesis are often too costly to perform as a large number of study subjects are required to capture tumor development. Bladder cancer patients at high risk of progression (based on histopathology and clinical parameters) undergo radical surgical interventions to prevent metastatic disease, limiting clinical sampling post that stage. However, an animal model of carcinogen induction can capture the timeline of early carcinogenesis as well as of the NMIBC and MIBC stages. A number of male and female OH-BBN induced Hgf-Cdk4 R24C animals displayed microhematuria already at day 24, but also throughout day 48 and the 10-week endpoint (Fig 5). The detection of microhematuria was used as a marker of early tumor lesion formation and urine collected at those time-points was used for proteomic analysis. Microhematuria or hematuria was not observed in day 0 samples and neither in health matched controls at any time-point. After the 10-week endpoint, animals displayed irregular hematuria until the progression to MIBC which resulted in weight loss, gross hematuria and palpable bladder tumors. Fig 5A-5H outlines the proteomic profile during the early tumor development phase in male and female mice. Early protein alterations were present already at day 24 (3.5 weeks) of carcinogen exposure. The number of elevated proteins in urine at day 24 was larger in females (17 proteins) (Fig 4A) compared to males (6 proteins) (Fig 5C). Functional enrichment of the upregulated proteins in males and females using STRING [58] revealed pathways and activities associated with immune and stress responses, inflammation and cytokine regulation, leukocyte homing, apoptotic processes and cancer development (S2-S5 Files). Interestingly, most upregulated proteins at day 24 re-surfaced at the MIBC time-point whereas they were not detected at day 48 and 10-week endpoint (Fig 5B, 5D, 5E and 5G). An intersection analysis

PLOS ONE
Single-cell RNAseq and longitudinal proteomic analysis of a novel semi-spontaneous urothelial cancer model

Effects of anti-PD1 and CpG ODN immunotherapy in Hgf-Cdk4 R24C autochthonous and subcutaneous tumors
To investigate whether or not anti-PD1 therapy can hinder progression to MIBC in Hgf-Cdk4 R24C mice with OH-BBN induced NMIBC we assessed it in this model (Fig 6A). No significant difference in survival was noted in response to anti-PD1 as compared to control survival measurements in the same strain (Median survival: Female untreated: 56 days; Male untreated: 43 days; Female anti-PD1: 77 days; Male anti-PD1: 54,5 days) (Fig 6C), although a trend of an increased median survival in both treated and untreated females was noted compared to male counterparts. Unexpectedly, many anti-PD1 treated animals died spontaneously before reaching the humane endpoint (included in the survival analysis shown), in contrast to the untreated group. Therefore, tumor measurement or other examinations of these animals were not performed. The repeated occurrence of spontaneous deaths could indicate other complications such as tumor spread, (uro)sepsis, anti-PD1 side effects, etc. Although there are no reports of direct immunotoxicity of anti-PD1 therapy in the Hgf-Cdk4 R24C mouse strain in previous studies [62], it is possible that anti-PD1 therapy affects the local tumor-microenvironment which may impact tumor growth and invasiveness. Therefore, to examine whether the unresponsiveness of OH-BBN induced Hgf-Cdk4 R24C tumors to anti-PD1 therapy was due to a progression-prone microenvironment and stroma, subcutaneously inoculated 74/7 tumors (growing on Hgf-Cdk4 R24C mice) were treated with anti-PD1 or CpG 1668 (Fig 6B). While anti-PD1 therapy did not impact tumor growth nor the well-being of the animals, subcutaneous tumor growth in Hgf-Cdk4 R24C mice was significantly reduced in CpG 1688 treated animals at day 22 (Fig 6D) (which was the latest time-point when all animals were still alive). When survival (Fig 6E) (Fig 6F-6K) were analyzed separately for each sex, a clear response of increased survival benefit (p = 0.05) ( Fig 6E) and decreased tumor growth in CpG 1688 treated females (p = 0.026) (Fig 6K), but not males, was observed.

Discussion
Recreating urothelial tumors in mouse models has been challenging. Gene expression signatures are not often consistent between transgenic models and human tumors, engraftment models do not develop relevant microenvironment, and carcinogen exposed models commonly do not reach the MIBC stage [18] or are not analyzed at that stage. Herein, we have induced urothelial cancer in transgenic Hgf-Cdk4 R24C mice whose germline genetic alterations lead to susceptibility towards development of tumors of epithelial origin spontaneously or upon carcinogen or UV radiation exposure [49,50]. Patients with bladder cancer present urothelial cell cycle dysregulation including aberrant activity of CDK4 [43,44], while overexpression of HGF is known to be involved in early stages of carcinogenesis [63]. Excessive expression of HGF by stromal cells promotes the survival and invasion of tumor cells [64], but also the polarization of immune microenvironment cells towards immunosuppressive phenotypes [65,66]. Thus, herein we present a model that harbors both pro-oncogenic factors leading to a dysregulated cell-cycle as well as invasiveness, along with a broad mutational burden by the induction of the tumorigenesis by OH-BBN.
The Hgf-Cdk4 R24C mice developed NMIBC tumors that progressed to MIBC over time suggesting that this model reflects a progressive rather than relapsing NMIBC. Additionally, the histopathology of Hgf-Cdk4 R24C bladder tumors was of urothelial carcinoma with squamous differentiation, which in humans is associated with progression, chemotherapy resistance and possibly disease progression after anti-PD1 therapy [67][68][69][70][71]. Also, the faster progression and shorter survival of Hgf-Cdk4 R24C mice compared to Cdk4 R24C with OH-BBN induced urothelial tumors indicates that the co-occurrence of both Hgf overexpression and overactivation of Cdk4 (R24C) led to a more progressive tumor phenotype.
Moreover, we observed a lack of metastases or tumor development in other sites than the bladder which were previously recorded in studies while using Hgf-Cdk4 R24C mice [49][50][51]. It is possible that this is a cancer specific property, as multiple metastases arise in very advanced stages in human patients usually after surgical removal of the bladder [72] or that circulating tumor cells did not evolve into tumor metastases due to the necessary humane endpoint reached when the animals reached the muscle invasive bladder tumor stage. We could not examine the metastatic potential, in relation to patients who present metastatic tumor growth after surgical removal of the bladder, due to lack of surgical measures as well as due to ethical constraints. Moreover, we attribute the lack of spontaneous tumors to the age of the mice used in our study (up to 7.5 months old), while it has been reported that mice of the Hgf-Cdk4 R24C mice develop tumors spontaneously after the average age of 8 months [52].
In this study, we generated and characterized six cell lines from single cell colonies of Hgf-Cdk4 R24C urothelial tumors, a strategy that might impact intratumoral heterogeneity and therefore limited tumor growth upon engraftment in vivo, as previously shown [73]. Indeed, five of six cell lines did not grow in vivo. Those were IFNγ responsive in vitro, a quality that is required of tumor cells for recognition and elimination by CD8 T cells. This could explain why 74/7, which did not respond to IFNγ by expression of immune recognition receptors, was the only cell line to form tumors upon engraftment.
The presence of genetically heterogeneous urothelial tumor foci was recently highlighted [74], while it became apparent that such heterogeneous tumor compartments are usually under-represented in patient tissue microarray biopsies [75,76]. Molecular subtypes are investigated by bulk sequencing on such biopsies, which probably limits the examination of intratumoral subtype heterogeneity. Recently, intratumoral heterogeneity in human urothelial bladder cancer was identified as clonal variety by DNA seq, where subclones of invasion were found in samples of progression, adding to the genetic mosaic theory in cancer [77]. In our study, the OH-BBN induced tumors displayed multifocal and heterogeneous histopathology with concomitant CIS similar to human bladder cancer [74][75][76]. To model a setup where this would be analyzed in a deeper level, we used single-cell RNAseq to decipher the transcriptomic heterogeneity of Hgf-Cdk4 R24C bladder tumor cells. Data analysis revealed distinct tumor and cancer stem cell-like cell clusters with different molecular subtype expression patterns in carcinogen induced tumors, contrary to homogeneous MB49 tumors. Our results altogether indicate that OH-BBN induced tumors in Hgf-Cdk4 R24C mice recapitulate human urothelial cancer induction, progression and heterogeneity. Moreover, we suggest that single cell analyses of clinical bladder cancer samples may shed light on whether coexisting heterogeneous tumor cell populations have impact on prognosis and/or therapy outcomes.
Despite the technical limitation of the PEA due to the lack of a reference protein to enable exact protein level comparison between mouse strains, proteomic analysis of liquid biopsies revealed that the novel mouse model has distinct profile from commonly used models, therefore suggesting an alternative shaping of tumor macro-and microenvironment in this strain. Repeated urine sampling of the transgenic strain over-time enabled us to capture the proteomic changes during early carcinogenesis, tumor establishment and at the MIBC stage. Those changes suggest that the tumor microenvironment was somewhat different between sexes, similar to humans [11] and that during carcinogenesis pre-tumoral urine proteomic alterations can be detected by highly sensitive PEA analysis.
Up to 75% of locally advanced MIBC or metastatic bladder cancer patients do not respond to anti-PD1 therapy [78,79]. This could be probably due to established immunosuppression in late stages that is beyond inversion or due to specific subtype related responses to immunotherapy as recently highlighted in a suggested classification system [80]. Despite the lack of a stratification method with high molecular precision in the clinic, anti-PD1 therapy is being tested in clinical trials to prevent recurrence and progression of high-risk NMIBC after tumor resection, as well as in conjunction with BCG therapy in a pre-surgical setting [81]. Recently, the FDA approved the use of anti-PD1 for BCG-unresponsive, high-risk NMIBC, after the therapy achieved responses in almost a third of patients [5]. In Hgf-Cdk4 R24C mice with OH-BBN induced tumors monotherapy with an anti-PD1 blocking antibody did not prevent the progression of non-invasive tumors. It is however plausible that tumor resection and combination therapies might be the key to achieve complete tumor regression for progressive disease or that the model could be a useful tool to further study mechanisms of immunotherapy resistance, since a recent study showed that anti-PD1 therapy of OH-BBN induced MIBC tumors in wild type C57BL/6 mice had a therapeutic effect [82]. Similar to OH-BBN induced bladder tumors, subcutaneous 74/7 tumors in Hgf-Cdk4 R24C mice did not respond to anti-PD1 monotherapy in vivo at a dose of 2-2.5 mg/kg. It is plausible that a higher dose anti-PD1 could achieve anti-tumor responses as a peritumoral dose of around 1.2-1.5 mg/kg did not impact MB49 tumor growth while an IV administration of 4-5 mg/kg was curative [83].
Interestingly, peritumoral CpG 1688 monotherapy significantly reduced tumor growth in female animals, while male mice did not respond to this therapy. The choice of this therapy was based on previous noted therapeutic impact on MB49 subcutaneous tumors inoculated in the C57BL/6 strain [15,83]. The systemic administration of CpG ODNs in oncology based clinical trials limited their potential, as CpG ODN infusion induced systemic cytokine release and narrowed the therapeutic window [84,85]. Therefore the administration route is key for success, as well as the presence of tumor antigens along with the adjuvant [86]. The use of these adjuvants as local immune-modulators in the neo-adjuvant setting in clinical trials is supported by publications [87]. Herein, we used peritumoral administration of CpG 1668 in vivo and observed a successful reduction of growth of 74/7 tumor graft in female hosts. Recently, a study elegantly showed that innate immune responses are more potent in females compared to males [12]. This suggests that in our study TLR stimulation by CpG 1688 possibly triggered a more potent anti-tumor response in female mice with subcutaneous 74/7 tumors. This is of particular interest as to date checkpoint blockade immunotherapy seems to yield a higher survival benefit for male rather than female patients, indicating that different therapeutic approaches should be considered among sexes [11,88,89]. This would require not only equal stratification of men and women to clinical trials, but also the use of pre-clinical models that are appropriate for drug target discovery and evaluation in terms of sex biology.
In conclusion the novel Hgf-Cdk4 R24C autochthonous urothelial cancer model and cell lines should be utilized to further understand the immuno-biology implications of tumor subtype and biological sex on therapeutic immune modulation in bladder cancer.