1. Computational and Systems Biology
  2. Stem Cells and Regenerative Medicine
Download icon

A protein phosphatase network controls the temporal and spatial dynamics of differentiation commitment in human epidermis

  1. Ajay Mishra
  2. Bénédicte Oulès
  3. Angela Oliveira Pisco
  4. Tony Ly
  5. Kifayathullah Liakath-Ali
  6. Gernot Walko
  7. Priyalakshmi Viswanathan
  8. Matthieu Tihy
  9. Jagdeesh Nijjher
  10. Sara-Jane Dunn
  11. Angus I Lamond
  12. Fiona M Watt  Is a corresponding author
  1. King’s College London, United Kingdom
  2. University of Cambridge, United Kingdom
  3. University of Dundee, United Kingdom
  4. University of Edinburgh, United Kingdom
  5. Université Paris Descartes, France
  6. Microsoft Research, United Kingdom
Research Article
Cite this article as: eLife 2017;6:e27356 doi: 10.7554/eLife.27356
4 figures, 2 data sets and 16 additional files

Figures

Figure 1 with 1 supplement
Genomic and proteomic analysis identifies protein dephosphorylation events that correlate with commitment.

(a) Schematic of experimental design. (b) Colony formation by cells harvested from suspension at different times. Representative dishes are shown together with % colony formation (n = 2 independent experiments, n = 3 dishes per condition per experiment; p-values calculated by Tukey’s multiple comparison test). (c) Cells isolated from suspension at different time points were labelled with anti-involucrin (IVL) antibody (green) and DAPI as nuclear counterstain (blue). IVL-positive cells were counted using ImageJ (n = 3 independent cultures; more than 300 cells counted per condition. p-value calculated by two-tailed t-test). Scale bars: 50 μm (d) RT-qPCR quantification of ITGα6, TP63, IVL and TGM1 mRNA levels (relative to 18 s expression) (n = 3 independent cultures). (e) t-SNE plot of genome-wide transcript expression by keratinocytes placed in suspension for different times. The t-SNE algorithm takes a set of points in a high-dimensional space and finds a faithful representation of those points in a lower-dimensional space, typically the 2D plane. (f) Heatmap representing hierarchical clustering of differentially expressed proteins (p<0.05). (g–i) Dot plots correlating expression of significantly differentially expressed peptides (p<0.05) that change twofold relative to 0 hr in at least one of the time points, with their corresponding differentially expressed transcripts. Pearson correlations (r) are indicated. (j) Histogram of normalised SILAC ratios corresponding to high confidence phosphorylation sites that differ between 0 and 4 hr. (k) Scatter plot correlating log2 normalised SILAC ratios for total protein changes (y-axis) with log2 phospho-peptide ratios (x-axis) between 0 and 4 hr. (b, c) *p<0.05; **p<0.01; ns = non-significant).

https://doi.org/10.7554/eLife.27356.002
Figure 1—figure supplement 1
Clonal growth, genomic and proteomic analysis of suspension-induced terminal differentiation.

(a) Keratinocytes harvested after 0, 1, 2, 4, 6, 8, 10, 12 or 24 hr in suspension in methylcellulose were seeded at 100, 500 or 1000 cells in 10 cm2 culture dish on mitotically inactivated J2 3T3 feeders and cultured for 12 days. Following staining, colonies were counted using ImageJ (n = 2 independent experiments with three replicate dishes per experiment; p-value calculated by two tailed t-test). (b) RT-qPCR quantification of IVL and TGM1 mRNA levels at different times in suspension (n = 3 independent cultures; p-value calculated by two tailedt-test). (c) RT-qPCR quantification of DLL1, Ki67 and Lrig1 mRNA levels at different times in suspension (n = 3 independent cultures; p-value calculated by one-way ANOVA). (d) Hierarchical clustering of significantly expressed transcripts at 0, 4, 8 and 12 hr in suspension (each time point represents the mean value of n = 3 independent experiments). (e) GO analysis of differentially expressed genes upregulated at 4, 8 and 12 hr relative to 0 hr. The bar plots represent –log10 of p-values of the identified GO terms. (f) Schematic of SILAC-Mass Spectrometry labelling protocol. (g) GO analysis of proteins ranked in the order of their expression level (fold increase relative to 0 hr) at 4, 8 and 12 hr. GO terms were fetched for individual proteins rather than for clusters of proteins. Bar plots represent –log10 of the p-values of the identified GO terms. (a–c) Error bars represent s.d. *p<0.05; **p<0.01; ***p<0.005; ns = non significant.

https://doi.org/10.7554/eLife.27356.003
Figure 2 with 4 supplements
Functional screen identifies candidate phosphatases that regulate commitment.

(a) Heatmap showing differential expression at 4, 8 and 12 hr relative to 0 hr of phosphatases that are upregulated at 4 hr in the microarray dataset. (b) Effect of knocking down the 22 phosphatases identified in (a) on clonal growth of keratinocytes. Values plotted are average % clonal growth in n = 3 independent screens with n = 3 independent cultures per screen. Green: siSCR control. Red, blue: phosphatases with statistically significant effects on colony formation are shown (red: increase; blue: decrease). Grey: no statistically significant effect. (c, d) Effect of knockdowns on clonal growth after 0, 4, 8 or 12 hr in suspension. (c) Representative dishes. (d) Quantitation of mean % clonal growth ± SD (n = 3 independent samples). p-values generated by unpaired T-test. (e, f) RT-qPCR quantification of TP63 (e) and TGM1 (f) mRNA levels (relative to 18 s expression) in the same conditions as in (c). n = 3 independent transfections. (g) Epidermal reconstitution assay following knockdown of DUSP6, PTPN1, PPP3CA or PTPN13. n = 2 independent transfections. Top row shows representative H and E images. Epidermal thickness was quantified in multiple fields from eight sections per replicate ± SD relative to scrambled control (siSCR). Middle row shows staining for TP63 (pink) with DAPI nuclear counterstain (blue). % DAPI-labelled nuclei that were TP63+ was quantified in n = 2–3 fields per replicate. Bottom row shows staining for Ki67 (brown) with haematoxylin counterstain (blue). % Ki67+ nuclei was quantified in n = 3–6 fields per replicate. Error bars represent mean ± s.d. p-values were calculated using one-way ANOVA with Dunnett's multiple comparisons test (*p<0.05; **p<0.01; ****p<0.0005; ****p<0.0001; ns = non significant).

https://doi.org/10.7554/eLife.27356.004
Figure 2—figure supplement 1
Expression in published datasets of the protein phosphatases identified by suspension-induced differentiation.

(a) Heatmap showing mRNA expression of candidate phosphatases during calcium-induced stratification of primary human keratinocytes (GSE38628). (b) Heatmap showing mRNA expression of candidate phosphatases during ex vivo human epidermal reconstitution (GSE52651).

https://doi.org/10.7554/eLife.27356.005
Figure 2—figure supplement 2
Effects of phosphatase knockdown on keratinocyte growth and differentiation.

(a) Effect of phosphatase knockdown on colony formation. Values plotted are average % colony formation in three independent screens with three replicate dishes per screen. p-values were generated using a two tailedt-test. (b) Following siRNA transfection keratinocytes were seeded in collagen-coated 96-well format dishes and cultured in KSFM medium for three days. 24 hr post-transfection the culture dishes were moved into an Incucyte microscope and scanned every hour for the next 48 hr to evaluate cell growth and proliferation. Each data point is the average of three replicate screens with eight independent well scans per condition per screen. Values at each time point were normalised to the respective 0 hr time point or the first scan point. A cumulative plot is shown with the relative frequency (%) of cell doublings. (c) RT-qPCR quantification of phosphatases mRNA levels relative to 18 s RNA as a function of time in suspension, showing efficiency of the different knockdowns in the starting populations. p-values were generated by the unpaired T-test. (d) Kinetics of change in the total area occupied by colonies formed by cells recovered from suspension. (n = 3 independent samples). p-values generated by unpaired T-test. (e) Growth rate of human keratinocytes after phosphatase knockdown, calculated by fitting the data in (b) to an exponential growth curve and averaging the rates. p-values generated by one-way ANOVA (f) Sections of reconstituted epidermis from the experiments shown in Figure 2g labelled with antibodies to ITGβ1 and IVL, with DAPI as nuclear counterstain. Merge and individual channel images are depicted. Scale bar: 50 μm. (g) Quantification of the experiments shown in Figure 2g, showing the number of TP63+ (top panel) or Ki67+ cells (bottom panel) per μm of basement membrane. DAPI-labelled nuclei that were TP63+ were quantified in n = 2–3 fields per replicate. % Ki67+ nuclei were quantified in n = 3–6 fields per replicate. Error bars represent mean ± s.d. p-values were calculated using one-way ANOVA. a,c,d,g. Error bars represent s.d. *p<0.05; **p<0.01; ***p<0.0005; ****p<0.0001; ns = non significant.

https://doi.org/10.7554/eLife.27356.006
Figure 2—figure supplement 3
Effect of shRNA-mediated knockdown of the phosphatases on epidermal reconstitution.

(a) RT-qPCR quantification of phosphatase mRNA levels relative to 18 s RNA upon knockdown by lentiviral-mediated shRNA one week after infection and selection. p-values were generated by one-way ANOVA. (b) Epidermal reconstitution following knockdown of DUSP6, PPTC7, PTPN1, PPP3CA or PTPN13. n = 2 independent infections. Representative H and E images. (c) Epidermal thickness was quantified in eight non-consecutive sections per biological replicate as explained in Figure 2—figure supplement 4 and normalised to shSCR values. p-values were generated by one-way ANOVA. (d) Representative images of epidermis reconstituted by keratinocytes expressing shSCR or sh1DUSP6 co-stained for TP63 (yellow) and Ki67 (pink) with DAPI nuclear counterstain (blue). White asterisk shows a TP63+ Ki67 cell. (e) % Ki67-labelled nuclei that were also TP63+ was quantified in n = 5 fields per replicate. p-values generated by one-way ANOVA were >0.05. (f) Representative images of epidermis reconstituted by keratinocytes expressing shSCR or sh1DUSP6 co-stained for TP63 (yellow) and IVL (pink) with DAPI nuclear counterstain (blue). White asterisk shows a suprabasal cell co-stained for IVL and TP63. a, c, e. Error bars represent s.d. *p<0.05; **p<0.01; ***p<0.0005; ****p<0.0001; ns = non-significant.

https://doi.org/10.7554/eLife.27356.007
Figure 2—figure supplement 4
Workflow for automated quantification of the epidermal thickness.

Sketch representing the automated pathway to extract epidermis thickness from tissue sections. The python script was divided into three steps. (1) Detection of tissue sections by processing from ‘.ndpi’ files. With intensity clustering, the background was suppressed and tissue sections were localised on the slide. Each section was extracted with the highest resolution and split into four vignettes. (2) The epidermis was isolated using information from H and E staining in each vignette. An automatic thresholding was applied to select only the epidermis. Potential artefacts that were small in size were suppressed by thresholding on each object’s size. (3) Epidermal area was measured by pixel counting. To measure epidermal thickness, area was divided by a 1-pixel wide reduced-topology.

https://doi.org/10.7554/eLife.27356.008
Figure 3 with 1 supplement
Pro-commitment phosphatases regulate MAPK signalling and AP1 transcription factors.

(a) Gene Ontology (GO) term enrichment analysis of ranked peptides dephosphorylated at 4 hr. (b) Venn diagram showing intersection of signalling pathways regulated at 4 hr. (c) Top 15 dephosphorylated peptide sites at 4 hr, showing ratio between change in phospho-peptides and change in total protein. Highlighted in orange are phosphorylations on MAPKs. (d, e) Representative western blot (d) and quantification (e) showing phospho-ERK1/2 and total ERK1/2 in cells harvested after 0, 4, 8 and 12 hr in suspension. Cyclophilin B: loading control. (f) Heatmap represents mRNA expression (relative to 18 s RNA) of AP1 transcription factor superfamily members during suspension-induced differentiation post-phosphatase knockdown (n = 3; values plotted are means of three independent transfections). See Supplementary file 6 for p-values generated by two-way ANOVA. (g–h) Representative western blot, with quantitation, of phospho-ERK1/2 and total ERK1/2 in suspended cells following DUSP10 knockdown. siSCR and loading controls are shown. (i) Heatmap showing mRNA expression (relative to 18 s mRNA) of JUN, FOS and MAF family members after DUSP6 and DUSP10 knockdown (values plotted are means of three independent transfections normalised against scrambled control; see Supplementary file 7 for p-values generated by two-way ANOVA). (j) Clonal growth (representative dishes and quantification) following doxycycline-induced over-expression of DUSP10, DUSP6 and mutant DUSP6C293S in primary keratinocytes (n = 3 independent cultures). p-values were calculated using one-way ANOVA with Dunn's multiple comparisons test (*p<0.05; ns = non significant).

https://doi.org/10.7554/eLife.27356.009
Figure 3—figure supplement 1
Effects of DUSP10 knockdown and DUSP6 and DUSP10 over-expression.

(a, b) Western blot (a) and RT-qPCR measurements (b) of DUSP10 knockdown by SMART pool siRNA. p-values were generated by one-tail t-test. (c, d) Western blot of phospho-p38 and total p38 in suspended cells following DUSP10 knockdown. siSCR, cyclophilin B loading control and quantitation are shown. (e) RT-qPCR quantification of doxycycline-induced over-expression of DUSP6, mutant DUSP6C293S and DUSP10 relative to 18 s mRNA. Cells were treated with 1 µg/ml doxycycline for 24 hr. p-values were generated by a one-tail t-test. (f) RT-qPCR quantification of Cumate-induced over-expression of DUSP6 and mutant DUSP6C293S relative to 18 s mRNA. p-values were generated by one-way ANOVA. (g) Representative dishes showing effects of overexpressing wild type and mutant DUSP6 on clonal growth (representative of n = 3 dishes). (b, e, f) Error bars represent s.d. *p<0.05; ***p<0.0005; ****p<0.0001; ns = non significant.

https://doi.org/10.7554/eLife.27356.010
Figure 4 with 3 supplements
An autoregulatory network of phosphatases controls commitment.

(a) Colours represent log2 fold-change in phosphatase mRNA expression relative to 0 hr (values plotted are means of three independent experiments normalised against siSCR control). (b) Experimental observations that are encoded as constraints on the Boolean network trajectories. We defined three experimental constraints (gene expression changes during suspension-induced differentiation in the absence of pharmacological inhibitors, as well as under TSA and PKCi treatments). Each constraint encodes discrete gene expression states at the indicated time steps. For cells in the absence of drugs we imposed a switching scheme, whereby the system must change the representative network in order to achieve the expression constraints. Yellow boxes indicate the network that represents the system at that step; blank boxes for a given step (column) indicate that we did not impose a specific network, so the system can remain in the current network or switch forwards. Active means the gene is considered to be ‘on’ at that step, inactive means to be ‘off’ at that step. The discretisation is available in Supplementary file 9. (c) Networks able to satisfy the model constraints of the Boolean formalism in (b) are depicted. Solid lines show interactions already calculated in (a), while dashed lines are possible interactions inferred from Figure 4—figure supplement 1. See also Supplementary file 8 and 9. (d) Heatmap represents mRNA expression (relative to 18 s mRNA) of individual phosphatases in cells treated in suspension for 12 hr with PKCi or TSA (values plotted are the means of three independent experiments normalised against vehicle-treated control). (e, f) Representation of commitment as two saddle-node bifurcations in a direction xi of the state space for control cells or cells treated with PKCi or TSA. In the control both stem and differentiated cell states are stable (attractors), while commitment is an unstable state. Since the 0 hr network is able to reach the expression constraint for PKCi at 12 hr, we hypothesise that on PKCi treatment, the only stable state is the stem state. On TSA treatment there is a mandatory switch from the 0 hr network but the 12 hr network cannot be reached at any time point; we therefore hypothesise that commitment becomes a stable state while the stem and differentiated cell states are unstable. (g) 3D-volume rendered confocal images of wholemounts of human epidermis labelled with antibodies against commitment phosphatases or ITGβ1 (green) and counterstained with DAPI (blue). The distribution of each phosphatase relative to ITGβ1 is also shown graphically. (h) 3D-volume rendered confocal images of primary keratinocytes cultured on PDMS topographic substrates and labelled with antibodies to DUSP10 and DUSP6. Scale bars: 100 μm.

https://doi.org/10.7554/eLife.27356.011
Figure 4—figure supplement 1
Simulated Boolean networks.

(a) Single ABN with all possible interactions among the six commitment phosphatases. (b) RE:SIN representation of the same networks as in Figure 4a. (c) RE:SIN representation of the networks as in b) (solid lines) and additional possible interactions (dashed lines) inferred from Figure 4—figure supplement 2d,e.

https://doi.org/10.7554/eLife.27356.012
Figure 4—figure supplement 2
Effects on phosphatase expression of knockdowns and treatment with TSA or PKCi.

(a) Heatmaps showing the effect of knocking down individual phosphatases on mRNA levels of other phosphatases with time in suspension (0, 4, 8, 12 hr). RT-qPCR is relative to 18 s mRNA (n = 3 independent transfections; see Supplementary file 9 for p-values generated by two-way ANOVA with Dunnett's multiple comparisons test). (b) Effects of TSA and PKCi on clonogenicity of keratinocytes recovered following suspension in methylcellulose for 12 hr. Representative dishes of n = 3. (c) mRNA levels of IVL, TGM, ITGα6, and TP63 in cells held in suspension for 12 hr. Cells were treated with TSA, PKCi or DMSO (vehicle control) (n = 3 independent treated cultures with two technical replicates each). p-values for the comparisons were generated by one-way ANOVA. (d) Western blots showing phosphatase levels in primary keratinocytes upon knockdown of scrambled control (siSCR), DUSP6, PPTC7, PTPN1, PTPN13 or PPP3CA and suspension for 0, 4, 8 or 12 hr. Cyclophilin B: loading control. (e) RT qPCR quantification of phosphatase mRNA levels (relative to 18 s mRNA) following doxycycline-induced over-expression of DUSP6, mutant DUSP6C293S and DUSP10. Cells were treated with 1 µg/ml doxycycline for 8, 24 or 48 hr (n = 3 independent cultures; see Supplementary file 10 for p-values generated by two-way ANOVA with Dunnett's multiple comparisons test).

https://doi.org/10.7554/eLife.27356.013
Figure 4—figure supplement 3
DUSP6 and DUSP10 distribution in human epidermal wholemounts.

(a, b) 3D-volume rendered confocal images of human epidermal wholemountslabelled with antibodies against DUSP6 (red, a; green, b) and ITGα6 (green, a) or DUSP10 (red, b). Wholemounts were counterstained with DAPI (blue).

https://doi.org/10.7554/eLife.27356.014

Data availability

The following data sets were generated
  1. 1
  2. 2

Additional files

Supplementary file 1

Log2 fold change of normalised gene expression for all pairwise comparisons of mRNA levels during suspension-induced terminal differentiation.

For each condition the mean of n = 3 independent replicates was used and the pairwise fold change comparison is between the means of both samples.

https://doi.org/10.7554/eLife.27356.015
Supplementary file 2

Proteomics data for all pairwise comparisons of protein levels at 4, 8 and 12 hr in suspension relative to the 0 hr control.

For each condition the mean of n = 3 independent replicates was used and the pairwise fold change comparison is between the means of both samples.

https://doi.org/10.7554/eLife.27356.016
Supplementary file 3

Phosphoproteomics data for pairwise comparisons at 4 and 8 hr in suspension relative to the 0 hr control.

https://doi.org/10.7554/eLife.27356.017
Supplementary file 4

p-values generated for RT qPCR of TP63 and TGM1 for each conditional time course relative to control time course (siSCR) by two-way ANOVA with Dunnett's multiple comparisons test (related to Figure 2e,f).

https://doi.org/10.7554/eLife.27356.018
Supplementary file 5

Log2ratio of phosphopeptides over total proteins at 4 hr.

https://doi.org/10.7554/eLife.27356.019
Supplementary file 6

Effect of phosphatase knockdown on AP1 transcription factor expression.

p-values generated for each conditional time course relative to control time course (SCR) by two-way ANOVA multiple comparisons (for AP1 superfamily factors). p-values generated for RT qPCR of AP1 factors for each conditional time course relative to control time course (siSCR) by two-way ANOVA with Dunnett's multiple comparisons test.

https://doi.org/10.7554/eLife.27356.020
Supplementary file 7

Effect of DUSP6 and DUSP10 knockdown on AP1 transcription factor expression.

p-values generated for RT qPCR of AP1 factors relative to control cells (siSCR) by two-way ANOVA.

https://doi.org/10.7554/eLife.27356.021
Supplementary file 8

Boolean expression patterns and phosphatases interactions used to generate Figure 4c,d.

https://doi.org/10.7554/eLife.27356.022
Supplementary file 9

p-values generated for RT-qPCR of phosphatases for each conditional time course relative to control time course (siSCR) by two-way ANOVA with Dunnett's multiple comparisons test.

https://doi.org/10.7554/eLife.27356.023
Supplementary file 10

One-way non-parametric ANOVA (Friedman test) with Dunn's multiple comparisons test for the effect of overexpressing DUSP6 and DUSP10 on mRNA levels of the pro-commitment phosphatases, determined by RT-qPCR.

https://doi.org/10.7554/eLife.27356.024
Supplementary file 11

siRNA library for phosphatase knockdown.

https://doi.org/10.7554/eLife.27356.025
Supplementary file 12

shRNA library for phosphatase knockdown.

https://doi.org/10.7554/eLife.27356.026
Supplementary file 13

List of qPCR primers.

https://doi.org/10.7554/eLife.27356.027
Supplementary file 14

Uncropped versions of the western blots presented in Figure 3d,g and Figure 3Figure 4—figure supplement 2c.

https://doi.org/10.7554/eLife.27356.028
Source code 1

Automated measurement of epidermal thickness.

https://doi.org/10.7554/eLife.27356.029
Transparent reporting form
https://doi.org/10.7554/eLife.27356.030

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)