The CD44s/ZEB1 feedback loop enables the miR-200/ZEB1 circuit to function as a three-way switch
The TCS model incorporates a direct ZEB1 self-activation link based on the contribution of ZEB1 in stabilizing the SMAD complexes; this link can reinforce ZEB1 levels that are raised via the TGF-β pathway (Hill et al., 2013). Here, we include in the model recently published evidence regarding a self-reinforcing feedback loop of ZEB1 via ESRP1 and CD44s (Preca et al., 2015) (Fig. 2a), and show that this feedback loop indeed enables the tristability of the miR-200/ZEB1 circuit (see Additional file 1: SI Section 3 for the model formulation of the CD44s/ZEB1 feedback loop and Additional file 1: Table S1 for relevant parameters).
CD44, a crucial stem cell marker, has two major isoforms – CD44s (the CD44 standard isoform) and CD44v (CD44 variant isoform) (Ponta et al., 2003). Whereas the total amount of CD44 is maintained largely unchanged during EMT, the isoform switch from CD44v to CD44s is essential for the progression of EMT (Brown et al., 2011; Zhao et al., 2016), and CD44s can upregulate ZEB1 (Preca et al., 2015). This isoform switch is regulated by a splicing factor - epithelial splicing regulatory protein 1 (ESRP1) - which promotes the splicing of CD44v and inhibits that of CD44s (Brown et al., 2011). The transcription of ESRP1 can be directly inhibited by ZEB1, and therefore ZEB1 can upregulate itself by promoting the production of CD44s (Preca et al., 2015; Brown et al., 2011).
To analyze the effect of the CD44s/ZEB1 feedback loop on the behavior of the miR-200/ZEB1 circuit, we calculated a bifurcation diagram (Fig. 2a) to illustrate the existence of and transitions among the different stable states. The CD44s/ZEB1 feedback loop enables the miR-200/ZEB1 circuit to acquire three stable states (phenotypes) – (low ZEB1, high ESRP1) corresponding to the E phenotype, (medium ZEB1, medium ESRP1) corresponding to the hybrid E/M phenotype and (high ZEB1, low ESRP1) corresponding to the M phenotype. The hybrid E/M phenotype is not seen upon removal of the CD44s/ZEB1 feedback loop (Additional file 1: Figure S1), demonstrating that self-activation in a toggle switch is critical for attaining more than two stable states (Huang et al., 2007; Zhou & Huang, 2011; Lu et al., 2013). To understand the diverse EMT-inducing results (e.g., by hypoxia, which can upregulate both SNAIL and ZEB1 or by TGF-β, which mainly upregulates SNAIL), we calculated the phase diagram where the miR-200/ZEB1 circuit is driven by two independent signals – a transcriptional inhibition signal S
1 on miR-200 and a transcriptional activation signal S2 on ZEB1 (Fig. 2b). Different combinations of S
1 and S
2 allow 7 possible phases – 3 monostable phases - {E}, {E/M} and {M}, 3 bistable phases - {E/M, E}, {E/M, M} {E, M} and 1 tristable phase {E, E/M, M} (Fig. 2b). The bistable and tristable phases imply the co-existence of multiple stable states. For example, the tristable phase {E, E/M, M} implies that all three phenotypes – E, E/M and M can exist in the given range of S
1 and S
2, i.e. cells can switch back and forth among these three phenotypes. Depending on the details of EMT induction, cells tend to follow different trajectories in the phase diagram and thereby undergo different phenotypic transitions. The tristability of the miR-200/ZEB1 circuit enabled by the CD44s/ZEB1 feedback loop is quite robust to parameter perturbation (Additional file 1: SI Section 4, Additional file 1: Figure S2).
Next, to investigate the gene expression levels of ZEB1 across E, hybrid E/M and M cells, we analyze the NCI-60 panel of cell lines that have been classified into E, E/M, and M phenotypes on a population level based on ensemble CDH1/VIM ratio (Park et al., 2008). CDH1 encodes the epithelial marker E-cadherin, loss of which is a hallmark of EMT. VIM encodes vimentin, which is a commonly used marker for the mesenchymal phenotype. Both ZEB1 and ESRP1 show three different expression levels across the E, E/M and M cell lines, and the difference among these levels is statistically significant (Fig. 2c). These results are consistent with the prediction from the TCS model, which proposes that ZEB1 (Fig. 2a, b) and ESRP1 levels (Additional file 1: Figure S3) in the hybrid E/M phenotype are different from that of the E phenotype. Further, across the entire NCI-60 panel, ESRP1 expression negatively correlates with ZEB1 and VIM, and positively correlates with CDH1 (Jia et al., 2015; Roca et al., 2013; Watanabe et al., 2014). These observations, together with the reported strong negative correlation between the expression of ZEB1 and ESRP1 in the lung, breast and pancreatic cancer patient samples (Preca et al., 2015), suggest that the ZEB1-ESRP1 axis is active across multiple cancer types (Fig. 2c).
In addition to the expression analysis of ZEB1 from the NCI-60 cell line panel, we further conducted experiments in three NSCLC cell lines – H820, H1975 and H1299 - that have been characterized as epithelial, hybrid E/M and mesenchymal respectively (Schliekelman et al., 2015). Notably, the hybrid E/M H1975 cells stain for both E-cadherin and vimentin at a single-cell level and display collective migration (Schliekelman et al., 2015; Jolly et al., 2016). The immunofluorescence staining of the H1975 cells clearly shows the expression of ZEB1 in the nucleus concomitantly with E-cadherin on the membrane, which demonstrates that ZEB1 is present in the hybrid E/M phenotype as compared to the lack of ZEB1 in the epithelial H820 cells (Fig. 2d, Additional file 1: Figure S4). The intermediate levels of ZEB1 in hybrid H1975 cells are further verified by RT-PCR (Fig. 2e) and Western blot (Fig. 2f, Additional file 1: Figure S5). Notably, the difference of SNAIL levels (3 fold) is much smaller compared with the difference of ZEB1 (20 fold) between epithelial cells - H820 and H1437 and hybrid E/M cells – H1975 (Fig. 2e, f), suggesting ZEB1 levels may correspond strongly with maintaining these phenotypes, at least in these cell lines. Interestingly, the levels of SNAIL in mesenchymal cells – H1299 and H2030 were lower than those in epithelial cells - H820 and H1437 (Fig. 2e, f), further arguing for a more critical role of ZEB1 in maintaining a mesenchymal phenotype, potentially via multiple feedback loops (Preca et al., 2015; Gregory et al., 2011).
These results indicate that ZEB1 self-activation (direct or indirect) should be considered an integral part of the core EMT circuit, thereby implying that TCS model captures the biological mechanisms more precisely than the CBS one, at least in the contexts discussed above. The parameter region of the bifurcation and phase diagram (Fig. 2a, b) here is quite consistent with that for the miR-200/ZEB1 circuit with direct ZEB1 self-activation (see Fig. 4 in (31) for comparison). For simplicity, the CD44s/ZEB1 feedback loop is represented by a direct ZEB1 self-activation in the following analysis.
Next, we focus on published experimental data that has been claimed to support the CBS model over the TCS model. These experimental results are – (a) during TGF-β treatment of MCF10A cells, the concentration of exogenous TGF-β required to increase the SNAIL abundance is lower than that needed for a comparable increase in ZEB1 abundance (Zhang et al., 2014), and (b) MCF10A cells that attain a hybrid E/M phenotype revert to being epithelial upon removal of exogenous TGF-β, but cells that attain a mesenchymal phenotype fail to do so (Zhang et al., 2014). In the following two sections, we show that both the CBS and the TCS models can recapitulate these two experimental observations.
Different responses of SNAIL and ZEB1 to the EMT-inducing signal- exogenous TGF-β
To test whether the TCS model can recapitulate the different TGF-β-dose responses of SNAIL and ZEB1, we apply an EMT-inducing signal S
A
on SNAIL to mimic the induction of exogenous TGF-β, and analyze the steady-state responses of SNAIL and ZEB1 to different levels of S
A
. To investigate how differently SNAIL and ZEB1 respond to TGF-β induction at a single-cell level compared with that at a cell population level, we design two different kinds of stochastic simulations. First, at a single-cell level (Fig. 3a), in absence of a detailed understanding of the gene expression noise in regulating EMT, we consider the effects of white noise. Secondly, to investigate population-averaged results (Fig. 3b), we include cell-cell variability, which is represented by randomly assigned parameters to each cell.
The TCS model results show that the levels of TGF-β (S
A
) required for the activation of SNAIL (S
A
= 20∗103 molecules) is indeed lower as compared to that for ZEB1 (S
A
= 40∗103 molecules) on both the single-cell (Fig. 3a-c) and population levels (Fig. 3d-f). The difference in the upregulation of SNAIL and ZEB1 at different levels of TGF-β might be attributed to the different numbers of binding sites on 3’UTR of mRNA for the corresponding microRNAs – 6 or more binding sites of miR-200 family on ZEB1 mRNA and 2 binding sites of miR-34 family on SNAIL mRNA (Kim et al., 2011; Gregory et al., 2008). The different endogenous levels of miR-34 and miR-200 family members may also affect the difference in ZEB1 and SNAIL levels. In addition, tristability is more apparent for the ZEB1 mRNA levels than for the ZEB1 protein levels. Again, this can be attributed to the microRNA-mediated repression. Therefore, with intermediate levels of ZEB1 mRNA and miR-200 present in the hybrid E/M phenotype, the ZEB1 protein does not necessarily reach the intermediate level, because of the strong repression by miR-200; this results in the relatively less separated E and E/M phenotypes for ZEB1 protein levels. The fraction of each group with the ZEB1 mRNA levels - low, medium, high - can be adjusted by the levels of TGF-β (S
A
) (Additional file 1: Figure S6).
Next, we compare the responses of ZEB1 on the single-cell and population-averaged levels; we find that ZEB1 mRNA levels display a more significant trimodal distribution in the single-cell simulation compared with that from the population-averaged result (Fig. 3c). The cell-to-cell variability, as reflected by parameter randomization, tends to smoothen the distributions of ZEB1 mRNA in the E and the hybrid E/M phenotypes, potentially because the difference of ZEB1 mRNA levels between E and E/M is smaller as compared to that between E/M and M, as reflected by the simulations of single cells undergoing EMT (Fig. 3a). The difference of ZEB1 mRNA levels in E, E/M and M phenotypes shown here should not be compared with the gene expression data from NCI-60 (Fig. 2c), which is not time-course data for individual cells undergoing EMT. Since the existing experimental data on population measurements only yield average values, this may help explain why ZEB1 mRNA appears to be bimodally distributed (Zhang et al., 2014). Consequently, the averaged ZEB1 levels may not be able to argue strongly in favor of one mathematical model over the other.
The inhibition of miR-34 by ZEB1 as well as autocrine TGF- β signaling stabilizes the mesenchymal phenotype
To decode possible mechanisms of the irreversible nature of the transition from hybrid E/M to M, we focus on the feedback loop from the miR-200/ZEB1 circuit to the miR-34/SNAIL circuit – the inhibition of miR-34 by ZEB1 (Fig. 4a). To understand the effect of this link on the transitions among E, E/M and M, we calculate the bifurcation diagram of the complete circuit – miR-34/SNAIL/miR-200/ZEB1 - driven by EMT-inducing signal S
A
on SNAIL for varying strengths of the inhibition of miR-34 by ZEB1 (Fig. 4b). The stronger the inhibition of miR-34 by ZEB1, the smaller the level of EMT-inducing signal S
A
required for the cells to transition into and maintain a M phenotype (Fig. 4b). In addition, a stronger inhibition of miR-34 by ZEB1 can decrease the duration of the hybrid E/M phenotype and can therefore promote a quicker transition from the hybrid E/M phenotype to the M phenotype during temporal dynamic simulation (Additional file 1: Figure S7).
To characterize the relative stability of E, E/M and M states, we calculate the effective landscape of EMT at different strengths of the inhibition of miR-34 by ZEB1. Here, the external activation signal S
A
is chosen at 50*103 molecules, which enables the existence of all three stable states – E, E/M and M in all cases with different strengths of the inhibition of miR-34 by ZEB1 (Fig. 4b). When there is no feedback from ZEB1 to miR-34, cells are mainly bimodally distributed in either the E or the E/M phenotype (left panel in Fig. 4b). An intermediate strength of the inhibition of miR-34 by ZEB1 enables all three phenotypes to be attained (middle panel in Fig. 4b). When the strength of the inhibition of miR-34 by ZEB1 is further increased, the M phenotype becomes the dominant one, thus cells in this case can maintain the mesenchymal phenotype without reverting (right panel in Fig. 4b). Therefore, a strong inhibitory feedback from ZEB1 to miR-34 stabilizes the mesenchymal phenotype and contributes to the irreversible transition from the hybrid E/M to the M phenotype.
Another possible mechanism accounting for the ‘irreversibility’ of the mesenchymal phenotype can be the autocrine TGF-β/miR-200 loop, as suggested by the CBS model. Here, we evaluate the effect of the TGF-β/miR-200 loop based on the TCS framework (Fig. 5a) and show that removal of exogenous TGF-β cannot induce MET but instead the circuit is able to maintain the mesenchymal phenotype (Phase {M} in Fig. 5b) as long as the innate production rate of endogenous TGF-β is high. This bifurcation diagrams (Fig. 5c, d) indicate that the ‘irreversibility’ of mesenchymal phenotype to switch back to a hybrid E/M phenotype depends directly on the production rate of endogenous TGF-β.
In summary, both the TCS and the CBS models can explain currently published experimental results regarding EMT – (a) EMT is a two-step process, from E to E/M to M. (b) The concentration of exogenous TGF-β required for the activation of SNAIL is lower than that for ZEB1, (c) The mesenchymal phenotype can be maintained without reverting when exogenous TGF-β is removed due to certain levels of endogenous TGF-β. In addition, we show that the inhibition of miR-34 by ZEB1 can be responsible for the irreversibility of a complete EMT. These EMT mathematical models are helpful to the extent that they can continue to explain and predict the regulatory effects of newly identified EMT players. To this end, we focus on adding recently identified EMT regulator – FOXC2, which plays a crucial role in EMT progression, to both the TCS and the CBS frameworks and test whether these two models can elucidate the function of FOXC2 in regulating EMT.
Overexpression of SNAIL is insufficient to initiate EMT in absence of ZEB1 and FOXC2
The transcription factor FOXC2 serves as a key mediator in regulating EMT and linking EMT with stem-like properties and with metastatic competence (Mani et al., 2007; Hollier et al., 2013; Paranjape et al., 2016; Werden et al., 2016). FOXC2 expression can be upregulated by multiple EMT-inducing signals, such as TGF- β1, Snail, Goosecoid and Twist. It has been shown that FOXC2 expression is required to maintain the mesenchymal phenotype, the invasive properties and the stem cell-enrichment of HMLE cells following EMT induction (Mani et al., 2007; Hollier et al., 2013).
FOXC2 regulates EMT through its interactions with core EMT components – SNAIL and ZEB1. Overexpression of SNAIL significantly upregulates the expression of FOXC2 while overexpression of FOXC2 does not affect SNAIL levels (Hollier et al., 2013; Tisza et al., n.d.) (Additional file 1: Figure S8). This suggests that SNAIL functions as an upstream regulator of FOXC2. In addition, FOXC2 directly upregulates the expression of ZEB1 by binding to its promoter region (Werden et al., 2016). Here, we have expanded both the TCS and CBS models to include transcriptional regulation by FOXC2 (Fig. 6A, Additional file 1: Figure S9A, see Additional file 1: SI Section 5&6 for model formulation and Additional file 1: Table S2&3 for parameters).
To analyze the EMT-inducing behaviors of the TCS model in the presence and absence of FOXC2, we study the effect of two external signals - S
A
representing an EMT-inducing signal on SNAIL, and S
I
representing an inhibitory signal on FOXC2. The presence of FOXC2 (S
I
= 0) enables tristability of EMT and accounts for the two-step transition - from E to E/M and from E/M to M (Fig. 6b, top panel). The absence of FOXC2 (S
I
= 2∗105 molecules) results in low levels of ZEB1 mRNA and consequently the maintenance of epithelial phenotype (Fig. 6b, bottom panel) irrespective of the high levels of SNAIL (Additional file 1: Figure S10).
To test the above-mentioned predictions from the TCS model, we examined the immortalized HMLE cells to assess the role of FOXC2 knockdown (FOXC2-KD) during EMT. Here we measure the protein levels of canonical EMT markers in the HMLE-Snail cells, that have already undergone a complete EMT via SNAIL overexpression, in the presence and absence of FOXC2. We find that FOXC2-KD in HMLE-Snail cells eliminates the expression of ZEB1, vimentin and fibronectin while it restores the expression of E-cadherin, thus inducing a complete MET irrespective of SNAIL overexpression (Fig. 6c). This experimental observation is consistent with the prediction from the TCS model, but not with the prediction from the CBS model that an EMT-inducing signal can still drive the transition from an epithelial state to a hybrid E/M state upon FOXC2-KD (Additional file 1: Figure S9).
Is it possible that FOXC2 acts together with SNAIL and together can induce (partial) EMT even without ZEB1? This would argue against TCS and would instead be consistent with a modified version of CBS. We do not think this is supported by existing data. Observations in LNCaP and DU145 cells show that ZEB1 mediates the effect of FOXC2 on tumor-initiating potential and drug resistance (Paranjape et al., 2016), traits that are often correlated with EMT (Jolly et al., 2014; Jolly et al., 2015; Mani et al., 2008; May et al., 2011; Morel et al., 2008). Decreased expression of ZEB1 in PANC-1 cells, which express both epithelial and mesenchymal markers (Gradiz et al., 2016) and are thus likely to be hybrid E/M cells, results in a complete MET despite upregulation of SLUG and SNAIL in TGF-β treated cells (Abshire et al., 2016). Moreover, knockdown of ZEB1, but not necessarily of SNAIL and SLUG, had pronounced effects in cells losing their EMT-like properties (Weitzenfeld et al., n.d.). Therefore, we believe that the most consistent interpretation of the data is that SNAIL (even with FOXC2) is insufficient without ZEB1 and conversely both ZEB1 and FOXC2 are needed (Fig. 6D) and contribute to different aspects of driving EMT such as repression of epithelial program (ZEB1 is a transcriptional repressor (Spaderna et al., 2008) and activation of mesenchymal program (FOXC2 is a transcriptional activator (Werden et al., 2016). However, further experiments such as overexpression of ZEB1 in FOXC2 knockdown and overexpression of FOXC2 in ZEB1 knockdown cells will be crucial in further delineating the mechanisms of EMT dynamics and distinguishing the principles underlying EMT tristability.