Integrative Single Cell Atlas Revealed Intratumoral Heterogeneity Generation from an Adaptive Epigenetic Cell State in Human Bladder Urothelial Carcinoma

Abstract Intratumor heterogeneity (ITH) of bladder cancer (BLCA) contributes to therapy resistance and immune evasion affecting clinical prognosis. The molecular and cellular mechanisms contributing to BLCA ITH generation remain elusive. It is found that a TM4SF1‐positive cancer subpopulation (TPCS) can generate ITH in BLCA, evidenced by integrative single cell atlas analysis. Extensive profiling of the epigenome and transcriptome of all stages of BLCA revealed their evolutionary trajectories. Distinct ancestor cells gave rise to low‐grade noninvasive and high‐grade invasive BLCA. Epigenome reprograming led to transcriptional heterogeneity in BLCA. During early oncogenesis, epithelial‐to‐mesenchymal transition generated TPCS. TPCS has stem‐cell‐like properties and exhibited transcriptional plasticity, priming the development of transcriptionally heterogeneous descendent cell lineages. Moreover, TPCS prevalence in tumor is associated with advanced stage cancer and poor prognosis. The results of this study suggested that bladder cancer interacts with its environment by acquiring a stem cell‐like epigenomic landscape, which might generate ITH without additional genetic diversification.

The distinctive features of Bs2/3 cells relative to BsP/Bs1 cells were the upregulated levels of PCSK5, SGK1, HES1, PRKCZ, HOXA3, and HOXA5 and the lack of HBEGF or GPX4 expression (Figure S9C).BsP cells exhibited upregulated expression of proliferation-related genes, such as MKI67 but did not exhibit HES4 expression.Additionally, the expression of Tm4sf1 was detected in BsP cells.Meanwhile, the expression levels of GATA3, uroplakin-encoding genes (UPK2/UPK3A/UPK3B), and PIGR were upregulated in intermediate/umbrella cells (Figure S9C).
CXCL6 was a marker of the umbrella cells.Diffusion map inference of the developmental trajectory revealed that the basal progenitor cell evolves into the following three branches: Bs1, Bs2/Bs3, and Im1/2/3/Um (Figure S9B).This is consistent with the findings of previous studies [1].Consistent with the results of trajectory analysis, analysis with the CellCycleScoring algorithm revealed that the proliferation of BsP cells was significantly higher than that of other urothelium cells, suggesting that BsP cells are the proliferating progenitors for uroepithelium (Figure S9D).Notch signaling is reported to be essential for urothelium differentiation [2].In the normal urothelium, the major source of Notch signaling may be from the stromal compartment.The expression levels of the Notch ligands DLL1 and JAG1 in fibroblasts were markedly higher than those in epithelial cells (Figure S9C and Figure S18).As Notch signaling is contact dependent, this expression pattern indicates that Notch signaling is spatially confined to the basal cell layer, which is close to the stromal cells.The specific expression of the Notch downstream targets HES1 and HES4 and the upregulation of the Notch receptors NOTCH1 and NOTCH3 in basal cells suggest that Notch activity is upregulated in basal cells (Figure S9C).Furthermore, cell-cell signaling analysis revealed mutual Notch signaling between basal and intermediate/umbrella cells (Figure S9E).These results suggest that the normal urothelium is derived from basal cells under the control of Notch.

DNA hypomethylation accompanies TPCS generation
Since epigenotypes are defined as cells with similar ChrAcc profiles on tumor-associated DMR, we first analyzed how these profiles change between epigenotypes.Two distinct patterns of ChrAcc shift were identified related to tumor-associated DMRs (Figure S10A).A subset of DMRs, mostly hypermethylated in BLCA (hyperDMR), was exclusively "open" and accessible in stromal cells but remained closed in all epithelial cells irrespective of their malignancy status (Figure S10A,   B).Meanwhile, a second class of DMRs, mostly hypomethylated in BLCA (hypoDMR), exhibited an increased accessibility in cancer cells compared with that in both healthy epithelial and stromal cells (Figure S10A, C).The TPCS-related genes KRT6A/KRT6B and KRT17 were associated with DNA hypomethylation and exhibited increased ChrAcc in cancer cells (Figure S10A).
As DNAm is typically accompanied by an inaccessible chromatin state, the findings of this study suggest that DNAm levels at many hyperDMRs remain unchanged during bladder cancer oncogenesis.Consequently, the apparent DNA hypermethylation of most of these regions in tumor tissue results from cancer clonal expansion.Conversely, active DNA demethylation occurs at hypoDMRs during BLCA oncogenesis.Such a model implies that epigenotype evolution is mainly accompanied by changes in the ChrAcc and DNAm states on hypoDMRs.On the other hand, ChrAcc and DNAm states on hyperDMRs can serve as a tool for distinguishing cancer cells originating from different cell types, as they are more likely to remain constant during tumorigenesis.The ChrAcc of hyperDMRs changed little between cancer and normal epithelial cells, whereas hypoDMRs were more accessible in cancer cells (Figure S10D).This result is consistent with that of previous studies, which reported predominant DNA hypomethylation in BLCA [3].

DNA hypermethylation marks the cell-of-origin in BLCA
We performed unsupervised hierarchical clustering analysis on the DNAm profiles of hyperDMR in BLCA.Interestingly, the papillary, low-grade, non-invasive Ta tumors and the high-grade invasive T1-T4 tumors, which include all muscle-invasive bladder cancer (MIBC) samples, could be classified by DNAm levels on a core subset of hyperDMR (Figure S10E).The correlation of DNAm levels at core hyperDMRs with their adjacent gene expression levels suggests that different BLCAs may arise from distinct cell-of-origin.In particular, hyperDMRs specific to NMIBC were detected adjacent to genes whose expression levels in basal cells were higher than those in intermediate cells (GPX4/PRKCZ/CD44) (Figure S9C and Figure S10E).Conversely, hyperDMRs specific to high-grade tumors were associated with genes upregulated in intermediate cells, such as YWHAZ and GATA3 (Figure S9C and Figure S10E).This study prospectively sampled BLCA tumors before the assay.Hence, the stage of the tumor may be directly associated with its aggressiveness.Tumors with a high malignant potential may be highly likely to be sampled at an advanced clinical stage.A subgroup of T3/T4 MIBC cases exhibited hypermethylation in all core hyperDMRs and exhibited a correlation with the BsP origin based on adjacent gene expression (Figure S9C and Figure S10E).This suggests that the most aggressive BLCA tumors directly originate from BsP or that aggressive tumors transform to adopt a BsP-like epigenome (Figure S10E).Together, our results suggest that human BLCA has a heterogeneous origin, with the papillary low-grade tumors arising from superficial urothelium cells and high-grade tumors developing from basal cells.Prostate), and lymph nodes (LN).For each tissue sample, we either conducted a single-cell assay (scRNA, scATAC, scTCR) or a combination of these assays.The specific count of each type of bulk assay performed is detailed in the lower panel.

Figure S1 . 1 .
Figure S1.Donors and tissues in single-cell study, related to Figure 1.The donor ID and tissue type are displayed on the X and Y axes, respectively, with each tile colored according to the single-cell experiment type, as indicated in the lower panel.The number of tissues for each type is listed on the right, and the number of donors is listed at the top.The single-cell study encompassed a total of 37 donors, which included 26 donors in this study and 11 donors from PRJNA662018 dataset.This analysis featured 67 tissue samples across various types.Specifically, from PRJNA662018 dataset, we analyzed 11 tissues, 7 cancer (Ca.) and 4 cancer adjacent (Ca.adj.) tissues, while the remaining samples were directly involved into this study.The tissue types analyzed in single-cell study comprised cancer tissues (Ca.), cancer adjacent tissues (Ca.adj.), normal bladder tissues from DCD donors (Norm.Bladder), normal prostate tissues

Figure S2 .
Figure S2.Donors and tissues in bulk study, related to Figure 1.The donor ID and tissue type are displayed on the X and Y axes, respectively, with each tile colored according to the bulk tissue experiment type, as indicated in the lower panel.The number of tissues for each type is listed on the right, and the number of donors is listed at the top.The bulk study encompassed a total of 79 tissues from 58 donors in this study.Bulk analysis examined a variety of tissue types: cancer tissues (Ca.), cancer adjacent tissues (Ca.adj.), normal bladder tissues from DCD donors (Norm.Bladder), normal prostate tissues (Norm.Prostate), peripheral blood mononuclear cells (PBMC) and metastatic prostate tissues (Met.Prostate).Each tissue sample was underwent either one bulk assay (Mutation analysis, bulk DNA mehtylation, bulk ATAC) or a combination of these assays.The specific count of each type of bulk assay performed is detailed in the lower panel.

Figure S3 .
Figure S3.scRNA cell features, related to Figure 1.(A) Tissue-of-origin of scRNA cells for each donor.Single or multiple samples could be taken from a single donor.Donor/patient ID on Y axis was designed by combing the tumor stage of a donor with a sequential identification number.(B) Cell cycle score (G2M.Score + S.Score) of scRNA cells visualized on UMAP by color.(C) Single cell gene expression correlation with survival (Scissor) analysis of scRNA cell types.Red: unfavorable clinical outcome.Orange: favorable clinical outcome.Gray: unclear correlation.Epithelial cells are classified as cancer (EpiM) or normal (EpiN).

Figure S5 . 3 .(
Figure S5.Donor origin of single cells in the multimodal integration dataset, related to Figure 3. (Left) UMAP projection of scATAC cells or (Right) scRNA cells in the multimodal integration dataset.Cells are colored by their donor origin.scATAC and scRNA cells from similar donor are spatially close to each other, indicating accuracy of the multimodal integration.Donor/patient ID was designed by combing the tumor stage of a donor with a sequential identification number.

Figure S6 . 3 .
Figure S6.Developmental trajectory between single cell clusters in the multimodal integration dataset, related to Figure 3. Single cells, irrespective of their origin (scATAC or scRNA), are clustered by multimodal integration (GLUE).Cell types of these single cell clusters are annotated by the majority of scRNA/scATAC cells.Developmental trajectory between these single cell clusters is inferred by PAGA.

Figure S8 .
Figure S8.Epigenotype bridges transcriptional transition in BLCA, related to Figure 3. (A) UMAP of scATAC epigenotypes (left) or scRNA cell types (right) in the multimodal integration dataset.Whilst scRNA cell types are clearly segregated from each other, scATAC epigenotypes are mixed, suggesting phenotypic variability within epigenotypes.(B) Transition graph between scRNA and scATAC clusters.Few direct RNA-RNA transitions (black dashed line) between scRNA clusters (triangle) are found.Most transitions are bridged by scATAC epigenotypes (circles).(C) Degree of connectivity for each scATAC epigenotype or scRNA cell type.(D) Degree of connectivity of different single cell clusters (TPCS, basal cancer, or luminal cancer).P-value: wilcox test.

Figure S9 .
Figure S9.Basal cells generate the normal urothelium, related to Supplementary Notes.

Figure S10 .
Figure S10.DNA hypomethylation and epigenome reprograming drives BLCA evolution, related to Supplementary Notes.(A) Single cell chromatin accessibility of non-hematopoietic cells on MIBC-specific differentially methylated regions (DMR).DMR are classified as hypomethylated or hypermethylated in bladder cancer.Single cells are classified by epigenotypes into stromal (2/3/1/4), normal epithelial (14), or cancer (5-13 and 15).MIBC hypermethylated DMR (hyperDMR) chromatin are only open in stromal cells and their chromatin accessibility does not change between normal epithelial cells and cancer.MIBC hypomethylated DMR (hypoDMR) chromatin are closed in not only stromal but also normal epithelial cells, and their chromatin accessibility varied between different cancer cells.

Figure S12 .
Figure S12.Joint scATAC/scRNA analysis of the developmental trajectory from BsP towards TPCS, related to Figure 6.(A) UMAP projection of scATAC normal (gray) or TPCS (purple) cells.(B) Single cells colored by clnical stage of their donor.T0: DCD healthy donor.(C) Single cells colored by their epigenotypes.C14: normal urothelium epigenotype; C5/6/7: TPCS epigenotype.(D-F, left) Developmental trajectory on UMAP.Single cells belonging to the trajectory are colored by their respective developmental pseudotime (early-to-late: purple to gold).(D-F, Middle) integrated scRNA transcription factor gene expression profile across the developmental trajectory.(D-F, right) integrated scATAC transcription factor activity profile, of the similar transcription factors as in middle panel, across the developmental trajectory.On the middle and right panels, early developmental stage-specific transcription factors are labelled by blue arrow, late developmental stage-specific transcription factors by red arrow, and lineage-specific transcription factors by green arrow.(D) Developmental trajectory 1, towards a GATA3 high , luminal-lineage priming TPCS cluster.(E) Developmental trajectory 2, towards a ETV4 high , basal-lineage priming TPCS cluster.(F) Developmental trajectory 3, towards a multilineage-potent TPCS cluster.

Figure S13 .
Figure S13.CD8 + T cells in BLCA, related to Figure 7. (A) UMAP projection of CD8 + T scRNA cells.(B) Differential expression of marker genes among the CD8 + T cell clusters in (A).(C-D) Single T cells of globally-expanded TCR clones (blue, left) or locally-expanded TCR clones (orange, right).(E) Developmental trajectory (arrows) and single cell pseudotime of the CD8 + T cells.(F) Differential gene expression between trTeff-IFNG lo and trTeff-CXCR6 clusters.(G) Detailed differential expression of genes in naive T progenitors and tumor-reactive, experienced Teff cells.(H) Cell type prevalence of each CD8 + T cell types in sampled tissue.(I) Relative cell frequency in donors of different clinical stages.(J) Survival probability in IMvigor210 cohort of patients, stratified by trTeff-CXCR6 signature gene (TGFBR3, PBX4, CXCR6) expression.(K) trTeff-CXCR6 gene signature is specific for trTeff-CXCR6 cell type.(L) Distribution of trTeff-CXCR6 gene signature expression (x axis) in patients stratified by immunocheckpoint-inhibitor therapy treatment response.(M) Distribution of trTeff-CXCR6 gene signature expression (x axis) in patients stratified by histology-classified immune infiltration phenotype.

Figure S14 .
Figure S14.CD4 + T cells in BLCA, related to Figure 7. (A) UMAP projection of CD4 + T scRNA cells.(B) Differential expression of marker genes among the CD4 + T cell clusters.Known CD4 + T cell marker gene annotations are labeled by the bottom strips.(C) Single T cells of globally-expanded TCR clones (blue) or (D) locally-expanded TCR clones (orange).(E) Diffusion map of effector regulatory T cells (eTreg and trTreg) developing from naive Treg.(F) Developmental trajectory (arrows) from naive Treg through cTreg towards trTreg or eTreg.(G) Differential gene expression between cTreg and trTreg.(H) Differential gene expression between cTreg and eTreg.(I) Cell type prevalence of each CD4 + T cell types in sampled tissue.(J) Relative cell frequency in donors of different clinical stages.(K) Survival probability in IMvigor210 cohort of patients, stratified by trTeff-CXCR6 signature gene expression and trTreg signature gene expression.(L) Survival probability in IMvigor210 cohort of patients, stratified by trTeff-CXCR6 signature gene expression and eTreg signature gene expression.Signature genes are shown on the right.

Figure S15 .
Figure S15.NK cells in BLCA.(A) UMAP projection of NK scRNA cells.(B) Differential expression of marker genes among the NK cell clusters.

Figure S16 .
Figure S16.B cells in BLCA.(A) UMAP projection of B scRNA cells.(B) Differential expression of marker genes among the B cell clusters.

Figure S18 .
Figure S18.Analysis of fibroblast cells in scRNA dataset, related with Supplementary Notes.(A) Classification of scRNA fibroblast (FB) cells.Cell types are denoted by their specific marker gene expression.R: RGS5-positive; P: PDGFRA-positive.(B) Diffusion map inferred trajectory of fibroblasts.(C) Differential expression of marker genes in each fibroblast cell cluster.