Abstract
Cristina Segovia1,2,3,*, Mónica Martínez-Fernández1,2,3,*, Marta Dueñas1,2,3,*, Carolina Rubio2,3, Fernando F. López-Calderón2, Clotilde Costa1,4, Cristina Saiz-Ladera1,5, María Fernández-Grajera1, José Duarte2,3, Huberto García Muñoz6, Federico de la Rosa2,3, Felipe Villacampa2,3, Daniel Castellano2,3, Jesús M. Paramio1,2,3
1Unidad de Oncología Molecular, CIEMAT (ed70A), Madrid, Spain
2Grupo de Oncología celular y Molecular, Hospital Universitario 12 de Octubre, Madrid, Spain
3Centro de Investigación Biomédica en Red de Cáncer (CIBER ONC), Spain
4Unidad Mixta Roche-CHUS, Hospital Universitario de Santiago de Compostela, Travesía de Choupana, s/n, Santiago de Compostela, A Coruña, Spain
5Departamento de Bioquímica y Biología Molecular I, Facultad de Biología, Universidad Complutense, Madrid, Spain
6Servicio de Anatomía Patológica, Hospital Universitario 12 de Octubre, Instituto de Investigación 12 de Octubre i+12, UCM, Av Cordoba s/n, Madrid, Spain
*These authors have contributed equally to this work
Correspondence to:
Jesus M. Paramio, email: [email protected]
Keywords: bladder cancer, PIK3CA, genomics, Ezh2, recurrence progression
Received: July 20, 2016 Accepted: December 12, 2016 Published: January 02, 2017
ABSTRACT
The high rates of tumor recurrence and progression represent a major clinical problem in non-muscle invasive bladder cancer. Previous data showed that EZH2-dependent signaling mediates these processes, whereas the frequent alterations of PIK3CA gene (copy gains and mutations) are predictive of reduced recurrence. Here we show, using clinical samples and bladder cancer cell lines, a functional interaction between EZH2- and PIK3CA-dependent signaling pathways. PIK3CA alterations mediated, on the one hand, the increased expression of two miRNAs, miR-101 and miR-138, which posttranscriptionally downregulate EZH2 expression. On the other hand, PIK3CA alterations facilitate the activation of Akt which phosphorylates EZH2 on Ser21, precluding the trimethylation of histone H3 in K27. Remarkably the increased expression of miR101 or miR138 and the expression of Ser21-phosphorylated EZH2 are good prognostic factors regarding non-muscle invasive bladder cancer recurrence and progression. Collectively, this study provides molecular evidences indicating that the gene expression rewiring occurring in primary bladder tumors, associated with increased EZH2 expression and activity and mediating the increased recurrence and progression risk, are prevented by PIK3CA-dependent signaling. This molecular process may have deep implications in the management of bladder cancer patients and in the design of novel molecularly targeted therapeutic approaches.
INTRODUCTION
Bladder cancer (BC) is a complex disease caused by both genetic and environmental factors [1]. At diagnosis, two major types of BC: non-muscle invasive (NMIBC), and muscle invasive tumors (MIBC) are characterized. This pathological classification also defines the possible therapeutic options. The NMIBCs represent around 70% of BC and are treated by transurethral resections, in some cases followed by intravesical therapy. The MIBC are regularly treated by cystectomy followed by chemotherapy. Although NMIBCs have a favorable prognosis, they also show one of the highest rates of recurrence, which in some cases can progress into muscle-invasive tumors. This makes necessary a regular surveillance by cystoscopy and urine cytology indefinitely (EAU guidelines) [2]. As a consequence, NMIBC represents one of the most costly malignancies to health care systems in developed countries [2].
Genomic studies of BC identified different subtypes, and revealed that MIBC and NMIBC are characterized by distinct molecular profiles, suggestive of possible distinct molecular origin [3]. A large fraction of bladder cancers displays alterations in PI3K pathway components, including PIK3CA gene mutations frequently associated to low grade and stage tumors [4]. Remarkably, different mutations of the PIK3-mTOR pathway occur in NMIBC and MIBC, being the PTEN gene alterations more frequent in MIBC and associated with poor outcome [4]. These differences remain poorly understood, but they might define better therapeutic options, such as the increased sensitivity to mTOR inhibitors displayed by BC patients bearing mTOR-activating mutations [5]. Using a highly sensitive PCR technique, we have previously demonstrated that PIK3CA gene alterations are extremely frequent in NMIBC being present in non-affected bladder tissue, and associate with low recurrence and progression [6]. Conversely, genomic studies in human samples and in mouse models of NMIBC revealed that the recurrence and progression are associated with increased expression of EZH2, which promotes global changes in gene expression [7], including various miRNAs, involved in the repression of the epithelial mesenchymal transition [8] and the aberrant expression of lncRNA HOTAIR [8]. The different clinical evolution of tumors characterized by PIK3CA gene alterations and those exhibiting increased expression and activity of EZH2, led us to hypothesize that these two pathways may exert opposite roles in NMIBC. Here, we report that PI3K-dependent signaling negatively regulates EZH2 expression and activity in NMIBC, thus providing a possible explanation for the observed contrasting roles of these two pathways in this disease.
RESULTS
Transcriptional changes in NMIBC samples bearing PIK3CA alterations
The alterations in PIK3CA gene (mutations and copy gains or amplifications) are very frequent in NMIBC. We have previously shown in a series of 87 BC samples, using high sensitive PCR-based approaches, that the activating mutations (E545K/D, E542K and H1047R) or PIK3CA gene copy gains or amplifications may reach up to 50% of tumor samples, and they are also present in non-tumoral tissue from BC patients, suggesting field cancerization processes [6]. Moreover, they also seemed to associate with reduced recurrence likelihood [6]. We confirmed this association using updated clinical data from the same patients (Supplementary Figure 2A), and we also observed that these alterations are also associated with reduced progression upon recurrence (Supplementary Figure 2B). This suggests that such PIK3CA gene alterations are indicatives of good clinical outcome in NMIBC. To monitor whether similar association also occurs in MIBC, we used TCGA database. We observed that PIK3CA gene alterations were associated with increased overall and disease specific survival (Supplementary Figure 2C and 2D).
In order to gain insights into the possible biological basis of this reduced recurrence and progression, we analyzed our previous global transcriptome study [7] to discern possible gene expression changes that discriminate tumors bearing or not PIK3CA gene alterations. This provided 457 transcripts (306 overexpressed and 151 underexpressed) that characterized tumors bearing PIK3CA gene alterations (Supplementary Table 2).
Gene Ontology (GOBP; Supplementary Table S3 and Supplementary Table S4) and GSEA (Supplementary Table S5 and Supplementary Table S6) revealed that the upregulated genes in tumors bearing PIK3CA gene alterations played a major role in various metabolic and ribosomal processes, whereas the downregulated genes displayed an association with cytoskeleton organization and interphase cell cycle control. In addition, the ChEA analysis revealed a significant primary binding by various polycomb members in upregulated and downregulated genes (Figure 1A and Supplementary Table S7 and Supplementary Table S8). Since, gene downregulation mediated by polycomb repressing complexes appears to play a major role mediating recurrence in NMIBC [7], we decided to explore this observation. We found a highly significant overlap between the downregulated genes characteristic of recurrent tumors with those upregulated in tumors with PIK3CA gene alterations (Figure 1B). Moreover, the upregulated genes also displayed overlap with genes downregulated in tumors that showed progression upon recurrence (Figure 1C), and those downregulated in tumors bearing mutant FGFR3 and wt PIK3CA (Figure 1C), which correspond to the group showing the earliest recurrence in our series [6]. Finally, we observed by ChEA analysis that the overlapping genes between downregulated in recurrent tumors and in those upregulated in tumors bearing PIK3CA alterations displayed significant enrichment in binding by polycomb members, and also in H3K27me3 marks (Figure 1D). Since the downregulated genes in recurrent tumors were predominantly associated with gene repression mediated by EZH2 [7], these observations might indicate that PIK3CA gene alterations would be acting opposite to EZH2-mediated gene repression. Moreover, GSEA analyses of PIK3CA-altered tumors identified various miRNAs, including those targeting EZH2 (Figure 1E), suggesting that PIK3CA alterations may affect the posttranscriptionally regulation of EZH2 expression through altered expression of miRNAs.
Figure 1: PIK3CA alterations oppose to EZH2 in NMIBC. Transcriptome studies of human NMIBC recurrence. A. Summary of ChEA analyses showing enrichment for Polycomb binding motifs in genes identified between tumors with or without PIK3CA gene alterations. Note the enrichment in gene upregulated in tumors bearing PIK3CA alterations, whereas the binding by E2F and Myc are significant in upregulated and downregulated genes. B. Overlap between genes upregulated in tumors with PIK3CA alterations and downregulated in primary tumors showing recurrence. Left panel shows a Heatmap showing the distribution of genes (rows) and samples (columns) following supervised clustering (Pearson correlation and Average Linkage method) of 28 tumor and 10 normal samples according the recurrent and non-recurrent tumors. Right panel shows the corresponding heatmap showing the distribution according the presence or not of PIK3CA gene alterations in tumors. A red (overexpressed) to blue (downregulated) scheme following the above scale limits (in log2 scale) is shown. Lower panel shows a Venn diagram showing the overlap between genes downregulated in tumors that subsequently recurred and upregulated in tumors with PIK3CA gene alterations. Statistic significance was estimated by F Fisher’s exact test. C. Venn diagrams showing the overlap between genes downregulated in tumors that showed progression upon recurrence and upregulated in tumors with PIK3CA gene alterations (right panel); and overlap between downregulated genes in tumors bearing mutant FGFR3 and wt PIK3CA and upregulated in tumors with PIK3CA gene alterations (right panel). Statistic significance was estimated by F Fisher’s exact test. D. Summary of ChEA analyses showing enrichment for Polycomb binding motifs in overlapping genes between those upregulated in tumors with PIK3CA gene alterations and downregulated in primary tumors showing recurrence. E. Summary of GSEA analyses showing enrichment for miRNA motifs in genes identified between tumors with or without PIK3CA gene alterations. NES: normalized enrichment factor. Negative values indicated enrichment in tumors bearing PIK3CA alterations.
mIR-101 and miR-138 levels define recurrence and progression in NMIBC
To explore the possible posttranscriptional regulation of EZH2, we monitored whether the expression of EZH2 protein (analyzed by immunohistochemistry) or mRNA (analyzed by RTqPCR) are equally predictive of clinical outcome. We found that tumors showing positive EZH2 staining dictated early recurrence (Figure 2A, 2B) and progression (Figure 2D), in agreement with our previous findings [7]. However, the increased EZH2 gene expression levels (Figure 2C) was not predictive of early recurrence in NMIBC (Figure 2C), reinforcing a possible post-transcriptional modulation of EZH2. Of note, the positive EZH2 staining was not discriminative of tumor grade, stage, number of implants, tumor size or previous smoking history of the patients (Supplementary Figure 3).
Figure 2: PIK3CA alterations oppose to EZH2 in NMIBC through miRNA expression. Transcriptome studies of human NMIBC recurrence. A. Representative immunohistochemistry images of EZH2 staining showing positive (left panel) and negative (right panel). Bar=150μm B. Kaplan-Meier analysis showing that NMIBC patients with positive EZH2 expression (assessed by immunohistochemistry) in the primary tumor showed an earlier recurrence (P value was obtained by the log-rank test). C. Kaplan-Meier analysis showing that EZH2 gene expression (mRNA values distributed according the median) in the primary tumor does not discriminate an earlier recurrence (P value was obtained by the log-rank test). D. Kaplan-Meier analysis showing that NMIBC patients with positive EZH2 expression (assessed by immunohistochemistry) in the primary tumor showed progression in recurrence (P value was obtained by the log-rank test). E. RTqPCR analyses showing the expression of miR-101 (left panel) and miR-138 (right panel) in recurrent and non-recurrent NMIBC samples Statistical significance was obtained by Mann–Whitney’s test. ** p value <0.01, * p value<0.05 E’) RTqPCR analyses showing the expression of miR-101 (left panel) and miR-138 (right panel) in NMIBC samples according the presence of alterations in PIK3CA gene in the primary tumors. Statistical significance was obtained by Mann–Whitney’s test. *** p value <0.005, ** p value <0.01. F, F’. RTqPCR analyses showing the expression of miR-101 (F) and miR-138 (F’) in NMIBC samples showing positive or negative EZH2 staining. Statistical significance was obtained by Mann–Whitney’s test. G. Kaplan-Meier analysis showing that NMIBC patients with low expression of miR-101 and miR-138 positive expression (according the median) in the primary tumor showed an earlier recurrence (p value was obtained by the log-rank test). H. Kaplan-Meier analysis showing that NMIBC patients with low expression of miR-101 and miR-138 positive expression (according the median) in the primary tumor showed progression in recurrence (p value was obtained by the log-rank test).
We next studied possible differences in the expression of the miRNAs identified by the GSEA analysis between recurrent and non-recurrent tumors (Figure 1E). We found that only miR-101 and miR-138 displayed significant downregulation in recurrent respect to non-recurrent NMIBC samples (Figure 2E) and increased expression in tumors bearing PIK3CA gene alterations (Figure 2E’). Moreover, we also observed that the reduced EZH2 staining, indicative of reduced protein levels, associated with increased levels of miR-101 and miR-138 (Figure 2F, 2F’).
Since the increased protein levels of EZH2 predict early recurrence and tumor progression in recurrence (Figure 2B; see also [7]), we sought to determine whether miR-101 and/or miR-138 could also be predictive biomarkers of clinical outcome in NMIBC. We found that, although the individual expression of miR-101 or miR-138 were unable to predict the clinical outcome, tumors displaying increased expression of miR-101 or miR-138 showed reduced recurrence and reduced tumor progression upon recurrence better than those showing reduced expression of both miRNAs (Figure 2G, 2H).
Expression of active Akt associates with Ser21 phosphorylated EZH2 in NMIBC
The EZH2 catalytic activity to generate H3K27me3 chromatin marks is regulated by Akt-mediated phosphorylation [9]. Since PIK3CA gene alterations are associated with increased Akt activity in NMIBC [6], we monitored the expression of Ser 473 phosphorylated Akt (Figure 3A, 3A’), EZH2 (Figure 3B, 3B’) and Ser21-phosphorylated EZH2 (Figure 3C, 3C’) in NMIBC by immunohistochemistry. We observed a significant association of phosphorylated Akt with PIK3CA gene alterations and with reduced EZH2 protein expression (Figure 3D).
Figure 3: EZH2 phosphorylation in Ser21 is associated with PIK3CA gene alterations in NMIBC. A-A’. Representative immunohistochemistry images of Ser473 phosphorylated Akt positive (A), and negative (A’), B, C. Representative immunohistochemistry images of primary tumors showing positive EZH2 (B) and negative Ser21 phosphorylated EZH2 (C) staining Bar=150μm B’, C’) Representative immunohistochemistry images of primary tumors showing positive EZH2 (B’) and positive Ser21 phosphorylated EZH2 (C’) staining Bar=150μm D. Contingency tables showing the significant association between phosphorylated Akt with tumors bearing PIK3CA alteration (upper panel), phosphorylated Akt (Ser473) with tumors showing negative EZH2 staining (mid panel) and phosphorylated Akt with tumors showing positive Ser21 phosphorylated EZH2 staining mutation (lower panel) E. Kaplan-Meier analysis showing that NMIBC patients with positive EZH2 and negative Ser21 phosphorylated EZH2 expression in the primary tumor showed an earlier recurrence (p value was obtained by the log-rank test). F. Kaplan-Meier analysis showing that NMIBC patients with positive EZH2 and negative Ser21 phosphorylated EZH2 expression in the primary tumor showed a tendency to display progression in recurrence (p value was obtained by the log-rank test).
Importantly, we observed a predominant expression of Ser21 phosphorylated EZH2 in phosphorylated-Akt positive tumors (Figure 3D). Since the Ser-21 phosphorylation precludes EZH2 to carry out the H3K27me3 chromatin marks [9], and the downregulated genes that predict early recurrence preferentially displayed this chromatin modification [7], we tested whether Ser21-phosphorylated EZH2 could help to refine the recurrence probability in NMIBC patients. Of note, the positive Ser21-phosphorylated EZH2 staining was not discriminative of tumor grade, stage, number of implants, tumor size or previous smoking history of the patients (Supplementary Figure 4). Moreover, double immunofluorescence staining showed that all tumors positive for Ser-21 phosphorylated Ezh2 were also positive for total Ezh2 (Supplementary Figure 5). The series was divided into two groups: 1) those patients displaying positive EZH2 and negative Ser21-P EZH2 staining, and 2) those negative for EZH2, or those positive for both EZH2 and Ser21-P EZH2 (Supplementary Figure 5). Kaplan Meyer analysis of these two groups revealed that a statistically significant early tumor recurrence was associated with patients showing high EZH2 and negative Ser21-P EZH2 expression (Figure 3E). We also observed that progression in recurrences was only observed in primary tumors positive for EZH2 and negative in our series (Figure 3F). In order to confirm these observations, we analyzed, in a semiquantitive manner, the correlation between Akt activation (by AKT Ser473P staining) and the expression of total EZH2 or the Ser21-P EZH2 staining by the histoscore (HS) approach (see Materials and Methods) for each sample. This analysis revealed a significant negative correlation between AKT Ser473P and total EZH2, and a significant positive correlation between AKT Ser473P and Ser21-P EZH2 (Supplementary Figure 6). Moreover, we also observed that samples corresponding to tumors showing recurrence and those showing progression upon recurrence were predominantly associated with low HS for AKT Ser473P and Ser21-P EZH2 and high HS for total EZH2 (Supplementary Figure 6).
PIK3CA opposes to EZH2 in BC cell lines
To further substantiate the potential interaction between PI3K- and EZH2-dependent signaling, we selected a series of BC cell lines of low malignant potential and with complete molecular and genomic characterization [10]. By immunoblotting, we observed that the presence of PIK3CA mutations or PTEN gene LOH were associated with increased Akt Ser 473 and Thr 308 phosphorylation., We also observed that these cells also showed a trend to display partial reduction of EZH2 and increased Ser 21 phosphorylated EZH2 levels (Figure 4A). We found no relationship with activated ERK1/2. Interestingly, the activated Akt was partially and inversely correlated with the lncRNA HOTAIR expression (Figure 4B). This might be related with our previous observation suggesting that HOTAIR expression is, at least in part, under EZH2 positive regulation in BC [8]. On the other hand, the expression of miR-101 and miR-138 were, in general, increased in those cell lines showing activated Akt (Figure 4C). These results confirmed that the observed cross talk between PI3K- and EZH2-dependent signaling was also present in BC cell lines.
Figure 4: Expression of active PIK3CA mediates Ezh2 degradation and catalytic inhibition BC cells. A. Immunoblot showing the expression of phosphorylated ERK1/2, phosphorylated Akt, total EZH2 and Ser 21 phosphorylated EZH2 in a collection of NMIBC cell lines, of known genomic charactersitics [10]. GAPDH was used to normalize loading. B. RTqPCR analyses showing the expression of the LncRNA HOTAIR in the BC cell lines. Asterisks denote those BC cell lines showing activated Akt according the immunoblot against Ser473 or Thr308 phosphorylated Akt (panel A). C. RTqPCR analyses showing the expression of miR-101 (upper panel) and miR-138 (lower panel) in the BC cell lines. Asterisks denote those BC cell lines showing activated Akt according the immunoblot against Ser473 or Thr308 phosphorylated Akt (panel A). D. Immunoblot showing the expression of p110α, Ser 473 and Thr308 phosphorylated Akt, EZH2, Ser21 phosphorylated EZH2, trimethylated K27 histone H3 and trimethylated K9 histone H3 in RT112 cells transfected with empty vector (control) or with plasmids coding for a myristoylated PIK3CA gene, or with E545K or H1047 PIK3CA gene. ACTIN was used to normalize the loading. E. qPCR analyses showing the expression of miR-101 and miR-138 (lower panel) in the RT112 cells transfected with the quoted plasmids. Statistical significance was obtained by Mann–Whitney’s test.*** p value <0.005, ** p value <0.01, * p value <0.05 F. Luciferase expression (normalized to Renilla values) of the construct containing the 3’UTR wt of EZH2 gene or a 3’UTR mutated in the seed sequences for miR-101 and miR-138 in RT112 cells transfected with the quoted plasmids. Statistical significance was obtained by Mann–Whitney’s test.** p value <0.01, * p value<0.05 G. qPCR analyses showing the expression of the LncRNA HOTAIR in the RT112 cells transfected with the quoted plasmids. Statistical significance was obtained by Mann–Whitney’s test.Test. *** p value <0.005, ** p value <0.01, * p value <0.05
To demonstrate the converse effects between PI3K and EZH2, we expressed E545K, H1047R mutants or a myristoylated (constitutively active) PIK3CA gene constructs in RT112 cells. Their expression resulted in increased Akt activity, increased Ser-21 phosphorylation of EZH2, whilst reduced the total EZH2 and the H3K27me3 chromatin marks (Figure 4D). Moreover, we also observed that the different PIK3CA constructs also induced significant increase in miR-101 and miR-138 levels (Figure 4E) with non-significant differences among the different mutants. To monitor whether this increase in miRNA expression could account for the overall reduction in EZH2 levels, we performed a co-transfection experiment in which the different PIK3CA constructs were co-expressed along with pSiCheck2 plasmid in which the 3’UTR corresponding to EZH2 were placed downstream of the Renilla luciferase translational stop codon (Figure 4F). As a control, we used a similar construct including a 3’UTR of EZH2 in which the seed sequence for miR-101 and mir-138 has been mutated. We observed that the inclusion of the wt 3’UTR of EZH2, but not the mutant, led to a significant reduction of the luciferase activity. Similar results were obtained using the MGHU3 cells (data not shown). Finally, since we have observed that the expression of lncRNA HOTAIR is partially modulated by EZH2 in BC cell lines [8], we also monitored HOTAIR expression in the transfected cells. We observed that the HOTAIR levels were significantly reduced upon expression of E545K or H1047R mutants or a myristoylated PIK3CA gene constructs in RT112 BC cells (Figure 4G), although minor non-significant differences were detected among the different mutants, similar to those observed in the case of miR-101 or miR-138, were detected.
These results supported our observations of clinical samples, providing a possible molecular mechanistic framework to explain the opposite clinical outcome observed in NMIBC samples with respect to PIK3CA activation and EZH2-mediated gene expression reprogramming.
DISCUSSION
The PI3K-Akt axis is a well-recognized oncogenic pathway in multiple tumor types, being considered a putative therapeutic target [11]. However, the significance of such alterations in relation to clinical outcome are not conclusive and may be related to the specific mechanism of activation (PIK3CA mutation or PTEN loss), or specific tumor subtypes. For instance, PIK3CA mutation associates good prognosis in breast cancer patients [12–14], and PIK3CA amplification predicts increased survival in colorectal cancer [15]. On the contrary, PIK3CA mutation favors malignant progression in head and neck squamous cell carcinoma [16], and poor prognosis in ovarian cancer [17], whereas PIK3CA amplification predicts reduced survival in gastric cancer [18]. To date, the relevance of PIK3CA alterations in bladder cancer recurrence remains poorly understood. Importantly, different types of alterations in PI3K-dependent signaling occur in NMIBC and in MIBC with potential different clinical consequences [4].
EZH2 aberrant expression is also frequent in multiple tumors thus becoming an attractive molecular target [19, 20]. The increased expression of EZH2 may be due to transcriptional and posttranscriptional mechanisms. In the first instance, EZH2 transcription is positively modulated by E2F transcription factors [21], and in particular by E2F3 [7, 22]. In this regard, we reported the existence of a pRb-E2F-Ezh2 loop as an essential mechanism in bladder cancer in vivo in transgenic mice, which is also associated with poor clinical outcome in human BC patients [7]. Regarding the posttranscriptional regulation, EZH2 expression is controlled by various miRNAs [23]. In particular, increased miR-101 or miR-138 leading to reduced EZH2 expression have been reported for various types of cancer including BC [24]. However, the possible mechanisms of deregulated miR-101 and miR-138, have not been totally elucidated. In particular, reduced miR-101 expression has been associated with deletion [25, 26], and epigenetic silencing, probably mediated by c-myc or EZH2 increased expression [27]. To our knowledge, our data are the first evidence of positive modulation of these miRNAs by PI3K-dependent signaling.
In addition, EZH2 phosphorylation promotes its activation or repression depending on the specific residues phosphorylated [28]. Interestingly, the specific phosphorylation of EZH2 in Ser21 mediated by Akt suppresses its histone methyltransferase activity [9]. This could induce a change in the specificity of EZH2, which can methylate and activate other targets, such as androgen receptor [29], or Stat3 [30], and can mediate the response to chemotherapy [31]. This change in specificity has important clinical consequences, being associated with poor clinical outcome in prostate cancer [29] and glioblastoma [30], and explaining why in these tumors the increased Akt activity also associates with bad prognosis. Our data indicate that Ser21 phosphorylation of EZH2 is associated with low recurrence in NMIBC. This highlights the relevance of histone modification in the BC pathogenesis previously described through genomic studies [7, 32]. Nonetheless, our data cannot discard the existence of other possible mechanisms in advanced invasive BC.
Our previous data have suggested that Pi3K- and Ezh2-dependent pathways display opposite correlation with clinical outcome in BC [6, 7]. Here we provide a possible molecular explanation of these contrasting activities. which can be supported by our present data providing a possible molecular explanation. Activating mutations in PIK3CA gene, through increased expression of specific miRNAs, lead to reduced EZH2 expression. In addition, increased Akt activity, mediated in part by these PIK3CA gene mutations, produced the phosphorylation of EZH2 in Ser21, which accounted for reduced H3K27me3 marks.
Our present work may have potential clinical relevance. In particular, there is a current need of biomarkers that may predict early recurrence and progression of superficial NMIBC. Our data indicate that the presence of high levels of Ezh2 non phosphorylated in Ser21 (probably together with reduced expression of miR-101 and miR-138 and/or determination of PIK3CA gene status and activation of Akt) could represent excellent biomarkers, allowing a better identification of high risk tumors than the current clinical assessments [33]. However, this would require a large cohort studies and, if possible, quantitative measurements. In this regard, ongoing work will determine whether Ezh2 levels, together with the determination of specific miRNAs, can be used for prognostic determination in liquid biopsies [34]. In addition, our present data, together with the previous findings on the epigenetic changes occurring in BC, may also allow the identification of possible new therapeutic strategies [20].
Future work aimed to determine the possible activation of other transcription factors by EZH2-mediated methylation, and to ascertain whether similar processes occur in MIBC, could help to define better clinical intervention protocols for BC management.
MATERIALS AND METHODS
Patients
Tumor samples and medical records were analyzed from 87 patients who had been consecutively evaluated and treated by transurethral resection at the Urology Department of the University Hospital “12 de Octubre” between January 2009 and October 2011. Informed consent was obtained from all patients and the study was approved by the Ethical Committee for Clinical Research of University Hospital “12 de Octubre” (2014/0362). The sample recollection and preservation procedures have been reported elsewhere [6, 7]. The pathologic and clinical data of patients are provided in Supplementary Table 1. The determination of PIK3CA gene mutations and copy number in primary tumor sample DNAs have been previously described [6]. Progression upon recurrence was considered when the tumor, at any recurrence event, displayed increased stage and/or grade as monitored by the same pathologist.
Tissue microarray (TMA) and immunohistochemistry
The aforementioned cases fixed in formalin (10%) and embedded in paraffin were included in two separate TMAs (1.5-mm core diameter), with at least two representative duplicate cores for each case [35], and constructed with a manual tissue arrayer (Beecher Instruments, Sun Prairie, WI) using a standard method [36]. The TMAs were stained with HE and were reviewed to confirm the presence of representative tumor tissue (at least 70% of tumor cells). Immunohistochemistry was performed essentially as previously described [6]. Antibodies used were anti phospho-Ser473-Akt (Cell Signaling), anti EZH2 (Abnova) and anti-phospho-Ser21-EZH2 (AbCam). Initially, the classification of samples was performed according the percentage of positive cells in each samples (Neg: 0-30%, Mid: 30-60% and High: 60-100%). However, in most cases no significant differences were observed between mid and high samples (Supplementary Figure 1). Accordingly, the mid and high samples were grouped for final determinations. Signal was amplified using avidin-peroxidase (ABC elite kit Vector) and peroxidase was visualized using diaminobenzidine as a substrate (DAB kit Vector). Negative control slides were obtained by replacing primary antibodies with PBS (data not shown). Double-blind scoring of the results and selection of the thresholds, internal controls for reactivity of each antibody, and tissue controls for the series were done according previously published methods [6]. Immunohistochemitry scoring (Histoscore) was calculated by a semiquantitative approach following a previously reported method [37]. Briefly, staining intensity (0, 1+, 2+, or 3+) and the percentage of cells for each intensity were determined in a high power fixed field. The histoscore, ranging from 0 to 300, was then calculated according the following formula:
[1 × (% cells 1+) + 2 × (% cells 2+) + 3 × (% cells 3+)]
RTqPCR
Total RNA was isolated using miRNeasy Mini Kit (Qiagen) according to the manufacturer’s instructions and DNA was eliminated (Rnase-Free Dnase Set Qiagen). Reverse transcription was performed using the Omniscript RT Kit (Qiagen), TBP gene was used as reference gene for normalization using 50ng of total RNA and specific primers (EZH2: 5’-CCTGTCGACATGTTTTGGTC-3’; TBP: 5’-GTGTTTAAAATCTACATA-3’). PCR was performed in a 7500 Fast Real Time PCR System using Go Taq PCR master mix (Promega) and 1 μl of cDNA as a template. Melting curves were performed to verify specificity and absence of primer dimers. Reaction efficiency was calculated for each primer combination. The sequences of the specific oligonucleotides used have been reported elsewhere [7, 8]. To measure quantitatively the expression of miRNAs, RNA was extracted using the same method as for the genes. Reverse transcription was carried out from 10 ng total RNA along with miR-specific primer using the TaqMan® MicroRNA Reverse Transcription Kit (Applied Biosystems). PCR assays were performed using TaqMan® Gene Expression Master Mix and 7500 Fast Real Time PCR System (Applied Biosystems) as reported [38]. For normalization, we used RNU6B.
Microarray data
To determine gene expression in samples bearing PIK3CA alterations, we used a previously reported whole dataset (GSE38264) [7]. The discrimination between tumors with or without PIK3CA gene alterations (copy gains, mutations or both) was obtained by a supervised analysis of differential gene expression using Student’s T test (Adjusted Bonferroni p-val≤0.01 using 1000 random permutations) in the Multiexperiment Viewer 4.5 (MeV 4.5) software [39]. Hierarchical clustering analysis was done using Pearson correlation and complete linkage method. To search for possible overlapping between the obtained upregulated gene signature and those genes downregulated in recurrent NMIBC samples previously reported [7], we used a Fisher's exact test, and was considered statistically significant for Odds Ratio ≥ 2 and p ≤ 0.05. Gene Ontology of Biological processes was performed using the DAVID webtool (https://david.ncifcrf.gov/) [40]. Gene Set Enrichment Analysis (GSEA) was performed using the MSignature and Motif databases [41]. Identification of transcription factor binding was performed using the ChiP enrichement analysis of the Enrich webtool (http://amp.pharm.mssm.edu/Enrichr/) [42, 43]
Cell lines and plasmids
Cells were kindly provided by Dr. FX Real (CNIO, Madrid, Spain) and routinely cultured in DMEM supplemented with 10% FBS and 1% of antibiotic-antimycotic. Transfections were performed as previously described [8, 38]. Only Mycoplasma-free cultures were used. Plasmid coding for myr-p110α, E545K and H1047R PIK3CA mutants (in pBABE backbone) were generously provided by Dr. A. Carnero (IBIS/HUVR/CSIC/Universidad de Sevilla, Seville, Spain). Luciferase assays were performed using pSiCheck2 (Promega) plasmid derivatives after cloning the wt 3′UTR sequence of EZH2 gene (CTCGAGCTGCCTTAGCTTCAGGAACCTCGA GTACTGTGGGCAATTTAGAAAAAGAACATGCAGT TTGAAATTCTGAATTTGCAAAGTACTGTAAGAATA ATTTATAGTAATGAGTTTAAAAATCAACTTTTTATT GCCTTCTCACCAGCTGCAAAGTGTTTTGTACCAGT GAAGCGGCCGC), or its mutated counterpart (CTCGAGCTGCCTTAGCTTCAGGAACCTCGcgcgcggcGGGC AATTTAGAAAAAGAACATGCAGTTTGAAATTCTG AATTTGCAAcgcgcggcAAGAATAATTTATAGTAATGA GTTTAAAAATCAACTTTTTATTGCCTTCTgcgatcgTG CAAAGTGTTTTGTACCAGTGAAGCGGCCGC), into the XhoI-NotI sites of the expression vector downstream of the Renilla luciferase translational stop codon. This plasmid includes a firefly luciferase for normalization. The luciferase (Renilla and Firefly) activity was measured following the manufacturer’s recommendations.
Western blot
Cell pellets were disrupted by freeze-thawing cycles in lysis buffer (200mM4-(2-hydroxyethyl)-1-piperazineethanesulfonic acid pH 7.9, 25% glycerol, 400mM NaCl, 1mM ethylenediaminetetraacetic acid, 1mM ethylene glycol tetraacetic acid, 1mg/ml aprotinin, 1mg/ml leupeptin, 1mM phenylmethanesulfonyl fluoride, 20mM NaF, 1mM NaPPi, 1mM Na3VO4, 2.5mM dithiothreitol, and centrifuged to obtain supernatant containing total protein. Thirty five micrograms of protein per sample were resolved by SDS PAGE (Invitrogen) and transferred to nitrocellulose membranes (Amersham). Membranes were blocked with 5% non-fat milk diluted in Tris-buffered saline (TBS) and incubated with the appropriate antibodies diluted in TBS-Tween 0.5% bovine serum albumin. Secondary antibodies were purchased from Jackson ImmunoResearch. Super Signal West Pico Chemiluminscence Substrate (Pierce) was used according to the manufacturer’s recommendations to visualize the bands. Antibodies used are mouse monoclonal anti EZH2 (Abnova), anti phosphoSer21-EZH2 (AbCam), anti phosphoSer473-Akt (Epitomics and Cell Signaling), anti phosphoThr308-Akt (Cell Signaling), anti K27 trimethylated (Millipore) or K9 trimethylated histone H3 (AbCam) and anti Phospho-ERK1/2 (Thr202/Tyr204-P; Cell Signaling). Loading was controlled by using an anti Actin or anti GAPDH antibody (Santa Cruz Biotechnology).
Statistical analysis
Comparisons were performed using the Wilcoxon–Mann–Whitney test (for unpaired samples without normal distribution) and Student’s t Test (for paired samples showing normal distribution). Survival analyses (recurrence free or tumor progression in recurrence) according to various variables were performed using the Kaplan–Meyer method and differences between the patient groups were tested by the log-rank test. Contingency analyses were performed using the F Fisher’s exact test. Discrimination between samples showing increased or decreased tumor/normal relative expression of either gene or miRNA expression was made using the median. SPSS 17.0 and Graph prism 6.0 software were used.
ACKNOWLEDGMENTS
We express our acknowledgement to Dr. FX Real (CNIO) and Dr A Carnero (IBIS/HUVR/CSIC/Universidad de Sevilla) for providing us cell lines and plasmids, respectively.
CONFLICTS OF INTEREST
The authors declare no conflicts of interest
FUNDING
FEDER cofounded MINECO grants SAF2012-34378 and SAF2015-66015-R, Comunidad Autónoma de Madrid grant S2010/BMD-2470 (Oncocycle Program), MSyC grant ISCIII-RETIC RD12/0036/0009, PIE 15/00076 and CB/16/00228 to JMP.
Author contributions
CS, MMF, MD, CR, FFLC, CC, CSL, MFG performed experimental work, JD, HGM, FV, FdR and DC obtained clinical data and samples. CS, MMF, MD and JMP performed the analysis and interpretation of data. JMP supervised the study, and all authors contributed to the writing, review, and revision of the manuscript.
REFERENCES
1. Wu X, Hildebrandt MA, Chang DW. Genome-wide association studies of bladder cancer risk: a field synopsis of progress and potential applications. Cancer Metastasis Rev. 2009; 28: 269-80.
2. Burger M, Oosterlinck W, Konety B, Chang S, Gudjonsson S, Pruthi R, Soloway M, Solsona E, Sved P, Babjuk M, Brausi MA, Cheng C, Comperat E, et al. ICUD-EAU International Consultation on Bladder Cancer 2012: Non-muscle-invasive urothelial carcinoma of the bladder. Eur Urol. 2013; 63: 36-44.
3. Knowles MA, Hurst CD. Molecular biology of bladder cancer: new insights into pathogenesis and clinical diversity. Nat Rev Cancer. 2015; 15: 25-41. doi: 10.1038/nrc3817.
4. Knowles MA, Platt FM, Ross RL, Hurst CD. Phosphatidylinositol 3-kinase (PI3K) pathway activation in bladder cancer. Cancer Metastasis Rev. 2009; 28: 305-16.
5. Wagle N, Grabiner BC, Van Allen EM, Hodis E, Jacobus S, Supko JG, Stewart M, Choueiri TK, Gandhi L, Cleary JM, Elfiky AA, Taplin ME, Stack EC, et al. Activating mTOR mutations in a patient with an extraordinary response on a phase I trial of everolimus and pazopanib. Cancer Discov. 2014; 4: 546-53. doi: 10.1158/2159-8290.CD-13-0353.
6. Duenas M, Martinez-Fernandez M, Garcia-Escudero R, Villacampa F, Marques M, Saiz-Ladera C, Duarte J, Martinez V, Gomez MJ, Martin ML, Fernandez M, Castellano D, Real FX, et al. PIK3CA gene alterations in bladder cancer are frequent and associate with reduced recurrence in non-muscle invasive tumors. Mol Carcinog. 2015; 54: 566-76. doi: 10.1002/mc.22125.
7. Santos M, Martinez-Fernandez M, Duenas M, Garcia-Escudero R, Alfaya B, Villacampa F, Saiz-Ladera C, Costa C, Oteo M, Duarte J, Martinez V, Gomez-Rodriguez MJ, Martin ML, et al. In vivo disruption of an Rb-E2F-Ezh2 signaling loop causes bladder cancer. Cancer Res. 2014; 74: 6565-77. doi: 10.1158/0008-5472.CAN-14-1218.
8. Martinez-Fernandez M, Feber A, Duenas M, Segovia C, Rubio C, Fernandez M, Villacampa F, Duarte J, Lopez-Calderon FF, Gomez-Rodriguez MJ, Castellano D, Rodriguez-Peralto JL, de la Rosa F, et al. Analysis of the Polycomb-related lncRNAs HOTAIR and ANRIL in bladder cancer. Clin Epigenetics. 2015; 7: 109. doi: 10.1186/s13148-015-0141-x.
9. Cha TL, Zhou BP, Xia W, Wu Y, Yang CC, Chen CT, Ping B, Otte AP, Hung MC. Akt-mediated phosphorylation of EZH2 suppresses methylation of lysine 27 in histone H3. Science. 2005; 310: 306-10.
10. Earl J, Rico D, Carrillo-de-Santa-Pau E, Rodriguez-Santiago B, Mendez-Pertuz M, Auer H, Gomez G, Grossman HB, Pisano DG, Schulz WA, Perez-Jurado LA, Carrato A, Theodorescu D, et al. The UBC-40 Urothelial Bladder Cancer cell line index: a genomic resource for functional studies. BMC Genomics. 2015; 16: 403. doi: 10.1186/s12864-015-1450-3.
11. Engelman JA. Targeting PI3K signalling in cancer: opportunities, challenges and limitations. Nat Rev Cancer. 2009; 9: 550-62.
12. Cizkova M, Susini A, Vacher S, Cizeron-Clairac G, Andrieu C, Driouch K, Fourme E, Lidereau R, Bieche I. PIK3CA mutation impact on survival in breast cancer patients and in ERalpha, PR and ERBB2-based subgroups. Breast Cancer Res. 2013; 14: R28.
13. Loi S, Haibe-Kains B, Majjaj S, Lallemand F, Durbecq V, Larsimont D, Gonzalez-Angulo AM, Pusztai L, Symmans WF, Bardelli A, Ellis P, Tutt AN, Gillett CE, et al. PIK3CA mutations associated with gene signature of low mTORC1 signaling and better outcomes in estrogen receptor-positive breast cancer. Proc Natl Acad Sci U S A. 2010; 107: 10208-13.
14. Kalinsky K, Jacks LM, Heguy A, Patil S, Drobnjak M, Bhanot UK, Hedvat CV, Traina TA, Solit D, Gerald W, Moynahan ME. PIK3CA mutation associates with improved outcome in breast cancer. Clin Cancer Res. 2009; 15: 5049-59.
15. Jehan Z, Bavi P, Sultana M, Abubaker J, Bu R, Hussain A, Alsbeih G, Al-Sanea N, Abduljabbar A, Ashari LH, Alhomoud S, Al-Dayel F, Uddin S, et al. Frequent PIK3CA gene amplification and its clinical significance in colorectal cancer. J Pathol. 2009; 219: 337-46.
16. Woenckhaus J, Steger K, Werner E, Fenic I, Gamerdinger U, Dreyer T, Stahl U. Genomic gain of PIK3CA and increased expression of p110alpha are associated with progression of dysplasia into invasive squamous cell carcinoma. J Pathol. 2002; 198: 335-42.
17. Woenckhaus J, Steger K, Sturm K, Munstedt K, Franke FE, Fenic I. Prognostic value of PIK3CA and phosphorylated AKT expression in ovarian cancer. Virchows Arch. 2007; 450: 387-95.
18. Shi J, Yao D, Liu W, Wang N, Lv H, Zhang G, Ji M, Xu L, He N, Shi B, Hou P. Highly frequent PIK3CA amplification is associated with poor prognosis in gastric cancer. BMC Cancer. 2012; 12: 50.
19. Kim KH, Roberts CW. Targeting EZH2 in cancer. Nat Med. 2016; 22: 128-34. doi: 10.1038/nm.4036.
20. Martinez-Fernandez M, Rubio C, Segovia C, Lopez-Calderon FF, Duenas M, Paramio JM. EZH2 in Bladder Cancer, a Promising Therapeutic Target. Int J Mol Sci. 2015; 16: 27107-32. doi: 10.3390/ijms161126000.
21. Bracken AP, Pasini D, Capra M, Prosperini E, Colli E, Helin K. EZH2 is downstream of the pRB-E2F pathway, essential for proliferation and amplified in cancer. Embo J. 2003; 22: 5323-35.
22. Foster CS, Falconer A, Dodson AR, Norman AR, Dennis N, Fletcher A, Southgate C, Dowe A, Dearnaley D, Jhavar S, Eeles R, Feber A, Cooper CS. Transcription factor E2F3 overexpressed in prostate cancer independently predicts clinical outcome. Oncogene. 2004; 23: 5871-9.
23. Cao Q, Mani RS, Ateeq B, Dhanasekaran SM, Asangani IA, Prensner JR, Kim JH, Brenner JC, Jing X, Cao X, Wang R, Li Y, Dahiya A, et al. Coordinated regulation of polycomb group complexes through microRNAs in cancer. Cancer Cell. 2011; 20: 187-99. doi: 10.1016/j.ccr.2011.06.016.
24. Friedman JM, Liang G, Liu CC, Wolff EM, Tsai YC, Ye W, Zhou X, Jones PA. The putative tumor suppressor microRNA-101 modulates the cancer epigenome by repressing the polycomb group protein EZH2. Cancer Res. 2009; 69: 2623-9.
25. Carvalho J, van Grieken NC, Pereira PM, Sousa S, Tijssen M, Buffart TE, Diosdado B, Grabsch H, Santos MA, Meijer G, Seruca R, Carvalho B, Oliveira C. Lack of microRNA-101 causes E-cadherin functional deregulation through EZH2 up-regulation in intestinal gastric cancer. J Pathol. 2012; 228: 31-44. doi: 10.1002/path.4032.
26. Varambally S, Cao Q, Mani RS, Shankar S, Wang X, Ateeq B, Laxman B, Cao X, Jing X, Ramnarayanan K, Brenner JC, Yu J, Kim JH, et al. Genomic loss of microRNA-101 leads to overexpression of histone methyltransferase EZH2 in cancer. Science. 2008; 322: 1695-9.
27. Kottakis F, Polytarchou C, Foltopoulou P, Sanidas I, Kampranis SC, Tsichlis PN. FGF-2 regulates cell proliferation, migration, and angiogenesis through an NDY1/KDM2B-miR-101-EZH2 pathway. Mol Cell. 2011; 43: 285-98.
28. Caretti G, Palacios D, Sartorelli V, Puri PL. Phosphoryl-EZH-ion. Cell Stem Cell. 2011; 8: 262-5. doi: 10.1016/j.stem.2011.02.012.
29. Xu K, Wu ZJ, Groner AC, He HH, Cai C, Lis RT, Wu X, Stack EC, Loda M, Liu T, Xu H, Cato L, Thornton JE, et al. EZH2 oncogenic activity in castration-resistant prostate cancer cells is Polycomb-independent. Science. 2012; 338: 1465-9. doi: 10.1126/science.1227604.
30. Kim E, Kim M, Woo DH, Shin Y, Shin J, Chang N, Oh YT, Kim H, Rheey J, Nakano I, Lee C, Joo KM, Rich JN, et al. Phosphorylation of EZH2 activates STAT3 signaling via STAT3 methylation and promotes tumorigenicity of glioblastoma stem-like cells. Cancer Cell. 2013; 23: 839-52. doi: 10.1016/j.ccr.2013.04.008.
31. Kikuchi J, Koyama D, Wada T, Izumi T, Hofgaard PO, Bogen B, Furukawa Y. Phosphorylation-mediated EZH2 inactivation promotes drug resistance in multiple myeloma. J Clin Invest. 2015; 2015. doi: 10.1172/JCI80325.
32. Cancer Genome Atlas Research N. Comprehensive molecular characterization of urothelial bladder carcinoma. Nature. 2014; 507: 315-22. doi: 10.1038/nature12965.
33. Babjuk M, Bohle A, Burger M, Capoun O, Cohen D, Comperat EM, Hernandez V, Kaasinen E, Palou J, Roupret M, van Rhijn BW, Shariat SF, Soukup V, et al. EAU Guidelines on Non-Muscle-invasive Urothelial Carcinoma of the Bladder: Update 2016. Eur Urol. 2016. doi: 10.1016/j.eururo.2016.05.041.
34. Martinez-Fernandez M, Paramio JM, Duenas M. RNA Detection in Urine: From RNA Extraction to Good Normalizer Molecules. J Mol Diagn. 2016; 18: 15-22. doi: 10.1016/j.jmoldx.2015.07.008.
35. Alonso SR, Ortiz P, Pollan M, Perez-Gomez B, Sanchez L, Acuna MJ, Pajares R, Martinez-Tello FJ, Hortelano CM, Piris MA, Rodriguez-Peralto JL. Progression in cutaneous malignant melanoma is associated with distinct expression profiles: a tissue microarray-based study. Am J Pathol. 2004; 164: 193-203.
36. Rimm DL, Camp RL, Charette LA, Costa J, Olsen DA, Reiss M. Tissue microarray: a new technology for amplification of tissue resources. Cancer J. 2001; 7: 24-31.
37. Nenutil R, Smardova J, Pavlova S, Hanzelkova Z, Muller P, Fabian P, Hrstka R, Janotova P, Radina M, Lane DP, Coates PJ, Vojtesek B. Discriminating functional and non-functional p53 in human tumours by p53 and MDM2 immunohistochemistry. J Pathol. 2005; 207: 251-9. doi: 10.1002/path.1838.
38. Martinez-Fernandez M, Duenas M, Feber A, Segovia C, Garcia-Escudero R, Rubio C, Lopez-Calderon FF, Diaz-Garcia C, Villacampa F, Duarte J, Gomez-Rodriguez MJ, Castellano D, Rodriguez-Peralto JL, et al. A Polycomb-mir200 loop regulates clinical outcome in bladder cancer. Oncotarget. 2015; 6: 42258-75. doi: 10.18632/oncotarget.5546.
39. Saeed AI, Sharov V, White J, Li J, Liang W, Bhagabati N, Braisted J, Klapa M, Currier T, Thiagarajan M, Sturn A, Snuffin M, Rezantsev A, et al. TM4: a free, open-source system for microarray data management and analysis. Biotechniques. 2003; 34: 374-8.
40. Huang da W, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009; 37: 1-13.
41. 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.
42. Chen EY, Tan CM, Kou Y, Duan Q, Wang Z, Meirelles GV, Clark NR, Ma'ayan A. Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinformatics. 2013; 14: 128. doi: 10.1186/1471-2105-14-128.
43. Kuleshov MV, Jones MR, Rouillard AD, Fernandez NF, Duan Q, Wang Z, Koplev S, Jenkins SL, Jagodnik KM, Lachmann A, McDermott MG, Monteiro CD, Gundersen GW, et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 2016. doi: 10.1093/nar/gkw377.