Molecular subtypes of urothelial carcinoma are defined by specific gene regulatory systems

Background Molecular stratification of bladder cancer has revealed gene signatures differentially expressed across tumor subtypes. While these signatures provide important insights into subtype biology, the transcriptional regulation that governs these signatures is not well characterized. Methods In this study, we use publically available ChIP-Seq data on regulatory factor binding in order to link transcription factors to gene signatures defining molecular subtypes of urothelial carcinoma. Results We identify PPARG and STAT3, as well as ADIRF, a novel regulator of fatty acid metabolism, as putative mediators of the SCC-like phenotype. We link the PLK1-FOXM1 axis to the rapidly proliferating Genomically Unstable and SCC-like subtypes and show that differentiation programs involving PPARG/RXRA, FOXA1/GATA3 and HOXA/HOXB are differentially expressed in UC molecular subtypes. We show that gene signatures and regulatory systems defined in urothelial carcinoma operate in breast cancer in a subtype specific manner, suggesting similarities at the gene regulatory level of these two tumor types. Conclusions At the gene regulatory level Urobasal, Genomically Unstable and SCC-like tumors represents three fundamentally different tumor types. Urobasal tumors maintain an apparent urothelial differentiation axis composed of PPARG/RXRA, FOXA1/GATA3 and anterior HOXA and HOXB genes. Genomically Unstable and SCC-like tumors differ from Urobasal tumors by a strong increase of proliferative activity through the PLK1-FOXM1 axis operating in both subtypes. However, whereas SCC-like tumors evade urothelial differentiation by a block in differentiation through strong downregulation of PPARG/RXRA, FOXA1/GATA3, our data indicates that Genomically Unstable tumors evade differentiation in a more dynamic manner. Electronic supplementary material The online version of this article (doi:10.1186/s12920-015-0101-5) contains supplementary material, which is available to authorized users.


TABLE OF CONTENTS SUPPLEMENTAL FIGURES
. Pathways indicated in urothelial differentiation. Figure S2. Promoter binding overlap in the SCC-like downregulated genes. Figure S3. Genes involved in retinoic acid synthesis and degradation. Figure S4. Identification of the SCC-like subtype in the TCGA bladder cancer dataset. Figure S5. Transcription factors involved in urothelial differentiation are retained in Urobasal A tumors and partially lost in the Genomically Unstable and SCC-like tumors. Figure S6. Sample ranking by HOX switch score. Figure S7. The HOXA/HOXB switch in urothelial and breast cancer. Table S1. Transcription factor binding enrichment in promoters of downregulated genes of the SCC-like subtypes in bladder, and basal-like subtype in breast. Table S2. Ranking of transcription factors within the list of downregulated genes in the SCC-like and basal-like subtypes.  Figure S1. Pathways indicated in urothelial differentiation. Varley and co-workers [23] examined which intermediary factors that act downstream of PPARG, as uroplakins and cytokeratin expression are late events following PPARG activation in NHU cell cultures and does not have any strong peroxisome proliferator response element (PPRE) binding sites for the PPARG-RXRA complex itself. IRF-1 and FOXA1 does contain PPRE binding sites, and their expression was rapidly induced through simultaneous PPARG activation and EGFR inhibition. siRNA knockdown of either FOXA1 or IRF-1 abrogated the PPARGinduced induction of multiple late urothelial differentiation-associated genes. They also find GATA3 upregulated by the simultaneous activation and inhibition of PPARG and EGFR respectively, although at later time point. Böck and co-workers [22] examined the expressionchanges occurring during induced differentiation of cultured NHU by two separate methodologies in two separate cell cultures derived from two non-cancerous human urothelial tissue samples, and found that ELF3 was one of the shared transcription factors induced by both differentiation strategies. ELF3 knockdown resulted in reduced expression of FOXA1 and GRHL3 in response to induced differentiation through PPARG activation. Pretreatment with the PPARG antagonist T0070907 prior to PPARG activation blocked ELF3 mRNA and protein induction, confirming that ELF3 is a downstream target of PPARG. GATA3 was found among the top 1000 genes overexpressed by both differentiation strategies. Mauney and co-workers [36] used all-trans retinoic acid (RA) to stimulate ESC which led to downregulation of the pluripotency factor OCT-5 and a strong upregulation of UP1A, UP1B, UP2, UP3A, and UP3B. The expression of the characteristic uroplakins was found to be dependent on GATA4 and GATA6, also upregulated by RA stimulation, highlighting the importance of GATA factors for the ESC specification towards an urothelial lineage. Boskovic and Niles [39] showed that TBX2 is a direct target of retinoic acid signaling. An increase in TBX2 mRNA was seen within 2 hours after RA stimulation, indicating that no new protein synthesis was required for TBX2 induction. They identified retinoic acid response elements (RARE) in the promoter region of TBX2, and showed that these were functional using a luciferase assay vector with the TBX2 promoter sequence. Mutating the RARE sites reduced or abrogated the ability to respond to RA. The TBX2 promoter has ChIP-Seq binding sites for GATA3 and FOXA1 according to ENCODE/CISTROME data, and for RXRA, PPARG, FOXA1, and RARA according to Kittler data. Ballim and co-workers [38] provided evidence that TBX3 is a direct RAR-RXR target by cotransfecting cells with a RAR and RXR expressing vector along with a luciferase vector containing the TBX3 promoter. Addition of RA to the cell culture resulted in a significant increase in luciferase reporter activity. Retinoic acid response elements were also found within the TBX3 promoter sequence, which when mutated lowered the response to RA. The TBX3 promoter has ChIP-Seq binding sites for GATA3, FOXA1, PPARG according to ENCODE/CISTROME data, and for RXRA, FOXA1, and RARA according to Kittler data. While the specific role of the TBX genes in bladder cancer is yet unclear, methylation of both TBX2 and TBX3 was found to be associated with tumor progression in bladder cancer [40].

Figure S2. Promoter binding overlap in the SCC-like downregulated genes.
Co-occurrence of FOXA1, RXRA, GATA3 and PPARG binding in the promoter regions of the 829 SAM significant genes downregulated in the SCC-like subtype. ENCODE/CISTROME dataset (left), and the Kittler dataset (right).

Figure S3. Genes involved in retinoic acid synthesis and degradation.
Expression levels of ALDH1A2 and CYP26B1 involved in retinoic acid synthesis and degradation is altered in the SCC-like subtype. Figure S4. Identification of the SCC-like subtype in the TCGA bladder cancer dataset. (A) Gene expression data (RNASeqV2) from 242 samples was downloaded from the TCGA bladder urothelial carcinoma data portal. The dataset included 223 MI tumor samples and 19 "normal urothelial samples". Rsem estimate values for each sample was multiplied with 1 000 000 and an offset of 1 was added to each value. The data was log2 transformed, and the 19 normal samples were removed from the data matrix. The standard deviation for each gene was calculated, and the genes with a SD>90th percentile (2053 genes) were used to cluster the samples using the ConsensusClusterPlus package in R [24]. ConsensusClusterPlus settings were; maxK=10, reps=1000, pItem=0.8, pFeature=1, clusterAlg="hc", innerLinkage="ward", finalLinkage="ward", distance="pearson".

Figure S5. Transcription factors involved in urothelial differentiation are retained in Urobasal A tumors and partially lost in the Genomically Unstable and SCC-like tumors.
Tumor adjacent urothelium show nuclear staining of RXRA, PPARG, FOXA1, and GATA3. KRT5 indicate the basal cell population and the identity of SCC-like tumors. Figure S6. Sample ranking by HOX switch score. The HOX switch score was calculated by adding the average expression of the early HOXA genes (HOXA2, HOXA3, HOXA4, HOXA5, and HOXA6) and the average expression of the HOXB genes (HOXB2, HOXB3, HOXB4, HOXB5, HOXB6, HOXB7, HOXB8) and subtracting the average expression values of the late HOXA genes (HOXA9, HOXA10, HOXA11, HOXA13). After sorting the samples on their HOX switch score, an arbitrary cutoff point of 1 and -1 was set to define early and late HOX expressers. This resulted in 73 samples being defined as HOXA2-A6/HOXB2-B8 expressers, 114 as intermediate, and 79 samples as HOXA9-A12 expressers. A SAM analysis between the HOXA2-A6/HOXB2-B8 expressers and the HOXA9-A12 expressers resulted in 672 upregulated genes.     In-silico ChIP-Seq enrichment of the Late Cell Cycle QTC. * p-values were calculated using a one-tailed Fisher's exact test on the number of promoters with binding in a given gene list. † %-presence indicates the percentage of genes in each list that had binding sites. Table S4.

Gene lists from SAM and QTC analyses.
Sheet 1: SAM two class unpaired between SCC-like bladder tumors VS rest. Lund dataset.
Sheet 3: SAM two class unpaired between TCGA breast basal-like tumors VS rest.