Cells interpret temporal information from TGF-β through a 1 nested relay mechanism 2

11 The detection and transmission of the strength and temporal quality of intracellular and 12 extracellullar signals is an essential cellular mechanism. While TGF-β signaling is one of the 13 most thoroughly studied signaling pathways, the mechanisms by which cells translate TGF-β 14 signals remain unclear. In this paper, through an integrated quantitative and computational 15 approach we demonstrate that crosstalk among multiple TGF-β activated pathways forms a relay 16 from SMAD to GLI1 that initializes and maintains SNAILl expression, respectively. This 17 transaction is smoothed and accelerated by another temporal switch from elevated cytosolic 18 GSK3 enzymatic activity to reduced nuclear GSK3 enzymatic activity. This nested relay 19 mechanism places SNAIL1 as a key integrator of information from TGF-β signaling 20 subsequently distributed through divergent pathways; essentially cells generate a transient or 21 peer-reviewed) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was not . http://dx.doi.org/10.1101/134106 doi: bioRxiv preprint first posted online May. 4, 2017;


Introduction
Cells live in as state of constant environmental flux and must reliably receive, decode, integrate 30 and transmit information from extracellular signals such that response is appropriate 1-4 . 31 Dysregulation of signal transduction networks leads to inappropriate transmission of signaling 32 information which may ultimately lead to pathologies such as cancer. Therefore a central  that each hormone has a single functional role. Since its discovery in 1980s 11 , researchers have 49 been puzzled by the enigma of "how cells read TGF-β signals", arising from seemingly 50 transient localization of pSMAD2/3 in the nucleus. pSMAD2/3 also promotes expression of 73 SNAIL1. Due to the double-negative feedback loops between SNAIL1 and miR-34 family 74 proteins, the expression of SNAIL1 is sustained after decay of nuclear pSMAD2/3. This 75 mechanism is supported by SMAD2/3 binding sites on the promoter region of snail1 18 . 76 Since quantifying SMAD proteins are the fundamental readouts of most current TGF-β signaling 77 studies, we set out to examine the downstream localization and abundance of SNAIL directly. 78 Human MCF10A cells were treated with recombinant TGF-β1 and multicolor 79 immunofluorescence (IF) was performed using antibodies directed against pSMAD2/3, SNAIL1. Importantly this is not cell type dependent as equivalent two-wave dynamics were seen for 89 SNAIL1 mRNA in MCF10A (Fig. 1D), MCF7 and A549 cells (Fig S1B), suggesting there may 90 be a secondary activator of SNAIL1. 91 To address the hypothesis that the activation of SNAIL1 is multifactorial and not solely 92 dependent on SMAD2/3 we treated MCF10A cells with a SMAD2/3 inhibitor LY2109761 in 93 addition to TGF-β treatment (Fig. 1D). When the inhibitor was added concurrently with TGF-β 94 treatment, the SNAIL1 mRNA level was reduced to ~ 9% of that of the control (no inhibitor) 95 network integrates multiple feed-forward loops that converge at the regulation of SNAIL1 119 transcription. We therefore hypothesized that at the point when the SMAD ceases to function as 120 a major transcription factor of SNAIL1, it has already induced GLI1 expression, and that the 121 concentration of GLI1 may be sustained through positive feedback and continued induction via 122 other impactful TGF-β signaling pathways. Essentially, a relay from pSMAD2/3 to GLI1 123 initiates and then maintains SNAIL1 expression. 124 To test this hypothesis, we performed microscopy studies of SNAIL1-GLI1 using MCF10A cells. 125 The distribution of SNAIL1 found in this study ( 2A also reproduced the temporal dynamics of pSAMD2/3 and SNAIL1 ( Figure S2C).

135
If GLI1 is involved in the maintenance of SNAIL1 expression subsequent to the drop in 136 pSMAD2/3 concentration, it is reasonable to predict (Fig. S2D) that inhibiting GLI1 activity, 137 either at the onset of or at some subsequent time after TGF-β treatment, would have minimal 138 effect on the initial wave of SNAIL1 expression since it is caused by pSMAD2/3. However, it 139 would eliminate the second wave of SNAIL1 expression. Indeed this was observed when the 140 GLI1 inhibitor GANT61 was added together with TGF-β at the beginning of the experiment 141 resulting in a reduction of the SNAIL1 mRNA level to be 55% (at 12 h and 24 h), 12% (at 48 h) 142 and 7% (at 72 h) compared to those without inhibition at the corresponding time points. In 143 another experiment adding the inhibitor 48 hours after TGF-β treatment also reduced the mRNA 144 level measured at 72 h to be 25% (Fig. 2D). These results are qualitatively different from those 145 with the SMAD inhibitor (Fig. 1D).

146
To confirm that GLI1 activation is not restricted to the MCF10A cell line, we also examined 147 MCF7 and A549 cells with TGF-β treatment, and observed similar increased and sustained GLI1 148 expression (Fig. S2E). Furthermore, early and late GLI1 inhibition lead to a reduction of the 149 SNAIL1 mRNA level to be 13% and 22% for MCF7 cells, and to a less extent of 57% and 66% 150 for A549 cells, respectively (Fig. S2F). Additionally, increased GLI1 expression after TGF-β 151 treatment has been found for multiple liver cancer cell lines 24 . In toto these results support the 152 potential role of GLI1 as a signaling relay from pSMAD2/3 to SNAIL1.

153
GSK3 in a phosphorylation form with augmented enzymatic activity accumulates at 154 endoplasmic reticulum and Golgi apparatus. 155 Next, we hypothesized that GSK3 is fundamental to the observed multi-phasic GLI1 dynamic 156 (Fig. 2). Most published studies suggest that GSK3 is constitutively active in untreated cells, 157 facilitating degradation of SNAIL1 and GLI1; TGF-β treatment leads to GSK3 phosphorylation 158 and inactivation, which leads to an accumulation of SNAIL1 and GLI1 25, 26 .

159
Initially we tested whether the above mechanism is sufficient to explain the multi-phasic GLI1 160 dynamics. We treated MCF10A cells in the absence of TGF-β with a GSK3 inhibitor that 161 suppresses GSK activity without interfering with its phosphorylation. Given the above 162 mechanism, one should expect the GSK3 inhibitor to promote both GLI1 and SNAIL1. In our 163 experiment, SNAIL1 did increase, but there was no change in GLI1 expression in either nucleus 164 or cytoplasm (Fig. 3A), suggesting additional signaling mechanisms may be involved. 165 Besides the inhibitory serine phosphorylation, the literature suggests that tyrosine (Y279 in 166 GSK-3α and Y216 in GSK-3β) phosphorylation leads to augmented enzymatic activity of GSK3 167 27 . As a convenience when discussing the three forms of GSK3, we refer the enzymatically active 168 unphosphorylated form and the more active tyrosine phosphorylated form as "GSK3 A " and 169 "GSK3 AA ", respectively, and the inactive serine phosphorylation form as "GSK3 D ". Also we 170 reserve "GSK3" for the total GSK3. As expected, microscopy studies showed an increased 171 concentration of GSK3 D peaking around 12 hours after TGF-β treatment (Fig. 3B, Fig. S3A).

172
Large cell-to-cell variations in the concentration of GSK3 D were observed, however, the 173 abundance of cytosolic and nuclear GSK3 D were essentially equivalent (the expression ratio was 174 close to one) for cells without TGF-β treatment (Fig. S3B). This observation corroborates earlier 175 report that the serine phosphorylation does not affect GSK3 nuclear location 28 . TGF-β treatment 2). Given that a function of active GSK3 is to modify target proteins post-translationally, our 183 observation suggests an unreported role for GSK3 AA accumulating at the ER and Golgi apparatus 184 is to modify newly synthesized proteins before their release to the cytosol. Specifically previous 185 studies showed that in mammalian cells a scaffold protein SUFU binds to GLI to form an 186 inhibitory complex, and SUFU phosphorylation by GSK3β prevents the complex formation, 187 exposes the GLI1 nuclear localization sequence 29 , which explains the observed increase of free 188 GLI1 in the cytosol followed by nuclear translocation (Fig. 2B). Since the two phosphorylation 189 forms, GSK3 AA and GSK3 D , coexist within single cells at defined time points, we performed co-190 immunoprecipitation and found that the probability of having the two GSK3 phosphorylation 191 forms in one molecule was undetectable (Fig. S3C).

192
Contrary to our observation that TGF-β regulates GSK3 AA dynamics, other studies posit that 193 GSK3 AA is not regulated by external cues 30 . To resolve this paradox, we measured the relative 194 amount of different GSK3 forms through silver staining (Fig. S3D). Among the three forms, the 195 overall percentage of GSK3 D increased from a basal level of 37% to 65% at 12 h after TGF-β 196 treatment. In contrast, only a small fraction of GSK3 molecules assumed the GSK3 AA form and 197 its overall abundance was stable over time (from ~10% basal level to ~13% at 8 h then back to 198 ~10% at 12 h after TGF-β treatment). Essentially levels of GSK3 AA did not change in abundance 199 but did change in localizations (homing to the ER and Golgi apparatus) to form a high local

233
To test the functional roles of GSK3 suggested above, we performed a series of GSK3 activity 234 inhibition experiments. First, we pretreated MCF10A cells with GSK3 inhibitor SB216763, 235 washed out the inhibitor then added TGF-β1 (Fig. S4B). We predicted that the treatment would 236 slow down GLI1 nuclear accumulation, and at later times decrease the overall increase of GLI1 237 and SNAIL1 compared to cells without GSK3 inhibitor. Indeed this was observed (Fig. 4C, adding TGF-β. In this case the inhibitor had opposite effects on the GLI1/SNAIL1 protein 245 concentration: it slowed down the initial release and translocation of GLI1 needed to accelerate 246 the GLI1 accumulation, but also decreased GLI1 and SNAIL1 degradation that becomes pre-247 eminent when the proteins were present at high levels. Compared to the samples grown in the 248 absence of the GSK3 inhibitor, we also observed slower and more scattered GLI1 nuclear 249 accumulation and SNAIL1 increase on day 2, but by day 3 the overall levels of GLI1 and 250 SNAIL1 were actually higher than the case without the inhibitor (Fig. 4C, TGF-β + GSK3_I).

251
The SMAD-GLI1 relay increases the network information capacity and leads to 252 differential response to TGF-β duration 253 Our results show that TGF-β1 signaling is effected through pSMAD2/3 directly with fast pulsed 254 dynamics concurrently with a relay through GLI1 which has a much slower dynamics. The 255 signaling ported by these two channels converges on SNAIL1 with a resultant two-wave 256 expression pattern. To further dissect the potential functional interactions between these two 257 pathways, we performed mathematical modeling and predicted that the two distinct dynamics 258 allows cells to respond to TGF-β differentially depending on stimulus duration (Fig. 5A). Short 259 pulses of TGF-β only activate pSMAD2/3 and the first wave of transient SNAIL1 expression.

260
When the signal duration is longer than a defined threshold value, activation of GLI1 will lead to 261 the observed second wave of SNAIL1 expression. We confirmed the predictions with MCF10A   The two-wave dynamic of TGF-β-induced SNAIL1 expression has been observed in several 298 cellular systems 21, 33 , supporting the underlying relay mechanism discovered in this work. The The functional switch from pSMAD2/3 to GLI1 relays information from TGF-β signaling 320 beyond the initial induction of SNAIL1, and this relay is facilitated by a second relay from the 321 active to the inactive phosphorylation form of GSK3 proteins. Active regulation of the 322 abundance and nuclear location of GSK3 AA form has been observed in neurons 36 . In contrast to 323 these earlier reports we observed an accumulation of GSK3 AA in the ER and Golgi apparatus.    global strategy was used to identify the nuclear shape, and the Otsu algorithm was used for further 452 calculation. Clumped objectives were identified by shape and divided by intensity. Next, using the shrank 453 nuclear shape as seed, cell shape was identified by the Watershed algorithm. For identifying the clusters 454 of GSK3 AA formed around a nucleus, the nuclear shape was shrank manually by 3 pixels and used as a 455 new seed to grown the boundary with the watershed method until reaching background intensity level. All 456 parameters were optimized through an iterative process of automatic segmentation and manual inspection. were quantified by the average pixel intensity within the nucleus or cytosol region of a cell. Next, the 461 quantified results were examined manually, and those cells with either cell area, nuclear area, or 462 fluorescent intensity beyond five folds of the 95% confidence range of samples from a given treatment 463 were discarded, which account for less than 1% of the cells analyzed. Immunofluorescence data were 464 further processed and plots were generated using customized R codes and Matlab codes. The R package 2 , 465 was also used in data analysis. were added into every 100 l of bead mixture respectively. The mixture was rotated at 4 °C for 3 hours. 480 Beads that were conjugated with antibodies were washed with PBS_Tween. An amount of 100 l of pre-481 cleaned lysis buffer was added into conjugated beads and rotated at 4 °C overnight. Targeted proteins 482 were eluded from beads by incubating with 40 l 1x Laemmli buffer with SDS at 70 °C for 10 minutes. 483 For the samples an amount of 5 l was used for western blot assay, and an amount of 30 l was loaded 484 for SDS-PAGE (Bio-Rad) and followed by silver staining (Fisher). 485

Network reconstruction and coarse-graining 486
The full network from TGF-1 to SNAIL1 (Fig. S2A) was generated with IPA (Qiagen®). Specifically, 487 all downstream regulators of TGF-1 and upstream regulators of SNAIL1 in human, mice and rat were 488 searched and added to the network. Then, direct or indirect relationships between every pair of regulators 489 were searched and added to the network. After obtaining the whole network, regulators that have been 490 reported to be activated later than SNAIL1 were removed. Examination of the network reveals that the 491 network can be further organized into three groups: the TGF--SMAD-SNAIL canonical pathway, the 492 TGF--GSK3--catenin pathway that has the most number of links, and others. We further noticed that 493 GLI1 is a central connector of TGF-, SMAD, GSK3 and SNAIL1. We performed western blot studies 494 on -CATENIN and found that neither its concentration nor its location changes significantly before day 495 3, therefore we removed -CATENIN from the network. In addition, previous studies report that the 496 SMAD-GLI axis plays important role in TGF- induced EMT 23 . Therefore we further grouped the 497 network as the SMAD module, the GLI module, and the GSK3 module, as well as the remaining ones that 498 we referred as "Others", and reached the network shown in Fig. 2A. Those molecular species not 499 explicitly specified in Fig. 2A either have their effects implicitly included in the links, for example the 500 link from TGF- to GSK3, or are included in the links of "Others". This treatment is justified since our 501 various inhibition experiments indeed showed that the three factors we identified affect SNAIL1 502 expression the most. These "other" species may contribute to snail1 activation at a time later than what 503 considered in this work. Therefore we emphasize the network in Fig. 2A is valid only within the time 504 window we examined, i.e., within three days after TGF-1 treatment for MCF10A cells.        Canonical TGF-β/pSMAD2/3/SNAILl pathway (Fig. 1A) 646 We used this ordinary differential equation (ODE) model in Fig. 1 and Fig. S1. 647 where [GLI] n is the concentration of nuclear GLI1. We used this ODE model to generate results in Fig.  669

671
Since the process involves many steps and a detailed model would require many parameters to determine, 672 instead we used two phenomenological time-dependent functions to qualitatively mimic the dynamics of 673 the enzyme activities of cytosol GSK3 and nuclear GSK3 we experimentally measured, 674 [GSK] c (t) = k GSKc * TGF * (1 − exp (− t a1 )) * exp (− t − b1 a1 ) , 675 [GSK] n (t) = 1 − k GSKn * TGF * (1 − exp (− t a2 )) * (exp (− t − b2 a2 )) . 676 Furthermore, the basal pool of cytosol GLI1 is considered, which is by sequestered in the cytosol by Sufu 677 but could translocate to the nuclear after Sufu is inactivated by the cytosol enzyme GSK3 activity. We 678 used a revised ODE of nuclear GLI1 concentration derived with the quasi-equilibrium approximation (see 679 below) 680  We used this ODE model to generate results in Fig. 5 and Fig. S5A. 726 Parameter space searching 727 Step 1: Calculate single cell distributions of experimental observables. We calculated histograms of the 728 distributions from the single cell experimental data. Suppose that we have N observables measured in M 729 time points, we have an N × M dimensional distribution of the data. Since we used fixed cells and we had 730 no information on the temporal correlation, we treated the distributions from different time points as 731 independent, i.e., = ∏ =1 . 732 Step 2: Define pseudo-potentials from the parameterized distribution. We defined a pseudo-scalar-