The homeodomain transcription factor PROX1 has been linked to several cancer types, including gliomas, but its functions remain to be further elucidated. Here we describe a functional role and the prognostic value of PROX1 in glioblastoma. Low expression of PROX1 correlated with poor overall survival and the mesenchymal glioblastoma subtype signature. The latter finding was recapitulated in vitro, where suppression or overexpression of PROX1 in glioma cell cultures transitioned cells to a mesenchymal or to a nonmesenchymal glioblastoma gene expression signature, respectively. PROX1 modulation affected proliferation rates that coincided with changes in protein levels of CCNA1 and CCNE1 as well as the cyclin inhibitors CDKN1A, CDKN1B, and CDKN1C. Overexpression of SOX2 increased PROX1 expression, but treatment with a CDK2 inhibitor subsequently decreased PROX1 expression, which was paralleled by decreased SOX2 levels. The THRAP3 protein was a novel binding partner for PROX1, and suppression of THRAP3 increased both transcript and protein levels of PROX1. Together, these findings highlight the prognostic value of PROX1 and its role as a regulator of glioblastoma gene expression subtypes, intratumoral heterogeneity, proliferation, and cell-cycle control.
Significance: These findings demonstrate the role and prognostic value of PROX1 in glioblastomas; low PROX1 levels correlate with a mesenchymal gene expression subtype and shorter survival in glioblastoma tumors. Cancer Res; 78(20); 5901–16. ©2018 AACR.
Glioblastoma represents the most common and aggressive primary brain tumor type in adults (1). Its intrusive growth and the plastic nature of the tumor make complete surgical resection impossible and the tumor cells prone to evade chemo- and radiation therapy (2). Stem cell regulatory pathways are shown activated in gliomas supporting self-renewal, tumor maintenance, and survival under stress (3). Furthermore, glioblastomas are very heterogeneous on an intratumoral level, composed of tumor cells displaying different gene expression signatures constituting the different glioblastoma tumor subtypes (4). Thus, a glioma stem-like phenotype, cell motility, and tumor cell heterogeneity are considered significant hurdles to overcome for developing new treatment against these tumors (5, 6).
Glioblastoma may arise from adult neural stem cells or multipotent neural progenitor cells that persist in proliferative niches in the human central nervous system, or alternatively from differentiated cells of different lineages as for example astrocytes or oligodendrocytes (7). The human glioblastoma transcriptome has been found to resemble normal outer radial glial cells and intermediate progenitors (8). In support of this, it was recently shown that glioblastoma initiation is associated to aberrant reactivation of a normal developmental program in the brain (9). Therefore, a better understanding of developmental pathways and their involvement in glioblastoma is thought to lead to new therapeutic possibilities.
Prospero-related homeobox 1 (PROX1) is a transcription factor that mediates cell fate decisions of neuroblasts, as reviewed in ref. 10. In several instances, PROX1 has been shown to play an active role in cancer. For example, PROX1 suppresses the growth of neuroblastoma (11), whereas it enhances colorectal cancer progression as a β-catenin/TCF/LEF target gene contributing to a transition from early to a dysplastic stage (12). Prox1 regulates the number of cancer stem cells, by promoting cell proliferation and thereby expanding the cancer stem cell population in intestinal adenomas and colorectal cancer after activation of the Wnt pathway (13). In human astrocytic brain tumors, the percentage of PROX1+ cells has been shown to increase with tumor grade (14). Furthermore, level of PROX1 predicted survival in grade II gliomas, where a percentage of PROX1+ cells over 10% correlated with worse outcome (15). Moreover, PROX1 has been proposed as a novel pathway-specific prognostic biomarker for high-grade astrocytomas (16). Here we sought further understanding of the role of PROX1 in glioblastoma by analyzing publicly available gene expression data from tumor samples and by performing in vitro experiments.
Materials and Methods
Glioblastoma cell lines U-343 MG, U-343 MGa, U-343 MGa Cl2:6, U118MG, U178MG, U373MG, and U1242MG have previously been characterized (17–21) and references therein. The high-grade glioma cultures 4, 11, and 18 were initially characterized based on gene expression and phenotypes (22), and have in other studies also been referred to as U2975, U2982, and U2987 (23), respectively. For unity, these cultures are referred to as U2975(4), U2982(11), and U2987(18) in this study. The cell lines and cultures above were retrieved from authors’ lab stocks (22) and subsequently cultured for less than 3 months. All cell lines were maintained in DMEM (Thermo Fisher Scientific) with 10% FBS, 2 mmol/L l-glutamine, and penicillin–streptomycin at 37°C with 5% CO2, and 5% O2 and >95% humidity. The other glioma cell lines (U87MG, U251MG, M059J, and M059K), as well as the osteosarcoma line U2-OS and colon cancer line SW480 were purchased from ATCC and cultured per standard guidelines for less than a month. SOX2 and YFP (control) overexpressing U-343 MG cultures were generated by lentiviral transduction with pLEX-Blast-V5-SOX2 and pLEX-Blast-V5-YFP by methods previously described (24). To confirm the relationship between the U-343 clones they were subjected to STR profiling at NGI-Uppsala, SciLifeLab, Uppsala University, using the AmpFISTR Identifiler PCR Amplification Kit (Thermo Fisher; Supplementary Table S1). The cultures used in this study were screened for Mycoplasma using the MycoAlert Mycoplasma Detection Kit (LT07-218; Lonza). Stock solution of CVT-313 (Santa Cruz) was prepared by dissolving powder in DMSO, stored at −20°C, and used by dissolving stock solution in cell culture media and incubating with cells at concentration and time indicated.
Lentiviral transduction and transfection of plasmids or siRNA constructs
To generate a PROX1 overexpressing cell culture, U-343 MG-PROX1, purified lentivirus (Lentifect) for PROX1 transduction into U-343 MG was purchased (PROX1 containing virus LPP-F0925-Lv105-100-S and negative control virus LPP-NEG-LV105-100-C; GeneCopoeia), and used according to manufacturer's instructions. To study PROX1 suppression in U-343 MGa, the cell line U-343 MGa-shPROX1 was generated using lentiviral vectors against PROX1 (TRCN0000232123; Sigma-Aldrich) or unrelated control shRNA targeting GFP (TRCN0000072178, clonetechGfp_438s1c1, previously described; ref. 24). Virus production was performed by calcium-phosphate-mediated cotransfection of HEK 293T cells with packaging plasmids (MISSION Lentiviral Packaging Mix, SHP001; Sigma-Aldrich). Twenty-four hours after transfection, the different supernatants were collected two times with 24-hour intervals, filtered and then used to infect U-343 MGa cells cultured under optimal conditions. Lentiviral transductions of glioblastoma cells were conducted at MOI of 1 to 2. The cells were analyzed and selected in the presence of 1 μg/mL puromycin for 72 hours or until all nontransduced control cells were killed. Transduced cells were kept in culture for 3 days, prior to experiments. For RNAi experiments SMARTpool mixes, ON-TARGETplus siRNA (GE Healthcare Dharmacon), were used to target PROX1 (L-016913-00-0005), SOX2 (L-011778-00-0005), or THRAP3 (L-019907-00-0010) or nontargeting control pool (D-001810-10-05) according to the manufacturer's instructions. Lipofectamine 2000 (Thermo Fisher Scientific), and RNAiMAX (Thermo Fisher Scientific) reagents were used for the transfection of plasmid DNA and siRNA, respectively, according to the manufacturer´s instructions.
Publicly available gene expression data
The Cancer Genome Atlas (TCGA) data from 539 glioblastoma samples, published in 2013 (25) was downloaded from UCSC University of California Santa Cruz Xena portal (www.xenabrowser.net). Data for glioblastoma gene expression based subtypes and patient survival were downloaded from cBioPortal for Cancer Genomics at http://cbioportal.org and merged with the expression data. Samples categorized as G-CIMP or taken from relapse tumors were removed yielding a final set of 375 tumors. Gene expression data for cell cultures was downloaded from the Cancer Cell Line Encyclopedia (CCLE) (26) database at www.broadinstitute.org/ccle. Cell cultures described to be derived from glioblastomas were selected, which generated a set of 45 cell cultures. Gene expression data for 48 glioblastoma cell cultures was downloaded from the human glioblastoma cell culture (HGCC) portal at www.hgcc.com (27).
Survival analysis of patients with glioblastoma
Analyses of overall and disease-free survival were conducted using data available for patients with glioblastoma, published in 2008 (28), extracted from http://cbioportal.org and patients with glioma from http://gliovis.bioinfo.cnio.es/. Statistical analysis of the data was performed using an unpaired two-tailed t-test assuming normal distributions (Prism 6.0; Graphpad software Inc.). Kaplan–Meier survival curves were plotted for patients grouped according to PROX1 expression levels above or below 1, and THRAP3 above or below 0.3. Log-rank (Mantel–Cox) tests were used to determine differences between survival curves.
Gene set enrichment analysis
For the gene set enrichment analysis (GSEA; ref. 29), we used the software on the GenePattern online server (www.genepattern.broadinstitute.org) module version 18 and the Molecular Signatures Database (MSigDB; ref. 30) version 6.0. Overlaps of PROX1 correlated gene sets were computed with the following categories in MSigDB collections using the default setting (FDR q value below 0.05). For GSEA of the CCLE data for 45 glioblastoma cell cultures, we divided the cultures into two groups based on the expression of PROX1. To determine at what cutoff to divide the sample data into two groups according to PROX1 expression, the cultures were first ranked per expression level of PROX1. The t tests were then performed between high and low PROX1 expression in the samples for all the genes. The t test was performed for all permutations starting with the two cultures with the highest PROX1 versus the rest, the three cultures with the highest PROX1 level versus the rest, and so on. For each permutation, the t-test significance was determined for all genes and internally compared with the P-value–based rank position of PROX1. The cutoff for high versus low PROX1 expressing cultures was set to where PROX1 had the highest relative P-value ranking compared with other genes. In other words, the cut-off was set to where PROX1 had the lowest relative P value compared with other genes. We found the optimal cutoff when the cultures were divided into 25 low versus 20 high PROX1 expressing cultures where PROX1 was the 11th most significant gene to distinguish the samples. We then used the 25 low versus 20 high PROX1 expressing culture division for further GSEA analyses. The analyses were run with 1000 permutations, phenotype as permutation type, collapsed data set, and the AFFYMETRIX chip platform file option. In this case, we investigated the gene set databases that included all subcategories including positional, curated, motif, computational, gene ontology, oncogenic signatures, immunologic signatures, and hallmarks.
RNA purification and quantitative real-time PCR
Total cellular RNA was extracted from 70% to 80% confluent cultures using PureLink Kits (Ambion RNA extraction products) according to manufacturer's instructions (Thermo Fisher Scientific). RNA concentrations were measured with a NanoDrop spectrophotometer and samples were stored at −80 °C. For gene expression analysis, real-time quantitative RT-PCR (qRT-PCR) was performed. The Power SYBR Green RNA-to-CT 1-Step Kit (Thermo Scientific) was used according to the manufacturer's instructions in conjunction with a 7500 StepOnePlus Real-Time PCR system (Thermo Scientific). GAPDH was used as the internal standard reference. Gene expression levels were calculated using the relative ΔΔCt method. The primers used for qRT-PCR are listed in Supplementary Table S2. Statistical analysis of the data was performed using an unpaired Student two-tailed t test (Prism 6.0; Graphpad software Inc.). P < 0.05 was considered statistically significant.
RNA sequencing and gene expression analysis
RNA quality was assessed by electrophoresis on the Tapestation instrument (Agilent Technologies) and samples with a RIN score above 8 were used in further analysis. RNA library preparation, sequencing, raw data processing, and quality control were performed at the National Genomics Infrastructure at Science for Life Laboratory, Stockholm, Sweden). Reads were mapped using Tophat (2.0.4.) and FPKM values were generated with Cufflinks (2.1.1.). For further gene expression based analyses, we filtered transcripts from the expression data (FPKM values) that were present for the CCLE gene-centric RMA-normalized Affymetrix data, retaining about 12,042 genes. Significant difference in expression was calculated in Prism applying a t test to two biological groups: (i) U-343 MG-control versus U-343 MG-PROX1 with two replicates in each group and (ii) U-343 MGa-control versus U-343 MGa-shPROX1 with three replicates in each group. Volcano plots were generated by using genes with FPKM values over 0.01, a two-fold change in expression and a P value below 0.05. Venn diagrams were generated using GeneVenn (http://genevenn.sourceforge.net) to identify PROX1 induced or suppressed genes in an extended list only with more than two-fold difference in expression between groups of replicates, where the average of either replicate group was not 0. The PROX1-induced and -repressed genes were investigated at MSigDB (30) to identify enriched gene sets. RNA-seq data have been deposited at EBI ArrayExpress database with the accession number E-MTAB-6991. To generate gene expression heat maps, the Heat map module at the GenePattern site was used (31).
Glioblastoma gene expression based subtype analysis
To assess changes of glioblastoma gene expression subtypes due to PROX1 overexpression in U-343 MG-PROX1 or suppression in U-343 MGa-shPROX1, we counted the number of subtype defining genes whose expression were altered in our RNA-seq experiment following modulation of PROX1. To assess the change of each gene expression subtype we used the subtype defining centroid genes previously described by (32), where each subtype was described by a set of 210 centroid genes with either high or low expression as compared with the others. For instance, in our analysis corresponding to U-343 MG-PROX1 versus U-343 MG for the mesenchymal subtype, 35 of 128 genes with available measurements were increased more than one-fold, whereas 12 of 33 genes decreased more than one-fold. This resulted in a rescaled index of −42 for the mesenchymal subtype, when −100 is full negative change of subtype, 0 is neutral, and 100 is full positive change.
Cell proliferation assay
Cell proliferation assay was performed by cell counting using Countess automated cell counter (Thermo Fisher Scientific). For proliferation analysis, cells were seeded in 12-well plates and the total cell number was determined at the indicated time points. The experiments were repeated three times using triplicate samples. Statistical analysis of the data was performed using an unpaired Student two-tailed t test. P values <0.05 were considered statistically significant.
Immunoblotting, coimmunoprecipitation, and mass spectrometry analysis
For detection of specific antigens the following antibodies were used; PROX1 (R&D Systems; AF2727), PROX1 (Abcam; ab37128), p53 (Santa Cruz; sc-6243), p53 (Santa Cruz; sc-6243-G), p57 Kip2 (Abcam; ab75974), P21 (Abcam; ab7960), P21 (Santa Cruz; sc-397), cyclin B1 (Santa Cruz; sc-245), cyclin D1 (Santa Cruz; sc-8396), cyclin E (Abcam; ab7959–1), cyclin A2 (Abcam; ab16726), Cdk2 (Santa Cruz; sc-6248), TRAP150 (Santa Cruz; sc-48779), β-actin (Sigma-Aldrich; clone AC15), SOX2 (EDM Millipore; AB5603), GFAP (Abcam; ab10062), and FN1 (BD Biosciences; 610078). For analysis of detergent-soluble proteins, we used 0.5% Nonidet P-40 lysis buffer with complete protease and phosphatase inhibitors (PhosSTOP, cOmplete ULTRA; Roche). Analysis of detergent-soluble proteins, coimmunoprecipitation (co-IP), and IB were conducted as described (33). For co-IP (“large-scale”) followed by mass spectrometry, the cells were harvested and prepared according to the nuclear complex co-IP kit instructions (Nuclear Complex Co-IP Kit; Nordic Biolabs), using the high stringency buffer option (Active Motif). For the initial co-IP in U2-OS cells, we used an antibody against the FLAG tag (Flag M2, F1804; Sigma-Aldrich). Immunoprecipitates were separated by SDS-PAGE (NuPAGE Novex 4–12% Bis-Tris Protein Gels; Life Technologies), the gels stained with colloidal CBB Kit (Life Technologies), bands unique to the PROX1-FLAG antibody lane were cut out, and identified by mass spectrometry (LC/MS-MS; Alphalyse A/S; Denmark). Protein samples were reduced and alkylated with iodoacetamide, that is carbamidomethylated, and trypsin digested. Peptides were concentrated on a ZipTip micropurification column and eluted onto an anchorchip target to be used for analysis on a Bruker Autoflex III MALDI TOF/TOF instrument. MALDI MS-MS was performed on 15 peptides for peptide fragmentation analysis. Peptide tolerance was set to 60 ppm (one miscleavage allowed). The MS and MS-MS spectra were combined and Mascot software version 2.2.03 was used. The database NRDB was used for protein identification.
Procedures for immunofluorescence staining and microscopy analysis have previously been published (33).
Low PROX1 mRNA level in glioblastoma tumors correlates with mesenchymal gene expression subtype and with shorter survival
Because PROX1 has been shown to have diagnostic and prognostic value in gliomas, we utilized gene expression datasets to further explore this. First, we investigated a combined low- and high-grade glioma dataset composed of 226 grade II, 244 grade III, and 150 grade IV gliomas from the GlioVis portal (http://gliovis.bioinfo.cnio.es/) to assess PROX1 expression levels (34). PROX1 transcripts were significantly less in grade IV tumors as compared with grade II (P < 0.0001) and III (P < 0.0001) tumors (Fig. 1A). In a survival analysis, patients with glioma with high PROX1 expression levels survived longer, with a median survival of 64.6 months versus 40.5 for the subgroup with low PROX1 expression levels (P < 0.0001; Fig. 1B). To investigate the prognostic value of PROX1 exclusively in glioblastomas, we extracted gene expression data available for 206 patients with glioblastoma from TCGA at http://cbioportal.org. Survival differences for different PROX1 expression z-score cutoffs (−2 to +2 thresholds) were analyzed and values between 0.8 and 1.3 for overall survival, and 0.9 and 1.3 for disease-free survival yielded significant differential survival (Fig. 1C). Specifically, the median overall survival for patients with a PROX1 expression higher than z-score 1 was 33.6 months as compared with 12.3 months for the subgroup with lower expression (16 vs. 190 patients; P = 0.0015), whereas the median disease-free survival for patients with a PROX1 expression higher than z-score 1 was 16 months compared with 6.5 (13 patients vs. 149 patients; P = 0.0118; Fig. 1D and E). To investigate the correlation between PROX1 expression and molecular subtypes in glioblastoma, we explored a TCGA dataset (25) containing 375 patients with glioblastoma with available subtype calling. We divided the patients into four equally sized groups based on PROX1 high to low expression values (Fig. 1F). We found that the group with highest PROX1 expression, patients 1 to 94 (group 1), was mostly represented by tumors of classical and proneural subtypes (Fig. 1G). On the contrary, in tumors with lowest PROX1 transcript levels, patients 283–375 (group 4), mesenchymal and neural subtypes were the most abundant. Specifically, significant lower PROX1 expression (z-score) was found in mesenchymal subtype compared with proneural (P < 0.0001) or classical subtype (P < 0.0001; Fig. 1H). Together, the proneural and classical have overall a higher PROX1 expression, whereas neural and mesenchymal cases have lower PROX1 expression levels. Thus, low PROX1 expression defines a subset of patients enriched for mesenchymal and neural glioblastoma subtypes that display shorter disease free and overall survival. Of note, the neural subgroup has been associated with contamination of normal neuroepithelial cells and nontumor-specific expression signature and should therefore be considered with caution (35).
Expression analysis of PROX1 in different glioblastoma cell lines identifies coexpressed genes and functional implications
To assess gene expression differences between glioblastoma cultures with high or low PROX1 expression levels, we performed GSEA for gene expression data from 45 glioblastoma cell lines in the CCLE (26). Twenty-five cultures were defined to contain high PROX1 levels versus 20 with low. A ranked list of genes was generated and queried against gene sets deposited at MSigDB (29). A total of 137 gene sets were significantly enriched for glioma cultures with relatively high PROX1 expression and 71 gene sets for low PROX1 expressing glioblastoma cultures (Supplementary Table S3).
Analysis of protein and gene expression in glioma cell cultures and generation of stable cell lines to assess PROX1 function
To evaluate PROX1 expression we compared protein levels in a set of glioma cancer cell cultures including the clonal model U-343 (19) by immunoblotting and immunofluorescence staining. We found that U-343 MGa 31L, U-343 MG, U87MG, U2982(11), U2975(4), M059J, and M059K had undetectable levels of PROX1 protein, whereas it was readily detected in U-343 MGa Cl2:6, U-343 MGa, U373MG, U251MG, and U2987(18) (Fig. 2A and B). PROX1 also displayed a nuclear localization as previously described in the positive control colorectal cancer cell line SW480 (Fig. 2B; ref. 12). U-343 MG and U-343 MGa (and its derivatives U-343 MGa Cl2:6 and U-343 MGa 31L) were originally derived from one glioblastoma tumor, and have previously been characterized (17, 19). The U-343 MG cells express fibronectin (FN1) but not glial fibrillary acidic protein (GFAP), and display a bipolar fibroblastic cell shape. The U-343 MGa cells however express GFAP but not FN1, and they display a polygonal morphology. Corresponding levels of glioma cell protein markers SOX2, GFAP, and FN1 were also determined in the panel of cell lines by immunoblotting. The astroglial stem cell related proteins SOX2 and GFAP were present in all the PROX1 expressing cells, whereas FN1 was expressed in all PROX1 negative cultures, except U-343 MGa 31L, and U87MG (Fig. 2A). Thus, PROX1 is differentially expressed in glioma cell lines and is often coexpressed with SOX2 and GFAP (Fig. 2A). The relatively high expression of FN1 in U-343 MG, and high expression of GFAP in U-343 MGa cultures was confirmed by qRT-PCR (Fig. 2C). Given that the two cell lines U-343 MG and U-343 MGa are derived from the same tumor but differ in levels of PROX1 protein (Fig. 2A), we chose these two for further analysis of PROX1 function. Using lentivirus, we generated stable cell lines with either PROX1 overexpression or shRNA-mediated PROX1 suppression. In U-343 MG-PROX1 or MGa-shPROX1, we achieved stable overexpression or suppression of PROX1 at the mRNA and protein level, respectively (Fig. 2D and E). Experiments using double-immunofluorescence staining have previously shown the existence of GFAP and PROX1 double positive cells in grade IV brain tumors (14). To investigate if PROX1 is a marker or driver of the U-343 MGa phenotype defined by high GFAP expression, we performed qRT-PCR analysis. We found increased GFAP mRNA levels upon overexpression of PROX1 in U-343 MG, and decreased levels in U-343 MGa upon suppression of PROX1 (Fig. 2F). Accordingly, immunofluorescence staining and immunoblotting of U-343 MGa-shPROX1 samples revealed reduced GFAP protein (Fig. 2G and H). Similar to U-343 MGa, the high-grade glioma culture U2987(18) expresses GFAP endogenously (22). We found that a transient suppression of PROX1 in this culture was sufficient to decrease GFAP transcript levels as determined by qRT-PCR (Fig. 2I). Taken together, PROX1 correlates with and regulates GFAP expression in human glioblastoma cells, thus could be a driver of the U-343 MGa phenotype.
Global gene expression analysis by RNA-seq upon PROX1 overexpression or suppression in glioblastoma cells
To gain a more comprehensive understanding of PROX1 as a transcription factor in glioblastoma, gene expression analysis was conducted by RNA sequencing. We analyzed differential gene expression upon PROX1 overexpression in U-343 MG or suppression in U-343 MGa cells (Fig. 3A). A total of 394 genes were significantly upregulated more than two-fold and 139 downregulated significantly more than two-fold after PROX1 overexpression. In the same way, 137 versus 93 genes were altered upon PROX1 suppression (Fig. 3A). For further analyses, we compared intersectional genes between PROX1 overexpression and suppression using an extended gene list applying only a cut-off of two-fold expression change. We selected 1,472 genes whose expression was increased more than two-fold upon PROX1 overexpression in U-343 MG and 733 that were decreased more than two-fold upon PROX1 suppression in U-343 MGa. From these, we identified 237 overlapping genes, which we defined as PROX1-induced genes (Fig. 3B). Conversely, the expression of 855 genes increased more than two-fold upon PROX1 suppression in U-343 MGa, and 1199 genes decreased more than two-fold upon PROX1 overexpression in U-343 MG. Out of these, 286 overlapped and are from here on referred to as PROX1-repressed genes (Fig. 3C).
From a database (MSigDB) analysis of the 237 PROX1-induced genes, we found enrichment of GO terms including neurogenesis and regulation of cell proliferation (Fig. 3D) and overlaps with hallmark gene categories (Fig. 3E). Similarly, results from the analysis performed for the PROX1-repressed genes are presented in Fig. 3F and G. Of note, both PROX1-induced and PROX1-repressed gene sets were enriched for hallmark genes implicated in epithelial to mesenchymal transition.
To assess the relevance of our identified PROX1-induced and -repressed genes in vivo, we performed further analysis of these by comparing them with PROX1 correlated genes in glioblastoma TCGA expression data (Fig. 3H and I). The Venn diagram in Fig. 3H presents that 16% (38 of 237) of PROX1-induced genes identified here overlap with the top 2,000 genes positively correlated to PROX1 expression in TCGA glioblastoma cases. Conversely, 32% of PROX1-repressed genes overlapped with the top 2,000 negatively correlated to PROX1 expression (Fig. 3I). GO terms and hallmark gene categories were also identified for a merged list of the intersecting genes from Fig. 3H and I, as illustrated in Fig. 3J. Collectively, our analysis of RNA-seq data facilitated the identification of PROX1-regulated functional gene networks, including those controlling neurogenesis, regulation of cell proliferation, and epithelial to mesenchymal transition.
PROX1 regulates gene expression-based subtypes in glioblastoma
Because PROX1 expression in our survey of TCGA data correlated with glioblastoma subtypes (Fig. 1G and H), we used our gene expression data to assess its regulation of these subtypes. We calculated a support index for each signature by counting the number of classifier genes (32) with altered expression for each corresponding subtype. We found that PROX1 overexpression decreased the support index for mesenchymal subtype and increased it for nonmesenchymal subtypes in U-343 MG cells. Upon a Chi-square test the change of genes supporting a mesenchymal subtype was found significant with P < 0.000095 (Fig. 4A; Supplementary Fig. S1). Conversely, PROX1 suppression significantly (Chi-square, P-value 0.017) increased the mesenchymal support index in U-343 MGa cells (Fig. 4B; Supplementary Fig. S1). Thus, PROX1 modulation altered the glioblastoma gene expression subtype in cell cultures, where overexpression shifted cells into a nonmesenchymal and suppression into a mesenchymal subtype.
To pinpoint a system independent core set of PROX1-regulated genes, we filtered genes that were correlated with PROX1 in TCGA, CCLE, and HGCC resource (Fig. 4C; ref. 27). In an extended exploration of these datasets, we found similarities between gene expression patterns in glioblastoma tumor material and cell cultures regarding PROX1 correlated genes, and identified GO terms suggesting involvement of PROX1 in several biological processes including regulation of nervous system development and chromatin modifications (Supplementary Fig. S2). The analysis revealed 82 intersecting PROX1 correlated genes between CCLE, HGCC, and TCGA (Fig. 4C). Strikingly, of these the top 10 PROX1-upregulated genes, ranging from a 203- to a 15-fold increase, all have been reported to be involved in neuroglial development, including SEMA6A, NES, PTPRZ1, FABP7, POU3F2, and GAP43 (Fig. 4C). Moreover, the intersect genes that were upregulated more than two-fold were enriched for gene sets including the “classical” and “proneural” subsets, “neural development” and “upregulated in glioblastoma cells with capacity to form neurospheres” (Supplementary Table S4). We also noted that FABP7, PTPRZ1, and GAP43 were among the top stemness signature genes in the single-cell gene expression data described by Patel and colleagues (4). Upon analysis of the Patel and colleagues stem cell signature genes, we found that PROX1 increased the expression level of 27 of the 44 genes (61%) with available expression values, highlighting the connection between PROX1 and stemness regulation in glioblastoma (Fig. 4D). To investigate PROX1 expression pattern at a single-cell level, we analyzed the single-cell RNA-seq data available for five glioblastomas (4). Because a PROX1 expression measurement was not available in the dataset, we developed a PROX1 proxy profile to assess PROX1-related expression patterns. Specifically, we calculated a proxy profile based on the average expression of the top five PROX1-regulated genes from the Patel and colleagues stemness signature, composed of SPARCL1, PTPRZ1, FABP7, GLUL, and GAP43. We filtered and ranked the top 100 positively correlated genes to the proxy profile in the Patel and colleagues dataset. In each of the tumors we observed cell populations with varying expression levels of genes positively correlated with the PROX1 proxy (Fig. 4E). By analyzing single-cell data, we demonstrate cell heterogeneity within tumors regarding genes correlated to the PROX1 proxy profile. Furthermore, we found that 71% of genes with a fold change more than 2 (30 of 42; among the top 100) could be upregulated by PROX1 overexpression, based on the RNA-seq analysis of U-343 MG-PROX1 (Fig. 4E). Together, the combined analysis of functional effects of PROX1 in cell cultures, coexpression in multiple gene expression datasets, including both bulk tumor and single-cell measurements, describes PROX1 as a regulator of glioblastoma tumor evolution that distinguishes subpopulations of cells in heterogeneous tumors and whose loss of expression enables a switch from a nonmesenchymal to a mesenchymal phenotype.
PROX1 overexpression induced cell-cycle regulators and increased cell proliferation
It has previously been shown that SOX2+/GFAP+ glioblastoma cultures have higher rates of proliferation and tumorigenicity, as compared with those with mesenchymal profile characterized by high expression of FN1 (22). We explored data from the HGCC biobank and found that high expression of PROX1 correlated with high tumorigenicity of the HGCC cultures as previously assessed by intracranial injection into NOD-SCID mice (Fig. 5A; ref. 27). We thus compared the proliferation rates between U-343 MG-control and U-343 MG-PROX1, and between U-343 MGa-control and U-343 MGa-shPROX1 cell lines by cell counting. Overexpression of PROX1 increased proliferation of U-343 MG-PROX1 cultures (Fig. 5B), whereas PROX1 suppression decreased the proliferation rate of U-343 MGa (Fig. 5C). In the RNA-seq experiment, we noticed changes in the levels of several cyclin and cyclin-dependent kinase (CDK) genes following PROX1 modulation (Fig. 5D). To investigate if the altered proliferation was paralleled by changes in levels of cell-cycle regulatory proteins, we performed immunoblotting analysis of cyclins. We found increased levels of cyclins A1 (CCNA1), B1 (CCNB1), D1 (CCND1), and E1 (CCNE1) proteins in U-343 MG-PROX1. PROX1 suppression resulted in decreased cyclins A1 and E1 in U-343 MGa-shPROX1, whereas no noticeable change was observed for cyclin D1 or cyclin B1 (Fig. 5E). In addition, we observed decreased p53 (TP53) and increased p21 (CDKN1A), p27 (CDKN1B), and p57 (CDKN1C) protein levels after PROX1 overexpression, whereas the levels remained unchanged in U-343 MGa-shPROX1 compared with the control (Fig. 5E). In summary, we found increased or decreased proliferation upon PROX1 overexpression or suppression, respectively, which coincided with changes of cell-cycle regulators.
PROX1 interacts directly and colocalizes with THRAP3 in the nucleus
The function of PROX1 and its regulation remains poorly understood in molecular details. We thus performed experiments to identify proteins that interact with PROX1. Initially, we used an osteosarcoma cell line that stably expresses PROX1-FLAG or empty-FLAG. A nuclear complex co-IP, followed by mass spectrometry analysis, identified proteins associated with PROX1, of which none has previously been described as a PROX1 interacting protein (Supplementary Fig. S3). One of the identified PROX1 associated proteins was THRAP3 (thyroid hormone receptor associated protein 3, aka TRAP150). THRAP3 peptides detected by MS/MS sequencing are shown (Fig. 6A). Subsequently, we confirmed the interaction of endogenous PROX1 with THRAP3 protein in the U-343 MGa cells via a reciprocal co-IP experiment (Fig. 6B). Given the important role of THRAP3 in several cellular processes related to transcriptional activities, we decided to investigate a functional link between THRAP3 and PROX1. We found that both proteins localize to the nucleus (Fig. 6C). Upon THRAP3 suppression following 5 days of siRNA treatment PROX1 levels were found to have increased (Fig. 6D). Similarly, 2-day THRAP3 siRNA treatment also led to increased PROX1 protein levels in U-343 MGa as well as in U-343 MG-PROX1 cells (Fig. 6E). In line with the result that PROX1 increased GFAP transcript and protein levels (Fig. 2F), we found that GFAP protein levels also increased upon THRAP3 suppression, although seemingly independent of PROX1 (Fig. 6F). To investigate if the increase in PROX1 and GFAP levels was due to transcriptional effects, we also assessed transcript levels and found an increase of PROX1 and GFAP expression upon THRAP3 suppression (Fig. 6G). In parallel, no transcript changes were observed for the oligodendroglial and neuronal makers OLIG2 and MAP2, respectively (Fig. 6G). To investigate the association between THRAP3 and PROX1 transcripts in glioma, we explored the TCGA data. THRAP3 expression was significantly less in glioma grade IV as compared with grade II (P < 0.0001) and grade III (P < 0.0001) tumors (Fig. 6H), displaying a similar pattern to PROX1 expression in the same samples (cf. Fig. 1A). Moreover, a positive correlation was found between THRAP3 and PROX1 transcript levels in glioblastoma (Pearson r = 0.58; Fig. 6I). Finally, patients with high or low THRAP3 transcript levels differed with regard to overall and disease-free survival (Fig. 6J, K, and L). A survey of genes with the highest correlation to both PROX1 and THRAP3 expression in glioblastoma samples revealed overlaps with hallmark gene sets including oxidative phosphorylation, biological processes including nucleoside triphosphate metabolic process, and cellular compartments such as respiratory chain and mitochondria (Supplementary Fig. S3). Taken together, these results identify THRAP3 as a novel interacting partner for PROX1, and suggest its involvement in the regulation of PROX1 and GFAP transcript and protein levels.
Correlation and regulatory connection of PROX1 and SOX2 in glioblastoma
SOX genes are key regulators of embryonic development of the CNS, and its maintenance in adults (36). Specifically, SOX2 has been reported to bind to the PROX1 promoter in a global ChIP-seq analysis (37), which suggests that SOX2 could be a determinant for PROX1 regulation. It was recently shown that SOX2 is stabilized through phosphorylation by cyclin E-CDK2 complexes in glioma (38). Given that PROX1 was found to affect expression of cyclins in glioma cells, we investigated the connection between PROX1, SOX2, and cyclins. As illustrated in Fig. 2A, we noted a positive correlation between PROX1 and SOX2 protein levels in a panel of glioma cell lines. In addition, we found a significant positive correlation between expression of PROX1 and SOX2 in 375 cases from TCGA clinical data (r = 0.43; Fig. 7A), in 48 glioblastoma cell cultures from HGCC (r = 0.38), and in 45 glioblastoma cell lines from CCLE (r = 0.55; Fig. 7B). In line with SOX2 being an upstream regulator of PROX1 we found that SOX2 overexpression in U-343 MG caused elevated PROX1 transcript and protein levels (Fig. 7C and D), and SOX2 suppression in U-343 MGa decreased PROX1 protein levels (Fig. 7E). However, suppression of PROX1 in U-343 MGa did not decrease SOX2 protein levels (Fig. 7D). Given that CDK2 activity can stabilize SOX2 (38), we investigated if treatment of U-343 MGa cells with the CDK2 inhibitor CVT-313 compound would result in decreased PROX1 protein levels and found that CVT-313 treatment decreased both SOX2 and PROX1 protein levels within 24 hours (Fig. 7F). It was confirmed that treatment of cells with 10 μmol/L of the compound for 24 hours also decreased cyclin E1 levels (Fig. 7F). Together, these data indicate that SOX2 controls PROX1 expression and that CDK2 inhibition leads to decreased levels of both SOX2, PROX1, and cyclin E1. It was recently reported that THRAP3 binds to SOX9 and regulates transcription (39). Also, we have previously reported a subset of glioblastoma cell cultures with high expression of SOX2 and SOX9 (22). We therefore investigated the correlation of expression between SOX2, SOX9, and THRAP3 in glioblastoma samples and found that THRAP3 was highly correlated with SOX2 and SOX9 (Fig. 7G and H). Including the results presented in Fig. 6i, this connects the expression levels of SOX2, THRAP3, PROX1, and SOX9 in glioblastoma samples.
Results presented in this work indicate that PROX1 mRNA level fluctuates with tumor grade of gliomas, and is expressed at a lower level in grade IV gliomas. In particular, for grade IV gliomas we found that the mesenchymal subtype was associated with relatively low PROX1 levels. In the light of recent findings about tumor heterogeneity and the intermixture of mesenchymal and nonmesenchymal gene expression subtype cells (4), the data present a mechanism of how these subtypes may be regulated by PROX1 during tumor cell evolution. In this analysis PROX1 is less expressed in the mesenchymal-like U-343 MG cell line as compared with its sibling line U-343 MGa, and the clonal derivatives thereof, U-343 MGa Cl2:6 and U-343 MGa 31L. According to a suggested model of tumor progression, from a nonmesenchymal to a mesenchymal subtype (40), the U-343 MG cells should then have arisen from a U-343 MGa ancestry clone where PROX1 expression has been lost. Per our gene expression analysis, the U-343 MGa culture became more mesenchymal-like upon PROX1 suppression, and conversely U-343 MG shifted to a nonmesenchymal subtype profile upon PROX1 overexpression. Thus, our data show that PROX1 can regulate the glioblastoma subtype in a bi-directional manner. In extension, the switch between glioblastoma transcriptional subtypes is reversible and thus allows glioblastoma cells to fluctuate between subtypes and give rise to heterogeneous tumor cell populations. The subtype gene expression signature is relatively stable for the U-343 cultures in vitro, suggesting an underlying genetic or epigenetic cause. Moreover, whether external queues from neighboring cancer cells and tumor microenvironment also may affect PROX1 levels remains to be studied. Beyond strictly clonal differences, maintained by genetic or epigenetic events, this could be an underlying explanation for previous observations where not all cancer cells within single tumors displayed similar PROX1 protein levels (14). Therefore, our findings underscore PROX1 function in gliomas where its decreased mRNA level is connected with tumor cell evolution and intratumoral heterogeneity.
PROX1 is not an established oncogene or tumor suppressor, but results presented here place it downstream of SOX2 that is a prominent glioblastoma cancer driver due to its focal amplification and tumor-initiating capacity (41, 42). The high correlation of PROX1 and SOX2 expression in tissue, and previous reports of SOX2 binding to the PROX1 promoter identified by ChIP-seq analysis of glioblastoma (37) highly suggest that PROX1 is regulated by SOX2. Of note, our GSEA of glioblastoma cell cultures from CCLE with high or low PROX1 expression identified a significant enrichment of SOX2 target genes. Thus, in addition to PROX1 being part of a normal CNS developmental program, we propose that abnormally regulated PROX1 is an important component of astrocytic tumor formation and growth. The decreased proliferation rate observed upon PROX1 suppression in U-343 MGa raises the question if PROX1 has a growth promoting effect in early stage tumors; whether it decreases with increasing grade and potentially regulates transcriptional subtypes during tumor progression. It has been reported that the tumor-initiating capacity in stem-like glioma cell cultures with high growth rate, as opposed to mesenchymal-like cultures with low proliferative rate and tumor-initiating capacity, is maintained by SOX2 (22). Consistent with this, in our analysis of HGCC data (27), PROX1 expression was also associated with higher tumor initiating capacity. However, we found that glioblastomas with lower PROX1 levels, mesenchymal subtype, have worse outcome. Similarly, Patel and colleagues demonstrated that glioblastomas with a high representation of proneural cells were associated with longer survival, compared with those with a high intermixture of mesenchymal cells (4). Furthermore, the glioblastoma cancer stem cell signature was found strongest in individual cells conforming to the proneural and classical subtypes, but underrepresented in cells of the mesenchymal subtype. This presents a conceptual contradiction, where tumor data connects the mesenchymal subtype with poor outcome, whereas in an experimental setting nonmesenchymal stem-like cells display higher tumorigenic potential. We speculate that as gliomas become more mesenchymal when they progress, mesenchymal subtype glioblastomas may represent more progressed tumors that per se have shorter survival. Many experimental models do not reflect this and thus may be the underlying cause for this discordance. Furthermore, the intermixture of nonmesenchymal and mesenchymal cells may have tumorigenic effects that are not observed in tumor models with low clonal heterogeneity.
The present results show that PROX1 and GFAP are positively correlated in glioblastomas and cell lines datasets, and GFAP protein is decreased upon PROX1 suppression in two independent experimental settings. IHC analysis by Elsir and colleagues identified up to 80% PROX1+/MAP2+ and around 30% PROX1+/GFAP+, but only 1% to 2% MAP2+/GFAP+ cells in high-grade gliomas, indicating that tumor cells represent either neuronal or glial differentiation pathways in vivo (14). Furthermore, Prox1 depletion has been shown to reduce the number of stem cells and the rate of cell proliferation in intestinal tumor growth (13). According to the present data analysis of glioblastoma samples, PROX1 expression correlates with GFAP and is frequent in the proneural and classical subtypes. Based on these observations, we speculate that PROX1 marks a transitory stage in the evolution of gliomas, possibly coinciding with the choice between neuronal and glial cell developmental paths. Similarly, studies on the role of Prox1 in mouse CNS development indicate a crucial function at this stage of CNS development (10, 43).
PROX1 is part of several complexes involved in development and tumorigenesis including NCoR and HDAC3 (44), as well as PGC-1α and ERRα, which together with BMAL1 play central roles in the transcriptional control of energy homeostasis and biological clocks in the mice liver (45, 46). Here we identified THRAP3 as a novel binding partner for ectopic and endogenous PROX1 in U2-OS and U-343 MGa cells, respectively. Suppression of THRAP3 in the glioblastoma culture U-343 MGa resulted in increased PROX1 and GFAP protein levels, findings that physically and functionally connects PROX1 and THRAP3. Moreover, THRAP3 suppression increased PROX1 and GFAP transcripts. Interestingly, Sono and colleagues recently demonstrated THRAP3 as an interactor with SOX9, where the two proteins together inhibited transcriptional activity during chondrogenesis (39). In the analysis of glioma we found high correlation between PROX1, THRAP3, SOX2, and SOX9. In addition to PROX1, SOX2 binds to the promoter region of THRAP3 and SOX9 (37). Together, these results imply that SOX2 acts upstream of transcriptional complexes composed of PROX1 or SOX9, and that THRAP3 is a modulator of their transcriptional activity. Further, THRAP3 is known as a clock-regulated component of CLOCK–BMAL1 complexes that promotes CLOCK–BMAL1 transcriptional activity and is important for circadian clock function (47). In the GSEA, we identified low PROX1 levels to be associated with increased levels of ERRα target genes, described in ref. 48. Based on our and others’ results this ties PROX1, through THRAP3, to ERRα and BMAL1. We thus speculate that THRAP3 interacts with PROX1 to modulate PROX1 transcriptional function and may take part in regulation of developmental pathways and of the cell-cycle machinery by direct or indirect control of metabolic pathways, for example via changes of estrogen signaling (49).
PROX1 has previously been shown to control CCNE1 by direct binding to its promoter (50). Adding to this, we show that PROX1 affects cyclin protein levels and growth rate in glioblastoma. Berger and colleagues unraveled a critical function for CCNE1 during neural cell fate specification in Drosophila, through interfering with Prospero localization to the nucleus and thereby regulating its function (50, 51). Based on mutant analysis in the CNS of flies, it has been proposed that cyclin E function in specifying a neuronal as opposed to glial cell fate is independent of its established role in G1–S phase transition (51). Furthermore, CCNE1/CDK2 has recently been shown to stabilize pluripotency factors, on a protein level, including SOX2, NANOG, and POU5F1 (OCT4; ref. 38). These observations together with data presented here suggest a feed forward loop where PROX1 enhances CCNE1 transcription, CCNE1 in complex with CDK2 stabilizes SOX2, and finally SOX2 enhances PROX1 transcription. This places PROX1 as a link between cell proliferation and stemness regulation in glioblastoma. The factors that might govern such a loop remain to be studied, but one could speculate that the regulation of such a feed forward loop is out of control in glioblastoma.
The characteristic intratumoral heterogeneity of glioblastoma is suggested to reflect neural development and is likely a key to understanding treatment failure (4). Based on the results in this paper PROX1 can be decreased by the CDK inhibitor CVT-313 and its SOX2 protein destabilizing effect (38). This places CVT-313 as a growth inhibitory drug that would target glioblastomas with high PROX1 expression. In contrast, a CVT-313 induced decrease of PROX1 would be predicted to transition cells from a nonmesenchymal to a mesenchymal subtype and thus potentially promote tumor progression. Further experiments are required to understand how the balance between the growth inhibitory and transition regulation of PROX1 translates in a tumor setting. In extension, the benefit of transitioning cells from a nonmesenchymal to a mesenchymal subtype or vice versa as a therapeutic approach remains to be studied.
In summary, this work highlights the utility of PROX1 as a prognostic marker and adds biological insight to its role in glioblastoma, where it maintains a neural stemness gene expression profile and a proliferative state with high levels of G1-cyclins. Furthermore, loss of PROX1 transitions proneural cells into a mesenchymal subtype, which during tumor evolution would yield heterogeneous cell populations. Finally, CDK2-inhibitors may pose an opportunity to target cell growth of glioblastoma cells with high PROX1 expression.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Conception and design: K.M. Goudarzi, J. Bartek, M. Nistér, M.S. Lindström, D. Hägerstrand
Development of methodology: K.M. Goudarzi
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): K.M. Goudarzi, M. Guo, M. Nistér, M.S. Lindström
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): K.M. Goudarzi, J.A. Espinoza, D. Hägerstrand
Writing, review, and/or revision of the manuscript: K.M. Goudarzi, J.A. Espinoza, J. Bartek, M. Nistér, M.S. Lindström, D. Hägerstrand
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): M. Guo
Study supervision: M. Nistér, M.S. Lindström, D. Hägerstrand
The authors appreciate discussions with the members of the Nistér laboratory, as well as the Bartek group and the SciLifeLab research community who helped make this article possible. This work was supported by KID-funding (to K.M. Goudarzi, 3595/2012), Åke Wiberg stiftelse (to M.S. Lindström, 390316483), Karolinska Institutet (to M.S. Lindström, 1885/12-226 and to D. Hägerstrand, 2014fobi41302), Magnus Bergvall's stiftelse (to M.S. Lindström and to D. Hägerstrand, 2014-00600), King Gustaf V's Jubilee Foundation (to M.S. Lindström, 164102), and the Swedish Research Council (to M.S. Lindström, K2012-99X-21969–01–3). M. Guo was supported by the Chinese Scholarship Council. M. Nistér was supported by grants from the Swedish Cancer Society (CAN 2014/836, contract 160334; CAN 2017/737, contract 170659), the Cancer Society in Stockholm (2015-151213), the Swedish Research Council (K2014-67X-15399-10-4), the Swedish Childhood Cancer Foundation (PR2014-0021; 2017-0002), Karolinska Institutet (2016fobi47658), and the Stockholm County Council (SLL), and J. Bartek by grants from the Swedish Cancer Society (contract 170176) and the Swedish Research Council (K2014-46602-117891-30).
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.