Research Papers:

Predictive molecular biomarkers for determining neoadjuvant chemosensitivity in muscle invasive bladder cancer

PDF |  Full Text  |  Supplementary Files  |  How to cite  |  Press Release

Oncotarget. 2022; 13:1188-1200. https://doi.org/10.18632/oncotarget.28302

Metrics: PDF 706 views  |   Full Text 1820 views  |   ?  

Neal Murphy _, Andrew J. Shih, Paras Shah, Oksana Yaskiv, Houman Khalili, Anthony Liew, Annette T. Lee _ and Xin-Hua Zhu _


Neal Murphy1,2,*, Andrew J. Shih3,*, Paras Shah4, Oksana Yaskiv1,5, Houman Khalili3, Anthony Liew3, Annette T. Lee1,3,# and Xin-Hua Zhu1,2,#

1 Donald and Barbara Zucker School of Medicine at Hofstra/Northwell, Hempstead, NY 11549, USA

2 Northwell Health Cancer Institute, Lake Success, NY 11042, USA

3 Feinstein Institutes for Medical Research, Manhasset, NY 11030, USA

4 Mayo Clinic, Rochester, MN 55902, USA

5 Northwell Health Department of Pathology, Greenvale, NY 11548, USA

* These authors contributed equally to this work and share first authorship

# These authors contributed equally to this work and share last authorship

Correspondence to:

Neal Murphy, email: [email protected]
Annette T. Lee, email: [email protected]
Xin-Hua Zhu, email: [email protected]

Keywords: muscle invasive bladder cancer; neoadjuvant chemotherapy; gene expression; molecular subtyping; canonical correlation analysis

Received: July 25, 2022     Accepted: October 12, 2022     Published: November 02, 2022

Copyright: © 2022 Murphy et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.


Introduction: Identifying neoadjuvant chemotherapy (NAC) response in patients with muscle invasive bladder cancer (MIBC) has had limited success based on clinicopathological features and molecular subtyping. Identification of chemotherapy responsive cohorts would facilitate delivery to those most likely to benefit.

Objective: Develop a molecular signature that can identify MIBC NAC responders (R) and non-responders (NR) using a cohort of known NAC response phenotypes, and better understand differences in molecular pathways and subtype classifications between NAC R and NR.

Materials and Methods: Presented are the messenger RNA (mRNA) and microRNA (miRNA) differential expression profiles from initial transurethral resection of bladder tumor (TURBT) specimens of a discovery cohort of MIBC patients consisting of 7 known NAC R and 11 NR, and a validation cohort consisting of 3 R and 5 NR. Pathological response at time of cystectomy after NAC was used to classify initial TURBT specimens as R (pT0) versus NR (≥pT2). RNA and miRNA from FFPE blocks were sequenced using RNAseq and qPCR, respectively.

Results: The discovery cohort had 2309 genes, while the validation cohort had 602 genes and 13 miRNA differentially expressed between R and NR. Gene set enrichment analysis identified mitochondrial gene expression, DNA replication initiation, DNA unwinding in the R discovery cohort and positive regulation of vascular associated smooth muscle cell proliferation in the NR discovery cohort. Canonical correlation (CC) analysis was applied to differentiate R versus NR. 3 CCs (CC13, CC16, and CC17) had an AUC >0.65 in the discovery and validation dataset. Gene ontology enrichment showed CC13 as nucleoside triphosphate metabolic process, CC16 as cell cycle and cellular response to DNA damage, CC17 as DNA packaging complex. All patients were classified using established molecular subtypes: Baylor, UNC, CIT, Lund, MD Anderson, TCGA, and Consensus Class. The MD Anderson p53-like subtype, CIT MC4 subtype and Consensus Class stroma rich subtype had the strongest correlation with a NR phenotype, while no subtype had a strong correlation with the R phenotype.

Conclusions: Our results identify molecular signatures that can be used to differentiate MIBC NAC R versus NR, salient molecular pathway differences, and highlight the utility of molecular subtyping in relation to NAC response.


In 2022, the American Cancer Society estimates there will be 81,180 new cases of bladder cancer diagnosed in the United States and 17,100 people will die from the disease [1]. About half of these initial cases are superficial to the muscularis propria and have a 5-year overall survival of 95%, while about one third of these initial cases are found to be muscle invasive bladder cancer (MIBC) with a significant decrease in 5-year survival rate to 69% [2, 3]. For MIBC, neoadjuvant chemotherapy (NAC) prior to cystectomy is standard of care and provides a 5–8% increase in 5-year overall survival compared to cystectomy alone [4, 5]. In addition, patients who achieve a pathologic response to NAC have a 5-year survival rate of about 80–90% compared to 40–50% in non-responders [6]. However, a significant number of patients do not respond to NAC, with reported pT0 rates of 38% to neoadjuvant methotrexate, vinblastine, doxorubicin, and cisplatin (MVAC) and accelerated MVAC [6, 7], as well as pT0 rates of 42% to dose-dense MVAC (ddMVAC) and 36% to gemcitabine, cisplatin (GC) [8]. The NAC non-responders suffer from unnecessary adverse effects and a delay in time to cystectomy leading to worse overall survival [9, 10]. Subsequently, there remains a critical need to understand the molecular biology behind NAC responsiveness, in order to better tailor individual NAC therapy.

Gene expression profiling of MIBC holds promise for future individualization of therapy, and several studies have attempted to predict chemosensitivity using the molecular profile of tumors. Somatic mutations in the following DNA damage repair genes have been found to correlate with response to chemotherapy: ERCC2, ATM, RB1, and FANCC [11, 12]. This discovery has led to ongoing clinical trials designed to provide bladder sparing approaches after NAC using these biomarkers (NCT03609216). However, these mutations are only one part of the complex biological pathway driving bladder cancer, with <15% of patients harboring mutations in ATM and ERCC2 in the 2017 TCGA analysis [13]. Smith et al. used a microarray analysis to identify genes from the National Cancer Institute’s Developmental Therapeutics Program (NCI-DTP) whose expression was related to in vitro drug sensitivity, and then determined which of these genes maintained concordant expression with in vitro chemosensitivity analysis from 40 commonly used bladder cancer lines [14]. The SWOG S1314 Phase II study was designed to evaluate this approach, named Coexpression Extrapolation (COXEN), in order to determine if a treatment-specific COXEN score can predict NAC pathologic response [15]. Results from the trial indicated that the COXEN scores were not significantly prognostic for chemotherapy response in the individual treatment arms [16].

Moving one step further, transcriptomic profiling and unsupervised class discovery performed in several studies have found prognostic molecular subtypes [13, 1723]. These classifications have been derived from non-overlapping datasets using different methods, which has generated conflicting results and no general consensus on a definitive subtype classification. Furthermore, none of the studies have attempted to develop a classification based on known NAC response, which in our study is defined as a pathological complete response pT0 for chemotherapy responders (R) and ≥pT2 for chemotherapy non-responders (NR).

Therefore, our MIBC patient population with its known chemotherapy response phenotype represents a unique cohort to understand both the molecular mechanisms driving NAC response and to identify a molecular signature that truly correlates with NAC response. Here we present the differential mRNA and miRNA expression analysis of a discovery cohort of known chemotherapy responders versus non-responders. Expression differences are used to identify NAC responsiveness, which was further validated in a validation cohort of eight patients, as well as identify relevant molecular pathways in NAC non-responders and responders, and examine the correlation between established molecular subtypes and NAC response.



Patients from 2008–2018 diagnosed with MIBC on a transurethral resection of bladder tumor (TURBT) were selected from the Northwell Health tumor bank. Patients received either 4 cycles of neoadjuvant GC or 3 or more cycles of ddMVAC chemotherapy. Pathological response at time of cystectomy was used to classify initial TURBT specimens as NAC responders with pathological complete response (pT0) versus non-responders with muscle invasive bladder cancer (≥pT2).

A total of 26 patients with known NAC response were identified, and inclusion in either the discovery or validation cohort was chosen at random. The discovery cohort consisted of 18 patients: 11 NR and 7 R. The average age of the entire cohort was 65, with about half receiving ddMVAC and half receiving GC. The majority of the patients were male. The NR group had an average time to recurrence of 13 months. The validation cohort consisted of 8 patients: 5 NR and 3 R. The average age of the entire cohort was 67 with the majority being male patients (see Table 1).

Table 1: Patient demographics

Discovery Cohort
Path StatusGenderAge at TURBTNACCys Path StageTime to Recur (Months)
NR6F72GCpT2bN0no recurrence
R1M64GCpTisN0no recurrence
R2F56GCpT0N0no recurrence
R3F57ddMVACpT0N0no recurrence
R4M59N/ApT0N0no recurrence
R5M65N/ApT0N0no recurrence
R6F70N/ApT0N0no recurrence
R7M75GCpTisN0no recurrence
Validation Cohort
NR13M52N/ApT2bN0no recurrence
R8M66ddMVACpT0N0no recurrence
R9M67GCpT0N0no recurrence
R10F57GCpT0N0no recurrence

RNAseq expression analysis

In the discovery cohort, 2039 genes had significant expression differences between the R vs. NR groups after accounting for multiple testing corrections (See Figure 1 for Top 20 genes in both cohorts and Supplementary Table 1 for the complete list of genes). Of those genes, a slight majority were seen in the R group (1,041 genes) compared to the NR group (998 genes). In the validation cohort, 602 genes had significant expression differences between the R vs. NR groups (p-value < 0.01, data included in Supplementary Table 2), with the majority seen in the R group (509 genes) compared to NR (93 genes).

Volcano plot of significant differentially expressed genes with the left volcano being the discovery cohort and validation cohort on the right.

Figure 1: Volcano plot of significant differentially expressed genes with the left volcano being the discovery cohort and validation cohort on the right. For each volcano plot, R are on the left side nand NR on the right side pf each volcano. X-axis is the log 2 fold change between groups while the y-axis is the -log 10 raw p-value. The top 20 genes by p-value are labeled.

Pathway enrichment analysis using GSEA (gene set enrichment analysis) identified 519 gene sets in the R discovery cohort (see Table 2 for top 12 biologically relevant gene sets for R and NR) and 155 gene sets in the R validation cohort with a p-value < 0.01. 33 of these gene sets overlapped. For the NR discovery cohort, there were 60 gene sets and 354 gene sets in the NR validation cohort with a p-value < 0.01. 18 of these gene sets overlapped (Supplementary Table 3).

Table 2: GSEA analysis, Top 12 pathways for the R discovery cohort on the top and NR discovery cohort on the bottom, based on p-value < 0.01


miRNA expression analysis

In the discovery cohort, no miRNA were found to meet statistical significance (Supplementary Table 4), however the validation cohort had 13 miRNAs that met a p-value threshold of .05. miR-18a was upregulated in the R group while the remaining 12 were seen in the NR group (Supplementary Table 5).

Canonical correlation analysis in the discovery and validation cohorts

Canonical correlation analysis (CCA) is used to determine and measure the association amongst two different sets of variables and redefine them as canonical components (CCs). Whereas Principal Component Analysis (PCA) focuses on finding linear combinations that contribute to the most variance in a data set, CCA focuses on finding linear combinations that correlate the most between two datasets. In essence, CCs are de novo combinations of genes and miRNA that have correlated expression patterns, and each CC can be seen as a co-expression network with its own molecular profile with the potential to influence a clinical outcome.

A total of 17 CCs were identified (Supplementary Table 6). To classify patients as R versus NR, univariate logistic regression was performed on the CCs. 3 CCs (CC13, CC16, and CC17) have an AUC >0.65 in both the discovery and validation datasets, with the top performer CC16 having an AUC of 0.81 in the discovery dataset (Table 3). Gene set enrichment showed CC13 as nucleoside triphosphate metabolic process and cell envelope, CC16 as cell cycle and cellular response to DNA damage and stress, CC17 as DNA packaging complex.

Table 3: Canonical components (CCs) with AUC >0.65 in both the discovery and validation cohorts, based on a p-value < 0.05

Canonical ComponentDiscovery AUCValidation AucGO Terms
CC 160.8180.667CELL CYCLE

Molecular subtype analysis

All 26 patients were classified using the following molecular subtyping: Baylor, UNC, CIT, Lund, MD Anderson, TCGA, and Consensus Class [13, 1720, 22]. The Baylor classification is divided into basal and differentiated. The UNC classification is split into basal-like and luminal. The CIT classification has MC1-7, while the Lund subtype breaks down into: urobasal A, genomically unstable, urobasal B, squamous cell carcinoma like, and infiltrated. MD Anderson (MDA) is divided into luminal, basal and p53-like. The TCGA has 5 subtypes: luminal infiltrated, luminal papillary, luminal, basal squamous, and neuronal. Lastly, the Consensus Classification was an international effort that reconciled six previously mentioned classification scheme into a new consensus classification that was divided into six subtypes: luminal papillary, luminal nonspecified, luminal unstable, stroma-rich, basal/squamous, and neuroendocrine-like.

Table 4 includes the subtype classification for each patient using the aforementioned 6 different molecular classification schemes. For the discovery cohort, 6 out of 11 NR were p53-like (MDA), which is known to be a chemo-resistant phenotype. In addition, the majority of the six patients who were classified as p53-like (MDA) in the NR discovery group were also defined as basal (UNC), MC4 (CIT), and stroma-rich (Consensus Class). Of the remaining patients in the NR discovery cohort, 3 were defined as the MDA basal subtype and 2 were defined as the luminal subtype. The R discovery cohort had one chemotherapy resistant p53-like patient (MDA), while the remaining 6 patients were evenly split between the MDA basal and luminal subtypes.

Table 4: Patient molecular subtype classifications, as defined by the baylor, UNC, CIT, lund, MD anderson, TCGA, and consensus class classifications

Discovery Cohort
Path StatusBaylorUNCCITLundMDATCGAConsensus Class
Validation Cohort

For the validation cohort, 3 out of the 5 NR patients were p53-like (MD Anderson), as well as differentiated (Baylor), basal (UNC), MC4 (CIT), luminal-infiltrated (TCGA), and stroma-rich (Consensus Class). The remaining 2 patients in the NR validation cohort were split between luminal and basal (MDA), while of the three patients in the R validation, 2 were defined as MDA luminal and 1 as basal.

When looking at NR for both discovery and validation combined, subtypes that correlated the most with being a NR include: basal (Baylor) 63% (7/11), differentiated (Baylor) 57% (8/14), basal (UNC) 68% (11/16), MC4 (CIT) 90% (10/11), p53-like (MDA) 90% (9/10), luminal-infiltrated (TCGA) 72% (8/11), basal-squamous (TCGA) 60% (6/10) and stroma rich (Consensus Class) 100% (8/8). Overall, the NRs from both cohorts have a significant portion of patients who are MC4 (CIT), p53-like (MD Anderson), and stroma-rich (Consensus Class). When examining the known chemotherapy resistant p53-like subtype within our entire cohort more closely, using Fisher’s exact test with a 2 × 2 contingency table (9 p53-like NRs and 1 p53-like R versus 7 non p53-like NRs and 9 non p53-like Rs), shows the p53-like subtype is significantly associated with being a NR (p-value of .03674, Odds Ratio 10.5). Lastly, for both the R discovery and validation cohorts, the subtypes for each classification were mixed without any definitive subtype making up a significant majority for each molecular classification.


The preceding analysis combines differential mRNA and miRNA expression in known MIBC phenotypes to better understand the molecular driver behind NAC response. The majority of molecular classification and subtyping to date has been performed with variations in gene expression technology, cohorts that have varied in size, and has primarily relied on unsupervised hierarchical clustering without factoring in NAC response. We report significant gene sets associated with NAC response phenotype, as well as three multigene and miRNA signatures generated by CCA that can be used to potentially classify NAC response.

To our knowledge, our study is the first to use combined differential mRNA and miRNA expression in MIBC to identify a NAC response signature. The mitochondrial and metabolic signature seen in CC13, the cell cycle signature and response to DNA damage seen in CC16, and the DNA packaging signature in CC17 all appear to be driven by the responders in the discovery cohort given the overlap seen when comparing these CC signatures to the GSEA results for this cohort (i.e., Gene Ontology Biological Process (GOBP)_mitochondrial gene expression, GOBP_DNA replication initiation, GOBP_DNA Unwinding involved in DNA Replication, seen in Supplementary Table 3).

Biologically, these CC and GOBP signatures make sense given MIBC’s historical success with cisplatin based chemotherapy regimens, in which cisplatin causes DNA damage, blocks cell division and triggers apoptotic cell death [24]. When looking at GOBP_mitochondrial gene expression, it is known that mitochondria take part in short and long patch base excision repair, which is achieved by DNA polymerase γ. This is critical for the DNA repair pathway to maintain whole genome stability, as well as being a primary source for reactive oxygen species (ROS) and regulating apoptosis [25, 26]. Within this pathway, several mammalian mitochondrial ribosomal small subunit (MRPS) genes are upregulated in the R discovery cohort, including MRPS12, MRPS34, MRPS28, MRPS14, and MRPS2. Upregulation of these genes in vitro may restore chemosensitivity in resistant bladder cancer cell lines, potentially by restoring the cells ability to undergo apoptosis.

Another intriguing pathway identified in the R discovery cohort was GOBP_Sulfur Compound Catabolic Process. Hydrogen sulfide inhibition in one study was shown to enhance cisplatin induced apoptosis both in vitro and in vivo, and exogenous administration in another study enhanced in vitro cell proliferation [27, 28]. Overall, there appears to be a metabolic shift in the R discovery cohort that likely contributed to chemosensitivity, as evident by other significant pathways expressed: GOBP_Glucosamine Containing Compound Metabolic Process, GOBP_Threonine Type Peptidase Activity, GOBP Isoprenoid Biosynthetic Process, and GOBP_ Cellular Response to pH.

When further examining GOBP_DNA replication initiation and GOBP_DNA Unwinding involved in DNA Replication , the minichromosome maintenance complex (MCM) genes play a crucial role in both DNA replication and unwinding. MCM2-7 forms a helicase complex that contributes to both the initiation and elongation phase of DNA replication [29]. MCM2-3 and 5–6 are significantly upregulated in the R discovery cohort. Although the knockdown of MCM genes has been linked to decreased cellular proliferation [30], increased expression could theoretically promote cells to undergo ineffective DNA replication due to platinum DNA intercalation that eventually leads to cell death. XPA is also substantially enhanced in the NR discovery cohort. It is a zinc finger protein that plays a central role in nucleotide excision repair (NER), which is responsible for repair of DNA adducts induced by cisplatin chemotherapy. The XPA protein binds to DNA and several NER proteins, acting as a scaffold to assemble the NER incision complex at sites of DNA damage [31]. Targeting MCM and XPA genes in vitro may be another potential candidate for inducing chemotherapy response.

Expanding upon the previous two pathways, DNA packaging plays a critical role in many nucleic acid processes, including transcription, DNA repair, and DNA replication. When DNA damage is induced by cisplatin, chromatin encounters dramatic changes in response to DNA damage both locally and genome-wide, which initiates the DNA repair process [32]. Genes implicated in the DNA packaging pathway may potentially change cisplatin response through change in their expression. ELK4 is markedly upregulated in the NR discovery cohort. This gene is a member of the Ets family of transcription factors, whose functions are DNA-binding transcription factor activity and chromatin binding [33]. FOXA3 is dramatically increased in the R discovery and validation cohort. This gene is a member of the forkhead class of DNA-binding proteins and also interacts with chromatin [34]. These chromatin interactive genes involved in DNA packaging will be further studied to predict cisplatin response.

Although no miRNA had statistically significant differential expression in the discovery cohort, several miRNAs significantly contribute to the three relevant CC signatures. The following miRNAs listed for each CC may be potential candidates for modulating chemotherapy sensitivity or serve as therapeutic targets. For CC13, miR-192 and mi-194 overexpression has been shown to inhibit cell proliferation in bladder cancer cells in vitro [35, 36]. In CC16, miR-15 has been shown to inhibit bladder cancer cell proliferation, migration and invasion in vitro by targeting BMI1 through the PI3K/AKT pathway [37]. Also, one study found that miR-34 upregulated PTEN and inhibited bladder cancer cell migration and invasion [38]. For CC17, miR-455-5p has been shown to induce cisplatin resistance in bladder cancer cells by regulating Notch1 [39].

Gene expression profiling of MIBC holds promise for future individualization of therapy, and several studies have attempted to predict chemosensitivity using the molecular profile of tumors. Two studies using whole-exome-sequencing and DNA sequencing targeting 287 cancer-related genes respectively, found that alterations in DNA repair genes may correlate with response to chemotherapy: ERCC2, ATM, RB1, and FANCC [11, 12]. Of those 4 genes, FANCC was the only one found to have significant expression, albeit in our NR discovery cohort.

Lastly, of note, two studies used microarray analysis of TURBTs to create a prediction scoring system that can potentially determine response to GC or MVAC NAC, respectively [40, 41]. Although these two studies created their prediction scoring system based on known NAC response, responders in these studies were defined as having had at least a partial response (either >60% tumor shrinkage or ≤pT1 after 2 cycles of chemotherapy) unlike our cohort where the majority of the responders were pT0 (we included 1 pTis) and response was determined after 4 cycles of chemotherapy (standard of care) [42]. Interestingly enough, 2 of the genes used in the prediction scoring system developed by Kato et al. had concordant significant expression in our discovery cohort: Spry1 was increased in NR and CRKL increased in R. PNPO was increased in our R discovery cohort, while it was increased in the NR cohort form Kato et al. [40]. None of the genes from the Takata et al. scoring system had significant expression in our discovery cohort.

In addition to the aforementioned molecular analysis, subtyping patients based on molecular expression is another promising approach that has potential for determining chemosensitivity. However, the TCGA and several other studies have identified molecular expression subtypes using hierarchical classification techniques that to date have served more as biological classifications rather than predictors of clinical outcome [13, 1722]. In the TCGA analysis, predicted chemotherapy response is hypothesis driven from patients with incomplete or unknown treatment history, and not derived from comparing molecular phenotypes between NAC responders and non-responders. The basal/squamous subtype is predicted to be the most chemotherapy sensitive, with luminal-infiltrated and luminal-papillary having lower response rates to NAC. However in the NR discovery cohort, 45% (5/11) were basal/squamous and 42% (3/7) of the R discovery cohort were either luminal infiltrated or luminal-papillary. Despite there being some basal/squamous patients in the NR discovery cohort, basal expression does appear to play a significant role in chemotherapy response as shown by the top two GSEA pathways in the R discovery cohort: GOBP_Keratinization and GOBP_Keratinocyte Differentiation. Significant genes within these two pathways are another candidate for further exploration to better characterize chemotherapy response.

Expanding upon established subtyping analysis, Seiler et al. used 149 genes from the subtypes identified by hierarchical clustering from four studies (TCGA, Lund, MDA, UNC), to develop a single-sample genomic subtyping classifier (GSC) that showed patient outcomes after NAC varies by four molecular subtypes: basal, claudin-low, luminal, or luminal infiltrated subtype [23]. The basal subtype benefited most from NAC, whereas the other three subtypes did not demonstrate a significant survival advantage with NAC. Although the four subtypes are associated with overall survival predictions after NAC, the subtypes could not establish which patients would have a pathological response to NAC. A lack of correlation between NAC pathological response and overall survival is a limitation that conflicts with previously mentioned studies, although the study may have been underpowered to find this correlation.

In contrast to the TCGA and Seiler et al. analysis, Choi et al. in their initial discovery cohort and in a subsequent phase 2 follow-up study, found that p53-like tumors were mostly resistant to NAC [21, 43]. 54% (6/11) of our NR discovery cohort were classified as p53-like, and when looking at all of the patients classified as p53-like in both the discovery and validation cohorts combined, 90% (9/10) were classified as NR. This correlation (being classified as p53-like and being a NR) was statistically significant in our cohort. Overall, this finding strengthens both the validity of the Choi et al. p53-like subtype and our own molecular expression analysis. Similar to what Choi et al. reported, the MDA basal and luminal subtypes for our two cohorts are mixed between NR and R. One patient in the R discovery cohort was defined as p53-like, which raises the point that relying on one classification system may have limited potential. However, given that MC4 (CIT) and stroma rich (Consensus Class) were also highly correlated with a NR phenotype, tumors that express these 2 subtype classifications, as well as p53-like may have the potential to be confidently classified as a NR. This may be worth assessing and validating in a larger follow-up study.

Limitations of our study include a small patient size in both the discovery and validation cohorts, as well as intratumoral heterogeneity that may limit reproducibility of results and subtype analysis. However, the overlap of a known chemoresistant molecular subtype (MDA’s p53-like subtype) and our cohort of NR, gives validity to our small sample size. Furthermore, CCA accounts for intratumor heterogeneity when creating de novo mRNA and miRNA modules across different samples and in effect theoretically limits heterogeneity. Another limitation is the decreased number of significant differentially expressed genes in the validation cohort as compared to the discovery cohort. We suspect that this is largely due to the fewer number of patients in the validation cohort.

In conclusion, our results identify molecular signatures that can be used to differentiate MIBC NAC responders versus non-responders. We have presented the salient molecular pathways and relevant genes, including mitochondrial response gene expression (MRPS12, MRPS34, MRPS28, MRPS14, and MRPS2), DNA replication initiation, and DNA unwinding and DNA damage (MCM2-3, MCM5-6 and XAP , ELK4, and FOXA3) that can be further analyzed to better understand NAC response. The above mentioned genes derived from their respective three pathways may be selected as part of a NAC response biomarker panel. In addition, we have highlighted the utility of molecular subtyping in relation to NAC response. If validated in a larger cohort, these findings may help deliver chemotherapy to those patients most likely to respond.

Materials and Methods

Patient selection

A total of 26 patients with known NAC response were identified, and inclusion in either the discovery or validation cohort was chosen at random. The discovery cohort consisted of 7 NAC responders and 11 non-responders, while the validation cohort consisted of 3 responders and 5 non-responders. The Northwell Health System Institutional Review Board and Regional Ethics Committee granted research approval for the study with waivers of HIPAA authorization and informed consent. Clinical investigation was conducted according to the principles outlined in the Declaration of Helsinki. TURBT specimens from the Northwell Health pathology department were received as formalin-fixed, paraffin-embedded (FFPE) tissue blocks. Pathologic response at the time of cystectomy was determined using the American Joint Committee on Cancer staging [44].

RNA extraction

A genitourinary oncology pathologist at Northwell Health reviewed all FFPE tissue blocks and their corresponding H&E slides, outlining neoplastic tissue on the H&E slide. Tissue blocks were matched up with their respective slides and approximately 35 mg of pre-identified tumor was cut away from the blocks using a scalpel. RNA was extracted from FFPE tissue using the RecoverAll™ Total Nucleic Acid Isolation Kit with quality control performed using an AB Bioanalyzer.

mRNA expression analysis

RNA was sequenced using the Illumina TruSeq RNA Access Library Prep Guide and NextSeq 500. Sequenced segments were aligned using the STAR2 aligner [45]. Gene counts were assessed using ht-seq counts [46]. Differential expression was calculated using DESeq2 [47]. Adjustment for multiple corrections was performed using Benjamini-Hochberg method. Pathway enrichment analysis was performed using GSEA [48, 49].

miRNA expression analysis

Differential expression of 754 miRNAs was analyzed using the TaqMan Open Array miRNA qPCR panel. A Ct threshold of 30 was used and miRNAs were filtered by removing miRNAs detected in <10% of samples. The top 10 most stable miRNAs from the NormPCR package were used as housekeeping miRNAs, which were used to normalize raw Ct values with deltaCt normalization [50, 51]. The limma package was used to calculate differential expression [52].

Canonical correlation analysis and validation

Sparse canonical correlation analysis was used to identify de novo mRNA and miRNA Canonical Components (CCs) with AUC to NAC response was determined using pROC [53, 54]. Data analysis was performed using R and tidyverse [55, 56]. Gene Ontology (GO) enrichment was performed using a Fisher’s test.

Subtype analysis

Subtyping was performed using the BCLAsubtyping R package created by Kamoun et al. [57].


MIBC: muscle invasive bladder cancer; NAC: neoadjuvant chemotherapy; mRNA: messenger RNA; miRNA: microRNA; NCI-DTP: National Cancer Institute’s Developmental Therapeutics Program; COXEN: Coexpression Extrapolation; R: responders; NR: non-responder; TURBT: transurethral resection of bladder tumor; GC: gemcitabine cisplatin; ddMVAC: dose-dense methotrexate vincristine adriamycin, cisplatin; CCA: canonical correlation analysis; CCs: canonical components; PCA: Principal Component Analysis; GSEA: gene set enrichment analysis; MDA: MD Anderson; GOBP: Gene Ontology Biological Process; ROS: reactive oxygen species; MRPS: mitochondrial ribosomal small subunit; MCM: minichromosome maintenance complex; NER: nucleotide excision repair; GSC: genomic subtyping classifier.

Author contributions

NM, AS, PS, ATL and XZ conceived and designed the study. OY, HK, NM, and AL acquired the data. AS did the statistical analysis. AS, NM and XZ analyzed and interpreted the data. NM, AS, and XZ drafted the manuscript. NM, AS, PS, XZ and ATL critically revised the manuscript for important intellectual content. OY, HK, AL, ATL, and XZ provided administrative, technical, and material support. XZ obtained funding and supervised the study.


The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflicts of interest.

Ethical statement and consent

The Northwell Health System Institutional Review Board and Regional Ethics Committee granted research approval for the study with waivers of HIPAA authorization and informed consent. Clinical investigation was conducted according to the principles outlined in the Declaration of Helsinki.


Funding for this project was done with institutional funds.


1. Key Statistics for Bladder Cancer. Available 2022 Feb 19, from https://www.cancer.org/cancer/bladder-cancer/about/key-statistics.html.

2. Cancer Facts & Figures 2019. American Cancer Society. Available 2019 Apr 28, from https://www.cancer.org/content/dam/cancer-org/research/cancer-facts-and-statistics/annual-cancer-facts-and-figures/2019/cancer-facts-and-figures-2019.pdf.

3. Bladder Cancer - Statistics. Cancer Net. 2012. Available 2019 May 5, from https://www.cancer.net/cancer-types/bladder-cancer/statistics.

4. Advanced Bladder Cancer (ABC) Meta-analysis Collaboration. Neoadjuvant chemotherapy in invasive bladder cancer: update of a systematic review and meta-analysis of individual patient data advanced bladder cancer (ABC) meta-analysis collaboration. Eur Urol. 2005; 48:202. https://doi.org/10.1016/j.eururo.2005.04.006. [PubMed].

5. Yin M, Joshi M, Meijer RP, Glantz M, Holder S, Harvey HA, Kaag M, Fransen van de Putte EE, Horenblas S, Drabick JJ. Neoadjuvant Chemotherapy for Muscle-Invasive Bladder Cancer: A Systematic Review and Two-Step Meta-Analysis. Oncologist. 2016; 21:708–15. https://doi.org/10.1634/theoncologist.2015-0440. [PubMed].

6. Grossman HB, Natale RB, Tangen CM, Speights VO, Vogelzang NJ, Trump DL, deVere White RW, Sarosdy MF, Wood DP Jr, Raghavan D, Crawford ED. Neoadjuvant chemotherapy plus cystectomy compared with cystectomy alone for locally advanced bladder cancer. N Engl J Med. 2003; 349:859–66. https://doi.org/10.1056/NEJMoa022148. [PubMed].

7. Plimack ER, Hoffman-Censits JH, Viterbo R, Trabulsi EJ, Ross EA, Greenberg RE, Chen DY, Lallas CD, Wong YN, Lin J, Kutikov A, Dotan E, Brennan TA, et al. Accelerated methotrexate, vinblastine, doxorubicin, and cisplatin is safe, effective, and efficient neoadjuvant treatment for muscle-invasive bladder cancer: results of a multicenter phase II study with molecular correlates of response and toxicity. J Clin Oncol. 2014; 32:1895–901. https://doi.org/10.1200/JCO.2013.53.2465. [PubMed].

8. Pfister C, Gravis G, Fléchon A, Soulié M, Guy L, Laguerre B, Mottet N, Joly F, Allory Y, Harter V, Culine S, and VESPER Trial Investigators. Randomized Phase III Trial of Dose-dense Methotrexate, Vinblastine, Doxorubicin, and Cisplatin, or Gemcitabine and Cisplatin as Perioperative Chemotherapy for Patients with Muscle-invasive Bladder Cancer. Analysis of the GETUG/AFU V05 VESPER Trial Secondary Endpoints: Chemotherapy Toxicity and Pathological Responses. Eur Urol. 2021; 79:214–21. https://doi.org/10.1016/j.eururo.2020.08.024. [PubMed].

9. Chu AT, Holt SK, Wright JL, Ramos JD, Grivas P, Yu EY, Gore JL. Delays in radical cystectomy for muscle-invasive bladder cancer. Cancer. 2019; 125:2011–17. https://doi.org/10.1002/cncr.32048. [PubMed].

10. Türk H, Ün S, Cinkaya A, Kodaz H, Parvizi M, Zorlu F. Effect of delayed radical cystectomy for invasive bladder tumors on lymph node positivity, cancer-specific survival and total survival. Tumori. 2018; 104:434–37. https://doi.org/10.5301/tj.5000626. [PubMed].

11. Van Allen EM, Mouw KW, Kim P, Iyer G, Wagle N, Al-Ahmadie H, Zhu C, Ostrovnaya I, Kryukov GV, O’Connor KW, Sfakianos J, Garcia-Grossman I, Kim J, et al. Somatic ERCC2 mutations correlate with cisplatin sensitivity in muscle-invasive urothelial carcinoma. Cancer Discov. 2014; 4:1140–53. https://doi.org/10.1158/2159-8290.CD-14-0623. [PubMed].

12. Plimack ER, Dunbrack RL, Brennan TA, Andrake MD, Zhou Y, Serebriiskii IG, Slifker M, Alpaugh K, Dulaimi E, Palma N, Hoffman-Censits J, Bilusic M, Wong YN, et al. Defects in DNA Repair Genes Predict Response to Neoadjuvant Cisplatin-based Chemotherapy in Muscle-invasive Bladder Cancer. Eur Urol. 2015; 68:959–67. https://doi.org/10.1016/j.eururo.2015.07.009. [PubMed].

13. Robertson AG, Kim J, Al-Ahmadie H, Bellmunt J, Guo G, Cherniack AD, Hinoue T, Laird PW, Hoadley KA, Akbani R, Castro MAA, Gibb EA, Kanchi RS, et al, and TCGA Research Network. Comprehensive Molecular Characterization of Muscle-Invasive Bladder Cancer. Cell. 2017; 171:540–56.e25. https://doi.org/10.1016/j.cell.2017.09.007. [PubMed].

14. Smith SC, Baras AS, Lee JK, Theodorescu D. The COXEN principle: translating signatures of in vitro chemosensitivity into tools for clinical outcome prediction and drug discovery in cancer. Cancer Res. 2010; 70:1753–58. https://doi.org/10.1158/0008-5472.CAN-09-3562. [PubMed].

15. S1314, Co-expression Extrapolation (COXEN) Program to Predict Chemotherapy Response in Patients With Bladder Cancer - Full Text View - ClinicalTrials.gov. Available 2019 May 5, from https://clinicaltrials.gov/ct2/show/NCT02177695.

16. Flaig TW, Tangen CM, Daneshmand S, Alva AS, Lerner SP, Lucia MS, McConkey DJ, Theodorescu D, Goldkorn A, Milowsky MI, Bangs RC, MacVicar GR, Bastos BR, et al. SWOG S1314: A randomized phase II study of co-expression extrapolation (COXEN) with neoadjuvant chemotherapy for localized, muscle-invasive bladder cancer. J Clin Oncol. 2019; 37:4506. https://doi.org/10.1200/JCO.2019.37.15_suppl.4506.

17. Mo Q, Nikolos F, Chen F, Tramel Z, Lee YC, Hayashi K, Xiao J, Shen J, Chan KS. Prognostic Power of a Tumor Differentiation Gene Signature for Bladder Urothelial Carcinomas. J Natl Cancer Inst. 2018; 110:448–59. https://doi.org/10.1093/jnci/djx243. [PubMed].

18. Damrauer JS, Hoadley KA, Chism DD, Fan C, Tiganelli CJ, Wobker SE, Yeh JJ, Milowsky MI, Iyer G, Parker JS, Kim WY. Intrinsic subtypes of high-grade bladder cancer reflect the hallmarks of breast cancer biology. Proc Natl Acad Sci U S A. 2014; 111:3110–15. https://doi.org/10.1073/pnas.1318376111. [PubMed].

19. Rebouissou S, Bernard-Pierrot I, de Reyniès A, Lepage ML, Krucker C, Chapeaublanc E, Hérault A, Kamoun A, Caillault A, Letouzé E, Elarouci N, Neuzillet Y, Denoux Y, et al. EGFR as a potential therapeutic target for a subset of muscle-invasive bladder cancers presenting a basal-like phenotype. Sci Transl Med. 2014; 6:244ra91. https://doi.org/10.1126/scitranslmed.3008970. [PubMed].

20. Sjödahl G, Lauss M, Lövgren K, Chebil G, Gudjonsson S, Veerla S, Patschan O, Aine M, Fernö M, Ringnér M, Månsson W, Liedberg F, Lindgren D, Höglund M. A molecular taxonomy for urothelial carcinoma. Clin Cancer Res. 2012; 18:3377–86. https://doi.org/10.1158/1078-0432.CCR-12-0077-T. [PubMed].

21. Choi W, Porten S, Kim S, Willis D, Plimack ER, Hoffman-Censits J, Roth B, Cheng T, Tran M, Lee IL, Melquist J, Bondaruk J, Majewski T, et al. Identification of distinct basal and luminal subtypes of muscle-invasive bladder cancer with different sensitivities to frontline chemotherapy. Cancer Cell. 2014; 25:152–65. https://doi.org/10.1016/j.ccr.2014.01.009. [PubMed].

22. Kamoun A, de Reyniès A, Allory Y, Sjödahl G, Robertson AG, Seiler R, Hoadley KA, Groeneveld CS, Al-Ahmadie H, Choi W, Castro MAA, Fontugne J, Eriksson P, et al, and Bladder Cancer Molecular Taxonomy Group. A Consensus Molecular Classification of Muscle-invasive Bladder Cancer. Eur Urol. 2020; 77:420–33. https://doi.org/10.1016/j.eururo.2019.09.006. [PubMed].

23. Seiler R, Ashab HAD, Erho N, van Rhijn BWG, Winters B, Douglas J, Van Kessel KE, Fransen van de Putte EE, Sommerlad M, Wang NQ, Choeurng V, Gibb EA, Palmer-Aronsten B, et al. Impact of Molecular Subtypes in Muscle-invasive Bladder Cancer on Predicting Response and Survival after Neoadjuvant Chemotherapy. Eur Urol. 2017; 72:544–54. https://doi.org/10.1016/j.eururo.2017.03.030. [PubMed].

24. Dasari S, Tchounwou PB. Cisplatin in cancer therapy: molecular mechanisms of action. Eur J Pharmacol. 2014; 740:364–78. https://doi.org/10.1016/j.ejphar.2014.07.025. [PubMed].

25. Chatterjee N, Walker GC. Mechanisms of DNA damage, repair, and mutagenesis. Environ Mol Mutagen. 2017; 58:235–63. https://doi.org/10.1002/em.22087. [PubMed].

26. Cormio A, Sanguedolce F, Musicco C, Pesce V, Calò G, Bufo P, Carrieri G, Cormio L. Mitochondrial dysfunctions in bladder cancer: Exploring their role as disease markers and potential therapeutic targets. Crit Rev Oncol Hematol. 2017; 117:67–72. https://doi.org/10.1016/j.critrevonc.2017.07.001. [PubMed].

27. Wahafu W, Gai J, Song L, Ping H, Wang M, Yang F, Niu Y, Xing N. Increased H2S and its synthases in urothelial cell carcinoma of the bladder, and enhanced cisplatin-induced apoptosis following H2S inhibition in EJ cells. Oncol Lett. 2018; 15:8484–90. https://doi.org/10.3892/ol.2018.8373. [PubMed].

28. Liu H, Chang J, Zhao Z, Li Y, Hou J. Effects of exogenous hydrogen sulfide on the proliferation and invasion of human Bladder cancer cells. J Cancer Res Ther. 2017; 13:829–32. https://doi.org/10.4103/jcrt.JCRT_423_17. [PubMed].

29. Gene group. Available 2022 Feb 20, from https://www.genenames.org/data/genegroup/#!/group/1085.

30. Yu S, Wang G, Shi Y, Xu H, Zheng Y, Chen Y. MCMs in Cancer: Prognostic Potential and Mechanisms. Anal Cell Pathol (Amst). 2020; 2020:3750294. https://doi.org/10.1155/2020/3750294. [PubMed].

31. Borszéková Pulzová L, Ward TA, Chovanec M. XPA: DNA Repair Protein of Significant Clinical Importance. Int J Mol Sci. 2020; 21:2182. https://doi.org/10.3390/ijms21062182. [PubMed].

32. Hauer MH, Gasser SM. Chromatin and nucleosome dynamics in DNA damage and repair. Genes Dev. 2017; 31:2204–21. https://doi.org/10.1101/gad.307702.117. [PubMed].

33. GeneCards Human Gene Database. Available 2022 Apr 1, from https://www.genecards.org/cgi-bin/carddisp.pl?gene=ELK4.

34. GeneCards Human Gene Database. FOXA3 gene - GeneCards. Available 2022 Apr 1, from https://www.genecards.org/cgi-bin/carddisp.pl?gene=FOXA3.

35. Jin Y, Lu J, Wen J, Shen Y, Wen X. Regulation of growth of human bladder cancer by miR-192. Tumour Biol. 2015; 36:3791–97. https://doi.org/10.1007/s13277-014-3020-8. [PubMed].

36. Zhang M, Zhuang Q, Cui L. MiR-194 inhibits cell proliferation and invasion via repression of RAP2B in bladder cancer. Biomed Pharmacother. 2016; 80:268–75. https://doi.org/10.1016/j.biopha.2016.03.026. [PubMed].

37. Zhang L, Wang CZ, Ma M, Shao GF. MiR-15 suppressed the progression of bladder cancer by targeting BMI1 oncogene via PI3K/AKT signaling pathway. Eur Rev Med Pharmacol Sci. 2019; 23:8813–22. https://doi.org/10.26355/eurrev_201910_19276. [PubMed].

38. Ding ZS, He YH, Deng YS, Peng PX, Wang JF, Chen X, Zhao PY, Zhou XF. MicroRNA-34a inhibits bladder cancer cell migration and invasion, and upregulates PTEN expression. Oncol Lett. 2019; 18:5549–54. https://doi.org/10.3892/ol.2019.10877. [PubMed].

39. Chen D, Xie S, Wu Y, Cui Y, Cai Y, Lan L, Yang H, Chen J, Chen W. Reduction of Bladder Cancer Chemosensitivity Induced by the Effect of HOXA-AS3 as a ceRNA for miR-455-5p That Upregulates Notch1. Front Oncol. 2020; 10:572672. https://doi.org/10.3389/fonc.2020.572672. [PubMed].

40. Kato Y, Zembutsu H, Takata R, Miya F, Tsunoda T, Obara W, Fujioka T, Nakamura Y. Predicting response of bladder cancers to gemcitabine and carboplatin neoadjuvant chemotherapy through genome-wide gene expression profiling. Exp Ther Med. 2011; 2:47–56. https://doi.org/10.3892/etm.2010.166. [PubMed].

41. Takata R, Katagiri T, Kanehira M, Tsunoda T, Shuin T, Miki T, Namiki M, Kohri K, Matsushita Y, Fujioka T, Nakamura Y. Predicting response to methotrexate, vinblastine, doxorubicin, and cisplatin neoadjuvant chemotherapy for bladder cancers through genome-wide gene expression profiling. Clin Cancer Res. 2005; 11:2625–36. https://doi.org/10.1158/1078-0432.CCR-04-1988. [PubMed].

42. bladder.pdf. Available from https://www.nccn.org/professionals/physician_gls/pdf/bladder.pdf.

43. McConkey DJ, Choi W, Shen Y, Lee IL, Porten S, Matin SF, Kamat AM, Corn P, Millikan RE, Dinney C, Czerniak B, Siefker-Radtke AO. A Prognostic Gene Expression Signature in the Molecular Classification of Chemotherapy-naïve Urothelial Cancer is Predictive of Clinical Outcomes from Neoadjuvant Chemotherapy: A Phase 2 Trial of Dose-dense Methotrexate, Vinblastine, Doxorubicin, and Cisplatin with Bevacizumab in Urothelial Cancer. Eur Urol. 2016; 69:855–62. https://doi.org/10.1016/j.eururo.2015.08.034. [PubMed].

44. Bochner BH, Hansel DE, Efstathiou JA, Konety B, Lee CT, McKiernan JM, Plimack ER, Reuter VE, Sridhar S, Vikram R. Urinary bladder. AJCC Cancer Staging Manual. 8th edn. Springer Cham: Chicago. 2017; 757–65.

45. Dobin A, Gingeras TR. Mapping RNA-seq Reads with STAR. Curr Protoc Bioinformatics. 2015; 51:11.14.1–11.14.19. https://doi.org/10.1002/0471250953.bi1114s51. [PubMed].

46. Anders S, Pyl PT, Huber W. HTSeq--a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015; 31:166–69. https://doi.org/10.1093/bioinformatics/btu638. [PubMed].

47. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014; 15:550. https://doi.org/10.1186/s13059-014-0550-8. [PubMed].

48. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005; 102:15545–50. https://doi.org/10.1073/pnas.0506580102. [PubMed].

49. Mootha VK, Lindgren CM, Eriksson KF, Subramanian A, Sihag S, Lehar J, Puigserver P, Carlsson E, Ridderstråle M, Laurila E, Houstis N, Daly MJ, Patterson N, et al. PGC-1alpha-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat Genet. 2003; 34:267–73. https://doi.org/10.1038/ng1180. [PubMed].

50. Perkins JR, Dawes JM, McMahon SB, Bennett DL, Orengo C, Kohl M. ReadqPCR and NormqPCR: R packages for the reading, quality checking and normalisation of RT-qPCR quantification cycle (Cq) data. BMC Genomics. 2012; 13:296. https://doi.org/10.1186/1471-2164-13-296. [PubMed].

51. Livak KJ, Schmittgen TD. Analysis of Relative Gene Expression Data Using Real-Time Quantitative PCR and the 2−ΔΔCT Method. Methods. 2001; 25:402–8. https://doi.org/10.1006/meth.2001.1262. [PubMed].

52. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015; 43:e47. https://doi.org/10.1093/nar/gkv007. [PubMed].

53. Witten DM, Tibshirani RJ. Extensions of sparse canonical correlation analysis with applications to genomic data. Stat Appl Genet Mol Biol. 2009; 8:Article28. https://doi.org/10.2202/1544-6115.1470. [PubMed].

54. Robin X, Turck N, Hainard A, Tiberti N, Lisacek F, Sanchez JC, Müller M. pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics. 2011; 12:77. https://doi.org/10.1186/1471-2105-12-77. [PubMed].

55. R Core Team. R: A language and environment for statistical computing. Vienna, Austria; 2013. Available from http://r.meteo.uni.wroc.pl/web/packages/dplR/vignettes/intro-dplR.pdf.

56. Wickham H, Averick M, Bryan J, Chang W, McGowan L, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen T, Miller E, et al. Welcome to the tidyverse. J Open Source Softw. 2019; 4:1686. https://doi.org/10.21105/joss.01686.

57. BLCAsubtyping: Transcriptomic tools to classify bladder tumours according to six published molecular classifications: Baylor, UNC, MDA, Lund, CIT-Curie, TCGA. Github; [cited 2022 Mar 22]. Available from https://github.com/cit-bioinfo/BLCAsubtyping.

Creative Commons License All site content, except where otherwise noted, is licensed under a Creative Commons Attribution 4.0 License.
PII: 28302