A five-gene expression signature to predict progression in T1G3 bladder cancer

European Journal of Cancer 64 (2016) 127e136

Available online at www.sciencedirect.com

ScienceDirect journal homepage: www.ejcancer.com

Original Research

A five-gene expression signature to predict progression in T1G3 bladder cancer Antoine G. van der Heijden a,1, Lourdes Mengual b,*,1, Juan J. Lozano c, Mercedes Ingelmo-Torres b, Maria J. Ribal b, Pedro L. Ferna´ndez d, Egbert Oosterwijk a, Jack A. Schalken a, Antonio Alcaraz b, J. Alfred Witjes a a

Department of Urology, Radboud University Medical Center, Nijmegen, The Netherlands Laboratory and Department of Urology, Hospital Clı´nic, Institut d’Investigacions Biome`diques August Pi i Sunyer (IDIBAPS), Universitat de Barcelona, Spain c CIBERehd, Plataforma de Bioinforma´tica, Centro de Investigacio´n Biome´dica en red de Enfermedades Hepa´ticas y Digestivas, Barcelona, Spain d Pathology Department, Hospital Clı´nic, Institut d’Investigacions Biome`diques August Pi i Sunyer (IDIBAPS), Universitat de Barcelona, Spain b

Received 2 May 2016; received in revised form 31 May 2016; accepted 10 June 2016

KEYWORDS (MESH) Biomarkers; Bladder cancer; Gene expression signature; Prognosis; Progressive disease

Abstract Objective: The aim of this study was to analyze tumour gene expression profiles of progressive and non-progressive T1G3 bladder cancer (BC) patients to develop a gene expression signature to predict tumour progression. Methods: Retrospective, multicenter study of 96 T1G3 BC patients without carcinoma in situ (CIS) who underwent a transurethral resection. Formalin-fixed paraffin-embedded tissue samples were collected. Global gene expression patterns were analyzed in 21 selected samples from progressive and non-progressive T1G3 BC patients using Illumina microarrays. Expression levels of 94 genes selected based on microarray data and based on literature were studied by quantitative polymerase chain reaction (qPCR) in an independent series of 75 progressive and non-progressive T1G3 BC patients. Univariate logistic regression was used to identify individual predictors. A variable selection method was used to develop a multiplex biomarker model. Discrimination of the model was measured by area under the receiver-operating

* Corresponding author: Laboratory of Urology, Hospital Clı´nic de Barcelona, Centre de Recerca Biome`dica CELLEX, Office B22, C/Casanova, 143, 08036 Barcelona, Spain. Tel.: þ34 93 227 54 00x4820. E-mail addresses: [email protected] (A.G. van der Heijden), [email protected] (L. Mengual), juanjo.lozano@ Ciberehd.org (J.J. Lozano), [email protected] (M. Ingelmo-Torres), [email protected] (M.J. Ribal), [email protected] (P.L. Ferna´ndez), [email protected] (E. Oosterwijk), [email protected] (J.A. Schalken), [email protected] (A. Alcaraz), [email protected] (J.A. Witjes). 1 Both authors contributed equally to this work. http://dx.doi.org/10.1016/j.ejca.2016.06.003 0959-8049/ª 2016 Elsevier Ltd. All rights reserved.

128

A.G. van der Heijden et al. / European Journal of Cancer 64 (2016) 127e136

characteristic curve. Interaction networks between the genes of the model were built by GeneMANIA Cytoscape plugin. Results: A total of 1294 genes were found differentially expressed between progressive and non-progressive patients. Differential expression of 15 genes was validated by qPCR in an additional set of samples. A five-gene expression signature (ANXA10, DAB2, HYAL2, SPOCD1, and MAP4K1) discriminated progressive from non-progressive T1G3 BC patients with a sensitivity of 79% and a specificity of 86% (AUC Z 0.83). Direct interactions between the five genes of the model were not found. Conclusions: Progressive and non-progressive T1G3 bladder tumours have shown different gene expression patterns. To identify T1G3 BC patients with a high risk of progression, a five-gene expression signature has been developed. ª 2016 Elsevier Ltd. All rights reserved.

1. Introduction Urothelial carcinoma is the most common bladder cancer (BC) with an incidence of approximately 90%. At diagnosis, two third of urothelial carcinomas are non-muscle invasive and one third are muscle invasive. Non-muscle invasive bladder cancer (NMIBC) is characterised by a high risk of recurrence (30e85% depending on stage and grade) after transurethral resection of the bladder tumour (TURBT). Moreover, up to 17% of all NMIBC will eventually progress to muscle-invasive bladder cancer (MIBC). Progression is mainly seen within 48 months after initial diagnosis and depends on stage, grade, and presence of concomitant carcinoma in situ (CIS) [1]. Progressive patients deserve careful attention, particularly because the treatment and pathogenesis of MIBC and NMIBC differ and conventional histopathological evaluation is inadequate to accurately predict the behaviour of high-risk NMIBC. Thus, there is a clear need for predictive biomarkers that can distinguish progressive from non-progressive NMIBC. Given the heterogeneity of BC, it is unlikely that a single marker can identify its progression. This underlines the importance of biomarker combinations. Gene expression analyses have had an important effect on the discovery of new molecular markers for diagnosis and prediction of disease outcome in BC [2e4]. In this work, differentially expressed genes between progressive and non-progressive T1G3 BC were studied. Subsequently, the recently introduced approach of Biomark Fluidigm Arrays based on the quantitative polymerase chain reaction (qPCR) allowed validation of identified differentially expressed genes in a larger series. A fivegene signature was identified to predict progression in T1G3 BC patients. 2. Materials and methods 2.1. Patients and samples A retrospective multicenter study in which a total of 96 patients (mean age 69 years, range 39e90 years; 24

females, 72 males) with T1G3 urothelial carcinoma of the bladder who underwent TURBT in two different centers (Hospital Clinic, Barcelona, Spain and Radboud University Medical Center, Nijmegen, The Netherlands) between 1993 and 2010 were included. Patients with concomitant CIS were excluded. Nearly all patients were at least treated with an induction cycle of BCG. Two groups of patients were formed; non-progressive patients consisting of recurrent non-progressive patients with at least two years of follow-up with one or more non-muscle invasive recurrences (N Z 56) and progressive patients consisting of individuals with a period of at least more than six months between the initial TURBT showing T1G3 BC and the progression to MIBC (N Z 40). Cytology and urethrocystoscopy 3 months after the initial TURBT had to be negative; otherwise, these patients were qualified as having residual and not progressive disease. The median time to progression was 1.0 (0.5e6.3) year. The median followup in non-progressive patients was 4.8 (2.0e10.9) years. Demographic and clinicopathological characteristics of enrolled patients are shown in Table 1 [5,6]. The clinical classifiers for progression, multiplicity and tumour size were analyzed in both progressive and nonprogressive patients. Multiplicity was found in 68% and 75% in recurrent and progressive patients, respectively. Tumour diameter of 3 cm was found in 9 of 56 (16%) recurrent and 6 of 40 (15%) progressive patients. Tissue samples were obtained under institutional review board-approved protocol. There was no central pathological review performed. 2.2. Tissue specimens and RNA isolation Upon obtainment, the tissue was fixed in 10% formalin within 24 h and subsequently embedded in paraffin. A slide of each specimen was stained with haematoxylineeosin to determine the presence of tumour. The tumour area was macro-dissected from slides (total thickness 80 mm) and total RNA was isolated using the RecoverAll Total Nucleic Acid Isolation Kit for FFPE (Ambion, Inc. Austin, TX, USA) according to the

A.G. van der Heijden et al. / European Journal of Cancer 64 (2016) 127e136

129

Table 1 Demographic and clinicopathological characteristics of enrolled patients

Progressive patients Male Female Mean age (years) Number of tumours Single Multiple Tumour size <3 3 cm Median time to progression (years) Subtotals Non-progressive patients Male Female Mean age (years) Number of tumours Single Multiple Tumour size <3 3 cm Median time follow up (years) Subtotal Total

Screening phase hospital clinic barcelona

Discovery phase hospital clinic barcelona

Discovery phase radboud UMC Nijmegen

11 1 74 (54e86)

8 4 74 (57e90)

12 4 66 (41e83)

2 10

4 8

4 12

9 3 1.0 (0.5e6.3) 12

11 1 0.9 (0.7e5.0) 12

14 2 0.9 (0.5e2.8) 16

8 1 67 (40e83)

18 8 72 (47e85)

15 6 61 (39e77)

0 9

4 22

14 7

9 0 4.9 (4.6e6.7) 9 21

22 4 4.7 (3.7e6.8) 26 38

16 5 5.9 (2.0e10.9) 21 37

manufacturer’s protocol. Total RNA was quantified by spectrophotometric analysis at 260 nm (NanoDrop Technologies, Wilmington, DE, USA). 2.3. Whole-genome gene expression microarray: global screening phase 2.3.1. Global mRNA profiling A flowchart of the entire study is shown in Fig. 1. Global mRNA profiling of 21 selected T1G3 BC cases from Hospital Clinic, Barcelona, Spain (12 progressive and 9 non-progressive), was performed using WholeGenome Gene Expression DASL HT Assay (Illumina, San Diego, CA, USA), according to manufacturer’s instructions [7]. All 21 selected samples had cycle quantification (Cq) values for RPL13A < 28, which is considered to be of acceptable RNA quality by the microarray manufacturer’s (data not shown). 2.3.2. Data analysis DASL gene expression data were processed employing quantile normalization using the Lumi bioconductor package. Next, a conservative probe-filtering step excluding those probes with a coefficient of variation of lower than 0.1 was conducted, which resulted in the selection of a total of 22,032 probes (corresponding to 16,653 genes) from the original set of 29,377. For the detection of differentially expressed probes, a linear model was fitted to the data and empirical Bayes moderated statistics were calculated using the Limma package

from Bioconductor [8]. Genes representing a change of 1.5-fold or greater and moderated p-value <0.05 were considered as differentially expressed. Further analysis was performed on the list of differentially expressed genes using Ingenuity Pathway Analysis (IPA) to identify networks and canonical pathways in which the genes were enriched. 2.4. Quantitative PCR using BioMark 96.96 Dynamic Arrays: classifier discovery phase 2.4.1. Validation of microarray expression data by quantitative PCR in an additional set of samples Differential expression of 94 genes was validated in an additional set of 75 tissue samples from two university medical centres (Table 1) using BioMark 96.96 Dynamic Arrays (Fluidigm, South San Francisco, CA, USA); 75 of 94 genes were selected from microarray expression experiments (p < 0.05). Additionally, 19 genes recently described in progressive NMIBC were selected from literature [9e11]. 2.4.2. cDNA synthesis and pre-amplification Complementary DNA (cDNA) was synthesised from 100 ng of total RNA, using the high-capacity cDNA reverse transcription kit (Applied Biosystems, Foster City, CA USA; hereafter referred to as AB) following manufacturer’s instructions, except that the final volume of the reaction was 25 ml. Each cDNA sample was used for the multiplex pre-amplification of the target cDNAs

130

A.G. van der Heijden et al. / European Journal of Cancer 64 (2016) 127e136

Global screening phase Whole genome gene expression microarray 12 prog; 9 non-prog. LIMMA 1294 genes differen ally expressed 75 genes from microarray (p<0.05) 19 genes from literature

Valida on independent set Quan ta ve PCR 28 prog; 47 non-prog. Rela ve quan fica on Differen al expression: 9 genes from microarray 6 genes from literature 15 genes differen ally expressed

Classifier discovery phase Variable selec on method

5-genes expression signature AUC=0.83

Fig. 1. Study outline. Tissue samples were obtained from a total of 40 progressive and 56 non-progressive high-risk NMIBC patients. Samples were split in a screening (21 samples) and validation phase (75 samples). Genes differentially expressed between progressive and non-progressive high-risk NMIBC patients were identified in the screening phase by using gene expression microarray. Genetic classifiers for NMIBC diagnosis and prediction of tumour aggressiveness were identified, which resulted in a five-gene expression signature. None of the samples from the screening set were employed for the validation process. In addition, an internal crossvalidation strategy was performed on the proposed signature. Abbreviations: AUC Z area under the curve; NMIBC Z Non-muscle invasive bladder cancer; prog Z progressive; non-prog Z nonprogressive.

using TaqMan PreAmp Master Mix kit following manufacturer’s instructions (AB), except that the final reaction volume of the reaction was 5 ml. A total of 94 target genes and the 2 endogenous controls (GUSB and PPIA) were used in the pre-amplification reaction. 2.4.3. Quantitative PCR using BioMark 96.96 Dynamic Arrays cDNA and TaqMan Gene Expression Assays 20X (AB) were loaded into the dynamic array following the

manufacturer’s instructions. The real-time quantitative PCR (qPCR) experiments were performed on the BioMark instrument. 2.4.4. qPCR data analysis The real-time qPCR analysis software was used to obtain cycle quantification (Cq) values. Threshold was manually calculated for each gene. Relative expression levels of target genes within a sample were expressed as DCq (DCq Z Cqendogenous control  Cqtarget gene). The mean Cq value of GUSB and PPIA were used as endogenous control. Genes with Cq values >35 were considered as lowly expressed, and their DCq were imputed to minimum DCq value for that gene. Fold-change values were generated from the median expression of the genes from the BioMark 96.96 Dynamic Arrays in the groups compared. Differences in gene expression levels between progressive and non-progressive patients were explored using the Student’s t test. Significance was defined as p < 0.05. R-software was used for all calculations and to construct the heat map. 2.4.5. Statistical analysis The association of each variable with NMIBC progression status was analyzed by univariate logistic regression. Significance was defined as p value < 0.05. Significant transcripts were subjected to a variable selection using elastic net [12] to identify the smallest number of transcripts to separate both groups. Since all the samples were used for the model generation, the performance of the model may be over-optimised. To correct this bias, we further performed a leave-one-out cross-validation (LOOCV). The optimal probability cut-off for the univariate study variables and logistic regression model was computed through a receiveroperating characteristic analysis. R-software was used for all calculations. Geneegene interaction networks for the genes of the model were built by the GeneMANIA Cytoscape 3.0.0 plugin [13]. Physical, co-expression, and pathway geneegene interactions were evaluated. 3. Results 3.1. Global gene expression profile analysis The analysis of gene expression profile of 12 progressive and 9 non-progressive T1G3 BC patients by Illumina microarray’s resulted in the identification of 1343 transcripts (1294 genes) differentially expressed between both groups; 620 down-regulated and 723 up-regulated transcripts in progressive compared with non-progressive cases. Heat map based on the 50 most differentially expressed genes shows a clear distinction between the progressive and the non-progressive groups (Fig. 2). Deregulated genes were evaluated using IPA. Several networks were significantly enriched by genes that had different expression between progressive and non-

A.G. van der Heijden et al. / European Journal of Cancer 64 (2016) 127e136

131

Fig. 2. Differentially expressed genes between progressive and non-progressive high-risk NMIBC. (A) Heat map displaying the 50 most differentially expressed genes between progressive and non-progressive samples. (B) Heat map displaying the 75 genes differentially expressed between progressive and non-progressive samples, to be further analyzed by RT-qPCR. Red pixels correspond to an increased abundance of mRNA in the samples, whereas green pixels indicate decreased mRNA levels. Rows represent genes and columns represent experimental samples.

progressive patients. The top five scoring networks for genes differentially expressed are shown in Supplementary Fig. 1S, being the network associated with Cell Death and Survival, Cancer, Organismal Injury and Abnormalities (score 30) illustrated in Fig. 3a. Significant IPA canonical pathways are depicted in figure Fig. 3b. Cancer was identified as the top disease associated with the gene set and nuclear factor-kappa B (complex) was found to be the top upstream regulator of our gene set. In addition, Cellular Development and Cellular Growth and Proliferation were identified as the two top molecular and cellular functions of our gene set (Supplementary Fig. 1S).

progression in T1G3 BC, we additionally analyzed expression of these genes in 28 progressive and 47 nonprogressive patients from two hospitals (Table 1) using BioMark 96.96 Dynamic Arrays. The relative gene expression values (DCq) for the 94 markers in progressive and non-progressive samples were analyzed by univariate analysis to see whether it is possible to separate both groups (Supplementary Table 1S). Fifteen genes were differentially expressed between progressive and non-progressive T1G3 BC samples (Table 2). Univariate analysis shows the performance of each individual biomarker as predictors of NMIBC progression (Table 3).

3.2. Prediction of NMIBC progression in an additional validation set

3.3. Classifier discovery phase

To test whether the expression profile of the 94 selected genes (Supplementary Table 1S) can be used to predict

To evaluate if a multiplex model could improve performance over single biomarkers, these 15 genes were subjected to a variable selection to identify the minimum

132

A.G. van der Heijden et al. / European Journal of Cancer 64 (2016) 127e136

A.G. van der Heijden et al. / European Journal of Cancer 64 (2016) 127e136

Table 2 Comparison of microarray and qPCR results for the 15 genes found differentially expressed in qPCR. Gene symbol

ANXA10 BNC2 BNIP3Lx C16orf74x DAB2 FGFR1 FLJ36031 HIST1H2AH HYAL2 MAP4K1x PFKFB4x PPARDx SERPINB5x SPOCD1 ZDHHC7

Microarrays

qPCR

Fold change

p

Fold change

p

3.09 2.38 1.27 1.07 2.85 12.07 2.16 2.78 2.65 1.44 1.39 1.07 2.18 8.00 2.40

0.041* 0.004* 0.600 0.868 0.003* 0.001* 0.002* 0.006* 0.003* 0.164 0.337 0.784 0.084 0.000* 0.004*

3.42 2.37 1.79 1.63 1.79 1.62 1.70 1.32 1.80 2.05 1.66 1.62 2.28 2.25 1.27

0.002* 0.006* 0.037* 0.028* 0.001* 0.019* 0.012* 0.035* 0.006* 0.003* 0.026* 0.026* 0.007* 0.003* 0.017*

Abbreviation: qPCR Z quantitative polymerase chain reaction. *p < 0.05. xGene selected from literature. Fold-change values were generated from the mean expression of the mRNAs in the groups compared.

Table 3 Univariate analysis of the individual biomarkers as predictors of NMIBC progression. Gene symbol

OR

p-value

AUC

SP

SN

NPV

PPV

ANXA10 BNC2 BNIP3L C16orf74 DAB2 FGFR1 FLJ36031 HIST1H2AH HYAL2 MAP4K1 PFKFB4 PPARD SERPINB5 SPOCD1 ZDHHC7

1.61 0.53 1.54 1.70 0.40 0.60 0.55 0.46 0.40 0.44 1.66 1.70 1.65 1.85 2.85

0.004 0.013 0.033 0.029 0.003 0.023 0.016 0.033 0.007 0.003 0.039 0.033 0.008 0.004 0.026

0.76 0.73 0.66 0.70 0.74 0.64 0.69 0.66 0.75 0.75 0.64 0.64 0.71 0.73 0.64

0.72 0.75 0.75 0.79 0.71 0.65 0.68 0.67 0.77 0.78 0.59 0.57 0.78 0.63 0.71

0.77 0.61 0.55 0.62 0.74 0.53 0.67 0.65 0.72 0.61 0.68 0.70 0.66 0.72 0.57

0.64 0.56 0.50 0.56 0.63 0.44 0.58 0.53 0.61 0.55 0.53 0.53 0.57 0.57 0.50

0.83 0.79 0.79 0.82 0.81 0.74 0.76 0.77 0.85 0.82 0.73 0.73 0.84 0.77 0.77

Abbreviations: AUC Z area under the curve; NMIBC Z Non-muscle invasive bladder cancer; NPV Z negative predictive value; PPV Z positive predictive value; OR Z odds ratio; SN Z sensitivity; SP Z specificity.

combination to accurately predict NMIBC progression. This analysis resulted in a final selection of a five-gene model that contains ANXA10, DAB2, HYAL2, SPOCD1 and MAP4K1. This model is capable to

133

discriminate progressive from non-progressive T1G3 BC cases with an overall sensitivity (SN) of 79% and a specificity (SP) of 86% (positive predictive value [PPV] Z 90%; negative predictive value [NPP] Z 71%; error rate (ER) Z 19%; area under the curve [AUC] Z 0.83) (Fig. 4). After applying the LOOCV analysis to this five-gene model, SN was 74% for discriminating between both groups of patients with a SP of 79% (PPV Z 85%; NPP Z 65%; ER Z 24%; AUC Z 0.751) (Fig. 4). The generated network by GeneMANIA shows that there are no direct interactions between the five genes of the model although two of them are in the same pathway (Fig. 5). 4. Discussion Prediction of disease outcome for patients diagnosed with high-risk NMIBC is an important clinical challenge. Histopathological and clinical tumour characterization is unable to accurately predict whether high-risk NMIBC will progress to MIBC. The currently used most important clinical classifiers for progression, stage T1, grade 3, and presence of CIS [1] were equal in our cohort of patients since all suffered from T1G3 BC without concomitant CIS. In our cohort, the less distinctive clinical classifiers, multiplicity and tumour size, did not distinguish progressive from nonprogressive patients. Since treatment regimens change when patients are likely to progress to MIBC, there is a critical need for objective methods to identify those patients. In this study, a set of genes related with progression of TIG3 BC was identified by microarray analysis. IPA analysis corroborates that the differentially expressed genes identified participate in canonical pathways related with carcinogenetic process. Validation of key genes in an independent series allowed developing a 5gene signature (ANXA10, DAB2, HYAL2, SPOCD1, and MAP4K1) to predict progression in T1G3 BC. Although LOOCV indicates a certain degree of overfitting, all data obtained after cross validation corroborate the SN and SP of the final model. The regulation of cancer progression by different mechanisms is shown by the fact that the genes of the model have no direct interactions and derive from different pathways. Other authors also investigated molecular

Fig. 3. Ingenuity pathway analysis (IPA). (a) Network associated with Cell Death and Survival, Cancer, Organismal Injury and Abnormalities derived from our set of genes differentially expressed (score 30). The Figure is a graphical representation between molecules identified in our data and their respective network. Green and red nodes indicate upregulatedregulated and down-regulated gene expression in progressive BC patients, respectively. Intensity of color is related to higher fold change. Nodes uncolored or gray indicate genes not differentially expressed in the data set but integrated into the computationally generated networks based on stored information in IPA knowledge memory, suggesting its relevance to the network. (b) Top five canonical pathways associated to our data set. The significance of the association between our differentially expressed genes and the canonical pathway was measured by using the ratio of the number of molecules from the data set that mapped to the pathway divided by the total number of molecules that map to the canonical pathway. Fisher’s exact test was used to calculate p values.

134

A.G. van der Heijden et al. / European Journal of Cancer 64 (2016) 127e136

Fig. 4. Receiver-operating characteristic curve analysis based on the predicted probabilities derived from the five-gene model (a) and from the LOOCV (b).

Fig. 5. Geneegene interaction network for the five genes of the prognostic model.

A.G. van der Heijden et al. / European Journal of Cancer 64 (2016) 127e136

markers for a better assessment of NMIBC prognosis [14e18]. Notably, Dyrskjot et al. constructed an 88gene classifier for Ta/T1 BC using microarray data to predict outcome [19]. This classifier achieved 66% sensitivity and specificity to predict progression in an international validation study [3]. Going more into detail, ANXA10, DAB2, HYAL2, and MAP4K1 had been previously associated with BC, but to the best of our knowledge, this is the first study that relates alterations in expression of SPOCD1 with BC. SPOCD1 encodes a protein that belongs to the Transcription factor S-II (TFIIS) family of transcription factors and the eightfold down-regulation that was observed may lead to a different transcription profile. ANXA10 participates in cell differentiation and tissue growth. ANXA10 was threefold down-regulated in progressive T1G3 BC patients, what is in concordance with previous reported data in other malignancies [20e22]. Munksgaard et al. showed that low ANXA10 protein expression was associated with concomitant CIS and a shorter progression-free survival [23]. Results about DAB2 expression levels in cancer appear to be conflicting as up- and down-regulation of DAB2 has been reported [24e27]. However, Karam et al. showed that DAB2 expression decreased in comparison to normal bladder and that decreased DAB2 expression was more prominent in advanced BC [28]. In our study, DAB2 was threefold up-regulated in progressive patients. Clearly, our patient population differs as we specifically studied T1G3 BC tumours and excluded patients with concomitant CIS, which may explain this seeming inconsistency. The exact role in tumour promotion or suppression of HYAL2 is not clear yet, but it plays a role in tumour cell adhesion and migration. In progressive T1G3 BC patients, HYAL2 was threefold up-regulated. This is in concordance with recently published data on other cancers [29,30]. MAP4K1 regulates cell cycle, cell adhesion, migration, and apoptosis. MAP4K1 was 1.5-fold up-regulated in our cohort of progressive patients. Wang and co-workers showed in BC cell lines also an up-regulated expression of MAP4K1 [31]. Given the genomic heterogeneity of BC, one limitation of this study is its sample size. Furthermore, our analysis may have underestimated the genetic variation present within individual tumours as we sampled only one segment of the tumour. In addition, although studies in two independent cohorts of patients showed comparable results, a final validation in an independent larger series is needed for clinical implementation. Another limitation is the retrospective non-randomised character of this study. However, prospective recruitment of patients is very difficult because the percentage of progressive T1G3 BC is low and the follow-up needs to be relatively long to observe progression; in our patient population, the longest time interval to progress to MIBC was 6.3 years.

135

An important benefit of this study is the fact that the methodology used is widely available. qPCR is reasonably simple to perform and is inexpensive, and thus, this five-gene model should be easily implemented in clinical practice. However, for molecular markers to be clinically useful, it is necessary to show that adding this molecular marker to an existing model based on the most important clinicopathological factors substantially improves the discrimination of the model. In this analysis, the most important clinical prognosticators (stage T1, grade 3 and concomitant CIS) according to the EORTC risk tables [1] were the same for progressive and non-progressive T1G3 BC patients. Thus, our signature helps selecting which T1G3 BC patients with the same clinicopathological risk factors have a higher chance of progression. In these patients, early cystectomy should be considered. Nevertheless, introducing this prognostic test for T1G3 BC patients in clinical practice requires further validation in a larger cohort of patients.

Acknowledgements The authors would like to thank the technical support from the staff of the Genome Analysis Platform at CIC bioGUNE and the Functional Genomic Platform at IDIBAPS. Appendix A. Supplementary data Supplementary data related to this article can be found at http://dx.doi.org/10.1016/j.ejca.2016.06.003. Grant information This work was supported by grants from the Dutch Cancer Society. Role of the funding source The study sponsor had no role in study design, in the collection, analysis, or interpretation of data neither in the writing of the report. Conflict of interest statement Antoine G. van der Heijden and Lourdes Mengual are inventors of patent EP15180037.2e1403 applied by Stichting Katholieke Universiteit on 12th January 2016. All other authors declared no conflicts of interest.

References [1] Sylvester RJ, van der Meijden AP, Oosterlinck W, Witjes JA, Bouffioux C, Denis L, et al. Predicting recurrence and progression in individual patients with stage Ta T1 bladder cancer using

136

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12] [13]

[14]

[15] [16]

A.G. van der Heijden et al. / European Journal of Cancer 64 (2016) 127e136 EORTC risk tables: a combined analysis of 2596 patients from seven EORTC trials. Eur Urol 2006;49(3):466e75. Modlich O, Prisack HB, Pitschke G, Ramp U, Ackermann R, Bojar H, et al. Identifying superficial, muscle-invasive, and metastasizing transitional cell carcinoma of the bladder: use of cDNA array analysis of gene expression profiles. Clin Cancer Res 2004;10(10):3410e21. Dyrskjot L, Zieger K, Real FX, Malats N, Carrato A, Hurst C, et al. Gene expression signatures predict outcome in non-muscleinvasive bladder carcinoma: a multicenter validation study. Clin Cancer Res 2007;13(12):3545e51. Mengual L, Burset M, Ars E, Lozano JJ, Villavicencio H, Ribal MJ, et al. DNA microarray expression profiling of bladder cancer allows identification of noninvasive diagnostic markers. J Urol 2009;182(2):741e8. Lopez-Beltran A, Sauter G, Gasser T, Hartmann A, SchmitzDra¨ger BJ, Helpap B, et al. Tumours of the urinary system. In: Eble JN, Sauter G, Epstein JI, Sesterhenn IA, editors. Pathology and genetics of tumours of the urinary system and male genital organs. World Health Organization Classification of Tumours. Lyon: IARC Press; 2004. p. 89e157. Sobin LH, Wittekind CH. TNM classification of malignant tumours. International Union Against Cancer. 6th ed. New York: Jonh Wiley & Sons; 2002. Fan JB, Gunderson KL, Bibikova M, Yeakley JM, Chen J, Wickham GE, et al. Illumina universal bead arrays. Methods Enzymol 2006;410:57e73. Smyth GK. Limma: linear models for microarray data. In: Bioinformatics and computational biology solutions using R and bioconductor. Springer; 2005. p. 397e420. Dyrskjot L, Zieger K, Kruhoffer M, Thykjaer T, Jensen JL, Primdahl H, et al. A molecular signature in superficial bladder carcinoma predicts clinical outcome. Clin Cancer Res 2005; 11(11):4029e36. Yun SJ, Jo SW, Ha YS, Lee OJ, Kim WT, Kim YJ, et al. PFKFB4 as a prognostic marker in non-muscle-invasive bladder cancer. Urol Oncol 2012;30(6):893e9. Kim WT, Yun SJ, Park C, Kim IY, Moon SK, Kwon TG, et al. Identification of C16orf74 as a marker of progression in primary non-muscle invasive bladder cancer. PLoS One 2010;5(12):e15260. Zou H, Hastie T. Regularization and variable selection via the elastic net. J R Statist Soc 2005;67:199. Montojo J, Zuberi K, Rodriguez H, Bader GD, Morris Q. GeneMANIA: fast gene network construction and function prediction for cytoscape. F1000Res 2014;3:153. Stein JP, Grossfeld GD, Ginsberg DA, Esrig D, Freeman JA, Figueroa AJ, et al. Prognostic markers in bladder cancer: a contemporary review of the literature. J Urol 1998;160(3 Pt 1): 645e59. Knowles MA. Molecular subtypes of bladder cancer: Jekyll and Hyde or chalk and cheese? Carcinogenesis 2006;27(3):361e73. van Rhijn BW, Vis AN, van der Kwast TH, Kirkels WJ, Radvanyi F, Ooms EC, et al. Molecular grading of urothelial cell carcinoma with fibroblast growth factor receptor 3 and MIB-1 is

[17]

[18] [19]

[20]

[21]

[22]

[23]

[24]

[25]

[26]

[27]

[28]

[29]

[30]

[31]

superior to pathologic grade for the prediction of clinical outcome. J Clin Oncol 2003;21(10):1912e21. Yates DR, Rehman I, Abbod MF, Meuth M, Cross SS, Linkens DA, et al. Promoter hypermethylation identifies progression risk in bladder cancer. Clin Cancer Res 2007;13(7): 2046e53. Karam JA, Shariat SF, Hsieh JT, Knowles MA. Genomics: a preview of genomic medicine. BJU Int 2008;102(9 Pt B):1221e7. Dyrskjot L, Thykjaer T, Kruhoffer M, Jensen JL, Marcussen N, Hamilton-Dutoit S, et al. Identifying distinct classes of bladder carcinoma using microarrays. Nat Genet 2003;33(1):90e6. Patsos G, Germann A, Gebert J, Dihlmann S. Restoration of absent in melanoma 2 (AIM2) induces G2/M cell cycle arrest and promotes invasion of colorectal cancer cells. Int J Cancer 2010; 126(8):1838e49. Matsubara D, Niki T, Ishikawa S, Goto A, Ohara E, Yokomizo T, et al. Differential expression of S100A2 and S100A4 in lung adenocarcinomas: clinicopathological significance, relationship to p53 and identification of their target genes. Cancer Sci 2005;96(12):844e57. Kim J, Kim MA, Jee CD, Jung EJ, Kim WH. Reduced expression and homozygous deletion of annexin A10 in gastric carcinoma. Int J Cancer 2009;125(8):1842e50. Munksgaard PP, Mansilla F, Brems Eskildsen AS, Fristrup N, Birkenkamp-Demtroder K, Ulhoi BP, et al. Low ANXA10 expression is associated with disease aggressiveness in bladder cancer. Br J Cancer 2011;105(9):1379e87. Zhang Z, Chen Y, Tang J, Xie X. Frequent loss expression of dab2 and promotor hypermethylation in human cancers: a metaanalysis and systematic review. Pak J Med Sci 2014;30(2):432e7. Xu HT, Yang LH, Li QC, Liu SL, Liu D, Xie XM, et al. Disabled-2 and Axin are concurrently colocalized and underexpressed in lung cancers. Hum Pathol 2011;42(10):1491e8. Anupam K, Tusharkant C, Gupta SD, Ranju R. Loss of disabled2 expression is an early event in esophageal squamous tumorigenesis. World J Gastroenterol 2006;12(37):6041e5. Zhou J, Hernandez G, Tu SW, Scholes J, Chen H, Tseng CP, et al. Synergistic induction of DOC-2/DAB2 gene expression in transitional cell carcinoma in the presence of GATA6 and histone deacetylase inhibitor. Cancer Res 2005;65(14):6089e96. Karam JA, Shariat SF, Huang HY, Pong RC, Ashfaq R, Shapiro E, et al. Decreased DOC-2/DAB2 expression in urothelial carcinoma of the bladder. Clin Cancer Res 2007;13(15 Pt 1): 4400e6. Bouga H, Tsouros I, Bounias D, Kyriakopoulou D, Stavropoulos MS, Papageorgakopoulou N, et al. Involvement of hyaluronidases in colorectal cancer. BMC Cancer 2010;10:499. Udabage L, Brownlee GR, Nilsson SK, Brown TJ. The overexpression of HAS2, Hyal-2 and CD44 is implicated in the invasiveness of breast cancer. Exp Cell Res 2005;310(1):205e17. Wang Y, Luo H, Li Y, Chen T, Wu S, Yang L. hsa-miR-96 upregulates MAP4K1 and IRS1 and may function as a promising diagnostic marker in human bladder urothelial carcinomas. Mol Med Rep 2012;5(1):260e5.