A Computational Systems Biology Approach Identifies SLUG as a Mediator of Partial Epithelial-Mesenchymal Transition (EMT)

Cancer metastasis remains the major cause of cancer-related deaths [Gupta and Massagué, 2006]. Metastasis is a dynamic process with extremely high (>99.5%) attrition rates; the minority of cells successful in establishing metastases are highly plastic, i.e., capable of adapting to their environmental changes [Celià-Terrassa and Kang, 2016; Jolly et al., 2017b]. An important axis of plasticity involved in metastasis is epithelial-mesenchymal plasticity (EMP), i.e., bidirectional transitions among epithelial, mesenchymal, and one or more hybrid epithelial/mesenchymal (E/M) phenotypes, that is, a complete or partial epithelial mesenchymal transition (EMT) and its reverse phenomena mesenchymal epithelial transition (MET) [Jolly and Celia-Terrassa, 2019]. Recent in vitro, in vivo, and in silico studies suggest that EMT or MET are not “all-or-none” processes; instead, cells can attain one or more hybrid E/M phenotypes which may be more plastic as compared to “fully” epithelial or mesenchymal ones [Font-Clos et al., 2018; Pastushenko et al., 2018; Tripathi et al., 2020]. Consistently, cells acquiring hybrid E/M phenotype(s) are crucial for tumorigenicity and/or metastasis-initiating properties across many cancer types [Pastushenko et al., 2018; Kröger et al., 2019], and undergoing a “full” EMT may lead to loss of such “stemness” traits [Celià-Terrassa et al., 2012; Bierie et al., 2017]. Thus, hybrid E/M cells are likely to be the “fittest” for metastasis [Jolly et al., 2018], a hypothesis which is strengthened by increasing reports suggesting an association of hybrid E/M phenotypes with worse clinicopathological features [Godin et al., 2020; Liao and Yang, 2020]. Therefore, it is critical to understand the molecular players regulating hybrid E/M phenotypes to decode their effect on aggravated metastasis and identify therapeutic approaches to restrict them.

Investigation of EMT/MET in developmental contexts have helped identify the modulators of cancer cell plasticity [Nieto et al., 2016]. One of the key molecules that has been associated with hybrid E/M phenotypes in various developmental contexts is SLUG (SNAIL2) [Jolly et al., 2015]; resonating observations have also recently been seen in clinical settings [Steinbichler et al., 2020]. However, this association has been mostly at a phenomenological level. SLUG belongs to a family of zinc-finger transcription factors SNAI. In mammals, this family consists of SNAIL (encoded by the SNAI1 gene; also SNAIL1), SLUG (encoded by the SNAI2 gene; also SNAIL2), and SMUC (encoded by the SNAI3 gene; also SNAIL3) [Barrallo-Gimeno and Nieto, 2005]. While SMUC is much less studied, SNAIL and SLUG are observed across species from Drosophila melanogaster to humans with functions varying from development to oncogenesis [Wu et al., 2005]. In Drosophila, they are involved in formation of mesoderm and development of the central nervous system [Hemavathy et al., 2000]. In vertebrates, they are involved in the formation of the neural crest [Hu et al., 2014] and EMT of mesodermal cells, and they can also participate in EMT-independent functions such as preventing apoptosis from gamma-irradiation in hematopoietic cells [Barrallo-Gimeno and Nieto, 2005]. The role of SLUG and SNAIL – the 2 more well-studied members of the SNAI family – has been found to be partially redundant [Stemmler et al., 2019; Hemavathy et al., 2000]. But a detailed comparison between them, at least in terms of their ability to induce a partial EMT, i.e., hybrid E/M phenotype(s), remains to be done to acquire mechanistic insights into their role in mediating the dynamics of EMP.

Here, we investigate the differences between the EMT-inducing potential of SLUG and SNAIL and investigate the contribution of SLUG in being able to maintain a hybrid E/M phenotype. Previous mathematical models developed for EMP signaling have often not distinguished between SNAIL and SLUG, and have taken them as a single node [Lu et al., 2013; Tian et al., 2013; Hong et al., 2015; Mooney et al., 2016]. With increasing knowledge about the different roles of SNAIL and SLUG in influencing EMP in different contexts, including their effect on one another [Chen and Gridley, 2013; Gras et al., 2014; Ganesan et al., 2015; Nakamura et al., 2018; Sundararajan et al., 2019; Swain et al., 2019; Hamidi et al., 2020; Varankar et al., 2020], it becomes imperative to investigate their differential roles from a detailed mechanistic perspective. Therefore, we have adopted a mechanism-based computational systems biology approach to elucidate these differences. Our results suggest that SLUG is a weaker inducer of EMT as compared to SNAIL; thus, SLUG, but not SNAIL, can often stabilize a hybrid E/M phenotype. Clinical analysis indicates high SLUG expression to be associated with worse patient survival across carcinomas, endorsing the growing evidence for a more aggressive behavior of hybrid E/M phenotype(s).

All computational analyses were performed using MATLAB (version 2018b), and microarray datasets were downloaded using the GEOquery R Bioconductor package [Davis and Meltzer, 2007]. Pre-processing of the microarray datasets was performed for each sample to obtain the gene-wise expression from probe-wise expression matrix using R (version 4.0.0).

The EMT scores for each dataset were calculated using 76 GS, KS, and MLR methods as described earlier [Chakraborty et al., 2020].

76 GS: A weighted sum of 76 gene expression levels was used to calculate the score for each sample. The mean across all tumor samples was subtracted to center the scores so that the grand mean of the score was zero. Negative scores can be interpreted as M phenotype, whereas the positive scores as E phenotype.

MLR: The EMT status by the MLR method is predicted based on order of the categorical structures and the principle that the hybrid E/M state falls in a region intermediary to E and M. The sample scores range between 0 and 2 where 0 corresponds to pure E and 2 corresponds to pure M, and a score of 1 indicates a maximally hybrid phenotype. These scores are calculated based on the probability of a given sample being assigned to the E, E/M, and M phenotypes.

KS: The scores were calculated by comparing cumulative distribution functions of epithelial (E) and mesenchymal (M) signatures. If the score is positive, then the sample is said to be M, and if it is negative, it is considered to be E.

The stemness score was calculated using CNCL as described [Akbar et al., 2020]. To identify the stem cell-like (CS/M) and non-stem cell like (NS/E) cells, the samples were clustered based on the expression of CSC/non-CSC gene list (CNCL). The quantitative stemness score was calculated by the correlation of CNCL gene expression values for each sample to the CS/M and NS/E matrices generated based on the median expression levels for each gene in the CNCL for CS/M and also for NS/E cells in the CCLE database. The 2 Pearson’s correlation values generated for each matrix, CS/M(r) and NS/E(r), are plotted.

As per the schematic shown in Figure 2a, the dynamics of all 4 molecular species (miR-200, SNAIL, ZEB, and SLUG) was described by a system of coupled ordinary differential equations (ODEs). The generic chemical rate equation given below describes the level of a protein, mRNA, or micro-RNA (X): dX/dt = g_XH^S(A,A⁰,n,λ) - k_XX

Where the first term g_X signifies the basal rate of production, the terms multiplied to g_X represent the transcriptional/translational/post-translational regulations due to interactions among the species in the system, as defined by the Hill function H^S(A,A⁰,n,λ). The rate of degradation of the species (X) is defined by the term k_XX based on first order kinetics. The complete set of equations and parameters are presented in the supplementary material (for all online suppl. material, see www.karger.com/doi/10.1159/512520).

As the degradation rate of ZEB mRNA and SLUG mRNA is much greater than that of ZEB protein and miR-200, and also the production rate of SLUG and SNAIL protein is much larger than that of ZEB protein and miR-200, we assumed that ZEB mRNA, SNAIL protein SLUG mRNA, and SLUG protein reach to the equilibrium much faster relatively, i.e., (dm_Z)/dt = 0, dS/dt = 0, (dm_Sl)/dt = 0 and dSl/dt = 0.This assumption reduces the equations given in the online supplementary material to 2 coupled ODEs of ZEB and miR-200. Then, the dynamical system was simulated to obtain the time evolution of ZEB protein and miR-200 in presence of external noise using Euler-Maruyama simulation. Dynamical states of the system were coarse grained as an itinerary of basins visited from the time evolution of ZEB and miR-200. MRT was calculated by multiplying the total number of successive states with it [Biswas et al., 2019].

We performed random circuit perturbation (RACIPE) analysis on the extended circuit (Fig. 4a) using the default settings for 10,000 parameter sets. We considered all RACIPE steady-state solutions up to tetra-stable parameter sets for principal component analysis (PCA) plots. PCA was performed on all the steady-state solutions from the RACIPE output data. The clusters were identified by performing hierarchical clustering on the RACIPE data by setting the number of clusters to 4, and the number of steady-state solutions ending up in each of these clusters were quantified. We also performed RACIPE on the same circuit but with either overexpressing SNAIL or SLUG by 10-fold. The same workflow as described above was used to analyze the overexpression-based analysis.

ProgGene [Goswami and Nakshatri, 2014] was used for conducting Kaplan-Meier analysis. The number of samples in SLUG-high versus SLUG-low categories are given below:

• GSE31210 (pathological stage I-II lung adenocarcinomas): n (high) = 113, n (low) = 113

• GSE9891 (285 serous and endometrioid tumors of ovary, peritoneum, and fallopian tube): n (high) = 138, n (low) = 138

• GSE30161 (Ovarian carcinoma): n (high) = 29, n (low) = 29

• GSE14333 (Primary colorectal cancer sample): n (high) = 94, n (low) = 93

• TCGA-LUAD (Lung cancer sample): n (high) = 75, n (low) = 75

• GSE17536 (Colorectal cancer sample): n (high) = 87, n (low) = 87

• GSE16125 (Colorectal cancer sample, stage [Duke’s I-IV]): n (high) = 16, n (low) = 16

• GSE50827 (103 primary pancreatic ductal adenocarcinoma samples): n (high) = 16, n (low) = 16

• GSE26712 (Advanced stage ovarian cancer sample): n (high) = 93, n (low) = 92

• GSE18520 (late-stage, high-grade papillary serous ovarian adenocarcinoma sample):

• n (high) = 27, n (low) = 26

• GSE13876 (157 advanced stage serous ovarian cancer sample): n (high) = 207, n (low) = 207

• GSE3143 (Breast cancer sample): n (high) = 80, n (low) = 77

• GSE2034 (Breast cancer sample): n (high) = 143, n (low) = 143

• GSE19615(Breast cancer sample): n (high) = 58, n (low) = 57

• GSE11121 (Node-negative breast cancer patient sample): n (high) = 103, n (low) = 97

• GSE2990 (189 invasive breast carcinoma sample): n (high) = 51, n (low) = 50

First, using publicly available transcriptomics datasets (GSE80042. GSE40690, GSE43495), we quantified the extent of EMT induced by SNAIL and SLUG individually in different cell lines. We used 3 different methods to quantify EMT – 76 GS, MLR, and KS. MLR and KS score the samples on a spectrum of [0, 2] and [-1, 1]; the higher the score, the more mesenchymal the sample is. On the other hand, there is no predefined range of scores calculated by 76 GS, and the higher the 76 GS score, the more epithelial the sample is [Chakraborty et al., 2020].

In LNCaP prostate cancer cells with inducible and reversible expression of SNAI1 or SNAI2 (GSE80042), those induced with SNAI1 showed a stronger extent of EMT as compared to those induced with SNAI2. This trend was consistently seen across the 3 scoring methods (Fig. 1a). Similarly, in immortalized human mammary epithelial cells transfected with retroviral constructs containing SNAI1, SNAI2 and SNAI3 (HMEC-hTERT), those induced with SNAI1 showed a significantly stronger EMT induction as compared to those induced with SNAI2 or SNAI3 (GSE40690); the degree of induction was the weakest for SNAI3. Again, this trend was consistent across scoring methods (Fig. 1b). Consistently, in human mammary epithelial cells, transduction of SNAI1 induced a stronger EMT than that of SNAI2 or TWIST (Fig. 1c) (GSE43495). In all these scenarios, the scores calculated across these methods correlated well with one another – while KS and MLR scores were positively correlated, 76 GS scores correlated negatively with KS and with MLR scores, as expected (online suppl. Fig. 1a–c).

Next, for GSE43495 and GSE40690 datasets, we also calculated the stemness scores using a gene list corresponding to EMT and stemness for breast cancer cells (CNCL) [Akbar et al., 2020]. Cells overexpressing SNAI1 were found to show a stronger enrichment of this gene signature as compared to those overexpressing SNAI2 and/or TWIST (Fig. 1d; left, top, bottom). When projected on the 2-dimensional representation of NS/E (non-cancer stem-like/epithelial) and CS/M (cancer stem cell-like/mesenchymal) axes, the cells induced with SNAI1 were found to be more of CS/M^Hi-NS/E^Low as compared to those induced with SNAI3 (CS/M^Low-NS/E^Hi) (Fig. 1d; right, top); the SNAI2 and/or TWIST-induced ones had a weaker enrichment of CS/M^Hi-NS/E^Low signatures (Fig. 1d; right, bottom). Put together, these results suggest that SNAIL is a stronger inducer of EMT than SLUG, and indicates the possibility of SLUG being associated with a partial EMT phenotype instead of a completely mesenchymal one.

To disentangle the different effects of SNAIL versus SLUG on EMT, we expanded the previous EMP circuitry we simulated earlier [Lu et al., 2013] by incorporating the interconnections of SLUG with miR-200 and SNAIL (Fig. 2a). SLUG and SNAIL can inhibit each other, albeit indirectly; also, SLUG and miR-200 have been known to form a mutually inhibitory feedback loop (online suppl. Table 2). These interconnections were represented by a set of coupled ODEs (online suppl. Section 1 and Table 1).

We started with first calculating the bifurcation diagram of ZEB mRNA levels for varying levels of an external EMT-inducing signal (Fig. 2b). As the value of EMT-inducing signal increases, the cells can switch from an epithelial state (ZEB mRNA < 100 molecules) to a hybrid E/M phenotype (100 < ZEB mRNA molecules <600) and finally to a mesenchymal state (ZEB mRNA > 600 molecules); this trend is seen for the circuit that does not include SLUG (Fig. 2b; top) as well as the one that includes SLUG (Fig. 2b; bottom). However, in presence of SLUG, an interesting feature was seen which was not present in the control case – a parameter region in which the only phenotype available to cells was the hybrid E/M one (Fig. 2b; bottom, gray region); thus, SLUG enabled the existence of a monostable hybrid E/M phenotype. Also, in the presence of SLUG, a higher level of EMT-inducing external signal I is required to push the cells to a completely mesenchymal state, representing the effect of SLUG in preventing cells from undergoing a full EMT.

To ascertain the robustness of this prediction, we performed sensitivity analysis where each kinetic parameter of the model was varied – one at a time – by ± 10%, and the corresponding change in the range of values of the external EMT-driving signal (I) enabling a hybrid E/M phenotype (i.e., the interval of x-axis between double-sided arrows shown at the bottom) was measured (online suppl. Fig. 2a). Except for a few parameter cases, the change in this interval of values of I was found to be less than 10% for a corresponding 10% change in the parameter values, indicating the robustness of the system. A change in few selected parameters such as those controlling ZEB self-activation or the effect of SLUG inhibiting SNAIL exhibited stronger sensitivity; but even in these cases, the decrease in the range of values of I for which a hybrid E/M existed was smaller as compared to the case in the absence of SLUG (dotted line in online suppl. Fig. 2a). Thus, SLUG can be thought of as a potential “phenotypic stability factor” for hybrid E/M phenotype.

We next sought to examine how this trait of enabling hybrid E/M phenotype changed when the strengths of interactions between SLUG-miR200 and SLUG-SNAIL were varied. While in the presence of SLUG, the hybrid state could exist in a monostable regime ({H}), it also coexisted with epithelial and mesenchymal states, thus giving rise to the multistable phases (i.e., combinations of coexisting phenotypes) {E,H}, {H,M} for varying strengths of external signal. On the other hand, neither the {E, H, M} phase nor the {E, M} one was seen in the circuit containing SLUG, as opposed to the control case (compare the top vs. bottom in Fig. 2b). First, we analyzed the effect of interaction between SLUG and miR200 on the dynamics (Fig. 2c). When the strength of repression of miR-200 by SLUG is reduced, 2 key changes are seen: (a) the {M} region expands while the {E} region shrinks, i.e., it is possible for cells to exit an epithelial phenotype and/or undergo a complete EMT at lower values of external signal I, and (b) the monostable hybrid region ({H}) disappears, and the {E, H, M} phase emerges (Fig. 2c; left). Similar trends, albeit to a weaker extent, were also seen when the threshold level of SLUG needed to repress miR200 is altered (online suppl. Fig. 2b). Further, to determine the effect of the reciprocal link (i.e., repression of SLUG by miR200), we varied the threshold levels of miR200 needed to repress SLUG. As this threshold value is increased (i.e., as the inhibition of SLUG by miR-200 is weakened), a larger {H} region is observed at the expense of shrinking {E} and {M} regions (Fig. 2c; right). Put together, these results suggests that as the influence of SLUG is reduced (i.e., weaker repression of miR-200 by SLUG, or stronger inhibition of SLUG by miR-200), and the existence of a monostable hybrid state ({H}) is disfavored. Conversely, when the effective impact of SLUG in modulating the EMP dynamics is increased (i.e., a stronger repression of miR-200 by SLUG, or weaker inhibition of SLUG by miR-200), the existence of a monostable {H} phase is facilitated.

Next, we quantified the effect of interactions between SLUG and SNAIL in altering the EMP dynamics (Fig. 2d). As the strength of repression of SNAIL by SLUG is decreased, the {E} phase shrinks and the {M} phase expands, i.e., when SNAIL is inhibited weakly by SLUG, a milder activation of SNAIL by the external signal I is sufficient to induce an EMT (Fig. 2d, left). However, no changes were observed when the threshold value of SLUG needed to repress SNAIL was altered (online suppl. Fig. 2c). To determine the effects of the reciprocal link (i.e., transcriptional inhibition of SLUG by SNAIL), we varied the strength of repression of the corresponding link, but did not observe any appreciable change in the phase boundaries (Fig 2d, right; online suppl. Fig. 2d).

SLUG is also known to self-activate in certain contexts [Kumar et al., 2015; Zhou et al., 2019]. Thus, we examined the changes in EMP dynamics upon including SLUG self-activation. We observed that irrespective of the link whose strength was being varied, including self-activation of SLUG increases the parameter region for the existence of monostable {H} region. For instance, the {H} phase was seen to persist even when SLUG repressed miR-200 weakly, and consequently, the emergence of the tri-stable {E,H,M} phase was not noticed in presence of SLUG self-activation (compare online suppl. Fig 2e with Fig. 2c, left). Similarly, transition into a completely mesenchymal state ({M}) required larger values of external signal I activating SNAIL (compare oline suppl. Fig 2f with Fig. 2d, left), leading to an expanded monostable hybrid ({H}) region. Therefore, SLUG self-activation may further enhance the likelihood of cells stably acquiring a partial EMT state and prevent them from transitioning to a fully mesenchymal state.

Finally, we replotted the abovementioned phase diagrams for a few varied values of other parameters pertaining to relative strengths of interaction of SLUG with miR-200 and/or SNAIL. Altering the degree of inhibition of SLUG by SNAIL did not significantly change the phase diagram (online suppl. Fig. 3a, b) as compared to the control case (Fig. 2d, left). However, increasing or decreasing the strength of repression of SNAIL by SLUG modified the phase diagram. When SLUG represses SNAIL strongly, a larger value of the external signal I (that activates SNAIL) is required to induce EMT (compare online suppl. Fig 3c with Fig. 2d, right). In contrast, when the repression of SNAIL by SLUG was weakened, the monostable epithelial and hybrid regions shrink, and a lower value of external signal (I) was sufficient to induce EMT (compare online suppl. Fig 3d with Fig. 2d, right). Because relatively major changes in the phase diagram were observed when the inhibitory action of SLUG was altered rather than the repression on SLUG being altered, we plotted these phase diagrams for 2 cases: (a) when SLUG inhibits both miR-200 and SNAIL strongly, and (b) when SLUG inhibits both miR-200 and SNAIL weakly. In the former case, the parameter region corresponding to a monostable hybrid E/M region ({H}) almost doubled, while the latter case led to a 4-fold decrease in the same (compare online suppl. Fig 3e, f with Fig. 2d, right).

Overall, these simulations suggest that the stronger the influence of SLUG in inhibiting SNAIL and/or miR-200 levels, the higher the likelihood of cells maintaining a hybrid E/M phenotype.

Previously identified “phenotypic stability factors” – GRHL2, OVOL1/2, and ΔNP63α – had another dynamical trait: their presence increased the mean residence time (MRT) of cells in a hybrid E/M phenotype. MRT is the mean (or average) time spent by the cells in a particular basin of attraction – E, M, and hybrid E/M here – as calculated via stochastic simulations [Biswas et al., 2019]. Thus, the larger the MRT of a phenotype, the higher the relative stability of the same. Hence, beyond enabling the parameter range of values of an external EMT-inducing signal for the existence of a hybrid E/M phenotype (as shown via bifurcation diagrams), an increase in MRT is usually thought of as a hallmark trait of a PSF. We investigated whether SLUG increases the MRT of cells in a hybrid E/M phenotype.

For the control case (i.e., circuit without SLUG), the epithelial state was found to be more stable in the {E,M} phase and the mesenchymal state is more stable in the {H,M} phase. All 3 states were found to be almost equally stable in the {E,H,M} phase (Fig. 3a). Cells were seen to be capable of transitioning among these different states under the influence of noise, depending on the level of the EMT-inducing external signal (I) (Fig. 3b). In the presence of SLUG, the hybrid state was found to be more stable than the epithelial state in the {E,H} phase and also in the {H,M} phase (Fig. 3c), and cells switched back and forth from hybrid E/M state to epithelial or mesenchymal state depending on the value of the external signal (I) (Fig. 3d). Thus, the presence of SLUG was capable of enabling a higher MRT for the hybrid E/M state.

However, in the {H, M} phase, as the levels of I increase, the MRT of mesenchymal state can be larger than that of a hybrid E/M state, as the system approaches a monostable mesenchymal state ({M}). These trends are also seen in potential landscape diagrams for corresponding cases (online suppl. Fig. 4).

Next, we deciphered the generic emergent properties of an extended gene regulatory network involving SNAIL, SLUG, ZEB1, miR-200, and E-cadherin (CDH1) (Fig. 4a). At a transcriptional level, ZEB1, SNAIL, and SLUG can all repress E-cadherin [Batlle et al., 2000; Mooney et al., 2016; Sterneck et al., 2020], and SNAIL and SLUG can activate ZEB1 [Wels et al., 2011; Wu et al., 2017]. Further, SLUG is known to self-activate [Kumar et al., 2015], but SNAIL self-inhibits [Peiró et al., 2006]. E-cadherin can also indirectly repress ZEB1 by regulating the localization of β-catenin [Mooney et al., 2016].

We used a computational framework RACIPE that enables mapping the general dynamic properties of gene regulatory networks by solving a set of coupled ODEs to obtain the possible steady-state solutions for a given network and corresponding ensemble of parameter sets. All the parameters and initial conditions required for solving the set of ODEs are randomly sampled from biologically relevant ranges [Huang et al., 2018]. We performed RACIPE analysis on this gene regulatory network (Fig. 4a) to find the robust steady-state solutions emerging from it, for an ensemble of parameter sets. We performed PCA on all the steady-state solutions produced by RACIPE (Fig. 4b). The PCA plot reveals the existence of 4 distinct clusters: epithelial (E), mesenchymal(M) and 2 hybrid clusters (H1 and H2), as confirmed by hierarchical clustering on all the steady-state solutions obtained via RACIPE (Fig. 4b). PC1 explains more than 65% of the variance in the data, and can be considered as the “canonical” EMT axis, given its composition that places CDH1 and miR-200 (epithelial markers) in one group and ZEB1, SLUG, and SNAIL (mesenchymal markers) in the other opposing one. Thus, E and M clusters are well-segregated along the PC1 axis, and 2 clusters potentially representing hybrid E/M states are sandwiched between the E and M ends. These 2 clusters (H1, H2) seem to be segregated along PC2 axis. Biologically, the hybrid cluster H1 can denote a SLUG + CDH1+ cell phenotype which has been documented recently experimentally as well (Fig. 4c, d) [Sterneck et al., 2020].

Next, we compared the effect of overexpression of SLUG versus SNAIL. Interestingly, these both scenarios showed opposite results: while SLUG overexpression enriched the SLUG⁺ CDH1⁺ phenotype, SNAIL overexpression depleted it (Fig. 4e–g), suggesting that SNAIL and SLUG may play different or even opposing roles in driving EMP, owing to the differences in underlying network topology. Furthermore, this enrichment of the SLUG⁺ CDH1⁺ population in the case of SLUG overexpression happens at the expense of the M cell population (Fig. 4g), lending credence to our previous results that SLUG is associated specifically with a partial EMT instead of a full EMT phenotype. Besides, this analysis endorses that network topology can be used to pinpoint potential PSFs for the hybrid E/M phenotypes [Jolly et al., 2016].

Hybrid E/M phenotype is thought of as more aggressive [Jolly et al., 2018; Pastushenko and Blanpain, 2019]. This hypothesis is further corroborated by the data showing poor patient survival for higher expression of PSFs such as GRHL2, OVOL1/2, and NRF2 [Jolly et al., 2016; Bocci et al., 2019]. To ascertain the effects of SLUG on patient prognosis, we investigated its role in patient outcomes. We observed that high SLUG expression associated consistently with poor relapse-free survival in lung, ovarian, and colorectal cancer (Fig. 5a–d), and with worse overall survival in lung, ovarian, colorectal, and pancreatic cancer (Fig. 5e–h; online suppl. Fig 5a-c), indicating that high SLUG levels can correlate with poor prognosis in many carcinomas.

In the case of breast cancer, the results were equivocal. High SLUG expression correlated with better overall patient survival (online suppl. Fig. 5d) and better metastasis-free survival (online suppl. Fig. 5e, f), but with poor relapse-free survival (online suppl. Fig. 5g, h). These context-specific differences of the association of SLUG with breast cancer prognosis traits remains to be better understood.

Cancer cell plasticity and metastasis are complex and emergent phenomena whose molecular underpinnings are yet to be comprehensively identified. An important axis of cellular plasticity is EMP. Various families of transcription factors have been known to affect EMP – TWIST1/2, ZEB1/2, and SNAIL1/2 can induce EMT [Lamouille et al., 2014], while OVOL1/2 and GRHL1/2/3 can induce MET [Saxena et al., 2020; Sundararajan et al., 2020]. Recent investigations have highlighted different and possibly even a contradictory role of different members of the same transcription factor family, for instance, the opposing roles of ZEB1 and ZEB2 in melanoma [Caramel et al., 2013]. Here, we examined the differences in the roles of SNAIL1 (encoded by the SNAI1 gene; also called SNAIL) and SNAIL2 (encoded by the SNAI2 gene; also called SLUG) in inducing EMT.

Our results suggest that SLUG, but not necessarily SNAIL, associates with a partial EMT or hybrid E/M phenotype. These observations are reminiscent of experiments suggesting that SNAIL is a stronger inducer of molecular and/or morphological features associated with EMT, as compared to SLUG and SNAIL3 [Gras et al., 2014]. Similar to other EMT-TFs, SLUG initiates EMT by repressing CDH1 by binding to the E-boxes [Hajra et al., 2002; Slabáková et al., 2011]; however, its binding affinity to the corresponding E-box DNA element was found to be two orders of magnitude less than that of SNAIL [Bolós et al., 2003]. Intriguingly, in TGFβ1-induced EMT in BPH-1 cells, SLUG expression was found to be decreased at later time points (72–96 h), but SNAIL and ZEB1 levels were consistently high [Slabáková et al., 2011], indicating a possible role of SLUG initiating rather than maintaining EMT [Savagner et al., 1997]. Consistently, SLUG levels were observed to peak in intermediate states of EMT in human mammary epithelial cells treated with TGFβ [van Dijk et al., 2018], endorsing the specific association of SLUG with a partially differentiated phenotype as witnessed in multiple scenarios involving EMT – embryonic development, wound healing, and carcinoma progression [Savagner et al., 2005; Côme et al., 2006; Leroy and Mostov, 2007; Hudson et al., 2009]. SLUG is also known to activate P-cadherin [Idoux-Gillet et al., 2018], another proposed marker for partial EMT [Ribeiro and Paredes, 2015].

The association of SLUG with partial EMT phenotype(s) also helps contextualize the growing literature of hybrid E/M phenotype being more stem-like as compared to a complete EMT [Jolly et al., 2019], and the functional role of SLUG in enabling stemness [Sterneck et al., 2020]. Mammary stem cells exhibit high SLUG levels, and SLUG-expressing normal mammary epithelial cells were seen to exhibit a partial EMT phenotype and generated 17 times more organoids than the control case [Guo et al., 2012; Ye et al., 2015]. Upon knockdown of SLUG, mammary stem cells are unable to transition to basal progenitor cells, leading to abnormal mammary architecture [Phillips et al., 2014]. The mammosphere-forming potential of basal progenitor cells also decreased upon SLUG knockdown, thus the primary mammospheres generated from SLUG-KO cells were unable to give rise to secondary mammospheres, highlighting the role of SLUG on self-renewal processes in mammary gland morphogenesis [Nassour et al., 2012]. Consistent observations have been made in breast cancer and prostate cancer, where SLUG degradation promotes the differentiation of breast cancer stem cells [Fraile et al., 2020], and the basal prostate cells that expressed SLUG exhibited a partial EMT phenotype and displayed increased stemness ability [Kahonouva et al., 2020]. The functional contribution of SLUG in vitro and/or in vivo enabling stemness is also observed in human epidermal progenitor cells [Mistry et al., 2014] as well as other cancers such as glioblastoma [Chesnelong et al., 2019], hepatocellular carcinoma [Sun et al., 2014; Tang et al., 2016], colorectal cancer [Kato et al., 2020], lung cancer [Kim et al., 2020], and squamous cell carcinomas [Yu et al., 2016; Moon et al., 2020]. Whether higher stemness is connected to a hybrid E/M phenotype in these cancers remains to be identified.

Our results indicate that the network topology that SLUG forms with SNAIL and miR-200 enables its role as a “phenotypic stability factor” for a hybrid E/M phenotype. It should be noted that the strengths of interactions of SLUG with SNAIL and/or miR-200 can be context dependent, thus altering its ability to be able to stabilize hybrid E/M phenotype(s). miR-200 family members have been shown to inhibit SLUG in esophageal cancer and prostate cancer [Liu et al., 2013; Zong et al., 2019; Basu et al., 2020], while SLUG can inhibit miR-200 in prostate cancer [Liu et al., 2013]. On the other hand, SLUG and SNAIL have been reported to inhibit each other in oral cancer [Nakamura et al., 2018] and breast cancer [Ganesan et al., 2015]. The effect of SLUG on EMP dynamics can also be influenced by upstream regulators of SLUG such as RNF8, YAP, C/EBPδ, miR-151, CBP, RCP, HDAC1, HDAC3, BRD4, and STAT3 [Daugaard et al., 2017; Hwang et al., 2017; Chesnelong et al., 2019; Cheng et al., 2020; Dai et al., 2020; Hu et al., 2020; Jury et al., 2020; Kuang et al., 2020; Shafran et al., 2020; Wang et al., 2020]. Moreover, SLUG may be involved in feedback loops with one or more p63 isoforms [Herfs et al., 2010; Dang et al., 2016; Srivastava et al., 2018] that can also enable a partial EMT [Jolly et al., 2017a; Westcott et al., 2020]. Thus, while the proposed role of SLUG in stabilizing hybrid E/M phenotype(s) is largely robust to parametric variations, including other links and/or nodes, it can alter the landscape of EMP dynamics. A better understanding of similarities and differences in terms of networks activated by SNAIL versus SLUG to influence the different states of EMT and associated traits such as drug resistance [Kurrey et al., 2009; Li et al., 2020], immune evasion [Tripathi et al., 2016], anoikis resistance [Huang et al., 2013], and autophagy [Xu et al., 2019] across different cancer types is therefore required.

This article has been uploaded on bioRxiv as a preprint: https://www.biorxiv.org/content/10.1101/2020.09.03.278085v1.

Not required for this type of study, as it does not involve any human subjects/animal models/identifiable human material and data. The study is purely computational in nature and any clinical data analyzed is publicly available in NCBI GEO.

The authors declare no conflicts of interest.

This work was supported by a Ramanujan Fellowship (SB/S2/RJN-049/2018) awarded to M.K.J. by the Science and Engineering Research Board (SERB), Department of Science and Technology (DST), Government of India.

M.K.J. designed and supervised research; A.R.S., S.S., and K.B. performed research and analyzed data. All authors contributed to writing and editing of the manuscript.