Diverse roles of guanine nucleotide exchange factors in regulating collective cell migration

Collective cell migration involves intercellular mechanical communication through adhesive contacts (Tambe et al., 2011; Weber et al., 2012; Zaritsky et al., 2015). In migrating monolayers, such communication is initiated by cells at the monolayer boundary and gradually transmitted to cells at the back of the group (Ng et al., 2012; Serra-Picamal et al., 2012; Zaritsky et al., 2014, 2015; Ladoux et al., 2016; Mayor and Etienne-Manneville, 2016). Effective cell–cell communication requires balanced control of contractility and cell–cell and cell–matrix adhesions (Hidalgo-Carcedo et al., 2011; Weber et al., 2012; Cai et al., 2014; Bazellières et al., 2015; Das et al., 2015; Hayer et al., 2016; Notbohm et al., 2016; Plutoni et al., 2016). Coordination between these processes is regulated, among several pathways, by signaling activities of the Rho-family GTPases (Wang et al., 2010; Hidalgo-Carcedo et al., 2011; Timpson et al., 2011; Omelchenko and Hall, 2012; Cai et al., 2014; Omelchenko et al., 2014; Reffay et al., 2014; Plutoni et al., 2016). Rho-family GTPases are spatially and temporally modulated by complex networks of upstream regulators, including 81 activating guanine nucleotide exchange factors (GEFs), 67 deactivating GTPase-activating proteins, and 3 guanine dissociation inhibitors (Jaffe and Hall, 2005; Omelchenko and Hall, 2012). The networks are composed of many-to-one and one-to-many interaction motifs; that is, individual GTPases are regulated by multiple GEFs, and one GEF often acts upon multiple GTPases. Moreover, some GEFs are effectors of GTPases, leading to nested feedback and feedforward interactions (Schmidt and Hall, 2002; Jaffe and Hall, 2005; Cherfils and Zeghouf, 2013; Hodge and Ridley, 2016). Such pathway “design” permits an enormous functional specialization of transient signaling events, at specific subcellular locations and with precise kinetics.

Our long-term goal is to disentangle these signaling cascades in the context of collective cell migration. Although the roles of GEFs and their interactions with Rho GTPases are widely studied for single-cell migration (Goicoechea et al., 2014; Pascual-Vargas et al., 2017), less is known about how they regulate collective migration (Hidalgo-Carcedo et al., 2011; Omelchenko et al., 2014; Plutoni et al., 2016). Here, we report a comprehensive and validated, image-based GEF screen that identified differential roles of GEFs. By design of quantitative measures that encode the collective dynamics in space and time, we were able to identify a surprising role of RHOA, RHOC, and a group of four upstream GEFs in modulating collective migration via efficient long-range communication.

Collective cell migration emerges from the individual motility of cells in an interacting group: an action of one cell affects its neighbor and can propagate over time to eventually coordinate distant cells (Zaritsky et al., 2015). To identify molecules implicated in this mechanism, we performed live-cell imaging of the wound-healing response of human bronchial epithelial cells from the 16HBE14o (16HBE) line (Fig. 1 A and Video 1). Cells formed apical junctions and maintained epithelial markers and group cohesiveness before scratching the monolayer, as assessed by the localization of E-cadherin and the tight-junction protein ZO1 at the lateral cell–cell contact areas (Fig. 1 B). Upon scratching, the monolayer transitioned over ∼2 h from a nonmotile phase to an acceleration phase to steady-state wound closure (Fig. 1 C). The acceleration phase was associated with a gradual transition of cells from unorganized local movements to a faster and more organized motility. Cells at the wound edge underwent this transition first, followed by a wave of coordinated motility propagating away from the wound edge (Fig. 1 A, insets). The propagation is thought to be driven by mechanical cell–cell communication (Matsubayashi et al., 2011; Ng et al., 2012; Serra-Picamal et al., 2012; Zaritsky et al., 2014, 2015; Notbohm et al., 2016).

To screen for GEFs implicated in the migration-initiating process, we devised an automated and robust pipeline quantifying the spatiotemporal dynamics of motility activation in the monolayer. The core of this pipeline relies on robust detection of the wound edge via a segmentation algorithm that considers both image texture and segmentation consistency over time (Materials and methods). We then applied particle image velocimetry at subcellular resolution to measure the speed of cells in the monolayer as a function of their distance to the edge over time. Spatiotemporal dynamics were quantified and visualized in kymographs, in which every bin records the mean speed in a band of constant distance from the wound edge at a particular time point (Fig. 1 D). As expected, the speed kymograph showed a backward-propagating wave (Fig. 1 E, top left).

Using this assay, we first measured the effects of depletion of the canonical Rho GTPases Cdc42, Rac1, and RhoA. shRNAs were designed to selectively knock down each GTPase (Fig. 1 F). As expected (Simpson et al., 2008; Vitorino and Meyer, 2008; Omelchenko and Hall, 2012; Plutoni et al., 2016), knockdown of Cdc42 and Rac1 reduced cell speeds throughout space and time (Fig. 1, E and H, wound-healing rate; confirmed by a second distinct hairpin, Fig. 1 F and Fig. S1 A). We also confirmed inhibition of cell speed and monolayer migration under shRNA-mediated knockdown of the Rac1- and Cdc42-GEF β-PIX (Omelchenko et al., 2014; Fig. S1, B–D).

Knockdown of the Rho-GTPase RhoA induced near concurrent transition of cells at the front and back of the monolayer from a nonmotile to a motile state (Fig. 1, E and G). Importantly, RhoA-depleted cell monolayers reached the same wound-healing rates as unperturbed monolayers (Fig. 1 H). These data suggest that RhoA is critically implicated in the communication chain required to synchronize cell behavior throughout the monolayer after wound infliction.

Given the differential effects observed for the knockdown of Rho GTPases, we expected that down-regulating GEFs might also differentially alter the dynamics of motility activation. To capture the shifts in these behaviors, we established a concise representation of a live wound-healing experiment between wound infliction and steady-state wound closure, which could be compared across many experiments. We defined a 12-dimentional feature vector to encode the information contained in a speed kymograph over the first 180 µm of the monolayer and 200 min after scratching by averaging the speed values in bins of 60 µm and 50 min, respectively (Fig. 2 A; Materials and methods). Accordingly, features 1–4 encode the acceleration of cells at the monolayer front, features 5–8 the acceleration of cells 60–120 µm behind the monolayer edge, and so on. Features 1, 5, and 9 encode spatial variations in migration shortly after wounding. In the example of Fig. 2 A, cells at the front have already developed motility (feature 1), but cells further back in the monolayer (feature 9) have not. Features 2, 6, and 10; features 3, 7, and 11; and features 4, 8, and 12 encode the same spatial variations at later times. The pattern of 12-dimensional feature vectors was replicated for 402 control experiments without shRNA treatment (Fig. 2 B, top). The broad variation in collective migration behavior, even for constant experimental conditions, was visualized by normalizing each feature across the population of all 402 experiments (Fig. 2 B, bottom; see Materials and methods).

To identify core features that would capture the relevant variation of these data over noise, we applied principal component analysis (PCA; Jolliffe, 2002) to orthogonally transform the set of possibly correlated 12-dimensional features to a smaller set of linearly uncorrelated features denoted as the principal components (PCs). Indeed, when using PCs, the first three components captured 92% of the variability in the population of control experiments, allowing a fourfold dimensionality reduction of the feature space. The variation of PC1 between experiments was highly correlated with the variation in wound-healing rate (Fig. 2 C). The same applied to a lesser extent to PC3, whereas PC2 was uncorrelated. We further analyzed the coefficients that map the original 12 features into the scalar value of the respective PC. For PC1, all 12 coefficients have a similar value. Hence the mapping corresponds to an averaging of the speed values across the kymograph (Fig. 2 D). For PC2, the coefficients are positive for features 1, 5, and 9 (i.e., the speed values in the first time window) and negative for features 4, 8, and 12 (i.e., the speed values in the last time window). Thus, PC2 captures the reverse temporal gradient of the speed values. For PC3, the coefficients are positive for features 1–4 (i.e., the speed values at the wound edge) and negative for features 9–12 (i.e., the speed values at the back of the monolayer). Thus, PC3 captures the spatial gradient of the speed values in the direction of wound closure.

Knockdown of Cdc42 and Rac1 decreased PC1, reflecting a significant overall reduction in speed (Fig. 2 E). PC2 was somewhat increased for shCdc42 but unaltered for shRac1, whereas PC3 strongly decreased for both. This shows that the global reduction in speed under these conditions is accompanied by flatter spatiotemporal gradients (Fig. 2 E; see Materials and methods). Knockdown of RhoA most prominently altered PC3, with no effect on PC1 (Fig. 2 E), reflecting an alteration in the spatial gradient, in accordance with Fig. 1 E.

When a monolayer migrates collectively, cells within the monolayer eventually have to migrate toward the monolayer edge. This directional cue is transmitted via backward propagation of a strain gradient (Zaritsky et al., 2014). Thus, the spatial propagation of migration directionality is a measure related to intercellular communication. To capture this aspect, we defined for every kymograph bin directionality as the absolute ratio between the mean velocity component perpendicular to the monolayer edge and the velocity component parallel to the monolayer edge. The directionality kymograph displayed a backward-propagating wave pattern, as previously observed for other cell systems (Ng et al., 2012; Zaritsky et al., 2014; Fig. S1 E). Analogous to the analysis of the speed kymographs, we binned directionality kymographs into a feature vector of 12 values (Materials and methods; Fig. 2 A) and applied PCA. The first three PCs resembled the PCs for speed (Fig. S1 F), suggesting a functional coupling between speed and directionality establishment in the process of activating collective migration. Similar to the speed analysis, PC1 of the directionality correlated strongly with the wound-healing rate, but PC2 and PC3 were uncorrelated (Fig. S1 G). Depletion of the GTPases Cdc42 and Rac1 primarily affected the overall magnitude of directionality, whereas RhoA depletion reduced the spatial gradient, indicative of rapid long-range communication between the monolayer front and back after wound infliction (Fig. S1 H). Altogether, our data establish the first three PCs of a feature space derived from kymographs as measures for the detection of alterations in magnitude and temporal or spatial gradients in speed and directionality and underline the roles RhoA plays in long-range communication across the monolayer.

Equipped with our image-based assay for collective migration, we targeted the 81 Rho-family GEFs known in the human genome by shRNA (Fig. 3 A). A custom library of pSUPER shRNA retroviruses (three hairpins per GEF) was used for the screening, targeting 80 of 81 GEFs. 75 of the remaining 80 GEFs were confirmed to be expressed in the 16HBE cell line by Western blotting or quantitative RT-PCR (qRT-PCR) in cases in which no reliable antibody was available (Materials and methods). Western blots and qRT-PCR were used to evaluate the knockdown efficiency of every hairpin, and a cutoff of 50% depletion was selected for analysis of wound-healing experiments (Fig. 3 B). 16 GEFs had knockdown efficiency of less than 50% for all three hairpins and thus were eliminated from the screen. Within the remaining group of 59 GEFs, the knockdown of 11 GEFs could not be validated by qRT-PCR, because of failures in the production of efficient primers. Nonetheless, we included those genes in our imaging-based screen and labeled them “unknown knockdown efficiency.” A total of more than 2,600 time-lapse movies from the screen and the follow-up experiments were analyzed (Materials and methods). Control videos showed broad day-to-day variability in wound-healing rates (Fig. 3 C).

To score the alteration of a knockdown well, we quantified in multiple image locations the differences from the control well on the same plate. For every location, we extracted PCs 1–3 for speed and directionality. For each of the six PCs, we quantified the separation of target well and control well using three different distance metrics that balance intra- and inter-experiment variability (Materials and methods). For each combination of PC and distance metric, we used 24 experiments with no measurable knockdown (KD efficiency = 0%) as a baseline to calculate z-scores (Materials and methods).

To identify high-confidence hits, we estimated an objective z-score threshold by balancing between false-positive and false-negative rates. As positive controls we grouped together known motility reducers in Cdc42, Rac1, and β-PIX. Experiments with no measurable knockdown of the target were used as off-target controls. Their mean and standard deviation were used to calculate z-scores for all experiments. Distributions of z-scores in the speed PC measures are shown in Fig. 3 D. As expected, the positive controls were most discriminative for PC1. Thus, we exploited the known motility reduction of the positive controls together with the corresponding distribution of off-target controls to estimate a threshold appropriate to identify hits (Materials and methods). The threshold was estimated as 9.8. Because the threshold was standardized relative to the off-target controls, it was transferrable to the other PC-measures.

In addition to β-PIX, which has previously been reported as required for collective cell migration (Omelchenko et al., 2014), 10 new GEFs were identified with significant knockdown effects. Among those, we focused on SOS1, ARHGEF18, ARHGEF28, ARHEG11, and ARHGEF3 for validation (Fig. 3 E, black). The other five (Fig. 3 E, magenta; and Fig. 3, F and G, blue crosses) were not followed up for reasons discussed in Fig. S2 (A–E).

We first investigated SOS1, a dual GEF for RAC1 and RAS (Nimnual et al., 1998) that was found to regulate epithelial tight-junction formation in 16HBE cells through the MEK/ERK pathway (Durgan et al., 2015). Knockdown of Sos1 by two different hairpins induced a reduction in wound-healing rate, cell speed, and directionality (Fig. 3, F and G; and Fig. S2 F), consistent with its role as a Rac1 activator. Following the experiments in Durgan et al. (2015), we then blocked MEK and ERK activity directly by small-molecule inhibitors and concluded that the SOS1-RAS pathway (Fig. S2 G; Materials and methods) is also required for collective migration (Fig. S2 H), in addition to epithelial tight-junction formation.

Four paralog RHOA-GEFs were identified as hits (Fig. 3 E). Depletion of these GEFs did not dampen wound-healing rate, cell speed, or directionality but enhanced front-to-back propagation of motility. Specifically, depletion of Arhgef18, similarly to depletion of RhoA, accelerated front-to-rear long-range communication (Fig. 4 A). Depletion of Arhgef3, Arhgef11, and Arhgef28 somewhat accelerated long-range communication in speed and directionality and enhanced directionality overall (Fig. 3 G and Fig. 4, B–E). GEF-H1 (also called Lcf or ARHGEF2), a mechanoresponsive RhoA-specific GEF and paralog of the four GEFs identified in the screen (Birkenfeld et al., 2007; Chang et al., 2008; Guilluy et al., 2011b), was excluded from the list, because all hairpins produced a knockdown <50%. Nonetheless, experiments using these relatively inefficient hairpins revealed a phenotype similar to depletion of RhoA (Fig. 4 F, compare with Fig. 2 E and Fig. S1 H for RhoA).

In a set of validation experiments, each RHOA-GEF hit was replicated with at least two different hairpins in 6–10 wells (a total of 23–53 locations per hit). When assessing multiple replicates, systematic alterations were recognized that were missed by the stringent criteria of the screen (Fig. 4 G). For example, faster long-range communication in directionality was a common outcome of depleting any RhoA-GEF, while increased directionality (Arhgef11, Arhgef28, and Arhgef3) and long-range communication in speed (Arhgef18) were GEF specific. This suggests differential roles among RhoA-GEFs in regulating functions downstream of RhoA signaling.

Transmission of motility guidance cues from cell to cell results in the formation of clusters of cells moving with coordinated trajectories (Zaritsky et al., 2014; Materials and methods). Cluster formation was quantified by recording the fraction of the monolayer in which cells migrated coordinately. In control experiments, the fraction increased steadily over time (Fig. S3 A), because of an expansion of clusters from the front into the monolayer (Fig. S3, B and C). In analogy to speed and directionality, the spatiotemporal dynamics of cluster formation was captured by three PCs, which encoded magnitude, temporal, and spatial gradients (Fig. S3 D). Importantly, these PCs explicitly capture the strength and propagation of short-range communication. PC1, which was most associated with the wound-healing rate (Fig. S3 E), was increased upon depletion of Arhgef28, while depletion of RhoA, Arhgef18, or Arhgef28 increased PC3 in coordination (Fig. S3 F), showing that under these manipulations, the front and rear of the monolayer are coordinated rapidly after wounding.

Altogether, these results establish roles for RhoA and four of its activating GEFs in regulation intercellular communication.

Previous work showed that partial down-regulation of actomyosin contractility enhances passive force transmission through cells, whereas high levels of myosin activity “scramble” mechanical signals (Ng et al., 2015). Hence we speculated that the enhancement of long-range communication induced by RhoA and RhoA-GEF depletion is related to reduced myosin activity under these perturbed conditions. To test this we inhibited myosin II directly using blebbistatin or via ROCK inhibition using Y27632. To assess the effect of short versus long-term treatment, the drugs were applied without or with 24 or 48 h preincubation to wounding experiment.

Treatment with low doses of Y27632 (15 or 20 µM) and blebbistatin (10 µM) in general increased cell speed and coordination (Fig. S3, G and H). Preincubation for 48 h led to increased long-range communication in directionality and/or increased coordination similar to the behavior observed upon depletion of RhoA/RhoA-GEFs (Fig. S3, G–I). Inhibition of the formin pathway, which is also downstream of RhoA (Higashida et al., 2004; Fig. S3 J), with low doses of SMIFH2 (5 or 10 µM; Fig. S3, J and K) increased speed, directionality, and coordination without changing long-range communication (Fig. S3 K). Treatment with higher doses of Y27632 (25 µM for 48 h) or SMIFH2 (25 µM) inhibited overall motility (Fig. S3, G and K).

The effect of a 48-h drug treatment is expected to be biologically more similar to the effect of protein depletion by shRNA than the effects of more acute treatments; thus, we next compared the effects of blebbistatin and Y27632 treatment over 48 h with the effects of RHOA/RHOA-GEF depletion (Fig. 5 A). We did not consider formin inhibition in this analysis, because of the absence of a communication phenotype that could match the RHOA-GEF depletion phenotypes. Similarity analysis was performed by representing each experiment by a nine-dimensional feature vector composed of the normalized PCs 1–3 for speed, directionality, and coordination and calculating the similarity between every pair of conditions with the L1 norm. The pairwise (symmetric) similarity matrix confirmed distinct but related functionalities for RHOA/RHOA-GEFs/contractility in control of intercellular communication (Fig. 5 B). ARHGEF11, ARHGEF28, and ARHGEF3 fall in a first cluster and ARHGEF18, RHOA, and blebbistatin treatment in a second cluster. Treatments with 15 or 20 µM Y27632, although similar among each other, differ in their effect from any of the two other phenotypes. This demonstrates the sensitivity of our analysis to distinguish perturbations of the RhoA-myosin pathway from perturbation of Rock, which in part affects the RhoA-myosin axis but appears to be implicated in additional pathways driving collective migration. More critically, the sensitivity of our assay predicts that ARHGEF18 (and likely also ARHGEF2) regulate specifically RhoA-mediated activation of contractility, while ARHGEF11, ARHGEF28, and ARHGEF3 coregulate RhoA-independent pathways that do not converge on myosin II promoted processes. These may include RHOB/C (but not CDC42 or RAC1) for ARHGEF11 (Rümenapp et al., 1999; Jaiswal et al., 2011), RHOC (but not CDC42 or RAC1) for ARHGEF28 (van Horck et al., 2001; Bravo-Cordero et al., 2011, 2013), and RHOB (but not RHOC) for ARHGEF3 (Arthur et al., 2002), although the abundance of RhoC expression in 16HBE is fourfold less than that of RhoA (Wallace et al., 2011). Acute inhibition of all Rho isoforms (A, B, and C) with the small-molecule inhibitor Rhosin (Shang et al., 2012) showed a general motility reduction in a dose-dependent manner, indicating that Rho isoforms are required for collective migration (Fig. 5 C). We also excluded the possibility of crosstalk between RhoA and RhoC at the expression level. Western blots verified that knockdown of RhoA did not reduce RhoC (Fig. 5 D).

To test the prediction of differential regulation of RhoA and RhoC, we performed a set of RhoC knockdown experiments. Depletion of RhoC increased cell speed and coordination (PC1), long-range communication in directionality, and coordination (PC3; Fig. 5 E). A pairwise similarity matrix indicated that RhoC had an intermediate phenotype between the RhoA-ARHGEF18 and the cluster of the other RhoA-GEFs (Fig. 5, F and G).

Together, these analyses identified ARHGEF18 as the only RHOA-GEF exclusively activating the RHOA isoform (Blomquist et al., 2000; Herder et al., 2013), although biochemically it has also been described as a GEF for RAC1 (but not CDC42; Niu et al., 2003). Arhgef18 activates RhoA at tight junctions, directly interacting with myosin IIA and regulating tight-junction assembly (Terry et al., 2011; Xu et al., 2013; Durgan et al., 2015; Kim et al., 2015). Thus, we interpret the similarity of defects in long-range communication induced by Arhgef18 and RhoA depletion and direct inhibition of myosin II contraction as an indication that ARHGEF18 locates specifically at the top of an ARHGEF18/RHOA/MYO-II pathway (Fig. 5 H). ARHGEF11, ARHGEF28, and ARHGEF3 on one hand target the RHOA/MYO-II pathway. On the other hand, these three GEFs also target migration directionality, which is unaffected by RhoA-mediated signals. The intermediate phenotype of RHOC could be explained by integration of multiple pathways with similar but differential function (Kafri et al., 2009) and/or dynamic interactions between the molecular components that can be regulated in time and space (Guilluy et al., 2011a). The functional similarity between the ARHGEF18 and RHOC phenotypes, notably in cell speed, could be explained by a similar coil domain structure (Cook et al., 2014), predictive of competition, binding, or indirect interactions through other proteins of ARHGEF18 with RHOC. Determining how the ARHGEF18/RhoA-, RhoC-, and ARHGEF3,11,28-mediated pathways mechanistically differentiate between long-range communication, speed, and directionality will require an analysis of the spatiotemporal activation patterns during migration of these upstream GEFs in conjunction with the targeted GTPases, for example by construction of GEF activation biosensors.

A second limitation of our approach is its inability to identify intercellular communication phenotypes upon overall reduction of motility because of inherent reduction of the spatial and temporal gradients in relation to controls (Fig. 2 E, PC3). For example, TRIO’s Rho-targeting GEF domain is activated by a noncanonical Notch signaling pathway (Le Gall et al., 2008; Song and Giniger, 2011), known to be important in intercellular communication (Lai, 2004; Grego-Bessa et al., 2007). On the basis of our findings of a RhoA axis in the regulation of intercellular communication, it would be obvious to hypothesize that this role of TRIO is modulated by its RhoA-targeting DHPH2 domains (Bellanger et al., 2000). Unfortunately, this hypothesis could not be tested by knockdown experiments as performed here. Trio knockdown primarily caused an overall reduction of motility, possibly through its activating interaction with Rac1/RhoG mediated by the DHPH1 domain. A differential analysis of domain-specific effects in GEFs with multiple GTPase interactions is beyond the scope of this work.

Consistent with several other studies (Vicente-Manzanares et al., 2009; Lim et al., 2010), our data suggest again that although necessary for the generation of motility forces, myosin II–mediated contractility acts as an inhibitor for cell-to-cell communication, likely because contractility intercepts passive mechanical force transduction through the cellular cortex (Ng et al., 2015). Hence, during collective migration myosin II contractility must be adjusted to balance two opposing objectives: generation of robust motile forces versus transmission of mechanical guidance cues. Here we show that signaling pathways upstream of myosin II, such as the RHOA/ROCK pathway, need to be balanced in the same manner.

The human bronchial epithelial cell line, 16HBE14o (16HBE), was provided by the laboratory of D.C. Gruenert (University of California, San Francisco, San Francisco, CA). 16HBE cells were cultured in MEM (Memorial Sloan-Kettering Cancer Center [MSKCC] core facility), supplemented with 10% FBS (lot number 169905; Omega Scientific), GlutaMAX (35050; Gibco), and a mixture of penicillin-streptomycin (100×, 10,000 U/ml; 15140; Gibco. Stable cell lines were selected with 1.5 µg/ml puromycin (P7255; Sigma-Aldrich). HEK293T cells were purchased from ATCC and grown in DME high-glucose + sodium pyruvate (MSKCC core facility) supplemented with 10% FBS (lot number 169905; Omega Scientific) and a mixture of penicillin-streptomycin (100×, 10,000 U/ml; 15140; Gibco). All cells were cultured in a 37°C incubator with 5% CO₂.

Primary antibodies used for Western blotting include ARHGEF18 (EB06163; Everest) at 1:500, BCR (N-20, sc-885; Santa Cruz Biotechnology, Inc.) at 1:1,000, Cdc42 (610929; BD) at 1:1,000, GAPDH (FL-335, sc-25778; Santa Cruz Biotechnology, Inc.) at 1:1,000, Intersectin 2 (H00050618-A01; Abnova) at 1:1,000, LARG (N-14, sc-15439; Santa Cruz Biotechnology, Inc.) at 1:1,000, α-PIX (4573S; Cell Signaling Technology) at 1:1,000, β-PIX (07-1450; EMD Millipore) at 1:2,000, Rac1 (23A8; Abcam) at 1:2,000, RhoA (sc-418; Santa Cruz Biotechnology, Inc.) at 1:500, SOS1 (C-23, sc-256; Santa Cruz Biotechnology, Inc.) at 1:1,000, SOS2 (C-19, sc-258; Santa Cruz Biotechnology, Inc.), Tiam1 (C-16, sc-872; Santa Cruz Biotechnology, Inc.), and α-tubulin (MCA77S; AbD Serotec) at 1:2,000. Secondary polyclonal antibodies conjugated with HRP for Western blot were from Dako and used at 1:5,000. Primary antibodies used for immunofluorescence include E-cadherin (13-1900; Invitrogen) at 1:100 and ZO-1 (61-7300; Invitrogen) at 1:100. Secondary antibodies conjugated with Alexa 488 or Alexa 568 (Invitrogen) were used at 1:400. Other reagents used include Hoechst 33342 (Sigma-Aldrich) at 1 µg/ml, Y27632 (stock in H₂O; HA139; Sigma-Aldrich) at indicated concentrations, blebbistatin (stock in DMSO; 203391; EMD Millipore) at indicated concentrations, ERKi (stock in DMSO; SCH772984; provided by N. Rosen laboratory, MSKCC, New York, NY) at 1 µM, GSK1120212 (stock in DMSO; S2673; Selleckchem) at 500 nM, PD0325901 (stock in DMSO; S1036; Selleckchem) at 500 nM, Rhosin (stock in DMSO; EMD Millipore) at indicated concentrations, and SMIFH2 (stock in DMSO; EMD Millipore) at indicated concentrations.

An shRNA library was constructed in pSUPERpuro, containing at least three hairpins per gene for 80 predicted human Rho GEFs. An empty pSUPERpuro vector was used as a negative control (termed “control” in the figures and “pSuper” in Table S5). ARHGEF3, Rac1, RhoA, RhoC, and TRIO shRNAs were obtained from The RNAi Consortium library collection, in the pLKO.1 vector (MSKCC RNAi core facility) and a pLKO.1 vector with a nontargeting sequence was used as a negative control (catalog number SHC002; Sigma-Aldrich) for these experiments. The RhoA, Rac1, Cdc42, and RhoC hairpins reported in the study were Cdc42 sh1 (5′-GGAGAACCATATACTCTTG-3′), Cdc42 sh2 (5′-GATGACCCCTCTACTATTG-3′), Rac1 sh1 (TRCN0000004871; 5′-TTAAGAACACATCTGTTTGCG-3′), Rac1 sh2 (TRCN0000004873; 5′-TAATTGTCAAAGACAGTAGGG-3′), RhoA sh1 (TRCN0000047710; 5′-GTACATGGAGTGTTCAGCAAA-3′), RhoA sh2 (TRCN0000047711; 5′-CGATGTTATACTGATGTGTTT-3′), RhoC sh1 (5′-CTACTGTCTTTGAGAACTATA-3′), RhoC sh2 (5′-GCGAACCGGATCAGTGCCTTT-3′), RhoC sh3 (5′-TGATGTCATCCTCATGTGCTT-3′), and RhoC sh4 (5′-GAATAAGAAGGACCTGAGGCA-3′). Hairpin sequences for the Rho GEFs are documented in Tables S2 and S3. Hairpins with unmeasurable knockdown were defined as “off-target” controls and were used to assess the extent of off-target effects and to define thresholds for hit identification in the screen that minimize false detection rates (Fig. 3 D).

For virus production, 90% confluent human embryonic kidney 293T cells cultured in a six-well plate were transfected with lentiviral or retroviral constructs using Lipofectamine 2000 (11668; Invitrogen) and Opti-MEM (31985; Invitrogen), according to the manufacturer’s instructions. The culture media were removed the day after transfection, and media were collected three times at 24-h intervals. To infect cells, 2 × 10⁵ 16HBE cells were seeded in each well of a six-well plate and incubated on the following day with 1.5 ml virus-containing media supplemented with 1.5 µl polybrene (8 µg/µl stock; Sigma-Aldrich). Spin infection was performed at 2,250 rpm for 30 min. Cells were selected starting 2 d after infection with 1.5 µg/ml puromycin (Sigma-Aldrich). Pooled selected cell lines were used for all experiments.

Primers were designed using National Center for Biotechnology Information Primer-Blast to target all transcription variants of the gene and to span exon-exon boundaries to avoid amplifying genomic DNA. For qRT-PCR, primers were selected to amplify a PCR product between 70 and 150 bp in length with melting temperatures between 57°C and 63°C. Two sets of primers were examined for each gene, and the primer pair with the highest efficiency by qRT-PCR was selected for quantifying GEF expression. For some genes, QuantiTect primer assays (QIAGEN) were used. All primers used to quantify GEF expression are shown in Table S1.

Total RNA was extracted using the RNeasy Plus Mini kit (74134; QIAGEN). cDNAs were prepared using Oligo dT or Random hexamer primers (IDT Technology). In brief, RNA mixtures were heated at 65°C for 5 min, then chilled on ice and mixed with 5× Reaction buffer (Thermo Fisher Scientific), RiboLock RNase inhibitor (EO0381; Thermo Fisher Scientific), dNTPs (Sigma-Aldrich), and RevertAid reverse transcription (EO0441; Thermo Fisher Scientific). Reverse transcription reactions were performed using the following PCR program: 25°C for 10 min, 42°C for 60 min, and 72°C for 10 min, followed by cooling at 4°C.

To examine gene expression, each PCR was performed in a total 20-µl reaction, with 100 ng cDNA as template, 0.5 µM forward and reverse primer, 2 µl 10× PCR buffer (without magnesium; Invitrogen), 1 µM deoxynucleotide triphosphates (Sigma-Aldrich), 1.5 mM MgCl₂, 0.2 µl Taq polymerase (10342; Invitrogen), and 9.2 µl H₂O. The PCR program used was 94°C for 3 min; followed by 30 cycles of 94°C for 45 s, 60°C for 30 s, and 72°C for 1 min; followed by 72°C for 10 min and 4°C. PCR products were electrophoresed on a 2% TAE gel and stained with ethidium bromide. For qRT-PCR, 1 µg cDNA was used per reaction, in 25-µl reactions containing 1.25 µl 5 µM forward and reverse primers and 12.5 µl Maxima SYBR Green/ROX qPCR Master Mix (2×; K0221; Thermo Fisher Scientific). qRT-PCRs were performed on a Bio-Rad Laboratories iQ5 Multicolor RT-PCR Detection System with the following conditions: 95°C for 10 min; 40 cycles of 95°C for 15 s and 60°C for 60 s for gene expression detection; followed by 71 cycles of 60°C for 30 s, with increase of 0.5°C per cycle for melting curve detection. Gene expressions were normalized by the expression of GAPDH and HPRT and triplicate measurements were used for each sample.

Cells were washed in ice-cold PBS and lysed in ice-cold RIPA buffer (50 mM Tris-HCl, pH 7.4, 150 mM NaCl, 2 mM EDTA, 1% NP-40, and 0.1% SDS) supplemented with 5 mM Na₃VO₄, 10 mM NaF, 25 mM β-glycerolphosphate, and 1 mM PMSF. Lysates were collected by cell scraping and cleared by centrifugation at 13,000 rpm for 1 min at 4°C and boiled in sample buffer (final concentration 50 mM Tris HCl, pH 6.8, 2% SDS, 10% glycerol, 0.1% bromophenol blue, and 100 mM DTT) for 5 min. Protein concentrations were determined by bicinchoninic acid assay (23225, Pierce BCA protein assay kit; Thermo Fisher Scientific). SDS-PAGE on 3%–8% Tris-acetate gels, 4%–12% Bis-Tris gels, or 12% Bis-Tris gels (NuPAGE; Thermo Fisher Scientific); transfer to polyvinylidene fluoride membranes (0.45 mm pore size; EMD Millipore); blocking; antibody binding; and ECL-mediated detection were performed as previously described (Durgan et al., 2015).

Western blots were used to assess protein expression levels for 9 GEFs that had validated antibodies available. qRT-PCR was applied for the remaining 71 GEFs to assess gene expression levels. The knockdown of 11 GEFs could not be assessed by qRT-PCR because of failure in the production of efficient primers.

16HBE cells grown on glass coverslips were washed with PBS and fixed in 3.7% (vol/vol) formaldehyde (Sigma-Aldrich) in PBS for 10 min at room temperature. For immunostaining, coverslips were washed three times in PBS and blocked with BTPA buffer (0.5% BSA, 0.02% sodium azide, and 0.25% Triton X-100 in PBS) for 30 min. After blocking, coverslips were incubated with primary antibodies against ZO-1 and E-cadherin (diluted in BTPA buffer) for 1 h at room temperate and washed three times in PBS for 5 min. Coverslips were then incubated with secondary antibodies and Hoechst (diluted in BTPA buffer) for 1 h at room temperature, washed three times in PBS for 5 min, and mounted onto microscope slides (Thermo Fisher Scientific) using fluorescent mounting media (Dako). Epifluorescence images were acquired with an upright Imager.A1 microscope (ZEISS), equipped with an EC-Plan-NEOFLUAR 40×/0.75 objective and a Hamamatsu Photonics Orca-ER 1394 C4742-80 camera, controlled by Axiovision software (ZEISS). Scale bars were added using ImageJ (National Institutes of Health).

For wound-healing assays, 3 × 10⁶ 16HBE cells were seeded in each well of a six-well tissue culture plate and incubated for 2 d before wounding. Wounding was performed by scratching with a P1000 pipette tip on the confluent monolayers, in the middle of each well, and a cell scraper was used to remove half of the cells from the plate. After washing with PBS several times to remove cell debris and adding fresh 16HBE media, plates were imaged on an inverted Axiovert 200M microscope (ZEISS) equipped with an EC Plan-NEOFLUAR 10×/0.3 Ph1 objective, a Hamamatsu Orca-ER 1394 C4742-80 camera, a 37°C incubator and a CO₂ controller, controlled by Axiovision software (ZEISS). The time-lapse image sequences were recorded for 16 h at 5-min intervals.

Velocity fields were computed using custom cross-correlation-based particle image velocimetry using nonoverlapping image patches of size 15 × 15 µm (Zaritsky et al., 2012). The frame-to-frame displacement of each patch was defined by the maximal cross-correlation of a given patch with the subsequent image in the time-lapse image sequence. The search radius was constrained to an instantaneous speed of 90 µm h⁻¹. At a frame rate of 5 min, this search radius corresponds to half the side length of a 15 × 15 µm patch.

Monolayer contours were calculated with a segmentation algorithm that classified the image regions as either “cellular foreground” or “background.” Our goal was to optimize the robustness of the segmentation to enable unsupervised analysis of thousands of movies, each containing tens of time points. Small inaccuracies in the segmentation have little effect on the resulting kymograph, every movie was manually validated, and fewer than 1% of the well locations were discarded on the grounds of failed segmentation. To achieve robustness, we introduced several priors to the algorithm: (1) the image contains one continuous region of “cellular foreground” and one continuous region of “background,” and (2) the contour advances monotonically over time. These priors allowed us to estimate the initial contour at time 0 and then use the segmentation at time t as a seed to expand the “cellular foreground” to time t + 1. The only pixels in question are those labeled “background” at time t and close enough to the “cellular foreground” region. The proximity threshold is calculated on the basis of the maximal cell velocity of 90 µm h⁻¹.

The segmentation was performed in super-pixels with a size equivalent to a 15 × 15 µm patch. After cross-correlation-based particle image velocimetry, which assigns to each super-pixel a displacement vector and a cross-correlation score defining the match of the corresponding patch signals between consecutive frames, we took advantage of the observation that super-pixels associated with “cellular foreground,” especially those at the front of the monolayer, had lower cross-correlation scores because of their textured and dynamic nature. We thus assigned to each super-pixel a pseudo-intensity value (1 – cross-correlation score) and applied the Rosin thresholding method (Rosin, 2001) to label the super-pixels as “cellular foreground” versus “background.” Background super-pixels enclosed by foreground regions were relabeled as “cellular foreground,” and isolated “cellular foreground” super-pixels, usually attributed to debris or textured plate patterns, were relabeled as “background.” We next calculated the union of the “cellular foreground” region in the previous frame with the new “cellular foreground” region. Morphological closing and hole filling defined the final segmentation. For the first image in the time-lapse sequence, we defined the “cellular foreground” as the union of the cross-correlation-based segmentation and a texture-based segmentation. The latter was implemented by first representing each super-pixel by the distribution of its local binary patterns, a widely used representation of local texture in images (Ojala et al., 2002), followed by unsupervised K-means clustering with K = 2. This heuristic was found to robustly complement the cross-correlation-based segmentation by identifying super-pixels mislabeled as “background” at larger distances from the wound edge, where cells do not sufficiently migrate and change appearance to generate a low enough cross-correlation score.

During the multilocation filming, the microscope stage exhibited repositioning errors that caused systematic shifts in all vector fields. We estimated the stage shift for every image in a time-lapse sequence and corrected the vector fields accordingly. To accomplish this, we exploited patches in the background region that are expected to stay in the same position and subtracted their robust mean velocity from the vector fields in the “cellular foreground.”

Speed kymographs were constructed by calculating the mean speed of all patches in spatial bands of 15 µm from the monolayer’s front through time. Cell directionality is defined as the absolute ratio between the velocity component perpendicular to the monolayer edge and the velocity component parallel to the monolayer edge. Each bin in the directionality kymograph was calculated as the ratio obtained by the two-component decomposition of the speed kymograph to a component normal and parallel to the monolayer front. These components were calculated by considering the orientation of the wound edge. Coordination (or short-range communication) was measured by detecting clusters of cells moving in coordination. Coordination kymographs were constructed by recording the fraction of patches that participated in these coordinated clusters for every spatial band and time frame. Explicit detection of coordinated clusters was performed for every time frame by applying region-growing spatial clustering of the image-patch grid on the basis of the correlation of their velocity fields, as previously described (Zaritsky et al., 2014). In brief, region-growing segmentation (Nock and Nielsen, 2004) started with regions containing a single image patch and iteratively merging spatially adjacent patches on the basis of their velocity vector similarity. Two regions are merged if their similarity is lower than a given threshold, and the merged region vector is updated to be the mean of all contributing patch vectors. Merging is performed in ascending order of the similarity between the adjacent patches. Kymographs were calculated for cells located up to 180 µm from the monolayer and for the first 200 min after wound infliction. We chose 180 µm because the cell monolayers captured in the field of view of the vast majority of experiments were equal or wider than this value. We chose 200 min to focus our analysis on the transient phase between wound infliction and steady-state collective migration. Extending the time window to 400 min did not change the conclusions drawn from kymograph analysis.

To obtain a compact representation of the kymographs, we averaged the kymographs in 3 × 4 bins at a resolution of 60 µm and 50 min (Zaritsky et al., 2012). To further reduce this 12-dimentional feature vector, we applied PCA (Jolliffe, 2002). The PCs are ranked by the spread of the data they capture. Mathematically, this is equivalent to the eigenvalues of the data covariance matrix. We used the 402 control experiments of the screen to calculate the PCA transformation for speed, directionality, and coordination and then applied the same transformation on control and knockdown experiments. First, the features were normalized (Fig. 2 B) to x′ = (x − μ)/σ, with feature x, mean μ, and standard deviation σ of the set of control experiments. The PCA transformation was calculated and applied to normalized features. Nearly identical PCA transformations were found when including the knockdown experiments (and for non-normalized features), indicating that the functional fluctuations between control and knockdown experiments are small compared with the day-to-day variation of control experiments. The coefficients of each PC defined the projection of the 12-dimensional feature vector into the direction of the PC. The coefficients of PC1 were similar for all 12 features, implying that PC1 encodes the mean of speed, directionality, or coordination (dependent on the considered kymograph) across time and space. The coefficients of PC2 and PC3 reflected the pattern of a temporal (PC2) or spatial (PC3) gradient (Figs. 2 D, S1 F, and S3 D; see text). It should be noted that conditions with overall reduced motility, such as with CDC42 and RAC1 knockdown, not only reduce PC1 but inherently flatten the spatial and temporal gradients in speed, thus reducing PC3 but increasing PC2 of this measure. This increase of PC2 is due to the positive coefficients for the (slower) speed at the onset of migration, which is minimally reduced under these perturbations, while the negative coefficients for reduced speed of perturbed experiment at later times increase overall values. In Figs. 2 E; 4, F–G; 5, C and E; S1 H; S2 H; and S3, F, G, and K, we depict the shifts induced by knockdown by subtraction of the PC of control experiments from the PC of the knockdown experiment. For PC1, negative values indicate reductions in speed, directionality, or coordination. For PC3, which reflects the spatial gradient, negative difference values imply more immediate front-to-back propagation and faster establishment of the steady state, while for PC2 this is reversed.

Time-lapse image sequences were acquired in at least three different locations in every knockdown and in every control well. To identify hits in the screen by comparison of knockdown and control conditions, we used metrics that take into account the variability within a well and between wells. Specifically, we applied three metrics: (1) the Davies-Bouldin index (Davies and Bouldin, 1979):

which divides the difference d(c₁, c₂) between the per-well mean values of a measurement extracted from multiple locations in the knockdown and control wells by the sum of the mean distances σ_i between the measurement from individual locations in the same well and the corresponding well mean value; (2) the Dunn index (Dunn, 1974):

which divides the difference between per-well mean values by the maximal distance in any one well between the measurement from individual locations and the corresponding well mean value; (3) the silhouette coefficient (Rousseeuw, 1987):

where n and m are the numbers of locations in control and knockdown wells, and

where a is the mean difference between the measurement extracted in location k and the measurements extracted in all other locations of the same well, and b is the mean difference between the measurement in location k and the measurements extracted in all locations of the other well.

Off-target control experiments (see text) were used to normalize the three metrics across different PC measurements and to define a z-score = (x − μ)/σ, with μ and σ denoting the mean and standard deviation of a particular PC metric in the off-target control population and x denoting the metrics of that PC for a knockdown experiment. We defined a single combined z-score scalar for every PC measure by accumulating the three metric-specific z-scores.

The 24 off-target and 18 positive controls were used to calculate a threshold for hit identification in the screen, which has the desired property of minimizing the number of false alarms to allow us to focus on real hits for follow-up experiments. Because the positive controls (CDC42, RAC1, and β-PIX) were selected on the basis of their known effect in reducing motility, the threshold was calculated to minimize false alarms in wound-healing rate. Calculating z-scores by standardization relative to the off-target controls made the threshold transferrable to the other PC measures. Specifically, we calculated the z-score threshold that maximizes the F-measure, defined as

with recall

and precision

derived from the true-positive (TP), false-negative (FN), and false-positive (FP) rates. The true-positive and false-positive rates were estimated on the basis the positive controls, whereas the false-negative rate was estimated on the basis of the off-target controls. Because of the small number of control experiments, the threshold had low confidence. To resolve this, we applied bootstrapping to estimate the distribution of the maximal F-measure thresholds. Off-target and positive controls were randomly selected with repetition from the original groups. The process was repeated 10,000 times to define the distribution of z-score values maximizing the F-measure, and the final threshold of z-score = 9.8 was selected with 99.9% confidence from this distribution. There were no false-positive occurrences in five of six PC measures, and one false positive occurred for the sixth PC measure (directionality PC2); the false-negative rate was 44.4%.

To validate the hits that were identified in the screen, we replicated experiments and statistically confirmed the phenotypes. For every experiment, we subtracted the mean PC measure of the daily control locations from the corresponding PC measure in every knockdown location. Thus, each data point represents the deviation of a location to the mean control. The null hypothesis was that the data come from a distribution whose median is zero. Statistical significance was inferred using the nonparametric Wilcoxon signed rank test, and a p-value of 0.01 was selected as the significance threshold. At least two different hairpins were required to validate a hit. SOS1 experiments were replicated with N (number of independent experiments) = 3–6 and n (total number of locations in knockdown wells) = 18–26; RHOA-GEFs were validated with n = 6–9 and n = 23–53 and contractility experiments (treated for 48 h) with n = 3 and n = 18. RhoGTPAses were validated with n = 6–9 and n = 24–47 and β-PIX with n = 3 and n = 11.

To calculate overall similarities between different experimental conditions, we included the three first PCs for speed, directionality, and coordination in a nine-dimensional vector for each well location of a knockdown experiment. As in the follow-up analysis, PC measures in the locations of a knockdown well were related to the mean of the corresponding PC measures in all locations of the daily control well. An experimental condition was represented as the mean nine-dimensional vector of all locations in all wells under the same treatment regimen. This defined a matrix with nine rows (one per measure) and k columns, where k is the number of different conditions examined for the analysis (e.g., eight for the analysis described in Fig. 5 B). Each of the nine measures was standardized across conditions by subtracting the mean and dividing by the standard deviation (z-scores) for each row in the matrix. This encoded the divergence of every specific condition from the mean across the examined experimental conditions. Similarity was calculated between every pair of conditions (columns) in nine dimensions of the standardized measures using the |L1| distance metric. High values reflect dissimilarity, and low values represent more similar phenotypes in this nine-dimensional space.

All following metadata are made available: (1) Primers used to quantify GEF expression are listed in Table S1. (2) shRNA hairpin sequences and knockdown efficiency (Western blots and qRT-PCR) are listed in Table S2 and summarized in Table S3. (3) Details of all screen experiments are provided in Table S4 and of all screen plus follow-up studies in Table S5: molecular perturbation, raw-data file name, file format (.zvi or .tif), physical pixel size in micrometers, measured knockdown efficiency (%), and experiment date. (4) Screen z-scores are provided in Table S6 and coordination z-scores from follow-up experiments in Table S7.

All raw and processed data are available upon request. These data include raw imaging data, processed kymographs (per location and mean per day), and quantification of phenotypes for genes that were followed up. The full data set size is 3 TB. Potential requesters are asked to send a drive of sufficient size.

The Matlab source code is available to the public at https://github.com/DanuserLab/MonolayerKymographs. It enables calculation of the velocity fields and segmentation of the monolayer front and the wound-healing rate for every time frame in the raw time-lapse images, and it generates kymographs for speed, directionality, and coordination for a full video. It will also produces the 12-dimentional feature vector representation and PC projection on the basis of the transformation that was calculated for this study. We also provide a script to perform PCA and compute the transformation given a set of high-dimensional experiments. Input is as follows: M = d × n matrix, where d is the number of dimensions (in our case, d = 12), and n is the number of experiments. Output is a PCA transformation (a wrapper for Matlab’s PCA code). This code includes means of visualization of the PCA weights (such as in Fig. 2 D) to enable interpretation. Documentation and a test example are included with the source code.

Fig. S1 depicts alterations of collective motility induced by depletion of Rho GTPases or the GEF β-PIX. Fig. S2 shows screen hits that were not followed up (A–E) and follows the SOS1-RAS pathway regulation of collective cell migration. Fig. S3 shows the effects of RhoA GEF depletion on short-range communication as quantified by the coordination parameter (A–F) and the effects of perturbation of actomyosin contractility on intercellular communication (G–I). Video 1 shows phase-contrast live imaging of the wound-healing response of 16HBE cells. Tables S1–S7 are available as Excel files. Primers used to quantify GEF expression are listed in Table S1. shRNA hairpin sequences and KD efficiency (Western blots and qRT-PCR) are listed in Table S2 and summarized in Table S3. Details of all screen experiments are provided in Table S4 and of all screen plus follow-up studies in Table S5. Screen z scores are provided in Table S6 and coordination z-scores from follow-up experiments in Table S7.

We thank A. Jamieson for helping package the source code and creating the repository. We thank C. Schaefer and N. Bendris for fruitful discussions and advice. We dedicate this work to our mentor and colleague A. Hall, who initiated this project.

This work was supported by the National Institutes of Health (P01 GM103723 to G. Danuser, A. Hall, and M. Overholtzer; P30 CA008748 to M. Overholtzer).

The authors declare no competing financial interests.

Author contributions: A. Hall and G. Danuser conceived the study. M. Overholtzer, A. Hall, and G. Danuser guided A. Zaritsky, Y.-Y. Tseng, MA. Rabadán, and S. Krishna. A. Zaritsky, Y.-Y. Tseng, and MA. Rabadán designed the experiments. Y.-Y. Tseng, MA. Rabadán, and S. Krishna performed all experiments. A. Zaritsky developed analytic tools and analyzed the data. A. Zaritsky and Y.-Y. Tseng interpreted the data. A. Zaritsky drafted the manuscript. All authors wrote and edited the manuscript and approve of its content.