Fluctuation of Rac1 activity is associated with the phenotypic and transcriptional heterogeneity of glioma cells

ABSTRACT Phenotypic heterogeneity of cancer cells is caused not only by genetic and epigenetic alterations but also by stochastic variation of intracellular signaling molecules. Using cells that stably express Förster resonance energy transfer (FRET) biosensors, we show here a correlation between a temporal fluctuation in the activity of Rac1 and the invasive properties of C6 glioma cells. By using long-term time-lapse imaging, we found that Rac1 activity in C6 glioma cells fluctuated over a timescale that was substantially longer than that of the replication cycle. Because the relative level of Rac1 activity in each cell was unaffected by a suspension–adhesion procedure, we were able to sort C6 glioma cells according to the levels of Rac1 activity, yielding Rac1high and Rac1low cells. The Rac1high cells invaded more efficiently than did Rac1low cells in a Matrigel invasion assay. We assessed the transcriptional profiles of Rac1high and Rac1low cells and performed gene ontology analysis. Among the 14 genes that were most associated with the term ‘membrane’ (membrane-related genes) in Rac1high cells, we identified four genes that were associated with glioma invasion and Rac1 activity by using siRNA knockdown experiments. Among the transcription factors upregulated in Rac1high cells, Egr2 was found to positively regulate expression of the four membrane-related invasion-associated genes. The identified signaling network might cause the fluctuations in Rac1 activity and the heterogeneity in the invasive capacity of glioma cells.


INTRODUCTION
Cancer cells that originate from a single cell acquire phenotypic heterogeneity because of genomic instability or heritable epigenetic changes (Lengauer et al., 1998;Shackleton et al., 2009). This heterogeneity is advantageous for the progression of cancer and promotes its highly invasive nature in tissues (Heppner, 1984;Rubin, 1990;Shackleton et al., 2009). Recently, however, it has been reported that the fate and behavior of mammalian cells, including cancer cells, can be determined by stochastic variation in gene expression (Brock et al., 2009). For example, heterogenous signaling patterns in monoclonal cancer cells can generate cells with diverse phenotypes, which have different drug sensitivities (Singh et al., 2010).
A typical example of a cancer that exhibits extensive heterogeneity is glioblastoma, which was previously termed glioblastoma 'multiforme', reflecting its histopathological divergence in size, shape, karyotype, etc. (Louis, 2006). Among the many experimental models of glioblastoma, the C6 glioma cell model has been frequently used to study the invasiveness of glioblastoma cells (Grobben et al., 2002). Glioblastomas formed from C6 glioma cells implanted into syngeneic Wistar rats share many histological hallmarks with human glioblastoma and C6 glioblastoma cells preferentially migrate along neuronal fibers and through the perivascular space, a pattern which resembles the spread of human glioblastoma.
Rho-family GTPases regulate cytoskeletal dynamics and, thereby, affect many cellular processes, including cell polarity, migration, vesicle trafficking and cytokinesis (Etienne-Manneville and Hall, 2002). In cancer cells, Rho-family GTPases play crucial roles in the acquisition of cancer-cellspecific behavior (Sahai and Marshall, 2002). Rac1, for example, accelerates tumorigenesis by regulating apoptosis, cell cycle progression, assembly and disassembly of tight junctions and adherens junctions, cell migration and cell invasion (Mack et al., 2011). Importantly, these pleiotropic functions of Rho-family GTPases have been characterized by comparing cancer cells with non-cancer cells. However, little is known about the heterogeneity and changes in the activity of Rho-family GTPases within cancer cell populations.
The ability to sort cells with respect to a property of interest is essential in order to study the heterogeneity of cancer cells by genetic, epigenetic, biochemical or cell biological analyses. Cell surface markers and cognitive antibodies have been routinely used for this purpose, but the methods used to sort cells depending on the intracellular activity of a signaling molecule are limited. Biosensors based on the principle of Förster resonance energy transfer (FRET) have been widely used to monitor the activity of the signaling molecules (Kiyokawa et al., 2006;Miyawaki, 2011); however, due to a lack of methods for the stable expression of the FRET biosensors, cell sorting with FRET-biosensor-expressing cells has been a difficult task.
Recently, we have reported the development of methods for the stable expression of FRET biosensors in cell lines and transgenic mice Komatsu et al., 2011). By using C6 rat glioma cells that stably expressed a FRET biosensor for Rac1, we found that C6 glioma cells penetrating the brain parenchyma show higher Rac1 and Cdc42 activities and lower RhoA activity than those advancing in the perivascular regions . This observation led us to investigate the mechanism by which the heterogeneity of Rho-family GTPase activity is generated, and the role of the heterogeneity in the invasion of glioma cells. For this purpose, we established a method to sort cells with respect to their levels of Rho-family GTPase activity. Then, by using next generation sequencers, we identified genes whose expression was correlated with Rac1 activity. Using this approach, we show here that fluctuations in the activity of Rac1 are associated with the heterogeneity of glioma invasion.

Distribution of Rac1 activity among C6 glioma cells
We have shown previously that C6 glioma cells that invade at the periphery of a tumor mass in the rat brain, or a three-dimensional (3D) spheroid, exhibit higher Rac1 activity than those trailing such leader cells . We speculated that such a distribution of Rac1 activity among C6 glioma cells might be autonomously generated during cell growth and spreading. To test this hypothesis, we prepared C6 glioma cells that stably express the FRET biosensor Raichu-Rac1, and visualized Rac1 activity on glass-bottomed dishes (Fig. 1A). We detected substantial variation in Rac1 activity, which exhibited a typical normal distribution (Fig. 1B,C). Correlation between the activity of Rac1 and the expression level of the biosensor was not observed (Fig. 1D).

Temporal fluctuations and resilience of Rac1 activity
The normal distribution of Rac1 activity probably reflected the noise of the system (Brock et al., 2009). To study the mechanism underlying the generation of the noise, we used time-lapse imaging to record Rac1 activity in C6 glioma cells for 5 days ( Fig. 2A; supplementary material Movie 1). C6 glioma cells that expressed Raichu-Rac1 were seeded onto a glass-bottomed dish that had a 282-mm-diameter spot, which prevented cells from moving out of the visual field. To maintain the cell density within the optimal range for cell growth, the Raichu-Rac1expressing cells were co-cultured with parental C6 glioma cells. We chose spots that had a single biosensor-expressing cell, and between one and several parental C6 cells for the tracking. Rac1 activity was averaged over the entire cell area and then plotted against time (Fig. 2B). Except for the rapid decline and increase during cell division, Rac1 activity exhibited changes over timescales that were longer than the duration of the cell cycle (.40 hours), which we refer to as fluctuations. Consequently, after 5 days, when the single cells had proliferated to give six to eight cells, Rac1 activity varied substantially among the daughter cells (supplementary material Fig. S1A). Analysis of plots that had been generated using power spectrum did not reveal any periodicity in the fluctuations in Rac1 activity (Fig. 2C). Of note, even after the substantial change during cell division, Rac1 activity returned to a level as before cell division, suggesting that the level of Rac1 activity is maintained by a mechanism that is unaffected by cell division (Fig. 2D). The range of Rac1 activity after 5 days (Fig. 2B) was similar to the range of Rac1 activity observed for the cell population (Fig. 1C); therefore, we concluded that the distribution of Rac1 activity was generated, primarily, by the fluctuations over timescales that occurred over periods longer than the duration of the cell cycle. To understand the biological relevance of the observed distribution of Rac1 activity, we examined the correlation between the cell area and Rac1 activity (Fig. 2E), and found a weak positive correlation. The positive, though weak, correlation between the cell area and Rac1 activity probably reflects the high level of Rac1 activity in lamellipodia (Itoh et al., 2002;Kraynov et al., 2000). We also examined the velocity of the migration (Fig. 2E) but could not observe any clear correlation with Rac1 activity.
Next, to examine the resilience of the level of Rac1 activity to suspension and re-adhesion procedures, we detached the cells by using trypsin and then induced cell adhesion by using trypsin inhibition in situ (Fig. 2F). Trypsinization induced cell rounding and decreased Rac1 activity, whereas trypsin inhibition induced cell adhesion and restored Rac1 activity. Notably, the relative Rac1 activity of each cell to other cells undergoing the suspension-adhesion protocol did not change before, during or after trypsinization (Fig. 2F). This observation agrees with the previous report that the suspension of adherent cells reduces Rac1 activity (del Pozo et al., 2000) and also suggests that the level of Rac1 activity in each cell is maintained by a mechanism that is unaffected by the suspension-adhesion procedure.  To understand the mechanisms underlying the fluctuations in Rac1 activity, we attempted to examine the transcriptomes of C6 glioma cells that had different levels of Rac1 activity. Given that the suspension-adhesion procedure did not alter the relative Rac1 activity of each cell, we sorted C6 glioma cells according to their Rac1 activity by using fluorescence-activated cell sorting (FACS). The FRET to cyan fluorescent protein (CFP) ratio was used as the index of FRET efficiency, as is used in microscopy. The FRET:CFP ratio was independent of the expression level of the biosensor (Fig. 3A) as observed under 2D conditions (Fig. 1D). C6 glioma cells in the highest and lowest decile, with respect to the FRET:CFP ratio, were named the Rac1 high and Rac1 low populations, respectively, and sorted (Fig. 3B).
Because Rac1 activity is closely associated with cell attachment (del Pozo et al., 2000), it was unknown whether the Rac1 activity monitored by the FRET:CFP ratio in FACS reflected the Rac1 activity of cells grown on the culture dishes. In Fig. 2F, we showed that the suspension-adhesion procedure did not alter the relative Rac1 activity of each cell in the cell population; here, we quantified the GTP that was bound to the endogenous Rac1 by using a pulldown assay. The amount of GTP-Rac1 in Rac1 high cells was larger than that in Rac1 low cells (Fig. 3C). Furthermore, we directly measured the GTP:GDP ratio by thin-layer chromatography (TLC) after cell sorting. Both Rac1 high and Rac1 low cells were plated on the culture dishes and radiolabeled with 32 P for 2 hours, followed by TLC analysis to measure the GTP:GDP ratio on the biosensor (Fig. 2D). Although the difference in the GTP:GDP ratio between Rac1 high and Rac1 low was smaller than the difference between wild-type and GTPase-deficient mutant Rac1 proteins, the GTP:GDP ratio of Rac1 high cells was constantly higher than that of Rac1 low cells, providing biochemical evidence that the FRET-based cell sorting reflects the Rac1 activity level. In addition, we found that cells in the G2/M phase were enriched in the Rac1 low population (supplementary material Fig. S1B), which agrees with the observation that Rac1 activity was transiently reduced during cell division (Fig. 2B).
We then examined the invasion of Rac1 high and Rac1 low cells into Matrigel by using a Boyden chamber assay. Cells that had reached the lower side of the filter were counted after 22 hours. Although the efficiency of invasion varied in each experiment, we constantly observed that Rac1 high cells invaded into Matrigel significantly faster than did Rac1 low cells (Fig. 3E). This observation agrees with the findings of our previous 3D spheroid assay, which showed that cells with higher Rac1 activity invaded into Matrigel first and guided cells with lower Rac1 activity ( Fig. 1B)  .
Finally, to investigate further our hypothesis that the fluctuations in the activity of Rac1 caused the observed distribution of Rac1 activity, the Rac1 high and Rac1 low cells were cultured for up to 9 days and re-analyzed by FACS. On the first day, the distribution of Rac1 activity within the sorted populations remained discrete, but after one week the distribution of Rac1 activity within each population was identical, supporting our hypothesis (Fig. 3F). We performed similar experiments with C6 glioma cells expressing a FRET biosensor for Cdc42. In agreement with the previous finding that glioma cells invading at the front exhibit high Rac1 activity and high Cdc42 activity , a similar result was obtained using Cdc42 high and Cdc42 low cells (Fig. 3F). Furthermore, we retrospectively analyzed long-term time-lapse images (supplementary material Fig. S1C), and found that, in agreement with the FACS data, cells in the highest and lowest decile, with respect to the FRET to CFP ratio, displayed similar profiles at 50 hours.

Transcriptional signatures of Rac1 high and Rac1 low cells
To investigate the association of the variation in Rac1 activity with transcriptional signatures, RNA sequencing (RNA-Seq) analysis was performed on C6 glioma cells that had been sorted by using FACS. mRNA was isolated from the Rac1 high and Rac1 low populations of C6 glioma cells and sequenced. The difference in the expression of genes between Rac1 high and Rac1 low cells was plotted against the calculated average expression (Fig. 4A). Similar analyses were performed for the cells expressing the biosensors for Cdc42 and RhoA to (A) C6 glioma cells expressing Raichu-Rac1 were analyzed by using FACS. The Rac1 activity (measured by using the FRET to CFP ratio, FRET/CFP) did not correlate with the FRET biosensor concentration (YFP) in each cell. (B) The top and bottom 10% of cells were sorted to obtain Rac1 high (red) and Rac1 low (blue) populations, respectively, (left panel). Small fractions of Rac1 high and Rac1 low cells were reanalyzed after sorting (right panel). (C) Rac1 high and Rac1 low cells were collected and analyzed by pulldown assay (left panel). Cells expressing CFP-Rac1V12 (constitutively active) and CFP-Rac1N17 (dominant negative) were used as positive and negative controls, respectively. White and black arrows indicate CFP-Rac1 and endogenous Rac1, respectively, quantification is shown in the right panel.
(D) Rac1 high and Rac1 low cells were plated on dishes and cultured for 3 hours. Cells were labeled with 32 P for 2 hours, lysed and immunoprecipitated by using an antibody against GFP, followed by TLC analysis (left panel) and quantification of GTP and GDP bound to the FRET biosensor (right panel). P-values were derived by using two-tailed Student's t-test. (E) Rac1 high and Rac1 low cells were used for the Matrigel invasion assay. (F) Rac1 high , Rac1 low cells, Cdc42 high and Cdc42 low cells were plated on dishes, cultured for the indicated periods and then the activity of Rac1 and Cdc42, respectively, was re-analyzed by FACS.
characterize the Rac1 high and Rac1 low populations (Fig. 4B). The difference in the expression of genes between Rac1 high and Rac1 low cells was positively correlated with that of Cdc42 high and Cdc42 low cells. However, there was no correlation between the expression of genes in Rac1 high and Rac1 low cells with RhoA high and RhoA low cells, or those of Cdc42 high and Cdc42 low and RhoA high and RhoA low cells. Again, this is in agreement with our previous observation that both the activity of Rac1 and Cdc42 were high in glioma cells migrating at the leading edge in rat brains and in 3D Matrigel experiments . To identify genes that were expressed differently, we used the weighted average difference method (WAD). The WAD method identified 713 differentially expressed genes using a cutoff of the top 5% of ranked genes. Gene ontology analysis, based on terms associated with biological processes, showed that the Rac1 high phenotype is associated most with the G-protein-coupled receptor signaling pathway, followed by cell-matrix adhesion and then the electron transport chain (Fig. 4C). The Rac1 low phenotype is associated with terms pertaining to cell division, the cell cycle and mitosis. In addition, analysis of terms relating to cellular components showed that genes that are associated with the respiratory chain, focal adhesion and mitochondrial respiratory chain complex I were upregulated in the Rac1 high population, and that the expression of genes relating to the cytoplasm and The relationship between the average expression [log 2 (Rac1 high 6Rac1 low )/2] and expression difference between Rac1 high and Rac1 low cells [log 2 (Rac1 high /Rac1 low )] is shown in a log ratio against abundance (M-A) plot. The WAD method identified 713 differentially expressed genes using cutoffs of the top 5% ranked genes. Genes enriched in Rac1 high and Rac1 low cells are depicted with pink and blue dots, respectively. The top 14 genes for the cellular component term 'membrane' are marked in orange, except for Pstpip2, Freq/NCS1, MMP15 and Tsn17, which are shown with red dots. Elmo1 and Egr2 are shown as green dots. (B) RNA-Seq analysis was similarly performed for the Cdc42 high , Cdc42 low , RhoA high and RhoA low cell populations. Scatter plots of the expression differences are shown. Correlation coefficient, r, between measures was calculated by using R software. (C,D) Gene ontology analysis with biological process terms (C) or cellular component terms (D) is shown for Rac1 high and Rac1 low cells (left and right panels, respectively). P-values were calculated by Pearson product-moment correlation coefficients.
nucleus, and that are integral to membrane, was increased in the Rac1 low population (Fig. 4D).

Identification of genes that regulate C6 glioma cell invasion
Based on the RNA-Seq data, we attempted to identify genes that regulate glioma invasion. For this purpose, we focused on the 23 genes that were most upregulated in the Rac1 high population and related to the cell component term 'membrane', which we term 'membrane-related genes' (supplementary material Table S1). Notably, 18 genes out of these 23 genes were upregulated in Cdc42 high cells in comparison with Cdc42 low cells, strongly suggesting that there is large overlap between gene expression in Rac1 high cells and Cdc42 high cells. The difference in gene expression between Rac1 high cells and Rac1 low cells was first confirmed by using quantitative (q)PCR, except for Ecop, for which we failed to prepare specific primers. Next, to characterize the membrane-related genes that were upregulated in Rac1 high cells, we knocked down the ten genes that were most upregulated in the Rac1 high population by using three different siRNAs, except for cathepsin L1 (Ctsl), against which we failed to prepare three effective siRNAs. From the remaining 12 genes, we arbitrarily chose and knocked down Ntrk2, Freq (also known as Ncs1), Il1rap and Pstpip2. C6 glioma cells were then examined for their invasive potential by using the Matrigel invasion assay. Of the 14 membrane-related genes, knockdown of matrix metalloproteinase 15 (MMP15, also known as MT2-MMP), tetraspanin 17 (TSN17), Pstpip2 and Freq, which we call membrane-related invasion-associated genes, significantly inhibited C6 glioma cell invasion (Fig. 5A). Except for MMP15, knockdown of the membrane-related invasionassociated genes suppressed Rac1 activity, suggesting that TSN17, Pstpip2 and Freq promote invasion through activation of Rac1 (Fig. 5B). Knockdown of the membrane-related but invasion-irrelevant genes Lgals3 and Rgs2 did not affect Rac1 activity. Notably, knockdown of TSN17 and Freq, but not Pstpip2, caused rounding of the cell shape (Fig. 5C). These two genes might be associated with Rac1-mediated membrane protrusion.
We next examined transcription factors that were upregulated in Rac1 high cells (supplementary material Table S2) and found that expression of Egr2 was reproducibly increased in Rac1 high cells and Cdc42 high cells. Similar to this, Elmo1 and PRex1, which are positive regulators of Rac1 activity, were found to be upregulated in Rac1 high cells (supplementary material Table S3). Knockdown of Egr2 and Elmo1, but not PRex1, suppressed C6 cell invasion and decreased Rac1 activity (Fig. 5D,E). Unlike the knockdown of the membrane-related genes TSN17 and Freq, knockdown of Egr2 or Elmo1 decreased Rac1 activity without affecting the cell shape (Fig. 5F).
Thus, we identified the genes of four membrane-related proteins, a transcription factor and a Rac1 activator as invasionassociated genes that are upregulated in Rac1 high cells. Notably, knockdown of Egr2, Elmo1, Pstpip2, Freq and TSN17 not only decreased the average Rac1 activity but also suppressed the fluctuations in the activity of Rac1 (supplementary material Fig.  S2), suggesting that the fluctuations in the activity of Rac1 might be associated with the basal level of Rac1 activity.

Hierarchy of the invasion-associated genes upregulated in Rac1 high cells
To delineate the signaling network of the invasion-associated genes that were upregulated in Rac1 high cells, we first transiently activated Rac1 by using the rapamycin-induced Rac1 activation system (Yagi et al., 2012). Upon rapamycin-induced membrane translocation of Tiam1, a guanine-nucleotide-exchange factor for Rac1, Rac1 activation was clearly detected by both FACS (Fig. 6A) and pulldown assay (Fig. 6B). Of the genes tested, only expression of Egr2 was rapidly and transiently induced (supplementary material Fig. S3). We next knocked down each gene and quantified the expression of the other invasionassociated genes by using qPCR (supplementary material Fig.  S4). The comprehensive knockdown experiments revealed intriguing features of the signaling network that regulates Rac1 activity (Fig. 6C). First, the membrane-related invasionassociated genes, Pstpip2, TSN17, Freq and MMP15 were clustered. Second, and more importantly, we could determine the hierarchy of the genes by assuming that the knockdown of a gene decreases the expression of the downstream genes and increases the expression of the upstream gene(s) by a negativefeedback loop. For example, Egr2 knockdown decreased expression of all four membrane-related invasion-associated genes but not the invasion-irrelevant genes, such as Rgs2 or Elmo1. By contrast, knockdown of the membrane-related invasion-associated genes increased Egr2 expression, suggesting that there is a negative-feedback loop to Egr2. Because the effect of Egr2 knockdown on the expression of Pstpip2 was not significant, this gene could be part of a different signaling pathway. From these data, we propose a model of the gene network that regulates Rac1 activity and the invasive capacity of C6 glioma cells (Fig. 6D).

DISCUSSION
The phenotypic heterogeneity of cancer cell populations is caused by genetic, epigenetic and non-genetic mechanisms. The nongenetic mechanisms that cause the variations in gene expression include transcriptional and translational noises (Brock et al., 2009). Although the precise nature of such noise remains largely elusive, we can speculate that the variation in gene expression reflects the intracellular signaling activities. Here, we have established a technique to sort cells according to the activities of intracellular signaling molecules and to examine the effect of the variation in signaling molecule activity on the biological or transcriptional heterogeneity of cancer cells.

FRET-based cell sorting
This technique is based on two assumptions: first, that the activity of the molecule of interest is maintained during FACS; second, that the transcriptome is not substantially perturbed during FACS. Given that it has been reported previously that the suspension of adherent cells reduces their Rac1 activity (del Pozo et al., 2000), it was of concern that the cell preparation process -i.e. trypsinization and suspension of adherent C6 glioma cellsmight mask the intercellular variation in the activities of the Rhofamily GTPases that are observed under both 2D and 3D conditions. Contrary to our expectation, the intercellular variation in Rac1 activity was reproduced after the cell preparation process. Re-analysis of Rac1 activity by FACS, TLC or pulldown assay demonstrated that Rac1 activity was conserved in both the Rac1 high and Rac1 low cell populations (Fig. 3A,B). Time-lapse imaging confirmed that the relative Rac1 activity of each cell was maintained before and after cell division or the suspensionadhesion procedure (Fig. 2D). These observations suggest that the activity of Rac1 in a single cell is dependent on both basal and external cues. The activation of external cues, which include growth factors and signaling mediated by integrin, rapidly changes Rac1 activity (Heasman and Ridley, 2008). Meanwhile, the basal Rac1 activity is determined by the status of intrinsic signals, which are unaffected by external cues and fluctuate over longer timescales.
The second assumption that requires further consideration is the effect of cell sorting on the transcriptome. The ontology analysis of genes that were upregulated in Rac1 high cells showed a close correlation with terms that are associated with biological processes that are linked to the function of Rac1 (Fig. 4C). The first and second most highly scored terms were the GPCR pathway and cell-matrix adhesion, respectively. Both of these are related to cell migration, with which the functions of Rac1 are most often associated (Sahai and Marshall, 2002;Sander and Collard, 1999). Another major function of Rac1 is the regulation of the NADH-mediated production of reactive oxygen species (Heasman and Ridley, 2008); therefore, it is not surprising that the electron transport chain was scored third. Furthermore, in agreement with the finding that Rac1 activity rapidly reduces during cell division (Yoshizaki et al., 2003), the first to third scores of genes upregulated in Rac1 low cells were assigned to the pathways of cell division, cell cycle and mitosis, respectively (Fig. 4C). Furthermore, among the 23 genes that were upregulated in the Rac1 high population and classified under the cell component term 'membrane', 13 genes are known to be involved in cancer cell invasion (supplementary material Table  S1). These observations support our assumption that the transcriptional profiles are reasonably conserved during FACS.
Genes associated with the Rac1 high phenotype Among the 14 genes that were classified under the cell component term 'membrane' and were upregulated in Rac1 high cells, knockdown of four genes inhibited C6 glioma cell invasion in the Matrigel invasion assay. Previous reports have indicated or suggested that proteins encoded by the four genes are more or less associated with invasion of cancer cells. MMP15 is expressed predominantly in glioblastoma (Lampert et al., 1998;Nakada et al., 1999), suggesting that MMP15 plays a major role in the degradation of the extracellular matrix during glioma invasion. TSN17 has been recently shown to regulate ADAM10, which has been shown to be involved in cancer progression (Dornier et al., 2012;Haining et al., 2012;Mochizuki and Okada, 2007). Pstpip2 regulates F-actin bundling and enhances the formation of filopodia (Chitu et al., 2005), which strongly points to a role of this protein in glioma invasion. Freq/NCS1 is a Ca 2+ -binding protein that is expressed predominantly in neuronal cells (Dason et al., 2012;Nakamura et al., 2006). In primary cultured adult cortical neurons, overexpression of NCS1 induces neurite sprouting; however, the role of NCS1 in glioma invasion has not been determined. Fig. 6. Identification of a signaling network comprising the genes that are upregulated in the Rac1 high cell population and associated with invasion. (A,B) C6 glioma cells expressing Raichu-Rac1 alone or Raichu-Rac1 plus plasma-membrane-targeted Lyn-FKBP12-rapamycin-binding domain (LDR) and FK506-binding protein (FKBP) fused with Tiam1 were stimulated with 10 mM rapamycin (solid line) or the solvent DMSO (dashed line) for 30 minutes. Rac1 activity was examined by using FACS (A) or pulldown assay (B). C6 glioma cells expressing CFP-Rac1N17 or CFP-Rac1V12 were used as negative and positive controls, respectively. White and black arrows indicate CFP-Rac1 and endogenous Rac1, respectively. Densitometry for GTP-bound Rac1 was normalized to the amount of the total Rac1 (right panel). (C) The genes listed in the left column were knocked down as described in Fig. 5, or Rac1 was activated as in A. qPCR was used to analyse the levels of the genes in the top row. Fold changes relative to the control siRNA-transfected cells are shown in the log2 scale. The genes were clustered by using the nearest neighbor method. The data on the invasion and Rac1 activity that is shown in Fig. 5 are also included.

Black cells indicate knockdown efficiency, yellow cells indicate upregulation (.50%), blue cells indicate downregulation (,250%), white cells indicate others. (D) A proposed model of Rac1 activity regulation in C6 glioma cells.
Genes responsive to the activation of Rac1 or Cdc42 have been identified by overexpressing constitutively active Q61L mutants of Rac1 or Cdc42 in NIH 3T3 cells (Teramoto et al., 2003). There are some similarities between this previous work and our present study. First, the expression profile of Rac1QL-expressing cells resembles that of Cdc42QL-expressing cells in the previous studies. We also found that the expression profile of Rac1 high cells resembles that of Cdc42 high cells ( Fig. 4B; supplementary material Table S1). Second, in cells expressing Rac1QL or Cdc42QL, genes related to the extracellular matrix and cell adhesion are upregulated and genes related to the cell cycle are suppressed, as in Rac1 high cells (Fig. 4C). However, there are also some discrepancies. For example, the collagen a1 chain precursor was upregulated 3.1-fold in Rac1QL-expressing cells but was suppressed to 56% in Rac1 high cells. This difference might have been caused by the lack of GTPase activity in cells expressing the Rac1QL mutant because the cycling between the GDP-bound and GTP-bound forms has been shown to play an important role in cell migration (Parrini et al., 2011). In another study, different expression levels of Rac1 have been used in colorectal cancer cells to identify the target genes of Rac1 by using microarray analysis (Gómez del Pulgar et al., 2007); however, we could not find any similarity to our data. In another study, C6 rat glioma cells were selected both in vitro and in vivo for phenotypes of high and low migratory capacity (Tatenhorst et al., 2005). By using microarray analysis, that study found 31 genes were differentially expressed in association with the migratory phenotypes. We could not detect a substantial resemblance between the gene expression profiles of that study and our present findings. Thus, the constitutive activation (or suppression) and intrinsic fluctuations in Rac1 activity might cause different transcriptional phenotypes. Alternatively, the effect of Rac1 on transcriptional profiles might be specific to the cell type. In any event, different approaches led to the identification of various genes related to glioma invasion. Further analyses will be required to find the cause of such divergence.
Hierarchy of invasion-associated genes that are upregulated in the Rac1 high population The comprehensive knockdown experiments strongly pointed to the role of Egr2 as a master regulator of C6 glioma invasion. Knockdown of Egr2 suppressed the expression of four membrane-related invasion-associated genes. By contrast, knockdown of the four membrane-related invasion-associated genes or Elmo1 increased the expression of Egr2, which suggests that there are negative-feedback loops between the invasive phenotype and Egr2 expression. Microarray analyses have revealed upregulation of Egr2 in metastatic squamous cell carcinomas (Kim et al., 2006;Liu et al., 2008). Furthermore, in fibroblasts infected with Kaposi sarcoma-associated herpesvirus, Egr2 induces MMPs and extracellular matrix metalloproteinase inducer (Emmprin) (Dai et al., 2012). These observations strongly support the proposal that Egr2 is a key regulator of glioma invasion.

Origin of the heterogeneity of Rac1 activity
What causes the heterogeneity of Rac1 activity among C6 glioma cells? The 5-day time-lapse imaging revealed that non-genetic fluctuations in the activity of Rac1, which occurred over timescales of more than 40 hours, caused the observed distribution of Rac1 activity ( Fig. 2A). This conclusion was also supported by the observation that culture of the isolated Rac1 high or Rac1 low cell population restored the original distribution of Rac1 within 1 week ( Fig. 3F; supplementary material Fig. S1C). Notably, our conclusion agrees with the variation of protein levels that has been seen previously in human H1299 lung carcinoma cells (Sigal et al., 2006), in which the expression levels of proteins have been shown to fluctuate over a timescale of more than 40 hours. By using knockdown experiments against genes upregulated in the Rac1 high population, we identified a gene network that regulates Rac1 activity (Fig. 6D). This network comprises both positive-and negativefeedback loops, which are sufficient to cause oscillation of the signaling network. Although we have not been able to confirm that the variation in Rac1 activity in vivo is also driven by the same mechanism, fluctuations in gene expression over prolonged periods of time, and the resulting fluctuations in Rac1 activity, could generate glioma cells that have different invasive capacities.

MATERIALS AND METHODS
Biosensors and cell lines C6 rat glioma cells were obtained from the American Type Culture Collection and cultured in Dulbecco's modified Eagle's medium containing 10% fetal bovine serum (FBS). The FRET biosensors for Rac1, Cdc42 and RhoA were Raichu-Rac1, Raichu-Cdc42 and Raichu-RhoA, respectively, and have been described previously (Itoh et al., 2002;Yoshizaki et al., 2003). To establish stable cell lines that expressed Raichu biosensors, we used two approaches: (1) we replaced CFP with teal fluorescent protein (TFP) and delivered the expression cassettes into C6 glioma cells by using a retroviral vector as described previously ; (2) we used piggyBac transposon-mediated gene transfer, which, more recently, has been used to stably express Raichu biosensors with higher sensitivity (Komatsu et al., 2011;Yusa et al., 2009). The cells were then single-cell cloned before further experiments unless described otherwise.

Time-lapse FRET imaging
FRET images were obtained and processed using the same conditions and procedures as previously reported (Aoki and Matsuda, 2009). Cells were plated on 35-mm-diameter glass-bottomed dishes (AGC Techno Glass, Shizuoka, Japan) or micro-patterned glass-bottomed dishes (CytoGraph; Dai Nippon Printing, Tokyo, Japan). Cells were imaged at 37˚C under 5% CO 2 by using an inverted microscope (IX81; Olympus, Tokyo, Japan) equipped with a 640 objective lens (UAPO/NA 1.35, Olympus), a 640 objective lens (UPLSAPO/NA 0.95, Olympus), a 660 objective lens (PlanApoPH/NA 1.40, Olympus), a cooled CCD camera (Cool SNAP-HQ or Cool SNAP-K4; Roper Scientific, Tucson, AZ), an LED illumination system (CoolLED precisExcite; Molecular Devices, Sunnyvale, CA), an IX2-ZDC laser-based autofocusing system (Olympus) and an MD-XY30100T-Meta automatically programmable xy stage (SIGMA KOKI, Tokyo, Japan). The following filters were used for the dual-emission imaging studies and were obtained from Omega Optical (Brattleboro, VT): an XF1071 (440AF21) excitation filter, an XF2034 (455DRLP) dichroic mirror, and two emission filters, XF3075 (480AF30) for CFP and XF3079 (535AF26) for yellow fluorescent protein (YFP). After background subtraction, FRET:CFP ratio images were created by using the MetaMorph software (Universal Imaging, West Chester, PA) and were represented by an intensity-modulated display mode. In the intensity-modulated display mode, eight colors from red through to blue are used to represent the FRET:CFP ratio, with the intensity of each color indicating the mean intensity of FRET and CFP. For the quantification, the FRET and CFP intensities were averaged over the whole cell area, and the results were exported to Excel software (Microsoft Corporation, Redmond, WA).
FRET-based cell sorting C6 glioma cells expressing Raichu-Rac1 were trypsinized, resuspended in PBS containing 3% FBS and analyzed and/or sorted with a FACSAria machine (Becton Dickinson, Franklin Lakes, NJ). We used the following combinations of lasers and emission filters for the detection of fluorescence from the biosensor: for the donor fluorescence of TFP and CFP, a 407-nm laser and a 480AF30 filter (Omega Optical); for the sensitized FRET from YFP, a 407-nm laser and a 535AF26 filter (Omega Optical); and for the acceptor fluorescence of YFP, a 475-nm laser and a 535AF26 filter (Omega Optical). Cells were first gated for size and granularity to exclude cell debris and aggregates. For cell sorting, C6 glioma cells in the highest and lowest decile with respect to the FRET to CFP (or TFP) ratios were sorted as Rac1 high and Rac1 low populations, respectively, into DMEM containing 10% FBS. Small fractions of Rac1 high and Rac1 low cells were re-analyzed for validation, and the remaining cells were snap-frozen and stored at 280˚C until RNA extraction. Detailed data analysis was performed using FlowJo. 7.6 software (Tree Star, Ashland, OR).

TLC of guanine nucleotides bound to GTPases
Guanine nucleotides bound to Raichu biosensors or GFP-tagged Rac1 proteins were quantified as described previously (Gotoh et al., 1995). Briefly, cells were sorted by FACS and plated on six-well dishes. After 3 hours, cells were metabolically labeled with 32 P for 2 hours and lysed with lysis buffer. The cell lysates were clarified by centrifugation and used to immunoprecipitate Raichu biosensors or GFP-tagged Rac1 by using an antiserum against GFP and Protein-A-Sepharose. Guanine nucleotides bound to the immuoprecipitates were separated by TLC and quantified with a BAS-1000 image analyzer.

Rac1 pulldown analysis
The Rac1 pulldown assay was performed according to the manufacturer's protocol (Cytoskelton, Denver, CO).

RNA extraction
Total RNA was isolated by using a Qiagen RNeasy Micro Kit (Qiagen, Hilden, Germany) or a Qiagen RNeasy Mini Kit according to the manufacturer's protocol. RNA preparations were confirmed to be free of proteins using a NanoDrop ND-1000 instrument (Thermo Fisher Scientific, Waltham, MA), and the integrity of these measurements was confirmed using a 2100 BioAnalyzer (Agilent Technologies, South Queensferry, UK). RNA that had an RNA integrity number of §8.6 was used for RNA-Seq.

Library preparation and sequencing
Total RNA was poly(A)-selected using poly(T) Dynabeads (Invitrogen, San Diego, CA). Sequencing libraries were prepared according to Illumina's mRNA-Seq protocol and sequenced at the Omics Science Center (OSC) RIKEN Yokohama Institute. Two independent libraries were analyzed for each data set. Sequence-read data have been submitted to the Sequence Read Archive at the DNA Data Bank of Japan (submission no. DRA000605).

Mapping and processing of RNA-seq reads
The reads of each dataset were aligned to the rat reference genome (rn4, Nov. 2004, version 3.4) using TopHat v1.3.0 (Trapnell et al., 2009). The resulting sequence alignment and map files in the BAM format were analyzed by using Cufflinks version 0.8.0 (Trapnell et al., 2010) to compute fragments per kilobase of transcript per million mapped reads (FPKM). Genomic annotations were obtained from Ensembl in gene transfer format (GTF). We used only reads mapped to 20 or fewer sites on the genome. The WAD method (Kadota et al., 2008) was then performed on the data of pairs of cells to generate expression differences. Differentially expressed genes were filtered for a WAD ranking cutoff of the top 5.0%. Gene Ontology (GO) annotations were used to assign biological functions to genes included in this study (Ashburner et al., 2000).

3D Matrigel invasion assay
3D Matrigel invasion assays were performed with a BD BioCoat Matrigel Invasion Chamber (Becton Drive, Franklin Lakes, NJ) according to the manufacturer's protocol. Briefly, 2610 4 cells were seeded on the membrane with or without Matrigel pre-coating. After 22 hours, cells were fixed, stained for nuclei with propidium iodide and imaged by using an epifluorescence microscope. The number of nuclei was counted with MetaMorph software (Universal Imaging). Data are expressed as the percentage invasion through the Matrigel matrix and membrane relative to the migration through the control membrane.

3D spheroid imaging
Organotypic culture was prepared as described previously (Gaggioli et al., 2007). In a 12-well plate coated with poly-(2-hydroxyethyl methacrylate) (Sigma-Aldrich, St Louis, MO) and containing 1 ml serum-free CO 2 -independent medium (Invitrogen), 10 6 cells were cultured overnight with slow agitation to form small aggregates. The aggregates were embedded in 6 mg/ml Matrigel, maintained in complete medium and observed under a two-photon microscope or a confocal microscope for up to 18 hours in an incubation chamber.

Quantitative PCR
RNA was reverse-transcribed by using a High Capacity cDNA Reverse Transcription kit (Applied Biosystems, Foster City, CA) according to the manufacturer's protocol. Then, the expression levels of each gene and GAPDH, used as a standard, were analyzed by using the Power SYBR Green PCR Master Mix (Applied Biosystems) and the ABI PRISM7300 Sequence Detection System (Applied Biosystems). The sequences of primers used for qPCR are shown in supplementary material Table S4.

siRNA knockdown experiments
Stealth RNAi Negative Control Duplex and Stealth RNAi against MMP15 were obtained from Invitrogen. Mission siRNAs against the other genes were purchased from Sigma-Aldrich (St Louis, MO). C6 cells that stably expressed Raichu-Rac1 were transfected with 20 mM siRNA by using Lipofectamine 2000 (Invitrogen). After 2 days of transfection, cells were used for invasion assay, qPCR or FRET imaging. The siRNA sequences are shown in supplementary material Tables S5.