Single-cell messenger RNA-sequencing analysis of basal epidermal stem cells

To explore the effects of non-cell autonomous Wnt signalling on epidermal cell state we sequenced the transcriptomes of single WT murine keratinocytes co-cultured with cells expressing an inducible form of stabilised beta-catenin (K14ΔNβ-cateninER) in a ratio of 9:1. We compared cells cultured in the absence of 4-hydroxy-tamoxifen (4OHT) with cells treated for 24 h with tamoxifen to induce beta-catenin. Cells were then disaggregated, loaded onto the C1 96-well microfluidic device (Fluidigm) and captured for sequencing. Owing to the single-cell capture method used, highly keratinized and terminally differentiated cells > 20 μm in diameter were excluded. We identified K14ΔNβ-cateninER cells by aligning reads to the transgene sequence and subsequently removed these cells from analysis (ten untreated cells and 14 activated cells). After quality control, we retained 125 WT control cells and 129 WT cells exposed to Wnt signalling neighbours. We recorded a median of 641,000 reads per cell equating to 4000–8000 genes expressed per cell (transcripts per million [TPM] > 1). Read alignment distribution was in line with other single-cell RNA-sequencing (RNA-seq) datasets with minimal ribosomal and intergenic reads (Additional file 1: Figure S1).

To explore cell-state heterogeneity in WT keratinocytes that had not been exposed to a neighbour in which beta-catenin was activated (untreated samples), we used reverse graph-embedding, a machine-learning technique. This enabled us to reconstruct a parsimonious tree connecting all observed epidermal cell states (DDRTree, Monocle 2) [12]. We applied the DDRTree algorithm to WT cells using expressed genes (TPM > 1) after removing cell-cycle associated genes. We identified five distinct in vitro cell states (Fig. 1a) forming three major branches that represent varying states of proliferation and differentiation.

Fig. 1 Molecular heterogeneity of epidermal cells in culture. a Epidermal cell transcriptomes and cell state relation visualised using DDRTree and coloured according to unsupervised clustering. Each data point is one cell and axes are dimensionally reduced components of the cell transcriptomes. Colours and numbers represent the five identified cell states. b Gene expression for four marker genes shown for each cell on the state map. Top left: Integrin beta-1 (Itgb1), a basal IFE marker. Top right: S100 differentiation associated genes. Bottom left: Mt2, a basal IFE marker. Bottom right: Keratin 10, a suprabasal IFE marker of commitment to differentiation. c Barplot showing mean expression per cell averaged over all keratin associated messenger RNAs (mRNA) (top) and all collagen mRNAs (bottom). d Distribution of the expression of cluster-specific markers coloured according to cell state. Two markers are shown per cell state. Black line represents the median expression of the cluster. e Median expression of the top three markers for each state per cell on the state map. f Heatmap showing the similarity between epidermal cells from this study and IFE cells from Joost et al. Similarity is measured by Pearson’s correlation coefficient. Cells from this study are coloured by cluster along the vertical axis. Cells from Joost et al. are coloured by differentiation status along the horizontal axis. Joost et al. IFE cell legend is shown in order of differentiation status. B basal IFE cells, DI/DII differentiated suprabasal IFE cells, KI/KII keratinized IFE cells. g Heatmap showing similarity between cluster average transcriptomes from this study and proliferative IFE cells from Zhang et al. [14] and activated vs quiescent IFE cells from Lien et al. [15] Full size image

States A and E showed highest expression of Mt2, a basal IFE marker, alongside markedly low expression of S100 epidermal differentiation complex genes in comparison with the remaining subpopulations. These expression patterns indicate states A and E represent transcriptomic signatures before commitment to differentiation (Fig. 1b; top left) [13]. Itgb1, another marker of basal IFE cells, showed variable expression in vitro in comparison to Mt2, but was expressed in all cells (Fig. 1b; top right). The differentiation marker Krt10 was variably expressed across all subgroups (Fig. 1b; bottom right) [4]. Separation between pre- and post-commitment cell states is further apparent when looking at pan-keratin and pan-collagen gene expression. In vivo, keratinocytes commit to differentiation upon detaching from the basement membrane, resulting in reduced collagen expression and increase in overall keratin content. States B, C and D express significantly more keratin messenger RNAs (mRNA) and conversely states A and E are characterised by higher collagen mRNA levels (Fig. 1c; Kolmogorov–Smirnov test, p < 0.05). For each cell state, we determined genes differentially expressed vs the remainder of the population and identified between six (state B) and 101 (state C) markers (Fig. 1d, Additional file 2: Table S1). We observed that the median expression of the top three markers for each state was sufficient to distinguish each state (Fig. 1e).

To determine the putative biological function of each cell state, we correlated the single-cell signatures with a comprehensive set of IFE subpopulation expression profiles identified by Joost et al. We further integrated two other bulk gene expression studies which identified signatures for proliferative and non-proliferative epidermal stem cells [14, 15]. From comparison with the Joost IFE subpopulations, all of our single cells correlate strongly with basal IFE stem cells, as expected since large (>20 μm) terminally differentiated cells were excluded from the analysis (Fig. 1f). States A, and to a lesser extent state E, showed high correlation with an isolated subpopulation of ‘non-proliferative self-renewing’ epidermal cells characterised by Lien et al., which we would describe, for simplicity, as non-proliferative stem cells. (Fig. 1g). We concluded that states A and E both exhibit a self-renewing cell gene expression signature but differ in proliferative state. States B, C and D formed a branch of the cell trajectory representing early commitment to differentiation, characterized by expression of proliferation associated genes from Roshan (e.g. MRPL33, YY1) and correlated strongly with the expression profile of proliferative keratinocytes from Zhang [3, 14]. This branch of the state tree shows expression of early differentiation markers such as MXD1, Dsc2, Dsg3 [16] and highest expression of S100 early differentiation-associated genes. Figure 1a summarises the state classification of cells as determined by our cluster and DDRTree analysis, highlighting the relationship between our three branches composed of five identified keratinocyte states.

Inducible Wnt signalling

In a previous study, we generated gene expression profiles from WT and beta-catenin activated adult mouse epidermis [17]. We reanalysed these data to estimate the relative proportion of cells in each of the cell states identified in vitro (Fig. 1a). We utilised CIBERSORT, a method for characterising the composition of tissue expression profiles resulting from mixtures of cells [18]. Our reanalysis indicated that epidermal beta-catenin signalling results in a depletion of cells in state A and increases the abundance of cells in state D (Fig. 2a). This is consistent with the in vivo observation that intrinsic activation of epidermal beta-catenin results in proliferation and expansion of the stem cell compartment [10].

Fig. 2 Induction of canonical Wnt signalling in a subpopulation of cells. a Estimated proportions of state A and state D cells in WT and beta-catenin activated epidermis samples from Collins et al. (GSE32966) [17]. b Immunofluorescence showing cytoplasmic and nuclear beta-catenin in WT and K14ΔNβ-cateninER keratinocytes after induction of canonical Wnt signalling using 4OHT. Scale bar, 20 μm. c Quantification of mean nuclear beta-catenin fluorescence intensity (n = 3 independent cultures). d Activation of canonical Wnt target genes in K14ΔNβ-cateninER cells upon induction with 4OHT quantified by quantitative reverse transcription polymerase chain reaction (n = 4). *P < 0.05, n.s. not significant. All data shown as mean ± SD Full size image

Next, to investigate whether epidermal cell states were altered by NCA Wnt signalling, we examined the treated sample, WT keratinocytes co-cultured with K14ΔNβ-cateninER keratinocytes in the presence 4OHT. We confirmed that K14ΔNβ-cateninER cells intrinsically activated canonical Wnt signalling in response to 4OHT by detecting beta-cateninER translocation into the nucleus (Fig. 2b, c; Lo Celso et al. [10]). We also validated upregulation of canonical downstream target genes such as Bmp4, Cyclin-D1 and Lef1 in K14ΔNβ-cateninER keratinocytes using quantitative reverse transcription polymerase chain reaction (qRT-PCR) (Fig. 2d).

Reconstruction of NCA Wnt induced state transition

Having identified several different states of WT keratinocytes and validated the intrinsic effects of beta-catenin activation, we used single-cell RNA-seq to deconvolve the effects of NCA Wnt signalling. Single-cell transcriptomes from WT keratinocytes co-cultured with 4OHT-activated K14ΔNβ-cateninER cells were compared with those of WT cells co-cultured with uninduced K14ΔNβ-cateninER cells and mapped onto the same dimensionally reduced space (Fig. 3a). To exclude the possibility of transcriptional changes resulting from 4OHT treatment alone, we screened bulk differential gene expression between the two cohorts of WT cells. We found no evidence of oestrogen receptor target genes among differentially expressed genes [19].

Fig. 3 Neighbouring Wnt activation reduces gene expression heterogeneity of protein-synthesis associated genes. a Keratinocytes exposed to Wnt signalling neighbours projected onto the WT epidermal cell state map. WT cells shown with grey circles, exposed cells shown coloured by cell cluster as triangles. Pie charts show relative proportions of WT vs signalling-exposed cells in each cell state cluster. States A and D are labelled as these states contain significantly different proportions of exposed vs unexposed cells. b Density of control vs NCA Wnt exposed cells along the reconstructed pseudotransition from states A to D. c Transcriptome coefficient of variation per-cell (TCOV), coloured by cell state and point size represents number of genes detected as expressed (TPM > 1) for that cell. Line depicts a loess-curve fit for the pseudotransition–TCOV relationship. d Scatterplot showing the log-ratio of coefficient of variation vs the log-ratio of gene expression between states A and D. Genes labelled depict the ten highest and lowest coefficient of variation ratios. Also see Additional file 1: Figure S2 for comparisons between states A and B and states B and D. e Violin plots depicting expression of cluster-specific genes showing differential dispersion between states. f Heatmap with symmetric gene expression scale (z-score normalised log2(TPM + 1)) showing reduction in gene expression heterogeneity of 281 genes between state A and the remaining transition states. g Gene rank enrichment analysis of log-fold change in gene coefficient of variation (CV) between cells in state D and the remainder of the cell population. Genes with lower variation in state 5 are enriched for translation and ribosome-associated gene ontology (GO) annotations (p value < 1e-6) Full size image

We observed the same five distinct transcriptional states in WT cells with or without NCA Wnt signalling (Fig. 3a), confirmed by independently analysing the treated cell population to reveal five equivalent subpopulations. However, exposure to non-cell autonomous Wnt signalling markedly changed the state distribution of keratinocytes (Fisher’s exact test, p < 0.05). Pie charts in Fig. 3a show the observed ratio of control and signalling-exposed cells. States A and D significantly deviated from the expected ratio (binomial test, p < 0.05). After exposure to NCA signalling, there was a depletion of cells in the non-proliferative stem cell state A and a higher than expected proportion of cells in state D, representing a transition towards a proliferative and more differentiated transcriptional state (Fig. 3a).

Taking the states with altered cell proportions and the transition states in between (states A, B, C and D), we reconstructed the state transition induced by neighbouring Wnt + keratinocytes using the Monocle pseudotime method [12]. WT and exposed cells were ordered from state A to state D to reconstruct the temporal order of gene expression changes for cells undergoing this transition, referred to as the pseudotransition. Figure 3b shows the proportion of control and Wnt signalling exposed cells along the reconstructed temporal transition from state A to state D; from this distribution it is clear that NCA Wnt exposed cells bias towards State D.

NCA Wnt signalling reduces heterogeneity in protein synthesis-associated transcripts

We next sought to understand why state A cells were uniquely depleted after neighbour Wnt signalling. Previous studies have shown that cell responses to extrinsic signalling are affected by intracellular and intercellular transcriptional noise [20,21,22]. We thus hypothesised that the response to NCA Wnt signalling involves changes in both the dynamic range of transcriptional variation (intracellular variation) and state-specific gene expression (intercellular variation).

We first examined whether there was a difference in intracellular transcriptomic heterogeneity between the three altered states and whether changes occurred along the pseudotransition. The resulting ordering of cells from state A to state D was used to examine the transcriptome coefficient of variation (TCOV) per cell (Fig. 3c). Here, TCOV is an intracellular measure of the spread of transcript abundance accounting for mean abundance. Notably, TCOV decreased over the state transition and was significantly higher in state A than states B, C and D (Kolmogorov–Smirnov test; p value < 0.05). This reduction in dynamic range of gene expression is consistent with previous studies that have shown that progenitor cells have a higher rate of stochastic multilineage gene expression that reduces upon cell-fate commitment [23, 24].

Next, we contrasted the heterogeneity of genes that do not change in expression level between the transition states (Additional file 1: Figure S2A–C). Figure 3d displays the relationship for the log-ratio of intercellular gene expression variation and expression level between the two extremes of the pseudotransition, states A and D, with the top ten differentially dispersed genes labelled. Of interest are genes that change in expression heterogeneity from state A to state D while remaining at constant expression levels. Notably, Baz2a, Sox2, Col7a1 and Calcrl were among the genes with reduced COV in state A without significant differential expression (Additional file 2: Table S2). Baz2a has been previously established as part of the nucleolar remodelling complex that is important for establishing epigenetic silencing and transcriptional repression of ribosomal RNA genes [25, 26]. Sox2 is an adult stem cell factor shown to be expressed in multiple epithelia [27]. Sox2 has been previously reported to be expressed in hair follicles but absent from the interfollicular epidermis [28]. Similarly, Col7a1 and Calcrl are significantly upregulated in hair follicle bulge stem cells [29]. Additional file 2: Table S2 lists the statistically significant and differentially heterogeneous genes characteristic of each state.

Figure 3e shows the expression variability for a selection of statistically significant (q < 0.05) differentially heterogeneous genes. We found multiple patterns of intercellular heterogeneity including genes that differed in variability across branches and genes that were selectively more variable in a single state. Peroxiredoxin-1 (Prdx1) is one such gene, which showed strongly variable expression in state A in comparison to the remaining states despite no statistically significant difference in median expression level.

When the four transition states were compared, we observed many more heterogeneously expressed genes in state A than states B, C or D. The contrast in number of differentially dispersed genes is demonstrated using the symmetric expression scale in Fig. 3f. We observed 281 significant and differentially heterogeneously expressed genes in state A, compared with only 19 genes in state D. A striking number of the heterogeneously expressed genes in state D are known regulators of stem cell identity such as Cdh11, Cav2 and Apc [30,31,32]. They are typically upregulated in basal stem cells; however, little is known about how their heterogeneity affects cell fate. Our analysis of intercellular variation suggests that in stem cells (state A) with low RNA and protein metabolism [33], transcriptional heterogeneity is lowest for stem-cell marker genes, emphasising the importance of transcriptional noise in addition to transcriptional amplitude. These observations are consistent with our hypothesis that state A cells are responsive to NCA Wnt signalling due to greater transcriptional variability. Exposure to this coordinated extracellular stimulus reduces transcriptional heterogeneity for these cells and biases their fate towards state D.

To determine genes essential for a cell to be receptive to neighbour Wnt activation, we analysed the fold-change in heterogeneity between states A and D, comprising the majority of genes with differential heterogeneity. We found strong enrichment for translation and ribosome-related genes, indicating a role for protein synthesis (p < 1e-6, Fig. 3g). We hypothesise that cells in state A exhibit a multilineage primed transcriptional programme with stochastic expression of metabolism associated genes. Upon fate commitment, cells in the IFE steadily increase their translational rate in a proliferation independent manner [33]. Hence, translation associated genes are subject to greater transcriptional regulation post commitment, independent of transcription level.

These data and our single-cell analysis identified an NCA Wnt-receptive subpopulation, state A, with greater dynamic range in gene expression (TCOV) and greater variation in the abundance of protein synthesis associated transcripts. Introduction of the NCA Wnt stimulus reduces variability in both aspects.

Transcription factors driving cell fate change

To understand drivers of the observed differential heterogeneity, we reconstructed transcriptional changes over time along the state trajectory. Expression of each gene was modelled as a non-linear function of pseudotransition time [12]. We found 632 genes that were dynamically regulated over the state transition (false discovery rate < 5%; Fig. 4a). Using hierarchical clustering, we grouped these genes into four patterns of dynamic expression. Group I genes, most highly expressed in state A, were enriched for methylation associated genes and histone modifiers such as Setd3 and Kdm7a. These genes represent a pre-transition transcriptional profile of state A cells without exposure to signalling beta-catenin induced cells. Group III genes show highest expression in state C, an intermediate transition state, with enrichment for desmosome genes such as Dsc2, Dsc3, Dsg2 and Dsp, which are most highly expressed in the suprabasal layers of murine epidermis indicative of early commitment to differentiation [2]. Group IV genes, predominantly expressed in state D, were enriched for protein synthesis associated genes and entry into mitosis, respectively. Notably, group III includes the transcription factor (TF) Klf5, a regulator of proliferation in intestinal epithelial cells [34], and E2f1, which leads to epidermal hyperplasia when overexpressed in mice [35].

Fig. 4 Induction of proliferation in Wnt + neighbours is regulated by 47 TFs. a Heatmap showing smoothed expression of pseudotransition-dependent genes (n = 632) ordered by hierarchical clustering and maximum expression. Top two enriched GO terms shown on left (all significant at q < 0.05). Genes (rows) are ordered by peak expression from state A to state D. b–d Heatmap showing expression of 47 TFs enriched as regulators of pseudotransition-dependent genes (b). Log-fold change in expression of TFs between state D and state A to show strongly directional TFs (c) and enrichment level of TFs (d). e Normalised proliferation index projected onto the cell state map. Arrow denotes the direction of pseudotransition. f–h Quantification of population proliferation by EdU assay in WT, K14ΔNβ-cateninER and co-cultured keratinocytes with and without 4OHT treatment (f). After stratification by Bcl3 nuclear abundance, cells in co-culture with activated K14ΔNβ-cateninER cells show a relative higher proliferation rate (g). Similarly, stratification by nuclear Smad4 shows higher proliferation in the treated co-culture condition (h) (n = 3 independent cultures). *P < 0.05. All data shown as mean ± SD Full size image

To gain insight into the regulation of the dynamically expressed genes induced by Wnt + cells, we performed a TF motif analysis (Fig. 4b–d). We calculated enrichment of TF binding sites from the ChEA database, removing TFs that were not expressed in any of the five cell states (log[TPM + 1] > 1). By analysing the promoters of the 632 dynamically expressed genes, we identified 47 TFs putatively regulating the state transition. TFs were separated into three groups according to the directionality of gene expression from state A to state D: positive; negative; and neutral. We noted that the activities of some identified TFs such as Smad3 and Smad4 are only partially dependent on expression level. Hence, we did not rule out TFs on the basis of expression.

From this analysis we predicted Smad3, Smad4, Kdm5b, E2f1 and E2f4 as previously unknown key regulators of the state transition, with Bcl3 as a likely regulator of the specific transition between states A and D. Known regulators of keratinocyte cell fate are shown in blue (Fig. 4d). Of note are Gata6 and Foxm1, two TFs upregulated in state D and previously shown to mark cells with multi-lineage differentiation potential and increased proliferative capacity, respectively [36,37,38]. We hypothesised that epidermal cells could be stratified based on the nuclear abundance of our identified novel state-regulating TFs, specifically Smad4 and Bcl3. Furthermore, our analysis on state D gene expression markers and transcriptome correlation indicated that this state shows a higher proliferation rate relative to states A and B. To investigate further, we calculated a cell proliferation index consisting of normalised expression of S-phase markers (Fig. 4e). This index demonstrated that states C and D comprise proliferative cells, with few cells in states A, B or E actively proliferating.

To confirm our findings, we used an EdU incorporation assay to distinguish proliferating cells and analysed whether keratinocytes positive for our predicted driver TFs (measured by nuclear intensity) were more likely to be proliferative. At the population level, there was no significant difference in proliferation when epidermal cells were co-cultured with 4OHT induced K14ΔNβ-cateninER cells (Fig. 4f). However, when cells were discriminated by nuclear intensity for Bcl3 or Smad4 we observed a significant difference between WT cells exposed to NCA Wnt signalling and WT or induced K14ΔNβ-cateninER cells alone (Fig. 4g and h). On average, 18% of Bcl3+ cells were positive for EdU uptake in WT or uninduced K14ΔNβ-cateninER cells; however, when WT cells were co-cultured with induced K14ΔNβ-cateninER cells the EdU+ fraction rose to 34%. The proportion of proliferative Smad4+ cells increased from 10–18% to 33% EdU+.

Taken together these results indicate that Bcl3 and Smad4 are specific markers of the epidermal state transition and mark cells moving along the trajectory between state A and state D (Fig. 1) during the first 24 h of exposure to a NCA Wnt signal.

NCA Wnt induced state transition is contact dependent

From our total population of 129 cells exposed to NCA Wnt signalling, cells in state E appeared to be unaffected. Our data suggest that state A comprises cells in a ‘responder’ state permissive to NCA Wnt signalling due to the presence of key TFs and a more heterogeneous gene expression programme. We sought to address whether the reduction in ribosome-related gene expression heterogeneity and the induced expression of transition TFs are contact- or distance-dependent. To answer this question, we labelled co-cultures of WT and K14ΔNβ-cateninER cells with a cell reporter of protein synthesis.

We measured global protein translation by assaying incorporation of O-propargyl-puromycin (OPP) and compared WT cells in contact with induced K14ΔNβ-cateninER cells to untreated control cells (Fig. 5a). We observed that WT cells showed higher translational activity when in contact with a Wnt + cell. To confirm this, we analysed the neighbours of over 10,000 Wnt + cells and compared the OPP fluorescence intensity distributions (Fig. 5b). We found a small but statistically significant increase in translation rate for both Wnt + cells and neighbour cells in the 4OHT treated condition. This suggested a contact-dependent mechanism for control of protein synthesis downstream of NCA Wnt signalling.

Fig. 5 Wnt + induced cell state transition is contact-dependent. a Appearance of typical WT cells and K14ΔNβ-cateninER co-cultures stained for OPP. b Left: Quantification of OPP in K14ΔNβ-cateninER cells and their neighbours in control and 4OHT treated conditions. Right: Similar quantification for cells co-cultured in low calcium medium (0.1 mM) to inhibit cell–cell contact (n = 3 independent cultures, n = 3 technical replicates, all pooled). c Appearance of typical WT cells and K14ΔNβ-cateninER co-cultures stained for Smad4. d Left: Quantification of Smad4 in K14ΔNβ-cateninER cells and their neighbours in control and 4OHT treated conditions. Right: Similar quantification for cells co-cultured in low calcium medium (0.1 mM) to inhibit cell–cell contact (n = 3 independent cultures, n = 3 technical replicates, all pooled). All scale bars, 20 μm, n.s. not significant Full size image

To confirm contact dependence and to rule out local diffusion of soluble factors, we repeated the assay in low calcium conditions. Keratinocytes cultured in low calcium medium do not form adherens or desmosomal cell contacts [39, 40]. Strikingly we observed no NCA Wnt effect under these conditions (Fig. 5b). Similarly, we observed no increase in nuclear abundance of Smad4, our predicted TF downstream of NCA Wnt signalling, in Wnt + cells, as predicted. However, in neighbouring cells there was a significant increase in nuclear Smad4 intensity, which is abrogated in low calcium conditions (Fig. 5c and d).

Taken together these data suggest that the Smad4-mediated cell state transition is downstream of non-cell autonomous Wnt signalling. Furthermore, the induction of this transition is contact dependent and does not occur under conditions where desmosomal adherens junctions are inhibited.