RND1 regulates migration of human glioblastoma stem-like cells according to their anatomical localization and defines a prognostic signature in glioblastoma

INTRODUCTION

Despite combined modality treatment, including surgery and radio-chemotherapy, the prognosis of patients with glioblastoma remains extremely poor [1]. Almost all the patients will die of a relapse in radiation fields or away from the radiation fields, in the brain parenchyma [2]. The failure of initial therapies, which is mainly dependent on tumor heterogeneity, is crucial in glioblastoma [3, 4]. In an attempt to classify diffuse gliomas, subgroups have been defined based on gene expression profiles and genetic and epigenetic alterations [5]. Moreover, high-resolution genome-wide studies have revealed that multiple clones harbouring a variety of genetic alterations coexist within the same tumor [6]. Molecular heterogeneity exists even at cellular level between cells that carry similar genetic alterations and is induced by environmental factors [3, 4]. Expanding single glioblastoma cells into clonal populations demonstrated unique properties including proliferation, differentiation, and different sensitivities to chemotherapeutic drugs [3]. Finally, among patients with proneural tumor, an increased tumor heterogeneity was correlated to a decreased survival [4].

A subset of glioma cells called glioma stem-like cells (GSCs) form heterogeneous glial tumors. They are responsible for the development and the maintenance of tumors [7] and have been proposed to be responsible for tumor recurrences. GSCs display higher resistance to conventional radio-chemotherapy treatments than non GSCs [8, 9]. After treating mice with temozolomide, an alkylating drug currently used in standard glioblastoma treatment, GSCs can still form new tumors [10]. GSCs preferentially reside in perivascular niches where they interact and communicate with tumor associated endothelial cells, via their basement membrane [11]. When orthotopically xenografted, GSCs form tumors that recapitulate the phenotype of patient tumors, notably the ability of glioblastoma cells to infiltrate diffusely [12]. Finally, GSCs are more invasive than their differentiated progeny cells [13].

To invade, glioma cells must initiate dynamic changes in the cytoskeleton organization, notably via integrins [14]. Integrins specifically bind to extracellular matrix proteins, connect their cytoplasmic domain to cytoskeleton and signalling proteins. In this way, integrins control RhoGTPases and actin reorganization that leads to cell migration. Among integrins expressed in glioblastoma cell lines, αvβ3 and αvβ5 integrins are involved in glioma invasion and progression to high-grade glioma [15]. The fibronectin receptor, α5β1 integrin, has also been shown as a promising therapeutic target for high-grade glioma [16]. Besides, elevated expression of α6 integrin from samples of glioblastoma patients is also correlated with a poor patient prognosis [17]. Two studies demonstrate that GSCs overexpressed α6β1 and a3β1 integrins in comparison to non-GSCs [17, 18] suggesting a role of these laminin receptors in GSC invasion. A better understanding of GSC migration toward their microenvironment, notably through integrin involvement, may influence the development of more effective therapies for glioblastoma.

In relation with tumor heterogeneity, clinical data, including ours, demonstrates that relapses in glioblastoma after radio-chemotherapy are more aggressive in tumors contacting the PVZ (PVZ+) than in PVZ– tumors. PVZ+ patients progress quicker and have a decreased overall survival compared to those with tumors not contacting the PVZ (PVZ–) [19–21]. Moreover, we have also recently shown that contact of glioblastoma with PVZ was an independent prognosis factor of shorter progressive free survival after re-irradiation [22]. We hypothesized that the high capacity of PVZ+ tumors to diffusely infiltrate could be explained by a highest migration capacity of GSCs of PVZ+ tumor in comparison to GSCs from the cortical zone (CT). To investigate this, we established GSC paired cultures from CT and PVZ tumor part for two patients diagnosed with glioblastoma. We demonstrated that gene expression hallmarks of GSCs are dependent on their anatomical origin and that PVZ GSCs possess higher migration capacities than CT GSCs. Among migration genes differentially expressed between CT and PVZ GSCs, we showed that overexpression of the RhoGTPase, RND1, dramatically inhibited PVZ GSC migration and, downregulation of RND1 increased CT GSC migration. In silico analysis showed that low expression of RND1 constitutes a bad prognosis factor for glioblastoma patients. Finally, based on signaling pathways activated in patients with low levels of RND1, we identified six genes that define an RND1^low signature that is an independent prognostic factor in glioblastoma.

RESULTS

Characterization and molecular heterogeneity of GSCs derived from CT and PVZ

Glioblastoma samples from two patients were removed from enhanced contrast regions on MRI in the CT and in the PVZ (Figure 1A). We established paired cultures of CT and PVZ GSCs (CT1, PVZ1, CT2, and PVZ2) and analyzed their stemness properties. CT and PVZ cells expressed neural tumor stem cell markers CD133, NESTIN, OLIG1, OLIG2, SOX2 and A2B5 [23] (Figure 1B, Supplementary Figures 1A and 1B). After culture in DMEM-F12 supplemented with 10% FCS (FCSM), CT and PVZ GSCs were able to differentiate into neuronal-like and astrocytic-like cells (Supplementary Figure 1C) and to express differentiation markers (GFAP, TUJ1, MAL and OMG) [23] (Figure 1C). We established with a classic limiting dilution assay that CT and PVZ neurospheres gave rise to secondary neurospheres (Figure 1D and Supplementary Table 1). Both CT and PVZ GSCs had the ability to form diffusely infiltrated tumors when xenografted in nude mice brain (Figure 1E). Altogether, our data shows that CT and PVZ cells derived from our patient samples possess GSC characteristics.

To determine whether PVZ GSCs possess a specific genomic signature, we performed a gene expression microarray. As shown on Figure 1F, different patterns of gene expression have been obtained according to the initial tumor location of the GSCs. We identified 108 genes differentially expressed between CT and PVZ cells (p < 0.001), associated with essential biological functions including cell adhesion, apoptosis, transcription and metabolic process. Up-regulated genes in PVZ GSC included RhoGTPase activating protein 18, transcription factor DP-2 and mannosidase alpha whereas down-regulated genes in PVZ cells included collagen type XI-alpha1, RhoGTPase, RND1 and protocadherin beta3 (Supplementary Table 2). Besides these protein coding genes, CT and PVZ cells differ in the expression of gene expression regulators (antisense RNA, miRNA, long intergenic RNA) and of regulators that guide chemical modifications of others RNAs (small nucleolar RNA…) (Supplementary Table 3). These results demonstrate the molecular heterogeneity of GSCs according to their brain tumor location, notably in the migration processes.

Invasion ability is increased in PVZ GSCs compared to CT GSCs

Cell spreading is the first step of cell invasion. GSCs reside preferentially in perivascular niches and interact with brain blood vessel basement membrane [24]. First, to compare CT and PVZ GSCs migration properties, we performed spreading assays on laminin, fibronectin and vitronectin, three extracellular matrix proteins found in the basement membrane of brain blood vessels and involved in glioma pathogenesis [15–17]. We showed that laminin is a critical extracellular matrix protein for CT and PVZ GSCs and that there was no difference of cell spreading on laminin according to the tumor location (Figure 2A and Supplementary Figure 2A). To further characterize migration properties of these GSCs, we performed directional migration assay in Transwells coated on their undersurface with fibronectin or laminin. None of the GSCs was able to migrate toward fibronectin whereas they all successfully migrated toward laminin (Figure 2B). No significant difference in GSC haptotaxis toward laminin was observed regardless of their initial location in the brain. To sharpen the characterization of GSC migration, we performed time-lapse videomicroscopy of single GSC seeded on laminin. Quantification of single cell migration revealed that PVZ GSCs migrated significantly faster than CT GSCs -as shown by the mean velocity determination- demonstrating differential migration capacities according to the original tumor location (Figure 2C). This result shows that PVZ GSCs have a higher capacity to explore their environment and to scatter into it than CT GSCs. All the GSCs migrated in different directions over the entire surface and no difference in directional persistence was observed according to the tumor location (Figure 2C). Finally, to assess the invasion ability of GSCs, we performed invasion assays in Transwells coating with growth factor reduced Matrigel. Figure 2D shows that PVZ cells invaded more than CT cells. These results showed that laminin is a permissive substrate for CT and PVZ GSC migration and that some GSC invasion properties are dependent on their location in the brain.

RND1 is down-regulated in PVZ GSCs in comparison to CT GSCs

To explain why PVZ GSCs migrate significantly faster than CT GSCs on laminin, we first determined the expression of the main receptors for laminin in glioma cells, ie α6β1, α6β4 and a3β1 integrins [18, 25]. By RT-qPCR, we showed that α6 integrin and β1 integrin mRNAs were strongly expressed in GSCs whereas α3 integrin and β4 integrin mRNAs were weakly expressed (Figure 3A). FACS analyses determined that CT and PVZ GSCs expressed α6 integrin and β1 integrin at high levels on their surfaces, and did not express β4 integrin (Figure 3B). Integrin expression may slightly differ between GSC lines, but independently of the initial tumor location of GSCs. In order to determine whether the integrin involved in laminin GSCs spreading could be different between CT and PVZ cells, we performed cell-spreading assays using functional blocking antibodies. As shown in Figure 3C, inhibition of α6 integrin significantly decreased CT1 and PVZ1 GSC spreading in a similar manner. Besides integrins, we highlighted in our gene expression microarray fourteen genes, including RND1, known to be involved in adhesion/migration that were differentially expressed between CT and PVZ cells (Figure 1F). Of the fourteen significantly altered genes, three were up-regulated and eleven were down-regulated in PVZ GSCs (Figure 3D). RND1 expression is significantly down-regulated in PVZ GSCs in comparison to CT GSCs (pPVZ1/CT1: 0.000019; pPVZ2/CT2: 0.00000497). RND1 is an atypical RhoGTPase that decreases cell adhesion via the inhibition of the formation of actin stress fibers [26]. As RhoGTPases are known regulators of cell migration and have been recently demonstrated as key elements of glioma pathogenesis [27], we focused our study on the role of RND1 in GSC migration. By RT-qPCR, we confirmed the significant down-regulation of RND1 mRNA in PVZ GSCs in comparison to CT GSCs (Figure 3E). To conclude, PVZ GSCs possess both lower levels of RND1 and higher migration properties in comparison to CT GSCs.

RND1 suppresses GSC spreading and migration towards laminin

To test a potential correlation between RND1 low-level expression and PVZ GSCs high migration, we investigated whether high levels of RND1 protein in PVZ cells could decrease their migration ability. PVZ1 cells were transfected with a plasmid encoding a fusion protein of EGFP and RND1 (PVZ1-RND1) or with a plasmid encoding EGFP (PVZ1-EGFP). GFP positive or GFP negative GSCs were selected by FACS (Supplementary Figure 3A). We confirmed RND1 overexpression in GFP positive PVZ1-RND1 GSCs by RT-qPCR (Figure 4A, Supplementary Figure 3B), which did not affect the cell viability (Supplementary Figures 3D and 3E). We then showed that, three hours after seeding on laminin, GFP positive PVZ1-RND1 cells remained round and stayed in suspension whereas GFP negative PVZ1-RND1 cells spread on laminin (Figure 4B) like the control cells (Supplementary Figure 3C). Furthermore, overexpression of RND1 in PVZ1 cells dramatically decreased their ability to migrate toward laminin after 24 h (Figure 4C). We next investigated the consequence of RND1 loss on migration in CT1 cells. CT1 cells were transduced with lentiviral particles expressing a shRNA directed RND1 (CT1 shRND1) or a control shRNA (CT1 shC). RND1 expression was down regulated in stably transduced CT1 shRND1 cells in comparison to CT1 shC (Figure 4D). RND1 loss has no effect on cell viability (Figure 4E) whereas it induces a slight but significant increase in mean velocity (Figure 4F) and in cell spreading (Figure 4G). These results demonstrated that RND1 could suppress spreading and migration abilities of GSCs. Together our data indicates that the differential RND1 expression level between PVZ and CT GSCs may explain, at least in part, their different migration profiles.

Lower expression of RND1 is correlated with a poor prognosis in patients with glioblastoma

As loss of RND1 is involved in GSCs migration from PVZ+ tumors, known to be more aggressive than PVZ– tumors [19–21], we tested whether RND1 gene expression could be correlated with glioblastoma prognosis. First, we analyzed RND1 expression by RT-qPCR in normal brain tissues, in several glioblastoma cell lines and in GSCs established in our laboratory. RND1 expression is significantly down regulated in glioblastoma cells compared to normal tissues (Figure 5A). To pursue our analysis, using gene expression databases in open access, a meta-analysis of RND1 expression revealed a significant down-regulation of RND1 in glioblastoma samples versus normal tissues (p < 0.05, Figure 5B). Using TCGA database, we next examined whether this down-regulation of RND1 was related to the prognosis of glioblastoma patients. Patients with a lower expression of RND1 (i.e <=4.8) showed a worse survival than those with a higher expression of RND1 (HR = 0.59, 95% CI:0.37–0.94, p = 0.028) (Figure 5C). Like patients with a higher expression of RND1, patients with a lower expression of RND1 are divided into the different biological clinical parameters (Supplementary Table 4). Interestingly, patients with a lower expression of RND1 mostly belong to the mesenchymal subtype (Supplementary Table 4). The mesenchymal subtype in glioma is defined by a genomic and transcriptomic profile, notably by higher expression levels of mesenchymal markers such as CHI3L1/YKL40 and MET [28, 29]. To determine whether a lower expression of RND1 correlates with an overexpression of mesenchymal genes, we performed Spearman’s correlation tests using the TCGA database. Supplementary Figure 4 shows that RND1 expression is inversely correlated with common mesenchymal genes [30–33] and specific glioma mesenchymal genes [28, 29, 34]. These results clearly demonstrated that a lower expression of RND1 correlates with a higher expression of mesenchymal genes and, is a factor of poor prognosis for glioblastoma patients.

RND1^low signature is an independent prognostic factor in glioblastoma

To determine which signaling pathways controlled by RND1 might be involved in the prognosis of glioblastoma recurrence and thus survival, we performed a functional enrichment of genes. In patients with low expression of RND1 (RND1^low), we identified thirteen signaling pathways that are activated (Figure 5D and Supplementary Table 5). The most significant being are the “extracellular matrix-receptor interaction” pathway (p = 1.32 e-09), the “focal adhesion” pathway (p = 1.05 e-07) and the “lysosome” pathway (p = 1.72 e-04). This raised the hypothesis that genes from these pathways could be involved in the worse survival prognosis of patients with low RND1 expression. To assess the relationship with overall survival of RND1 and genes from “extracellular matrix-receptor interaction” and “focal adhesion” pathways, a penalized cox regression model with lasso selection was used. We identified six prognostic genes -ITGA5, COL3A1, COL5A1, MET, COL1A2 and, LAMC1- (Supplementary Table 6) that we gathered under the name of RND1^low signature. These genes were all significantly overexpressed in patients with low RND1 expression compared to patients with high RND1 expression (Figure 5E). We then calculated the signature risk score for each patient in the TCGA database and, divided them into a high-risk group and a low-risk group by taking the median value of risk score. The median overall survival in the low-risk group is 17.8 months versus 13.8 months for the high-risk group (Figure 5F, top). The RND1^low signature was significantly associated with overall survival (p < 0.001). To validate our prognostic signature, the training model was applied to glioblastoma patients from REMBRANDT database. Consistent with TCGA results, our RND1^low signature predicts survival of glioblastoma patients (Figure 5F, bottom). Using TCGA, a multivariate cox regression analysis with clinical parameters was carried out to test the strength of RND1^low signature in its ability to predict survival. This analysis showed that the RND1^low signature remains a strong prognostic factor, independently of clinical parameters (Table 1, p = 0.0042). To conclude, we identified an RND1^low signature that is an independent prognostic factor in glioblastoma.

DISCUSSION

The aim of this work was to establish whether GSC migration heterogeneity exists according to the initial location of these cells within the tumor (PVZ+ or PVZ–). By using an original model of GSCs isolated from CT and PVZ, we demonstrated that PVZ GSCs migrated faster and invaded more than CT GSCs and, that their migration may be controlled by RND1. Moreover, we demonstrated that low-expression of RND1 in glioblastoma patient samples was correlated with a worse prognosis for patients. Finally, we identified an RND1^low signature that predicts outcome for glioblastoma patients.

A key finding of our pilot study is that the migration/invasion properties of GSCs depend on their anatomical location. This work will be the foundation for further studies, with a larger number of patients, but shows the increased migration ability of PVZ GSCs compared to CT-derived cells. The increased invasion ability of PVZ GSCs could explain the worse clinical outcome of PVZ+ patients. Only a few studies have previously described different GSC migration abilities according to their location. In fact, a GSC line derived from the PVZ, injected into a mouse brain was shown to invade the corpus callosum and the contralateral hemisphere whereas a GSC line derived from CT was not able to invade these sites [35]. More recently, it has been shown that the GSCs from peritumoral parenchyma are much more invasive than the GSCs from the tumor mass [36]. As CT and PVZ samples come from the same patient, our study illustrates the importance of the intratumoral heterogeneity on tumor behavior. In this study, we investigated and identified genes that may distinguish CT and PVZ GSCs. These genes are involved in cell migration, metabolism, transcription, translation, intracellular traffic, regulation of apoptosis and cell survival (Supplementary Table 2). A recent study analyzed proteins that are significantly altered in PVZ+/PVZ– tissues of glioblastoma patients [37]. In accordance to our results, they found that pathways involved in metabolism (notably in oxidative phosphorylation), extracellular matrix receptor interaction and, migration are affected in PVZ+ glioblastoma. All these functions are hallmarks of cancer and open new lines of research to explain the higher resistance to treatment and the poor overall survival of patients with a glioblastoma contacting the PVZ.

Besides this, we demonstrated in our GSC model that the higher migration ability of PVZ GSCs is potentially associated to lower levels of RND1. Only four recent studies explored the role of RND1 in cell migration. Consistent with our present data, inactivation of RND1 induces the invasion of immortalized breast cells in 3D matrigel and, overexpression of RND1 diminishes lung colonization in mice xenografted with breast cancer cells [38]. On the contrary, overexpression of RND1 in esophageal carcinoma cells promotes their migration [39]. This discrepancy concerning the role of RND1 in invasion is correlated with the difference of RND1 misregulation in these cancers. In fact, RND1 expression is down-regulated in glioblastoma patients and in the most aggressive subtypes of breast cancers [38] but it is up-regulated in esophageal squamous cell carcinoma [39]. Overexpression of RND1 suppresses focal adhesion sites [26] whose formation and turnover are crucial for cell migration [40]. It is known that there is a bi-phasic migration response to cell adhesion since both too weak and too strong adhesion can reduce cell migration [41]. Our functional enrichment of genes in patient tumors revealed that low expression of RND1 in glioblastoma induces an overexpression of focal adhesion proteins like extracellular matrix proteins (COL1A1 and LAMB1); integrins (ITGA5 and ITGB1); actin-binding proteins (FLNA; ACTN1) and vinculin (Supplementary Table 5). Moreover, it was previously shown that overexpression of RND1 in fibroblasts decreases the expression of vinculin [26]. We could hypothesize that misregulation of RND1 expression in glioblastoma cells leads to an optimal formation and turnover of focal adhesion sites and therefore, increased migration.

Recurrence of glioblastoma is caused by the combination of local invasion and therapy resistance. Using a data-driven approach, it has been recently demonstrated that the expression of members of RhoGTPases family is a key marker of glioma progression [27]. In fact, the overexpression of RND3, another member of RND subfamily, enhances the invasion of glioblastoma and is correlated with a poor prognosis. For its parts, in response to a protein complex containing pleiotrophin, secreted by neural precursor from the PVZ, RhoA signaling is activated in glioma cells and induces the migration of glioma cells to the PVZ [42]. In this study, we showed that low levels of RND1, involved in GSC migration, are also related to a decreased overall survival in patients. Moreover, we found that in tumors, a lower expression of RND1 correlates with an overexpression of mesenchymal genes and that tumors with a lower expression of RND1 mostly belong to the mesenchymal subtype. The mesenchymal subtype is associated with a poorer survival than other groups of glioblastoma patients [34, 43]. Interestingly, in breast, depletion of RND1 in immortalized mammary epithelial cells promotes epithelial to mesenchymal transition associated notably with disruption of adherens junctions, downregulation of E-cadherin and up-regulation of fibronectin [38]. Altogether, these results suggest that loss of RND1 in glioblastoma tumors could drive a mesenchymal subtype and hence, could decrease the overall survival in patients.

Based on functional enrichment of genes in glioblastoma patients, we identified six genes whose expression is inversely correlated to RND1 and that predict the survival of glioblastoma patients. The RND1^low signature gathered three qualities: it was discovered from a homogeneous population of glioblastoma patients treated with standard radio-chemotherapy; it involves a short list of genes; and it remains a prognostic factor by itself, independently of clinical parameters. Moreover, the predictive power of the RND1^low signature remains significant for both the training (TCGA) and validation sets (REMBRANDT). Thanks to these qualities, the RND1^low signature could be useful in clinical practice to predict the survival of glioblastoma patients. The RND1^low signature could also lead to clinical application to improve glioblastoma treatment through the targeting of genes involved in this signature. Indeed, ITGA5 and MET were found to be key contributors to the RND1^low signature with their high BSS (Supplementary Table 6). Integrin α5b1 has recently been described as a fine regulator of glioblastoma cell migration [44]. MET and its ligand HGF create an autocrine signaling loop that promotes GSC invasion [45]. In consequence, targeting ITGA5 or MET genes could inhibit the invasive capacity of glioblastoma cells induced by low RND1 expression and especially the one of PVZ+ cells. In conclusion, our study contributes to explain the shorter time to progression of patients with PVZ involvement and highlights RND1 as a gene involved in glioblastoma heterogeneity. Our study suggests that genes establishing the RND1^low signature could be interesting targets for optimizing glioblastoma treatment.

MATERIALS AND METHODS

Extracellular matrix proteins, antibodies and primers

Extracellular matrix proteins, antibodies and primers are respectively depicted in Supplementary Tables 7, 8 and 9. As previously described [38], commercial antibodies did not show sufficient affinity to allow the detection of endogenous RND1.

Tumor samples

Before any therapy, glioblastoma samples were obtained after informed consent from patients admitted to the neurosurgery department at Toulouse University Hospital. Tumors were histologically diagnosed as glioblastoma according to WHO criteria. For patients 1 and 2, two tumor samples were removed from the cortical area (CT1, CT2) and from the periventricular zone (PVZ1 and PVZ2) by utilizing stereotactic image-guided sampling. These patients had a large tumor that was in contact with both CT and PVZ. After mechanical dissociation of tumor tissues, cells were seeded at 37° C in a humid atmosphere of 5% CO₂ in glioblastoma stem cell medium (GSM) composed of DMEM-F12 (Lonza) supplemented with B27 and N2 additives (Invitrogen), EGF (20 ng/mL) and basic FGF (20 ng/mL) (Peprotech). When neurospheres were formed, they were isolated, dissociated with trypsin and cultured as previously described [23]. The percentage of GSCs into the neurosphere was evaluated by flow cytometry [23] (Supplementary Table 10). For other patients (A1, G, I, K, SC1, SC3), only one tumor sample was removed from different brain zones.

Cell culture and limiting dilution assays

Limiting dilution assays were performed on GSCs [23]. For cell differentiation, GSCs were grown in DMEM-F12 supplemented with 10% FCS (FCSM) for two weeks [23].

Human LN18, U87, U138, U251, SF763 and SF767 glioblastoma cells were maintained in DMEM (Lonza) supplemented with 10% FCS.

Cell transfection, cell sorting and cell transduction

Five hundred thousand PVZ1 cells per well were seeded in 6-well plates and then transfected using Fugene HD (Promega) with three μg of p-EGFP-RND1 (Addgene) or of p-EGFP (Clontech). A second transfection was realized seven days after the first transfection to improve gene expression. One week after the second transfection, GFP-positive or GFP-negative GSCs were sorted by FACS. Sorted GSCs were immediately seeded on culture plates to study their ability to migrate or pelleted to quantify RND1 mRNA expression levels.

Twenty five hundred thousand CT1 cells were transduced with lentiviral particles (MOI of 10:1) containing the pLK0.1-neo-CMVtGFP-shRNA plasmid with a sequence directed against RND1 mRNA or a control sequence (Sigma-Aldrich). Five days after transduction, transduced cells were selected with G418.

To establish PVZ1-RND1 cell lines, 25,000 PVZ1 cells were transduced with lentiviral particles (MOI of 10:1) containing the pLX317-puromycin-RND1 (which contains the cDNA of RND1; Sigma-Aldrich) or a control sequence (tGFP; Sigma-Aldrich). Five days after transduction, cells were selected with puromycin.

RT–qPCR and differential expression analysis

RNA from normal human cortex and white matter were obtained from Biochain, Origene, Clontech, and Agilent. Total RNA extraction, RT-qPCR protocol and ΔCt analysis were previously described [23]. Beta2 microglobulin or actin was used as endogenous control in the ΔCt analysis. Results of RT-qPCR are expressed either in 1/ΔCt values -used to express the levels of mRNA without comparison to a standard- or fold induction -used when we compared the levels of mRNA under a given condition to a standard.

After RT on GSCs RNA (at least three samples per cell line) with biotinylated desoxyribonucleotides, cDNA were hybridized on an Affymetrix Human Gene 2.0 ST array. Then, the DNA complexes were revealed by fluorescent streptavidin. Images were analyzed by Command Console and normalized with RMA method (data were normalized per cell line). Two lists of differentially expressed genes (CT1 vs PVZ1 and CT2 vs PVZ2) were established using a criteria based on adjusted p-value cut off of 0.001 and log2 fold change >0.65 or <-0.65. In the heatmap, we only showed the differentially expressed genes that were common to both lists and had the same directional change.

Orthotopic xenograft generation and immunohistochemistry

In accordance with ARRIVE guidelines, the French Institution animal ethics committee approval was obtained for the protocols used on animals. Orthotopic human glioblastoma xenografts were established in 4-6 weeks-old female nude mice (Janvier) with 2.5 × 10⁵ cells as previously described [23]. Each GSC line was xenografted at least in three mice. Mice were sacrificed at the appearance of neurological signs. Immunohistochemistry analysis was performed on the excised brains on paraffin-embedded sections (5 μm) [23].

Immunofluorescence

Immunofluorescence was performed as previously described [46]. For nestin and sox2 stainings, isolated GSCs were seeded on laminin-coated Labtek slides for 24 hours. For β3-tubulin and GFAP stainings, GSC neurospheres were seeded on laminin-coated Labtek slides and were grown in FCSM for five days.

Cell spreading assays

Fifteen thousand GSCs were seeded on pre-coated wells with extracellular matrix proteins at 1.5 μg/cm² and incubated at 37° C. When mentioned, GSCs were pre-incubated or not with 20 μg/ml of function-blocking antibodies for 30 min at 37° C. Three random fields per well from duplicate wells were pictured under a 10× objective. Cells were manually delineated. Cell surface (A), perimeter (P) and circularity (C = 4π(A/P²)) of at least 30 cells per experiment were calculated using the NIS-Elements Advanced Research 3.0 software (Nikon). Cells were classified into two groups: rounded cells and polarized cells as previously described [47].

WST-1 cell viability assays

CT1 and PVZ1 cells were seeded in triplicate into 96-well microplates at a density of 3,000 cells per well and allowed to form neurospheres. At indicated time, the WST-1 reagent (Roche Diagnostics) was applied for 1 h at 37° C. The formazan dye was quantified at 450 nm using a plate reader (FLUOstar Optima, BMG Labtech). The proliferation rate (PR) was calculated in the exponential phase of cells growth during seven 7 consecutive days: PR = (Absorbance_Day+8/Absorbance_Day+1) × 100.

Directional migration assays

Directional migration assays were performed as previously described [46] except for the following: twenty thousand GSCs per well were seeded in Transwells (24 wells, BD Biosciences) pre-coated on their undersurface with 1.5 μg/cm² of fibronectin or laminin and incubated for 24 h at 37° C.

Non-directional migration assays

GSCs (0.75 × 10⁴ cells/cm²) were seeded on laminin-coated wells (2 duplicate wells per condition) and were allowed to migrate for 4 h at 37° C, 5% CO₂. Using an inverted Nikon microscope at 10× magnification, two fields per well were imaged and followed at 3 min intervals with a Coolsnap HQ camera (Photometrics, Tucson, AZ). Manual tracking of the nucleus was performed to follow individual cell migration (at least 30 individual cells per condition per experiment) using NIS-Elements AR 3.0 software as previously described [48]. Directional persistence and mean cell velocity were calculated from time-lapse movies as previously described [48].

Cellular invasion

GSCs were trypsinized, washed with pure DMEM-F12 medium and seeded (100,000 cells/insert) in pure DMEM-F12 on growth factor reduced Matrigel Transwell chambers (8 μm, Corning). The lower wells contain complete GSM. Cells were allowed to migrate for 48 h, non invading cells were removed from upper wells using a cotton swab and invading cells adherent to the bottom of the insert were fixed with 100% methanol and stained with amido black. Invading cells were enumerated by counting the number of cells in 3 distinct fields for each insert under a 10× objective using an optical microscope. Assays were performed in duplicate. The invasion index was calculated with the following formula: (the mean number of invading cells per field for PVZ/ the mean number of invading cells per field for CT) × 100.

Flow cytometry

Flow cytometry was performed as previously described [49]. To specifically determine integrin expression in GSCs, the gated strategy was based on the previously described protocol [23].

RND1 gene expression meta-analysis

For the meta-analysis of RND1 gene expression, we selected the thirteen gene expression datasets (depicted in Supplementary Table 11) comparing glioblastoma samples to normal tissues available on the NextBio research browser (https://www.nextbio.com/) [50].

Survival analysis and functional enrichment

For survival analysis, using the glioblastoma database of TCGA (https://genome-cancer.ucsc.edu/), we focused on patients treated with standard radio-chemotherapy for primary glioblastoma, excluding patients with prior glioma history (n = 184 patients). Minimum p-value approach was used to dichotomize RND1 expression (if RND1>4.8). Overall survival rates were estimated using Kaplan-Meier method and univariate analyses were performed using logrank test. Two-sided p-values of less than 0.05 were considered statistically significant.

Using TCGA, two thousand genes that are the most differentially expressed in patients with low and high RND1 expression were identified by Student’s t-test on all available TCGA patients. The p-values were adjusted using the Benjamini-Hochberg procedure for multiple testing [51]. These genes (p-value cut off of 0.01 and log2 fold change>1.1 or <0.909) were then analyzed for functional enrichment using the Cytoscape (version 3.4.0) plugin ClueGO (version 2.2.5) [52] compared to the KEGG term. A lasso penalized cox regression was used to identify correlations between overall survival and RND1, genes from (KEGG_ECM_RECEPTOR) and (KEGG_FOCAL_ADHESION) pathways [53]. A 10-fold cross validation was realized to select the best penalty parameter lambda. Using a resampling approach, bootstrap selection stability (BSS) was computed for each parameter. From genes selected by the lasso procedure, a risk score prediction was created. It is based on the linear predictor given by the model. This score was then dichotomized by taking the median value of the risk score (threshold = 2.28). Thus, two groups were established (poor versus good prognostic) and corresponded to the signature for this dataset. The high-risk group of the RND1^low signature corresponds to a group with a value superior to the median value of the risk score and, the low-risk group corresponds to a group with a value inferior to the median value of the risk score. To validate our signature, the training model obtained from TCGA was then applied on glioblastoma patients (n = 178) from REMBRANDT database (http://www.betastasis.com/glioma/rembrandt/) and a new risk score was obtained. The risk score was dichotomized like mentioned above, by taking the median value of the new risk score (threshold = 3.15).

Statistical analysis

To compare the average from different experiments, Student’s test was used. Differences were considered statistically significant at p < 0.05; risk of 5%.

To examine whether an underexpression of RND1 correlates with an overexpression of mesenchymal genes, Spearman’s correlation tests were performed on patients treated with standard radio-chemotherapy for primary glioblastoma, excluding patients with prior glioma history.

Abbreviations

FCSM: medium supplemented with 10% of FCS; GSCs: glioblastoma stem-like cells; GSM: glioblastoma stem cell medium; CT: cortical area; PVZ: periventricular zone; PVZ+: tumors contacting the PVZ; PVZ–: tumors not in contact with the PVZ.

ACKNOWLEDGMENTS

The ARC foundation, Fondation de France and ITMO Cancer; ITMO Cancer and Region Midi-Pyrenees/INSERM supported Sabrina Boyrie, Laure Malric and Laetitia Mouly respectively. We thank Judith Martinez Gala, Manon Farce and Anne-Laure Iscache for their kind technical assistance. We also thank the department of pathology of Claudius Regaud Institute; Jean-José Maoret; Lidwine Trouilh and Sophie Lamarre from GenoToul for their kind help in the analysis of microarrays.

CONFLICTS OF INTEREST

ECJM has served as advisory board member for Merck-Serono. The other cited authors, which all concurred with the present article, declare no potential conflicts of interest.

FUNDING

This work was supported by ARC fondation (M2R20100600500), Fondation de France 201300038423, ITMO Cancer (ASC15031BSA and ASC14074BSA), Region Midi-Pyrenees (INSERM R14075BB), INSERM (U1037SE13BA, U1037SE14BA, U1037SE15BA), RITC fondation (RPH11003BBA), the Tumor Heterogeneity and Ecosytem (THE) program (Moglimaging project) and GRICR (R080).

REFERENCES

1. Stupp R, Hegi ME, Mason WP, van den Bent MJ, Taphoorn MJ, Janzer RC, Ludwin SK, Allgeier A, Fisher B, Belanger K, Hau P, Brandes AA, Gijtenbeek J, et al. Effects of radiotherapy with concomitant and adjuvant temozolomide versus radiotherapy alone on survival in glioblastoma in a randomised phase III study: 5-year analysis of the EORTC-NCIC trial. Lancet Oncol. 2009; 10:459–66. https://doi.org/10.1016/S1470-2045(09)70025-7.

2. Konishi Y, Muragaki Y, Iseki H, Mitsuhashi N, Okada Y. Patterns of intracranial glioblastoma recurrence after aggressive surgical resection and adjuvant management: retrospective analysis of 43 cases. Neurol Med Chir (Tokyo). 2012; 52:577–86.

3. Meyer M, Reimand J, Lan X, Head R, Zhu X, Kushida M, Bayani J, Pressey JC, Lionel AC, Clarke ID, Cusimano M, Squire JA, Scherer SW, et al. Single cell-derived clonal analysis of human glioblastoma links functional and genomic heterogeneity. Proc Natl Acad Sci U S A. 2015; 112:851–6. https://doi.org/10.1073/pnas.1320611111.

4. Patel AP, Tirosh I, Trombetta JJ, Shalek AK, Gillespie SM, Wakimoto H, Cahill DP, Nahed BV, Curry WT, Martuza RL, Louis DN, Rozenblatt-Rosen O, Suva ML, et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science. 2014; 344:1396–401. https://doi.org/10.1126/science.1254257.

5. Ceccarelli M, Barthel FP, Malta TM, Sabedot TS, Salama SR, Murray BA, Morozova O, Newton Y, Radenbaugh A, Pagnotta SM, Anjum S, Wang J, Manyam G, et al. Molecular Profiling Reveals Biologically Discrete Subsets and Pathways of Progression in Diffuse Glioma. Cell. 2016; 164:550–63. https://doi.org/10.1016/j.cell.2015.12.028.

6. Sottoriva A, Spiteri I, Piccirillo SG, Touloumis A, Collins VP, Marioni JC, Curtis C, Watts C, Tavare S. Intratumor heterogeneity in human glioblastoma reflects cancer evolutionary dynamics. Proc Natl Acad Sci U S A. 2013; 110:4009–14. https://doi.org/10.1073/pnas.1219747110.

7. Singh SK, Clarke ID, Terasaki M, Bonn VE, Hawkins C, Squire J, Dirks PB. Identification of a cancer stem cell in human brain tumors. Cancer Res. 2003; 63:5821–8.

8. Bao S, Wu Q, McLendon RE, Hao Y, Shi Q, Hjelmeland AB, Dewhirst MW, Bigner DD, Rich JN. Glioma stem cells promote radioresistance by preferential activation of the DNA damage response. Nature. 2006; 444:756–60. https://doi.org/10.1038/nature05236.

9. Liu G, Yuan X, Zeng Z, Tunici P, Ng H, Abdulkadir IR, Lu L, Irvin D, Black KL, Yu JS. Analysis of gene expression and chemoresistance of CD133+ cancer stem cells in glioblastoma. Mol Cancer. 2006; 5:67. https://doi.org/10.1186/1476-4598-5-67.

10. Chen J, Li Y, Yu TS, McKay RM, Burns DK, Kernie SG, Parada LF. A restricted cell population propagates glioblastoma growth after chemotherapy. Nature. 2012; 488:522–6. https://doi.org/10.1038/nature11287.

11. Schonberg DL, Lubelski D, Miller TE, Rich JN. Brain tumor stem cells: Molecular characteristics and their impact on therapy. Mol Aspects Med. 2013; 39:82–101. https://doi.org/10.1016/j.mam.2013.06.004.

12. Singh SK, Hawkins C, Clarke ID, Squire JA, Bayani J, Hide T, Henkelman RM, Cusimano MD, Dirks PB. Identification of human brain tumour initiating cells. Nature. 2004; 432:396–401. https://doi.org/10.1038/nature03128.

13. Qiu B, Zhang D, Tao J, Tie X, Wu A, Wang Y. Human brain glioma stem cells are more invasive than their differentiated progeny cells in vitro. J Clin Neurosci. 2011; 19:130–4. https://doi.org/10.1016/j.jocn.2011.06.014.

14. Guo W, Giancotti FG. Integrin signalling during tumour progression. Nat Rev Mol Cell Biol. 2004; 5:816–26. https://doi.org/10.1038/nrm1490.

15. Uhm JH, Gladson CL, Rao JS. The role of integrins in the malignant phenotype of gliomas. Front Biosci. 1999; 4:D188-99.

16. Janouskova H, Maglott A, Leger DY, Bossert C, Noulet F, Guerin E, Guenot D, Pinel S, Chastagner P, Plenat F, Entz-Werle N, Lehmann-Che J, Godet J, et al. Integrin alpha5beta1 plays a critical role in resistance to temozolomide by interfering with the p53 pathway in high-grade glioma. Cancer Res. 2012; 72:3463–70. https://doi.org/10.1158/0008-5472.CAN-11-4199.

17. Lathia JD, Gallagher J, Heddleston JM, Wang J, Eyler CE, Macswords J, Wu Q, Vasanji A, McLendon RE, Hjelmeland AB, Rich JN. Integrin alpha 6 regulates glioblastoma stem cells. Cell Stem Cell. 2010; 6:421–32.

18. Nakada M, Nambu E, Furuyama N, Yoshida Y, Takino T, Hayashi Y, Sato H, Sai Y, Tsuji T, Miyamoto KI, Hirao A, Hamada JI. Integrin alpha3 is overexpressed in glioma stem-like cells and promotes invasion. Br J Cancer. 2013; 108:2516–24. https://doi.org/10.1038/bjc.2013.218.

19. Adeberg S, Konig L, Bostel T, Harrabi S, Welzel T, Debus J, Combs SE. Glioblastoma recurrence patterns after radiation therapy with regard to the subventricular zone. Int J Radiat Oncol Biol Phys. 2014; 90:886–93. https://doi.org/10.1016/j.ijrobp.2014.07.027.

20. Jafri NF, Clarke JL, Weinberg V, Barani IJ, Cha S. Relationship of glioblastoma multiforme to the subventricular zone is associated with survival. Neuro Oncol. 2013; 15:91–6. https://doi.org/10.1093/neuonc/nos268.

21. Khalifa J. TF, Lusque A, Plas B, Lotterie JA, Benouaich-Amiel A, Uro-Coste E, Lubrano V, Cohen-Jonathan Moyal E. Subventricular zones: new key targets for glioblastoma treatment. Radiation Oncology. 2017; 12:67.

22. Attal J, Chaltiel L, Lubrano V, Sol JC, Lanaspeze C, Vieillevigne L, Latorzeff I, Cohen-Jonathan Moyal E. Subventricular zone involvement at recurrence is a strong predictive factor of outcome following high grade glioma reirradiation. J Neurooncol. 2018; 136:413–9. https://doi.org/10.1007/s11060-017-2669-4.

23. Dahan P, Martinez Gala J, Delmas C, Monferran S, Malric L, Zentkowski D, Lubrano V, Toulas C, Cohen-Jonathan Moyal E, Lemarie A. Ionizing radiations sustain glioblastoma cell dedifferentiation to a stem-like phenotype through survivin: possible involvement in radioresistance. Cell Death Dis. 2014; 5:e1543. https://doi.org/10.1038/cddis.2014.509.

24. Bellail AC, Hunter SB, Brat DJ, Tan C, Van Meir EG. Microregional extracellular matrix heterogeneity in brain modulates glioma cell invasion. Int J Biochem Cell Biol. 2004; 36:1046–69. https://doi.org/10.1016/j.biocel.2004.01.013.

25. Delamarre E, Taboubi S, Mathieu S, Berenguer C, Rigot V, Lissitzky JC, Figarella-Branger D, Ouafik L, Luis J. Expression of integrin alpha6beta1 enhances tumorigenesis in glioma cells. Am J Pathol. 2009; 175:844–55. https://doi.org/10.2353/ajpath.2009.080920.

26. Nobes CD, Lauritzen I, Mattei MG, Paris S, Hall A, Chardin P. A new member of the Rho family, Rnd1, promotes disassembly of actin filament structures and loss of cell adhesion. J Cell Biol. 1998; 141:187–97.

27. Clarke K, Daubon T, Turan N, Soulet F, Mohd Zahari M, Ryan KR, Durant S, He S, Herbert J, Ankers J, Heath JK, Bjerkvig R, Bicknell R, et al. Inference of Low and High-Grade Glioma Gene Regulatory Networks Delineates the Role of Rnd3 in Establishing Multiple Hallmarks of Cancer. PLoS Genet. 2015; 11:e1005325. https://doi.org/10.1371/journal.pgen.1005325.

28. Phillips HS, Kharbanda S, Chen R, Forrest WF, Soriano RH, Wu TD, Misra A, Nigro JM, Colman H, Soroceanu L, Williams PM, Modrusan Z, Feuerstein BG, et al. Molecular subclasses of high-grade glioma predict prognosis, delineate a pattern of disease progression, and resemble stages in neurogenesis. Cancer Cell. 2006; 9:157–73. https://doi.org/10.1016/j.ccr.2006.02.019.

29. Verhaak RG, Hoadley KA, Purdom E, Wang V, Qi Y, Wilkerson MD, Miller CR, Ding L, Golub T, Mesirov JP, Alexe G, Lawrence M, O’Kelly M, et al. Integrated genomic analysis identifies clinically relevant subtypes of glioblastoma characterized by abnormalities in PDGFRA, IDH1, EGFR, and NF1. Cancer Cell. 2010; 17:98–110. https://doi.org/10.1016/j.ccr.2009.12.020.

30. Gibbons DL, Creighton CJ. Pan-cancer survey of epithelial-mesenchymal transition markers across the Cancer Genome Atlas. Dev Dyn. 2018; 247:555–564. https://doi.org/10.1002/dvdy.24485.

31. Smith BN, Bhowmick NA. Role of EMT in Metastasis and Therapy Resistance. J Clin Med. 2016; 5. https://doi.org/10.3390/jcm5020017.

32. Zeisberg M, Neilson EG. Biomarkers for epithelial-mesenchymal transitions. J Clin Invest. 2009; 119:1429–37. https://doi.org/10.1172/JCI36183.

33. Zhao M, Kong L, Liu Y, Qu H. dbEMT: an epithelial-mesenchymal transition associated gene resource. Sci Rep. 2015; 5:11459. https://doi.org/10.1038/srep11459.

34. Wang Q, Hu B, Hu X, Kim H, Squatrito M, Scarpace L, deCarvalho AC, Lyu S, Li P, Li Y, Barthel F, Cho HJ, Lin YH, et al. Tumor Evolution of Glioma-Intrinsic Gene Expression Subtypes Associates with Immunological Changes in the Microenvironment. Cancer Cell. 2017; 32:42–56.e6. https://doi.org/10.1016/j.ccell.2017.06.003.

35. Tchoghandjian A, Baeza-Kallee N, Beclin C, Metellus P, Colin C, Ducray F, Adelaide J, Rougon G, Figarella-Branger D. Cortical and subventricular zone glioblastoma-derived stem-like cells display different molecular profiles and differential in vitro and in vivo properties. Ann Surg Oncol. 2012; 19:S608-19. https://doi.org/10.1245/s10434-011-2093-5.

36. Ruiz-Ontanon P, Orgaz JL, Aldaz B, Elosegui-Artola A, Martino J, Berciano MT, Montero JA, Grande L, Nogueira L, Diaz-Moralli S, Esparis-Ogando A, Vazquez-Barquero A, Lafarga M, et al. Cellular plasticity confers migratory and invasive advantages to a population of glioblastoma-initiating cells that infiltrate peritumoral tissue. Stem Cells. 2013; 31:1075–85. https://doi.org/10.1002/stem.1349.

37. Gollapalli K, Ghantasala S, Kumar S, Srivastava R, Rapole S, Moiyadi A, Epari S, Srivastava S. Subventricular zone involvement in Glioblastoma - A proteomic evaluation and clinicoradiological correlation. Sci Rep. 2017; 7:1449. https://doi.org/10.1038/s41598-017-01202-8.

38. Okada T, Sinha S, Esposito I, Schiavon G, Lopez-Lago MA, Su W, Pratilas CA, Abele C, Hernandez JM, Ohara M, Okada M, Viale A, Heguy A, et al. The Rho GTPase Rnd1 suppresses mammary tumorigenesis and EMT by restraining Ras-MAPK signalling. Nat Cell Biol. 2014; 17:81–94. https://doi.org/10.1038/ncb3082.

39. Xiang G, Yi Y, Weiwei H, Weiming W. RND1 is up-regulated in esophageal squamous cell carcinoma and promotes the growth and migration of cancer cells. Tumour Biol. 2015; 37:773–9. https://doi.org/10.1007/s13277-015-3855-7.

40. O’Neill GM. The coordination between actin filaments and adhesion in mesenchymal migration. Cell Adh Migr. 2009; 3:355–7.

41. Huttenlocher A, Ginsberg MH, Horwitz AF. Modulation of cell migration by integrin-mediated cytoskeletal linkages and ligand-binding affinity. J Cell Biol. 1996; 134:1551–62.

42. Qin EY, Cooper DD, Abbott KL, Lennon J, Nagaraja S, Mackay A, Jones C, Vogel H, Jackson PK, Monje M. Neural Precursor-Derived Pleiotrophin Mediates Subventricular Zone Invasion by Glioma. Cell. 2017; 170:845–59.e19. https://doi.org/10.1016/j.cell.2017.07.016.

43. Bhat KPL, Balasubramaniyan V, Vaillant B, Ezhilarasan R, Hummelink K, Hollingsworth F, Wani K, Heathcock L, James JD, Goodman LD, Conroy S, Long L, Lelic N, et al. Mesenchymal differentiation mediated by NF-kappaB promotes radiation resistance in glioblastoma. Cancer Cell. 2013; 24:331–46. https://doi.org/10.1016/j.ccr.2013.08.001.

44. Blandin AF, Noulet F, Renner G, Mercier MC, Choulier L, Vauchelles R, Ronde P, Carreiras F, Etienne-Selloum N, Vereb G, Lelong-Rebel I, Martin S, Dontenwill M, et al. Glioma cell dispersion is driven by alpha5 integrin-mediated cell-matrix and cell-cell interactions. Cancer Lett. 2016; 376:328–38. https://doi.org/10.1016/j.canlet.2016.04.007.

45. Joo KM, Jin J, Kim E, Ho Kim K, Kim Y, Gu Kang B, Kang YJ, Lathia JD, Cheong KH, Song PH, Kim H, Seol HJ, Kong DS, et al. MET signaling regulates glioblastoma stem cells. Cancer Res. 2012; 72:3828–38. https://doi.org/10.1158/0008-5472.CAN-11-3760.

46. Monferran S, Paupert J, Dauvillier S, Salles B, Muller C. The membrane form of the DNA repair protein Ku interacts at the cell surface with metalloproteinase 9. EMBO J. 2004; 23:3758–68. https://doi.org/10.1038/sj.emboj.7600403.

47. Sadok A, Bourgarel-Rey V, Gattacceca F, Penel C, Lehmann M, Kovacic H. Nox1-dependent superoxide production controls colon adenocarcinoma cell migration. Biochim Biophys Acta. 2008; 1783:23–33. https://doi.org/10.1016/j.bbamcr.2007.10.010.

48. Taboubi S, Garrouste F, Parat F, Pommier G, Faure E, Monferran S, Kovacic H, Lehmann M. Gq-coupled purinergic receptors inhibit insulin-like growth factor-I/phosphoinositide 3-kinase pathway-dependent keratinocyte migration. Mol Biol Cell. 2010; 21:946–55. https://doi.org/10.1091/mbc.E09-06-0497.

49. Monferran S, Skuli N, Delmas C, Favre G, Bonnet J, Cohen-Jonathan-Moyal E, Toulas C. Alphavbeta3 and alphavbeta5 integrins control glioma cell response to ionising radiation through ILK and RhoB. Int J Cancer. 2008; 123:357–64. https://doi.org/10.1002/ijc.23498.

50. Kupershmidt I, Su QJ, Grewal A, Sundaresh S, Halperin I, Flynn J, Shekar M, Wang H, Park J, Cui W, Wall GD, Wisotzkey R, Alag S, et al. Ontology-based meta-analysis of global collections of high-throughput public data. PLoS One. 2010; 5. https://doi.org/10.1371/journal.pone.0013066.

51. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Series B Stat Methodol. 1995; 57:289–300.

52. Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, Fridman WH, Pages F, Trajanoski Z, Galon J. ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics. 2009; 25:1091–3. https://doi.org/10.1093/bioinformatics/btp101.

53. Tibshira R. The lasso method for selection in the cox model. Stat Med. 1997; 16:385–95.