Since keratinocytes increase in size as they differentiate (Adams and Watt, 1989), we enriched for undifferentiated cells by filtration prior to placing them in suspension for up to 12 hr (Figure 1a). By determining when cells recovered from suspension could no longer resume clonal growth on replating (Figure 1b; Figure 1—figure supplement 1a), we confirmed that there is a marked drop in colony forming ability between 4 and 8 hr. This correlates with an increase in the proportion of cells expressing the terminal differentiation markers involucrin (IVL) and transglutaminase 1 (TGM1) (Figure 1c,d; Figure 1—figure supplement 1b) and downregulation of genes that are expressed in the basal layer of the epidermis, including integrin α6 (ITGα6) and TP63 (Figure 1d) and the stem cell markers DLL1 and Lrig1 (Tan et al., 2013; Figure 1—figure supplement 1c). As expected, Ki67 expression was also reduced in suspension, reflecting the drop in proliferation upon differentiation (Figure 1—figure supplement 1c).

Figure 1 with 1 supplement see all

Download asset Open asset

To define commitment at the molecular level, we next collected keratinocytes after 4, 8 and 12 hr in suspension and performed genome-wide transcriptomics, using Illumina-based microarrays, and proteome-wide peptide analysis, by SILAC-Mass-Spectrometry (MS) (Figure 1e,f; Figure 1—figure supplement 1d–g; Supplementary file 1, 2). Keratinocytes collected immediately after trypsinisation served as the 0 hr control. When comparing the starting cell population (0 hr) with cells suspended for 4, 8 or 12 hr, t-SNE plots of differentially expressed genes (Figure 1e) and unsupervised hierarchical clustering of genes and proteins (Figure 1f; Figure 1—figure supplement 1d) indicated that the 4 hr samples clustered separately from the 8 and 12 hr samples. GO term enrichment analysis of differentially expressed transcripts or proteins showed enrichment of terms associated with epidermal differentiation at 8 and 12 hr (Figure 1—figure supplement 1e,g) (Sen et al., 2010; Mulder et al., 2012), which is consistent with the drop in clonogenicity seen at these time points.

When we compared the significantly differentially expressed proteins (p-value<0.05) that changed ≥2 fold at one or more time points with their corresponding transcripts (Figure 1g–i), there was a moderately positive correlation at 8 and 12 hr (Pearson correlations of 0.51 and 0.68, respectively), consistent with the correlation between bulk mRNA and protein levels seen in previous studies of mammalian cells (Schwanhäusser et al., 2011; Ly et al., 2014). However, at 4 hr transcripts and proteins were only weakly correlated (Pearson correlation 0.19, Figure 1g).

The poor correlation between protein and transcript levels at 4 hr suggested a potential role for post-transcriptional mechanisms in regulating commitment. To investigate this, we performed unbiased SILAC-MS-based phospho-proteomic analysis. SILAC-labelled peptides isolated from cells at 0, 4 and 8 hr time points were enriched for phosphopeptides using HILIC pre-fractionation and titanium dioxide affinity chromatography (Supplementary file 3). Over 3500 high confidence phosphorylation sites were identified with an Andromeda search score >= 30 and quantified at both the 4 and 8 hr time points. At 4 hr approximately two thirds of the changes involved dephosphorylation (Figure 1j), but these dephosphorylation events could not be attributed to decreases in protein abundance, as shown by the discordance between total protein abundance and changes in protein phosphorylation (Figure 1k).

We next identified protein phosphatases that were differentially expressed between 0 and 4 hr. We intersected the proteomic and genomic datasets, revealing 47 differentially expressed phosphatases, 22 of which were upregulated (Figure 2a). Interrogation of published datasets revealed that these phosphatases were also dynamically expressed during calcium-induced stratification of human keratinocytes (Hopkin et al., 2012) and differentiation of reconstituted human epidermis (Lopez-Pajares et al., 2015) (Figure 2—figure supplement 1a,b).

Figure 2 with 4 supplements see all

Download asset Open asset

To examine the effects on keratinocyte self-renewal of knocking down each of the 22 phosphatases identified in the screen, we transfected primary human keratinocytes with SMART pool siRNAs and measured colony formation in culture (Figure 2b; Figure 2—figure supplement 2a). Live-cell imaging was used to monitor cell growth for 3 days post-transfection (Figure 2—figure supplement 2b). Knocking down seven of the phosphatases significantly increased clonal growth (Figure 2b), with five phosphatases – DUSP6, PPTC7, PTPN1, PTPN13 and PPP3CA - having the most pronounced effect (p-value<0.001). Silencing of these five phosphatases also increased colony size and delayed the loss of colony forming ability in suspension (Figure 2c,d; Figure 2—figure supplement 2c–d). Consistent with these findings, phosphatase knockdown delayed the decline in TP63 and increase in TGM1 levels during suspension-induced differentiation (Figure 2e,f; Supplementary file 4). Knocking down each phosphatase did not have a major effect on the growth rate of keratinocytes and live-cell imaging did not reveal any apoptosis (Figure 2—figure supplement 2e). Cumulatively, the effects on keratinocyte self-renewal and differentiation suggest that DUSP6, PPTC7, PTPN1, PTPN13 and PPP3CA are pro-commitment protein phosphatases.

The effects of knocking down DUSP10 at 0 hr differed from those of the other five phosphatases. DUSP10 knockdown significantly reduced clonogenicity (p-value<0.001; Figure 2b; Figure 2—figure supplement 2a) and, in contrast to the other phosphatases, decreased the growth rate of keratinocytes (Figure 2—figure supplement 2b,e). In addition, whereas expression of DUSP6, PPTC7, PTPN1, PTPN13 and PPP3CA declined by 8 hr in suspension, DUSP10 expression remained high (Figure 2a). Taken together, their effects on keratinocyte self-renewal and differentiation suggest that all six are commitment-associated protein phosphatases, with DUSP10 differing from the other phosphatases in potentially antagonising commitment.

To examine the effect of knocking down the pro-commitment phosphatases on the ability of keratinocytes to reconstitute a multi-layered epithelium, we seeded cells on de-epidermised human dermis and cultured them at the air-medium interface for three weeks (Figure 2g; Figure 2—figure supplement 2f,g). Knockdown of DUSP6, PTPN1, PPP3CA and PTPN13 did not prevent cells from undergoing terminal differentiation, as evidenced by suprabasal expression of involucrin and accumulation of cornified cells (Figure 2g; Figure 2—figure supplement 2f). However, the number of TP63-positive cells was increased, whether measured as a proportion of the total number of cells (Figure 2g) or per length of basement membrane (Figure 2—figure supplement 2g). Knockdown of DUSP6, PTPN1 or PPP3CA increased epidermal thickness without increasing the proportion of proliferative, Ki67-positive cells (Figure 2g; Figure 2—figure supplement 2g). Conversely, PTPN13 knockdown led to an increase in the percentage of Ki67-positive cells and a slight reduction in epidermal thickness (Figure 2g), which could reflect an increased rate of transit through the epidermal layers. In addition, Ki67 and TP63-positive cells were no longer confined to the basal cell layer of epidermis reconstituted following phosphatase knockdown but were also present throughout the viable suprabasal layers.

To complement our studies with siRNAs, which achieve transient knockdown, we performed stable knockdown of the pro-commitment phosphatases with shRNA lentiviral vectors (two different shRNAs per target; Figure 2—figure supplement 3a) and performed epidermal reconstitution experiments (Figure 2—figure supplement 3b,c). We also developed a tool for automated, unbiased measurements of epidermal thickness (Figure 2—figure supplement 4, Source code 1). No apoptotic cells were observed in reconstituted epidermis, as evaluated by lack of Caspase-3 labelling. Consistent with the effects of siRNA-mediated knockdown, stable knockdown of DUSP6, PPTC7, PTPN1 and PPP3CA increased epidermal thickness, whereas PTPN13 knockdown did not (Figure 2—figure supplement 3b,c). We also confirmed expansion of the TP63+ compartment (Figure 2—figure supplement 3d). Virtually all Ki67+ cells co-expressed TP63 (Figure 2—figure supplement 3d,e). Furthermore, many suprabasal cells co-expressed TP63 and IVL (Figure 2—figure supplement 3f). Thus, on both transient and stable knockdown, the transition from the stem to the differentiated cell compartment was disturbed (Figure 2—figure supplement 3f).

To identify the signalling networks affected by upregulation of phosphatases during commitment, we performed GO analysis of ranked peptides that were dephosphorylated at 4 hr. The top enriched pathways were ErbB1 signalling, adherens junctions, insulin signalling and MAPK signalling (Figure 3a). Several of the proteins we identified are components of more than one pathway (Figure 3b) and all have been reported previously to regulate epidermal differentiation (Connelly et al., 2010; Haase et al., 2001; Kolev et al., 2008; Trappmann et al., 2012; Scholl et al., 2007). In particular, constitutive activation of ERK delays suspension-induced differentiation (Haase et al., 2001).

Figure 3 with 1 supplement see all

Download asset Open asset

We next ranked protein phosphorylation sites according to the log₂-fold decrease at 4 hr, plotting the ratio between the change in phosphorylation sites and the change in total protein (Supplementary file 5). To specifically identify dephosphorylation events, we excluded from the ranking phosphorylations that remained constant while total protien abundance increased by more than 0.5 in log₂. Consistent with the predicted dynamic interactions between signalling pathways (Figure 3b), phosphorylation sites on MAPK1 (ERK2) and MAPK3 (ERK1) were identified in the top 15 most decreased sites (Figure 3c). Other proteins in the top 15 included components or regulators of the cytoskeleton (FLNA), Rho signalling (DOCK5, ARHGEF16, CIT) and EGFR signalling (EPS8), again consistent with the GO terminology analysis.

We performed western blotting to confirm that ERK1/2 activity was indeed modulated by suspension-induced terminal differentiation and by the candidate pro-commitment phosphatases (Figure 3d,e). As reported previously, the level of phosphorylated ERK1/2 diminished with time in suspension (Janes et al., 2009) (Figure 3d). siRNA-mediated knockdown of DUSP6, PPTC7, PTPN1, PTPN13 or PPP3CA resulted in higher levels of phosphorylated ERK1/2 relative to the scrambled control, both at 0 hr and at some later time points (Figure 3d,e). These effects are consistent with the known requirement for ERK MAPK activity to maintain keratinocytes in the stem cell compartment (Trappmann et al., 2012).

Transcriptional regulation of epidermal differentiation is mediated by the Activator Protein 1 (AP1) family of transcription factors (Eckert et al., 2013), which is the main effector of the MAPK and ErbB signalling cascades (Chang and Karin, 2001). Quantification of the levels of AP1 transcripts during suspension-induced terminal differentiation revealed that different AP1 factors changed with different kinetics, as reported previously (Gandarillas and Watt, 1995) (Figure 3f). Notably, the level of several members of the MAF subfamily of AP1 factors (MAF, MAFB and MAFG) significantly increased during differentiation, consistent with recent evidence that they mediate the terminal differentiation program in human keratinocytes (Lopez-Pajares et al., 2015). In line with these observations, knockdown of individual pro-commitment phosphatases reduced the induction of MAF AP1 factors in suspension (Figure 3f; Supplementary file 6). These experiments are consistent with a model whereby induction of phosphatases in committed keratinocytes causes dephosphorylation of ERK MAPK and prevents the increase in expression of AP1 transcription factors that execute the terminal differentiation program.

Given that DUSP10 knockdown differed from the other protein phosphatases in decreasing the clonogenicity (Figure 2b; Figure 2—figure supplement 2a) and growth rate of keratinocytes (Figure 2—figure supplement 2b,e), we examined the effect of DUSP10 knockdown on ERK1/2 activity and AP1 factor expression. Unlike knockdown of DUSP6, PPTC7, PTPN1, PTPN13 or PPP3CA, knockdown of DUSP10 (Figure 3—figure supplement 1a,b) did not increase ERK1/2 activity (Figure 3g,h). Consistent with its known regulation of p38 MAPK (Caunt and Keyse, 2013), DUSP10 knockdown increased phospho-p38 at 0 hr; however, this effect was not observed at later times in suspension (Figure 3—figure supplement 1c,d). When the effects of DUSP10 knockdown on AP1 transcription factors were compared with those of DUSP6, it was clear that these phosphatases differentially regulated several AP1 factors, including members of the JUN, FOS and MAF subfamilies (Figure 3i; Supplementary file 7).

We confirmed the differing effects of DUSP10 and DUSP6 on cell fate by using either Cumate or Doxycycline inducible overexpression in human keratinocytes (Figure 3j; Figure 3—figure supplement 1e–g). Whereas overexpression of DUSP10 increased colony formation, overexpression of DUSP6 reduced clonal growth (Figure 3j). A dominant negative mutant of DUSP6 (C293S) lacking phosphatase activity (Okudela et al., 2009) had no effect (Figure 3j; Figure 3—figure supplement 1f,g).

We next examined whether the six phosphatases we identified experimentally could form a single phosphatase-phosphatase interaction network governing commitment and differentiation. To do so, we employed a Boolean network abstraction and constructed a model in which all phosphatases, DUSP6, PPTC7, PTPN1, PTPN13, PPP3CA and DUSP10, were allowed to interact either positively or negatively (Figure 4—figure supplement 1a). For this we utilised the concept of an Abstract Boolean Network (ABN) (Yordanov et al., 2016), which by assigning definite or possible interactions, implicitly defines approximately 10¹⁷ concrete different Boolean network models (Figure 4—figure supplement 1a). In accordance with this modelling formalism, we abstracted discretised gene expression levels as ‘on’ or ‘off’ if their mean expression was respectively higher or similar/lower than the average of all genes (Supplementary file 8). From the set of all possible Boolean networks (the ABN), we sought to test whether there would be a network consistent with the discretised gene expression states over the differentiation time course. To this end, we employed the automated reasoning approach encapsulated in the tool RE:IN, as described by Dunn et al., 2014. We encoded a set of experimental constraints comprising the different time points analysed during suspension-induced differentiation (Figure 1a). These constraints define whether each gene should be on or off at each time step along the differentiation trajectory (Materials and methods). Using this approach, we formally deduced that a single Boolean network was unable to recapitulate the measured gene expression dynamics (Figure 4—figure supplement 1a).

As we were not able to find a single interaction network sufficient to model the discretised gene expression states, we further investigated each individual time point to identify whether the network topology might differ over the time course in suspension. To obtain the necessary experimental data to identify critical network interactions, we performed individual knockdowns of the phosphatases after 0, 4, 8 or 12 hr in suspension. Next, we measured by RT-qPCR the effect that each individual knockdown had on the mRNA level of any of the other phosphatases (Figure 4—figure supplement 2a) at each time point. We used the fold-change of the interaction to inform its directionality (activating or inhibitory) only if it was statistically significant (p-value<0.05, Supplementary file 9). This led to the inferred mechanistic networks depicted in Figure 4a, where the node colours show fold-change in each phosphatase with time in suspension, relative to the 0 hr time-point. Arrows indicate positive effects on expression and T-bars show inhibitory effects. These networks reveal different interactions at play at each of the time points, suggesting that the network topology reconfigures during commitment and differentiation.

Figure 4 with 3 supplements see all

Download asset Open asset

In light of these inferred mechanistic interactions, we again employed a Boolean network abstraction to explore the dynamics of epidermal stem cell commitment. Using the discretised gene expression levels defined above (Supplementary file 8), we sought to test whether the network topologies derived for each time point (Figure 4a,b) were consistent with the data. For this we used an extension of RE:IN, called RE:SIN, which allows the user to explore switching networks, as described by Shavit et al. (2016). RE:SIN allowed us to determine whether the inferred networks could recapitulate the dynamic changes in gene expression by similarly encoding these as expected states on the commitment trajectory (Figure 4—figure supplement 1b). We first found that the networks inferred by the genetic knockdown experiments could not satisfy the experimental constraints. This suggested that there were additional interactions between the phosphatases. To identify such interactions, we analysed the effect that each individual knockdown had on the protein levels of the other phosphatases (Figure 4—figure supplement 2d; Supplementary file 10) and the effect of DUSP6 and DUSP10 overexpression on phosphatase mRNA levels (Figure 4—figure supplement 2e). We included those as possible interactions, again using the concept of an ABN at each time step (Figure 4—figure supplement 1c; Supplementary file 8).

By defining an ABN at each time point that incorporates the possible phosphatase interactions deduced from the knockdown experiments, we found that the experimental constraints could be satisfied. Furthermore, our analysis revealed which of the possible interactions were required to meet these constraints. These could therefore be considered as definite interactions (Figure 4c, compare with Figure 4—figure supplement 1c). We thus used the approach encapsulated in RE:SIN to uncover the networks governing differentiation at each time point.

The Boolean network analysis indicates a key role for DUSP6 at commitment (4 hr), when DUSP6 expression is positively regulated by all the other phosphatases in the network and by a self-activating loop (Figure 4a). In contrast, at all other time points the interactions between individual phosphatases were predominantly negative. Several phosphatases, including PTPN1 and PTPN13, showed positive autoregulation at one or more time point. It is also notable that DUSP10 was predominantly negatively regulated by other phosphatases, except at 4 hr, when it was positively regulated by PTPN1 and by an autoregulatory loop (Figure 4a).

To test the robustness of the network, we examined the effects of treating keratinocytes with the histone deacetylase inhibitor Trichostatin A (TSA) and a Protein Kinase C inhibitor (PKCi). Both drugs blocked suspension-induced terminal differentiation, as measured by expression of IVL and TGM1, and prevented downregulation of ITGα6 (Figure 4—figure supplement 2b,c). However, cells treated with TSA still underwent commitment, as evaluated by loss of colony-forming ability and downregulation of TP63, whereas those treated with the PKC inhibitor did not. At 12 hr in suspension, both drugs reduced expression of DUSP10 and PTPN13 and increased expression of DUSP6 and PPTC7 relative to untreated cells; however, they differentially affected PPP3CA and PTPN1 (Figure 4d).

We next used our modelling approach to test whether network switching still occurs under PKCi and TSA treatment using the constrained ABNs that we derived for each stage. As before, we discretised gene expression levels from these experiments as ‘on’ or ‘off’ (Supplementary file 8). Under both treatments, we found that we could not transit through the networks while respecting the expected expression states, corroborating the experimental findings that terminal differentiation is blocked in these conditions. The PKCi phosphatase expression pattern at 12 hr was compatible with the phosphatase interaction network derived at 0 hr, supporting the conclusion that PKCi arrests cells in the stem cell compartment. In TSA-treated cells the network must switch from that derived at 0 hr to that derived at 4 hr, and subsequently to that derived at 8 hr, or straight from 0 hr to 8 hr. Neither PKCi nor TSA treatment resulted in an expression pattern compatible with the network derived at 12 hr for untreated cells, consistent with the inhibition of differentiation.

The above findings can be summarised using dynamical systems terminology. Epidermal differentiation can be described by two saddle-node bifurcations (Figure 4e,f). We start with a single minimum that corresponds to a stable state: the stem cell. Then we pass through a first saddle-node bifurcation where we have two minima, so the system can switch from one state to the other and this corresponds to commitment. Next the global minimum changes and we finally have a second saddle-node bifurcation that leads again to a single steady state, which corresponds to the differentiated cells (Zhang et al., 2012). From here, the stem cell state and the terminally differentiated state emerge as stable states, while commitment is inherently unstable and serves as a biological switch.

Given the role of the protein phosphatases in regulating commitment of keratinocytes in suspension, we examined whether phosphatase expression was also subject to spatial regulation, since spatiotemporal coordination of stem cell behaviour contributes to epidermal homeostasis (Doupé and Jones, 2012). By wholemount labelling the basal layer of sheets of human epidermis (Jensen et al., 1999), we found that DUSP6, PTPN1 and PPP3CA were most highly expressed in cells with the highest levels of β1 integrin or α6 integrin, which correspond to stem cell clusters (Jensen et al., 1999) (Figure 4g; Figure 4—figure supplement 3a,b). In contrast, PPTC7 was enriched in the integrin-low regions, while PTPN13 and DUSP10 were more uniformly expressed throughout the basal layer (Figure 4g; Figure 4—figure supplement 3a,b). The inverse relationship between the patterns of PPP3CA and PPTC7 expression is in good agreement with the network analysis demonstrating negative regulation of PPP3CA by PPTC7 in undifferentiated keratinocytes (0 hr; Figure 4a).

The patterned distribution of phosphatases could be recreated in vitro by culturing keratinocytes on collagen-coated PDMS elastomer substrates that mimic the topographical features of the human epidermal-dermal interface (Viswanathan et al., 2016). We observed clusters of cells with high levels of DUPS6 on the tops of the features, where stem cells expressing high levels of β1 integrins accumulate (Viswanathan et al., 2016) (Figure 4h). In contrast, DUSP10 was uniformly expressed regardless of cell position (Figure 4h). These results indicate that the phosphatases are subject to spatial regulation that is independent of signals from cells in the underlying dermis.

Primary human keratinocytes (strain km) were isolated from neonatal foreskin and cultured on mitotically inactivated 3T3-J2 cells in complete FAD medium, containing one part Ham's F12, three parts Dulbecco's modified eagle medium (DMEM), 10⁻⁴ M adenine, 10% (v/v) FBS, 0.5 μg/ml hydrocortisone, 5 μg/ml insulin, 10⁻¹⁰ M cholera toxin and 10 ng/ml EGF, as described previously (Gandarillas and Watt, 1995). Prior to suspension-induced differentiation in methylcellulose (Adams and Watt, 1989) we enriched for stem cells by filtering the disaggregated keratinocytes twice through a 40-µm sterile membrane. For knockdown or overexpression experiments, keratinocytes were grown in Keratinocyte-SFM medium (Gibco) supplemented with 0.15 ng/ml EGF and 30 μg/ml BPE. All isoforms of PKC were inhibited using 5 μM GF 109203X (Tocris; at lower concentrations GF190203X preferentially inhibits the -α and -β isoforms) and histone deacetylase was inhibited with TSA (Sigma Aldrich, UK). PDMS substrates that mimic the topography of the epidermal-dermal junction were generated as described previously (Viswanathan et al., 2016).

Genome-wide expression profiling was performed using the Illumina BeadArray platform and standard protocols. Data were processed using R (http://www.r-project.org/) or Genespring GX13.1 software. The gene expression data are deposited in the GEO databank (GSE73147).

Microarray initial processing and normalisation were performed with BeadStudio software. BeadChip internal p-values (technical bead replicates) were used to identify genes significantly expressed above the background noise. To filter out genes with signals that were not significant, a p-value of 0.05 was used as the cut-off value and only genes with a p-value<0.05 in at least one sample passed the filter. From the original set of 34,685 gene targets, 23,356 targets met this criterion. The data were imported into GeneSpring v13, normalised using quantile normalisation and the biological replicates averaged for subsequent analysis. We performed pairwise comparison between 0 hr and 4, 8 and 12 hr. Genes showing a fold change higher than two compared to control (and both p-values were significant) were subjected to GO analysis with GeneSpring v13 (152 genes between 4 hr and 0 hr, 553 between 8 and 0 hr, 1136 genes between 12 and 0 hr). Hierarchical clustering based on Pearson’s uncentered distance was performed on time course gene expression data and the results presented as a heatmap. Additionally, t-SNE plots were generated for the transcriptomics data using the t-SNE package for R.

A two-way ANOVA with multiple comparisons corrected using the Holm-Sidak test was used to identify the statistically significant effects of single phosphatase knock-down on the expression of the other phosphatases (p_ij, Supplementary file 9). The weight of each edge was calculated as the inverse p-value for the respective interaction (w_ij = 1/p_ij). For simplicity, we kept only the significant links (p-value<0.05). The networks depicted in Figure 4a were drawn with Cytoscape (cytoscape.org, V3.2.1).

The Boolean network analysis was performed using RE:IN and RE:SIN software, which is designed to encode and test whether a single or set of Boolean networks are consistent with observed changes in gene expression states. The methodology has been described previously (Dunn et al., 2014; Yordanov et al., 2016; Shavit et al., 2016). Briefly, a set of possible and definite interactions comprises an Abstract Boolean Network (ABN), which implicitly defines a set of concrete Boolean networks. Each concrete network has a unique topology defined by the different combinations of the possible interactions, which can be present or absent. This allows the user to explore alternative model topologies, and to test the requirement for different interactions that are suggested by experimental data. The software allows the user to seek the set of concrete Boolean network models that are consistent with a set of experimental constraints, and then to use this constrained set of models to make predictions of untested behaviour. Experimental constraints are defined as expected states along network trajectories (Figure 4b), which are constructed from discretised gene expression patterns. The software is freely available to use and can be accessed together with tutorials and FAQ at research.microsoft.com/rein.

We considered Boolean networks with synchronous updates, such that each gene in the network updates at each step in the network trajectory. Logical update functions for each gene, which define how the gene switches on or off in response to its regulators, were as defined by Yordanov et al., 2016. We first used RE:IN to test whether a single network (one that does not switch along the differentiation trajectory) is consistent with the observed changes in gene expression. We tested all possible single networks by allowing possible positive and negative interactions between each pair of genes. None of these networks was found to be consistent with the gene expression changes. We next constructed a switching problem in which we defined four ABNs, which individually correspond to the time points measured in our knockdown experiments. We defined the allowable network switching to be monotonic: networks can progress forward in time (i.e. 0 hr to 4 hr, 0 hr to 8 hr etc.) but not backwards (i.e. 8 hr cannot switch to 4 hr or 0 hr). We then used RE:SIN to test whether it was possible to find a sequence of networks that allowed the gene expression changes to be satisfied. We found that this was possible, and, furthermore, identified interactions that were required/disallowed to meet the constraints (Figure 4c).

The network governing differentiation must be consistent with observed changes in gene expression patterns over time. To investigate whether a given network could recapitulate the observed behaviour, we defined the expected states along a differentiation trajectory starting from the initial gene expression state at 0 hr. Three experimental constraints were defined, corresponding to measured gene expression changes at specific time points during suspension-induced differentiation, PKC inhibition and TSA treatment (Figure 4b). Individual gene expression measurements were discretised as ‘on’ or ‘off’ if their mean expression was respectively higher or similar/lower than the average of all genes. We set one step in the network trajectory to correspond to 1 hr in the time course, such that we defined expected gene expression states at 0, 4, 8 and 12 hr.

A pre-confluent keratinocyte culture was split into two separate cultures with FAD^{-lysine-arginine} medium (Sigma) differentially supplemented either with K⁰R⁰ or with K⁸R¹⁰ (stable isotopes of amino acids Lysine and Arginine; Cambridge Isotope Laboratories) (Ly et al., 2014). FCS used for SILAC medium was also depleted of Lysine and Arginine (Sigma). Cells were grown in SILAC medium for 5–6 days to reach 70–80% confluence and later harvested for downstream assays. Light labelled cells served as the 0 hr sample whereas heavy labelled cells were suspended in methylcellulose and harvested at 4, 8 and 12 hr. Cell extracts prepared from 0 h cells were mixed individually with 4, 8 and 12 hr samples at a 1:1 ratio of total protein. The mixed samples were then subjected to mass spectrometry (MS).

For measurement of protein abundance changes by SILAC, lysates containing ~80 µg protein were prepared in LDS sample buffer containing reducing agent (TCEP). Proteins were separated by electrophoresis on a 4–12% Bis-Tris NuPAGE gel, which was Coomassie-stained and cut into eight equally sized gel slices. Gel-embedded proteins were reduced with TCEP, alkylated with iodoacetamide, and trypsin-digested to release peptides. Peptides were extracted from gels using 50% acetonitrile (ACN) containing 5% formic acid (FA), dried and resuspended in 5% FA. Peptides were then analysed by LC-MS/MS on a Dionex RSLCnano coupled to a Q-Exactive Orbitrap classic instrument. Specifically, peptides were loaded onto a PepMap100 75 µm x 2 cm trap column, which was then brought in-line with a 75 µm x 50 cm PepMap-C18 column and eluted using a linear gradient over 220 min at a constant flow rate of 200 nl/min. The gradient composition was 5% to 35% B, where solvent A = 2% ACN+0.1% FA and solvent B = 80% ACN+0.1% FA. Peptides were eluted with a linear elution gradient (5% B to 35% B) over 220 min with a constant flow rate of 200 nl/min. An initial MS scan of 70,000 resolution was acquired, followed by data-dependent MS/MS by HCD on the top 10 most intense ions of ≥2+ charge state at 17,500 resolution.

For phosphoenrichment analysis, lysates containing ~2 mg of protein were chloroform:methanol-precipitated. Pellets were resuspended in 8 M urea in digest buffer (0.1 M Tris pH 8 + 1 mM CaCl₂), diluted to 4 M urea with additional digest buffer, and digested with Lys-C for 4 hr at 37°C. The digests were diluted to 1 M urea and trypsin-digested overnight at 37°C. Digests were then acidified and desalted using SepPak-C18 vacuum cartridges. Desalted peptides were resuspended in mobile phase A for HILIC (80% ACN + 0.1% TFA). HILIC chromatography was performed using a Dionex Ultimate 3000 with a TSK Biosciences Amide-80 column (250 × 4.6 mm) (Di Palma et al., 2013; McNulty and Annan, 2008; Navarro et al., 2011). Peptides were eluted using an exponential gradient (80% B to 60% B) composed of A (above) and B (0.1% TFA) at 0.4 ml/min over 60 min. 16 fractions were collected from 25 to 60 min. These fractions are enriched for hydrophilic peptides, including phosphopeptides. The fractions were dried before further phosphoenrichment by titanium dioxide (TiO₂).

For TiO₂ enrichment, HILIC fractions were resuspended in loading buffer (70% ACN + 0.3% lactic acid + 3% TFA). 1.25 mg of TiO₂ (GL Sciences) was added to each fraction and incubated for 10 min to bind phosphopeptides. Beads were washed with loading buffer, and two wash buffers, composed of (1) 70% ACN + 3% TFA and (2) 20% ACN + 0.5% TFA. Phosphopeptides were eluted in two steps, first with 4% of ammonium hydroxide solution (28% w/w NH₃) in water for 1 hr and again with 2.6% of ammonium hydroxide solution in 50% ACN overnight. Elutions were collected, dried, resuspended in 5% FA and analysed by LC-MS/MS. The LC-MS analysis was performed similarly to above with the following modifications: peptides were chromatographed on an EasySpray PepMap 75 µm x 50 cm column, and a ‘Top30’ method was used where the top 30 most intense ions were chosen for MS/MS fragmentation. Raw data were then processed using MaxQuant, which implements the Andromeda search engine, for peptide and protein identification and SILAC quantitation.

For the proteomics dataset, 1155 out of 2024 proteins had a significant p-value<0.05 (calculated in MaxQuant). To identify differentially expressed proteins between time 0 and 4, 8 and 12 hr, we collected per time point the proteins whose absolute ratios were >0.5. For each time point, the proteins were then split according to whether the ratio was positive or negative. Six lists were annotated for GO terminology using GeneSpring. The extensive list of GO-terms was submitted to Revigo (Supek et al., 2011) to reduce complexity and the resulting GO-categories depicted in Figure 1—figure supplement 1g Figure 1f. The mass spectrometry proteomics data have been deposited with the ProteomeXchange Consortium via the PRIDE partner repository with the dataset identifier PXD003281.

siRNA nucleofection was performed with the Amaxa 16-well shuttle system (Lonza). Pre-confluent keratinocytes were disaggregated and resuspended in cell line buffer SF. For each 20 μl transfection (program FF-113), 2 × 10⁵ cells were mixed with 1–2 μM siRNA duplexes as described previously (Mulder et al., 2012). Transfected cells were incubated at ambient temperature for 5–10 min and subsequently replated in pre-warmed Keratinocyte-SFM medium until required for the downstream assay. SMART pool ON-TARGET plus siRNAs (Ambion/GE Healthcare) were used for gene knockdowns. Each SMART pool was a mix of four sets of RNAi oligos. The sequences of the siRNA oligos are provided in the Supplementary file 11.

shRNA lentiviral vectors were purchased from Sigma-Aldrich (Mission shRNA). Two shRNAs per target were chosen. Their sequences are listed in Supplementary file 12. Primary keratinocytes were transduced with lentiviral particles for 48 hr, before being subjected to 2 μg/ml Puromycin selection for 72 hr. Cells were then expanded for experimental analysis.

For the Doxycycline inducible system we used the pCW57-GFP-2A-MCS lentiviral vector (gift from Adam Karpf; Addgene plasmid # 71783). Primary keratinocytes were transduced with lentiviral particles containing protein expression vector encoding genes for wild-type DUSP6, mutant DUSP6^C293S and DUSP10. 2 days post-transduction, cells were subjected to 2 μg/ml Puromycin selection for 3 days and then 1 µg/ml Doxycycline was added to the growth medium.

For Cumate induction, we used the lentiviral QM812B-1 expression vector (System Biosciences). Cells were transduced and selected as described for the Doxycycline system and protein expression was induced by addition of 30 µg/ml Cumate solution to the growth medium.

One-hundred, 500 or 1000 keratinocytes were plated on a 3T3 feeder layer per 10 cm diameter dish or per well of a six-well dish. After 12 days, feeders were removed and keratinocyte colonies were fixed in 10% formalin (Sigma) for 10 min then stained with 1% Rhodanile Blue (1:1 mixture of Rhodamine B and Nile Blue A [Acros Organics]). Colonies were scored using ImageJ and clonogenicity was calculated as the percentage of plated cells that formed colonies.

Pre-confluent keratinocyte cultures (km) were disaggregated and transfected with SMART pool siRNAs or non-targeting control siRNAs, or with lentiviral shRNAs. 24 hr post-siRNA transfection or >7 days post-shRNA infection, keratinocytes were collected and reseeded on irradiated de-epidermised human dermis (Sen et al., 2010) in six-well Transwell plates with feeders and cultured at the air–liquid interface for 3 weeks. The cultures were then fixed in 10% neutral buffered formalin (overnight), paraffin embedded and sectioned for histological staining. 6 μm thick sections were labelled with haematoxylin and eosin or appropriate antibodies. H and E stained images were acquired with a Hamamatsu slide scanner and analysed using NanoZoomer software (Hamamatsu).

H and E stained sections were scanned with a Hamamatsu slide scanner and analysed in Python. Briefly, tissue sections were automatically detected based on k-means segmentation, and colours were rebalanced. The epidermis was then isolated by colour separation (extraction of haematoxylin color channels by a deconvolution algorithm) (Ruifrok and Johnston, 2001). After performing an automatic thresholding, epidermal area was measured by pixel counting. To measure epidermal thickness, the area was divided by a 1-pixel wide reduced-topology (van der Walt et al., 2014). Full script is available as Supplementary Source Code 1.

For IVL staining of cells in suspension, cells were recovered from methylcellulose, spun down onto a slide using a Cytospin and then fixed with 4% paraformaldehyde and permeabilised. Sections of paraffin embedded reconstituted skin were dewaxed, rehydrated and subjected to heat-mediated antigen retrieval. Immunofluorescence labelling of cells and sections was performed with the appropriate antibodies. Wherever indicated, DAPI was used to label nuclei. Image acquisition was performed using a Nikon A1 confocal microscope (Nikon Instruments Inc.). To determine % antibody-labelled cells in sections, the total numbers of labelled cells and DAPI+ cells were assessed in Image J using the Cell Counter plugin. To determine the number of labelled cells per μm of basement membrane, the number of labelled cells was divided by the length of the basement membrane.

Cells were lysed on ice in 1x RIPA buffer (Bio-Rad) supplemented with protease and phosphatase inhibitor cocktails (Roche). Total protein was quantified using a BCA kit (Pierce). Soluble proteins were resolved by SDS-PAGE on 4–15% Mini-PROTEAN TGX gels (Bio-Rad) and transferred onto Trans-Blot 0.2 µm PVDF membranes (Bio-Rad) using the Trans-Blot Turbo transfer system (Bio-Rad). Primary antibody probed blots were visualised with appropriate horseradish peroxidase-coupled secondary antibodies using enhanced chemiluminescence (ECL; Amersham). The ChemiDoc Touch Imaging System (Bio-Rad) was used to image the blots. Quantification of detected bands was performed using Image Lab software (Bio-Rad). Uncropped versions of the western blot presented in that manuscript are available in Supplementary file 14.

The procedure was modified from previous reports (Jensen et al., 1999). Skin samples from either breast or abdomen were obtained as surgical waste with appropriate ethical approval and treated with Dispase (Corning) overnight on ice at 4°C. The epidermis was peeled off as an intact sheet and immediately fixed in 4% paraformaldehyde for 1 hr. Fixed epidermal sheets were washed and stored in PBS containing 0.2% sodium azide at 4°C. Sheets were permeabilised and stained with specific antibodies in a 24-well tissue culture plate. Image acquisition was performed using a Nikon A1 confocal microscope. 3D maximal projection (1024 × 1024 dpi), volume rendering and deconvolution on stacked images were generated using NIS Elements version 4.00.04 (Nikon Instruments Inc.).

Antibodies against the following proteins were used: P-ERK (Cell Signaling # 9101; western blot – 1:1000 dilution), ERK (Cell Signaling # 9102; western blot – 1:1000), P-p38 (Cell Signaling # 9211; western blot – 1:1000), p38 (Cell Signaling # 9212; western blot – 1:1000), Cyclophilin B (R and D # MAB5410; western blot – 1:4000), MKP-3/DUSP6 (Abcam # ab76310; western blot - 1:1000 and R and D Systems # MAB3576-SP; immunostaining – 1:200), PPTC7 (Abcam # ab122548; western-blot 1:250 and Sigma # HPA039335; immunostaining – 1:200), PTPN1/PTP1B (Sigma # HPA012542; western-blot - 1:500; R and D Systems # AF1366-SP; immunostaining – 1:200), PTPN13 (R and D Systems # AF3577; western-blot - 1:300 and immunostaining – 1:200), PPP3CA/Calcineurin A (Sigma # HPA012778; western-blot 1/1000 and R and D Systems # MAB2839-SP; immunostaining – 1:200), DUSP10 (Abcam # 140123; western-blot - 1:1000 and immunostaining – 1:200), TP63 (SCBT # sc-8431 or Abcam # ab735; immunostaining – 1:100), Involucrin (SY5, in-house; immunostaining – 1:500, or DH1, in-house; immunostaining – 1:200), Ki67 (Abcam # ab16667; immunostaining – 1:100) and β1-Integrin (clone P5D2, in-house; immunostaining – 1:200). Species-specific secondary antibodies conjugated to Alexa 488 or Alexa 594 were purchased from Molecular Probes, and HRP-conjugated second antibodies were purchased from Amersham and Jackson ImmunoResearch.

Total RNA was isolated using the RNeasy kit (Qiagen). Complementary DNA was generated using the QuantiTect Reverse Transcription kit (Qiagen). qPCR analysis of cDNA was performed using qPCR primers (published or designed with Primer3) and Fast SYBR green Master Mix (Life Technologies). RT-qPCR reactions were run on a CFX384 Real-Time System (Bio-Rad). Heatmaps of RT-qPCR data were generated by Multiple Expression Viewer (MeV_4_8) or GraphPad Prism 7. Sequences of qPCR primers are provided in Supplementary file 13.

Statistical analysis of the quantifications presented in the Figure legends was performed using GraphPad Prism 7.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

This research was funded by the Medical Research Council and the Wellcome Trust. The authors also gratefully acknowledge the use of Core Facilities provided by the financial support from the Department of Health via the National Institute for Health Research (NIHR) comprehensive Biomedical Research Centre award to Guy's and St Thomas' NHS Foundation Trust in partnership with King's College London and King’s College Hospital NHS Foundation Trust. Bénédicte Oulès is the recipient of fellowships from Societé Francaise de Dermatologie, Collège des Enseignants en Dermatologie de France and Fondation d’Entreprise Groupe Pasteur Mutualité. We thank Giacomo Donati, Arsham Ghahramani and Klaas Mulder for helpful advice and discussions and are grateful to Fredrik Pontén for providing anti-phosphatase antibodies.

© 2017, Mishra et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.