ERK basal and pulsatile activity are differentially regulated in mammalian epidermis to control proliferation and exit from the stem cell compartment

Fluctuation in signal transduction pathways is frequently observed during mammalian development. However, its role in regulating stem cells has not been explored. Here we tracked spatiotemporal ERK MAPK dynamics in human epidermal stem cells. While stem cells and differentiated cells were distinguished by high and low stable basal ERK activity, respectively, we also found cells with pulsatile ERK activity. Transitions from Basalhi-Pulselo (stem) to Basalhi-Pulsehi, Basalmid-Pulsehi, and Basallo-Pulselo (differentiated) cells occurred in expanding keratinocyte colonies and in response to a range of differentiation stimuli. Pharmacological inhibition of ERK induced differentiation only when cells were in the Basalmid-Pulsehi state. Basal ERK activity and pulses were differentially regulated by DUSP10 and DUSP6, leading us to speculate that DUSP6-mediated ERK pulse downregulation promotes initiation of differentiation whereas DUSP10-mediated downregulation of mean ERK activity promotes and stabilizes post-commitment differentiation. Quantification of MAPK1/3, DUSP6 and DUSP10 transcripts in individual cells demonstrated that ERK activity is controlled both transcriptionally and post-transcriptionally. When cells were cultured on a topography that mimics the epidermal-dermal interface, spatial segregation of mean ERK activity and pulses was observed. In vivo imaging of mouse epidermis revealed a patterned distribution of basal cells with pulsatile ERK activity and downregulation was linked to the onset of differentiation. Our findings demonstrate that ERK MAPK signal fluctuations link kinase activity to stem cell dynamics. Significance Understanding how intracellular signaling cascades control cell fate is a key issue in stem cell biology. Here we show that exit from the stem cell compartment in mammalian epidermis is characterised by pulsatile ERK MAPK activity. Basal activity and pulses are differentially regulated by DUSP10 and DUSP6, two phosphatases that have been shown previously to regulate differentiation commitment in the epidermis. ERK activity is controlled both transcriptionally and post-transcriptionally. Spatial segregation of mean ERK activity and pulses is observed both in reconstituted human epidermis and in mouse epidermis. Our findings demonstrate the tight spatial and temporal regulation of ERK MAPK expression and activity in mammalian epidermis.

showed significantly lower ERK pulse frequency and higher basal activity compared ERK activity profiles (Fig. S2 A and B). As expected, during NHK colony growth, the 163 number of cells expressing Involucrin-mCherry increased and they tended to have 164 lower ERK basal activity and pulse frequencies ( Fig. 2E and Fig. S2C). The reduction 165 in pulse frequency, however, did not show a strong correlation with Involucrin-mCherry 166 expression levels in individual cells. (Fig. 2E). 167 To further dissect the change in ERK activity associated with differentiation we 168 analysed the trajectories of ERK activity (mean and variance) and Involucrin 169 expression over time (Fig. 2F). We observed co-evolution of ERK activity variance and 170 Involucrin expression, showing that there is a strong tendency for undifferentiated cells 171 to downregulate pulse frequencies coupled with differentiation (Fig. 2G). In contrast 172 there was a gradual convergence of basal ERK activity towards a low level as 173 differentiation proceeded (Fig. 2H). This difference suggests that ERK pulses and 174 mean levels are subject to distinct regulatory mechanisms, and that the 175 downregulation in ERK pulses has a role in switching stem cells to differentiated cells. 176 The trajectory of ERK dynamics in Basal hi -Pulse hi (Day3 after plating) cells showed 177 that the cells tend to decrease pulse frequency and increase mean activity, which will 178 lead to Basal hi -Pulse lo ERK activity (Fig. S2D). This suggests that the ERK dynamics 179 is partially reversible in the early stage while it is irreversible once they are committed 180 to differentiation. 181 Epidermal keratinocytes not only transition from the stem cell compartment to the 182 differentiation compartment but can also transition between cell division and 183 quiescence(14). To determine whether ERK pulse activations were associated with 184 cell division, we followed the fate of Basal hi -Pulse lo cells (Day2 after plating) and 185 Basal hi -Pulse hi (Day3 after plating) that did not subsequently express Involucrin (Fig. 186 1E). Basal hi -Pulse hi cells had a high probability of dividing, whereas Basal hi -Pulse lo 187 cells did not (Fig. 2I). 188 We also recorded the ERK pulse frequency of individual cells and whether or not they 189 subsequently divided (within 48 hours). We found that 70% of cells with low ERK pulse 190 frequency (< 1.5 pulse/h) remained as single cells (Fig. 2J, left), whereas 60% of cells 191 with high ERK pulse frequency (> 1.5 pulse/h) proliferated, giving rise to 2 or more 192 cells during the recording period (Fig. 2J,right). This indicates that cells with a pulsatile 193 ERK profile are more likely to divide than cells with a stable-high ERK profile. 194 These results, combined with our observations of cells expressing  suggest that stem cells transition from Basal hi -Pulse lo to Basal hi -Pulse hi when 196 proliferating, and that a subsequent transition to Basal mid -Pulse hi is associated with 197 differentiation, leading to Basal lo -Pulse lo ERK activity in differentiated cells (Fig. 2K). 199 We next tested the effect of different differentiation stimuli on ERK activity dynamics.  12-O-Tetradecanoylphorbol-13-acetate (TPA) is known to stimulate Involucrin 211 expression(31) and also increases ERK activity(32). When cells were stimulated with 212 10 ng/ml TPA, they exhibited highly pulsatile ERK activity (Fig. 3 C-F). Overall ERK 213 pulse frequency increased in a dose dependent manner up to 40 ng/ml TPA (Fig. 3G), 214 while mean ERK levels peaked at 10 ng/ml (Fig. 3H). The time-course analysis 215 revealed that TPA induction of ERK pulses was transient, peaking at 9 hr after the 216 start of treatment ( Fig 3I). As before (Fig. 2 G and K), the onset of Involucrin expression 217 coincided with the downregulation of ERK pulses (Fig. 3J).

218
Like TPA, epidermal growth factor (EGF) transiently enhanced ERK pulse levels (Fig. 219 3 K and L), although, in contrast to TPA, EGF decreased overall ERK pulse frequency 220 in a dose dependent manner without affecting ERK mean levels ( Fig. 3 L and M). ERK 221 pulses peaked at 11 hr after EGF treatment and then showed a significant decrease, 222 again coincident with the onset of Involucrin expression (Fig. 3 O and P). Thus TPA and EGF had different effects on overall ERK pulse frequency (Fig. 3 G and M) but 224 the time course of changes is similar and in both cases ERK was downregulated when 225 cells expressed Involucrin (Fig. 3 I and J and O and P). This leads us to speculate that 226 the downregulation of ERK pulses triggers differentiation (Fig. 2K).

227
When human keratinocytes were treated with a MEK inhibitor PD0325901 (MEKi), 228 they exhibited a dose-dependent reduction in both ERK pulse frequency and basal 229 levels ( Fig. 3 Q and R and Fig. S4 A-D). We tested the inhibitor on different days after   (Fig 4F). This suggests the intriguing possibility that DUSP6-mediated ERK pulse 257 downregulation promotes the initiation of differentiation whereas DUSP10-mediated 258 downregulation of mean ERK activity promotes and stabilizes post-commitment 259 differentiation. This is consistent with the finding that DUSP6 is transiently upregulated 260 on commitment, while DUSP10 upregulation is more sustained(23). 261 We also tested siRNA-mediated knockdown of DUSP6 and DUSP10(23). Reduced 262 DUSP6 expression did not have a significant effect on ERK pulse frequency or basal 263 activity ( Fig. 4 G and H). However, reduction in DUSP10 expression led to increased 264 basal ERK activity while maintaining pulse frequency (Fig. 4 G and H). This confirms 265 that DUSP10 controls basal ERK activity.

266
As a further means of perturbing ERK activity we lentivirally transfected HNK cells with Together, β1-integrin, EGF and its downstream effectors underlie the ERK dynamics 287 transitions to achieve different cellular outcomes (Fig. 4Q)  and C). 309 We observed a patterned distribution of ERK activity on the substrates. Cells on the     to that in suprabasal differentiated cells (Fig. 6N). This rules out the possibility that the Institute of Kanazawa, Japan, and was subcloned into pCSII vector. To construct pLenti-430 Involucrin-mCherry, the sequence comprising the full length 3.7kb Involucrin promoter, 431 Involucrin intron, and SV40 splice donor and acceptor (SDSA) was subcloned from a 432 previously-reported beta-galactosidase reporter(54) into pLenti backbone, and mCherry was 433 tagged to the carboxy (C) terminus of the promoter sequence. Doxycycline-inducible DUSP 434 expression plasmids (pCW57-DUSP6 and pCW57-DUSP10) were previously described(23) 435 except that GFP was removed from the vector for compatibility with EKAREV. Doxycycline-436 inducible MEK EE expression plasmid (pCW57-MEK EE ) was constructed by subcloning 437 constitutive MEK1 mutant, MEK EE into pCW57 vector (55) . 438

Mice 439
Transgenic mice expressing EKAR-EVnls were obtained from Laboratory Animal Resource Patterned PDMS substrates were pre-coated with rat tail collagen type I (Corning) at 20 µg/ml 464 for 3 hours. Cells were seeded in complete FAD medium at a density of 75,000 cells cm -1 for 465 45 minutes at 37°C. Substrates were rinsed gently once with FAD medium to remove non-466 adherent cells and transferred to 6 cm dishes containing inactivated J2-3T3 cells seeded at a 467 density of 20,000 cells/cm 2 . 468 EGF, TPA, and PD0525901 were added to the medium 6 hours prior to live imaging at final 469 concentrations of 1 µg/ml and 1 µM, respectively. 470

Live imaging of human keratinocytes 507
2D culture images were acquired on a Nikon A1R laser scanning confocal microscope with 508 GaAsp detectors using a 20x Plan Apo VC 0.75 NA objective (Nikon) and NIS-Elements 509 (Nikon). Live cells were imaged in a temperature-controlled chamber (37°C) at 5% CO2. For 510 nuclear staining, Hoechst 33342 was added to the culture medium 30 min prior to imaging at 511 a final concentration of 5 µg/ml. Images were acquired every 5 min for up to 24 hours. For 512 FRET imaging, the EKAREV biosensor was excited by a 445 nm laser, and 482/35 nm and 513 525/50 nm emission filters were used to acquire CFP and FRET images, respectively. 514 Involucrin-mCherry was excited by a 560 nm laser and detected by a 595/50 nm emission 515 filter. The Hoechst 33342 signal was excited by a 405 nm laser and detected by a 450/50 nm 516 emission filter. 517 Images for cells on PDMS substrates were acquired on a Nikon A1R confocal/multiphoton 518 laser scanning microscope 25x Apo LWD 1.1 NA objective. Mineral oil was used to cover the 519 surface of the medium and prevent evaporation. Live cells were imaged in a temperature-520 controlled chamber (37°C) and pH was maintained with 15 mM HEPES buffer. Images were 521 acquired every 15 min. 522

Single cell proliferation assay 523
Mitotically inactivated J2-3T2 feeder cells were plated on collagen-coated 384 well glass 524 bottom plates in complete FAD medium at the density of 20,000 cells/cm 2 . Single NHKs were 525 plated onto each well. 6 hours after plating, only wells that accommodate single cell were 526 subjected to live imaging to measure ERK pulse levels for 6 hours. Cells were incubated for 527 another 48 hours and the same wells were revisited by microscope to observe cell numbers. 528

In vivo imaging of mouse epidermis 529
Methods for live imaging of mouse epidermis were previously described (11). Briefly, the ear 530 skin of anaesthetized mice was depilated and stabilized between a cover glass and thermal 531 conductive silicon gum sheet. In vivo imaging of tail epidermis was performed similarly. 532 Depilated tail was flanked with two silicon gum sheets and a cover glass was placed on the 533 top. In vivo live imaging was performed by ZEISS 7MP multi-photon microscope, equipped 534 with W Plan-Apochromat 20x/1.0 DIC VIS-IR M27 75mm water-immersion objective lens and 535 Coherent Chameleon Ti:Sapphire laser. EKAR-EVnls signal was detected by BP 500-550 and 536 BP 575-610 for CFP and FRET, respectively.

FRET analysis 538
Single cell ERK activity was measured by ratiometry of CFP and FRET signals (FRET/CFP) 539 since EKAREV functions by intramolecular FRET and the molecular number of the two 540 fluorescence proteins are considered to be equal. Each signal level was measured by mean 541 pixel intensity in individual cell areas. 542

Automated cell tracking 543
Tracking was performed by script-based operation of a FIJI plugin, Trackmate 544

Semi-automated single-cell tracking of cells expressing EKAR-EVnes 558
Cells expressing EKAR-EVnes cultured on the patterned PDMS substrate (Fig. 4D) were 559 tracked in a semi-automated manner with a custom-made program for FIJI/ImageJ. Cellular 560 centre locations were tracked by eye based on the lack of nuclear signal of EKAR-EVnes. 12-561 pixel square regions were automatically created around each centre and Huang's fuzzy 562 thresholding was applied to obtain cytoplasmic regions expressing EKAR-EVnes. The mean 563 CFP was FRET signals were obtained for each region and used for ERK activity (FRET/CFP). 564

Segmentation of tip and base areas on the patterned PDMS substrate 565
Tip and base regions were demarcated by circles with a diameter of 200 μm centred at each 566 tip (Fig. 4D, white dotted circles). ERK high and ERK low cells in base areas (troughs) were gated 567 by the 1.2 value of FRET/CFP (Fig. 4 E-G).

Instantaneous variance as a measure of population ERK pulse level 569
To quantify the level of ERK activity pulses of a population of cells at a given time point we 570 measured the variance of ERK activity (FRET/CFP) at that time point among all cells 571 (instantaneous variance). An increase in the instantaneous variance indicates a higher 572 variability of ERK activity in the population at a specific time point. Pulses were observed to 573 behave as stochastic events, as suggested by the exponential distribution of interpulse 574 intervals (Fig. S1D). This, combined with the fact that the mean ERK activity did not change 575 between conditions (Fig. 4A, Control/DUSP6) justifies the use of the instantaneous variance 576 as a measure of the level of ERK activity pulses of a population at a specific time point. When 577 the instantaneous variance remains unchanged between conditions (Fig. 4B, 578 Control/DUSP10) we say that the level of ERK activity pulses are the same for both 579 populations, regardless of changes in the mean level. 580

Moving variance as a measure of ERK pulse level in a time window 581
In order to study the change in ERK activity pulses over time in individual cells, we analysed 582 overlapping moving time windows of 50 minutes. Each time window was small enough for the 583 mean ERK activity, within the window, to be considered fixed, but long enough to 584 accommodate an ERK activity pulse (typical pulse ~0.25 hr). For each window we computed 585 the variance of ERK activity. The variance is a measure of dispersion of the measurements in 586 the window, and quantifies the extent of the deviations of the signal from its mean value. These 587 deviations can occur due to pulses, or noise (which is of much smaller amplitude than pulses). 588 The variance captures both the amplitude and number of pulses in the time window, giving a 589 quantitative measure of the pulsing level of ERK activity in the time window. The minimum 590 value for the variance is zero, which corresponds to a signal without pulses or fluctuations; 591 larger values of the variance indicate a higher level of pulsation. 592 This method allows a continuous assessment of the pulse level over time. Other methods, 593 such as peak detection and pulse count, amount to a discrete measurement of pulses, which 594 does not lend itself naturally to a detailed quantitative analysis of the temporal evolution of 595 ERK activity pulses. 596

Phase diagram 597
To construct the phase diagram of ERK activity variance and Involucrin mean level, cells co-598 expressing EKAREV-nls, and the Involucrin reporter, Involucrin-mCherry were considered. The values of variance of ERK activity and mean Involucrin level were plotted for every time 603 window, providing a trajectory in the plane spanned by ERK activity variance and mean 604 Involucrin level. This was repeated for all cells. 605 The ERK activity variance vs Involucrin mean level plane was then divided into regular blocks. 606 The trajectories that lay within each block were then averaged to obtain a mean direction for 607 each block. This procedure resulted in the phase diagram of ERK activity variance vs mean 608 Involucrin level. Every block in the phase diagram corresponds to a pair of ERK activity 609 variance/Involucrin mean level values, while the arrow in the block indicates the mean 610 direction to which these values changed in time. 611 The same procedure can be followed to construct other phase diagrams, such as mean ERK 612 activity vs mean Involucrin level. 613

Phase diagram normalization and transition probabilities 614
Each arrow in the phase diagram was decomposed into its x and y components. The x 615 components were rescaled by the maximum value of the Involucrin mean level, i.e. 616 x'=x/max(Inv. Mean level), while the y components were rescaled by the maximum value of 617 the ERK activity variance, i.e. y'=y/max(ERK activity variance). This rescaling amounted to 618 normalising both axes of the phase diagram to the range [0,1], and allowed the comparison 619 between the x' and y' components of the arrows. 620 The transition probabilities between neighbouring blocks corresponded to rx=x'/(|x'|+|y'|) and

Cluster analysis 628
To analyse the spatial organisation of cells in ear and tail epidermis we measured the radial 629 distribution function (RDF) g(r), which measure the deviations of the density of cells from that 630 of a random distribution, as a function of the distance r from a cell of reference (57). To 631 construct the RDF we constructed rings of radius r and width Dr around every cell of interest, 632 and counted total number N(r, Dr) of cells that lie within the rings. We then constructed the 633 reference number Nref(r, Dr), which was computed by measuring N(r, Dr) for a set of randomly 634 distributed point in the field of interest. The field of interest was constructed by considering 635 only regions of space where cells were observed in the experiment, as indicated by the black 636 outline in Fig. 6 G and M, and S5 C and D. This was done to prevent artefacts due to boundary 637 effects that might bias the clustering results. Finally, we constructed the RDF as the ratio g(r) 638 = N(r, Dr)/Nref(r, Dr). 639 We constructed the null RDF gnull(r)= Nnull(r, Dr)/Nref(r, Dr), against which the experimental 640 measurements were be compared. In this case, both Nnull(r, Dr) and Nref(r, Dr) were measured 641 for random distributions of points. In the null model, all particles are independent, randomly 642 distributed points, hence the gnull(r) is equal to unity. From 50 realisations of gnull(r), we 643 computed the 95% confidence intervals. When the experimental RDF lied above (below) the 644 95% c.i. the cells were considered clustered (dispersed) for that particular distance r, around 645

cells. 646
The RDF was also used to study the clustering between two distinct populations, high and low 647 pulsing cells. For this, we considered the cells of one group (high pulsing) as the reference 648 cells around which the rings are constructed, while we counted the number of cells of the 649 second group (low pulsing) that lie within the rings. In this case, when the experimental 650 observations lie below (above) the 95% c.i. both cell populations are considered segregated 651 (grouped). 652

Statistics 653
Statistical analyses were performed using MS Excel or MATLAB R2018a (Mathworks) 654 Software. We made use of the two-tailed Student's t-test for unpaired data to quantify 655 differences between experimental groups. Kolmogorov-Smirnov test was used to compare