表皮干细胞动态中的性别二态性

Main

The skin epidermis is a compartmentalized tissue that is subdivided into the hair follicle, interfollicular epidermis (IFE), infundibulum and sebaceous gland. All of these cell types are maintained by distinct populations of stem cells (SCs)^11,12. The IFE is a stratified epithelium comprising a proliferative basal layer and multiple layers of differentiated suprabasal cells¹¹. Following terminal differentiation, post-mitotic basal cells delaminate into the suprabasal cell compartment before shedding from the skin surface¹¹.

Several models have been proposed to explain the cellular mechanisms that maintain homeostasis of squamous epithelia¹². From studies of rats, it was initially proposed that oesophageal epithelium is maintained by a single equipotent population in which cell loss through differentiation is compensated by the duplication of neighbours¹³. Subsequently, on the basis of histological and proliferation kinetic analyses of skin IFE in mice, others proposed that the IFE comprises a mosaic of epidermal proliferative units in which individual slow-cycling SCs support a pool of transit-amplifying (TA) cells and their differentiated progeny¹⁴. However, clonal lineage tracing challenged the idea of epidermal proliferative units. Instead, it was suggested that the IFE is maintained by a single proliferative compartment in which most (around 80%) of divisions result in an asymmetric fate outcome. The remaining compartment leads to balanced stochastic SC loss through differentiation or replacement through duplication^1,2,3. Further short-term intravital imaging of the epidermis supported a single proliferative compartment model, with spatiotemporal coordination between sister cells⁴ and local coordination of cell division and loss⁵.

However, several observations have challenged the one-compartment model. First, lineage tracing revealed differences in the long-term survival and growth potential of clonal populations using CreER lines linked to the expression of different targeted promoters^6,7,8,9,10. Second, fate-mapping studies of back skin epidermis under different experimental conditions showed evidence of an enrichment of clones with even cell numbers⁹. This feature is difficult to explain in the framework of a one-compartment model. Third, fate mapping of the ear skin epidermis based on a keratin 10 (Krt10) promoter showed evidence of proliferative cells with limited proliferative potential. This result suggests that subsequent commitment to differentiation cells maintains the capacity to divide, which requires multiple proliferative populations or a continuum of SC states¹⁵. Finally, single-cell RNA sequencing (scRNA-seq) of both mouse^{9,15,16,17,18,19,20,21,22} and human tissue^23,24 revealed heterogeneity in the basal cell population. However, the functional identity and significance of transcriptional states remained unclear.

Besides reproductive organs, the sex of an animal has been associated with differences in organ morphology and tissue composition^25,26,27. In skin, morphological differences between male and female animals have been reported^28,29, including impacts of hormones on wound healing^30,31,32. However, little is known about the role of sex or sex hormones on epidermal SC fate in homeostasis and repair.

To determine the functional identity of cell populations in the skin IFE, emphasis has been placed on clonal data obtained from fixed samples, in which temporal information on cell lineages can only be inferred indirectly. Here we examine the fate of individual lineage-labelled cells and their progeny in the ear skin epidermis by tracking the dynamics of clones in time lapse over weeks to months. By combining single-cell gene expression profiling with biophysical modelling, we show that IFE maintenance relies on a locally coordinated pattern of stochastic SC renewal in which commitment to differentiation is not direct but occurs through an intermediate TA-cell-like state with limited proliferative potential. Furthermore, detailed features of clonal dynamics reveal temporal coordination of sister cell events that is strongly enriched in female mice. Functional experiments link this coordination with the role of sex hormones in the regulation of basal cell dynamics.

Long-term imaging of IFE cells

To determine the repertoire of proliferative cell states, their lineage relationships and pattern of fate in the ear skin IFE, we used intravital microscopy to quantify the size, composition and spatial organization of clones over 5 months of time-lapse imaging. Low-dose tamoxifen administration to Krt14^creER;iChr-control-mosaic mice induced stochastic expression of one of three nuclear fluorescent reporters (mCherry, GFP or Cerulean), which were sufficiently spaced to enable unambiguous clone tracking by time-lapse imaging. By recording the spatial coordinates of clones, their dynamics could be monitored in a precise manner over weeks and months. Notably, it was possible to score the frequencies of events leading to clone loss (‘extinction’). This process enabled the discrimination of clones that become lost after one or two divisions from those that expand further before incidental loss (Fig. 1a–d and Extended Data Fig. 1a–c).

To generate a statistically meaningful dataset, we tracked the temporal evolution of 1,171 clones from 10 mice (5 females and 5 males) over a 140-day time course, with time points extending from once a week to once every 2 weeks then once a month (Fig. 1a–d and Supplementary Fig. 5a–j). After combining clonal data from all animals, irrespective of sex, the mean basal cell number of surviving clones steadily increased from 1.47 ± 0.03 cells at day 7 (D7) to 6.5 ± 0.6 cells at D84. By contrast, the clonal persistence—the probability that a clone retains at least one basal cell—decreased to 27 ± 2% at D84 (Extended Data Fig. 2a,b). Notably, the majority of clones visible at D7 (58 ± 4%, n = 653) did not expand beyond 3 cells over the 5-month time course. Of those that did expand beyond this threshold, around 23 ± 3% (n = 292) were subsequently lost through differentiation over the time course, whereas 19 ± 2% (n = 226) were still present at the final time point (Fig. 1d). The average labelled cell density—the product of mean clone size and clonal persistence—rapidly converged to a constant value, which indicated that Krt14^creER targets renewing cells in a representative manner³³ (Extended Data Fig. 2c).

Sex-specific epidermal cell fate

After grouping by sex, we observed differences in the average properties of clones between male and female mice, with the average basal clone size consistently higher in female mice (Fig. 1e–g and Supplementary Note). Moreover, the frequency of clones that stayed within a maximum of 3 cells was higher in male mice (65 ± 3%, n = 268) than in female mice (52 ± 6%, n = 437; z-test z = 3.4, P < 0.001). Both the basal and suprabasal cell number of clones showed a marked enrichment in even sizes in female mice. This enrichment persisted throughout the entire 5-month time course and was apparent even among the largest clones (Fig. 1h,k and Extended Data Fig. 2d–f). A similar enrichment pattern was observed in the back and tail skin IFE of female mice at D28 after labelling. This finding indicates that this pattern is a general feature of epidermal cells (Extended Data Fig. 2i–m).

To further validate these findings, we performed unbiased short-term lineage tracing using 5-ethynyl-2′-deoxyuridine (EdU) pulse–chase experiments. EdU is incorporated during S phase and diluted through subsequent rounds of division, which enables fate mapping of dividing cells. After a 6-day chase, the pseudo-clone size distribution revealed a clear enrichment of even-sized basal cell numbers, a result that mirrored the long-term time-lapse imaging results (Extended Data Fig. 2g,h).

Modelling reveals a proliferative hierarchy

To determine the functional heterogeneity and fate behaviour of basal cells, we questioned the origin of the persistent and robust even-odd signature of clone sizes and why the strength might differ between male and female animals. In anticipation of the role of local correlations in the coordination of cell division and delamination events, we considered a minimal spatial model of the basal layer based on a 2D Voronoi-type framework³⁴ (for details, see Supplementary Note and ref. ²²). Noting the cohesive nature of clones, we positioned the model in a solid-like regime. Here in the absence of turnover, cells are jammed so that neighbour exchange is suppressed³⁵.

To model the clonal data, we first considered a one-proliferative compartment model. To ensure homeostasis of the SC density, the frequency of symmetric divisions was subject to mechanical feedback based on an area-based rule. That is, when the local cell density (as indexed through the cell area) is low, SC divisions are tilted towards duplication. Conversly, when the cell density is high, SC divisions are biased towards terminal differentiation. Notably, although this model could provide a satisfactory fit to the observed dependencies on the average clone size, the hallmark even-odd clone-size correlations in female mice could not be reproduced. To capture the abrupt increase followed by the shallow linear-like growth in the average basal clone size (Fig. 1e), the frequency of SC divisions that result in asymmetric fate must be high (around 70–90%, as reported in other squamous tissues)^1,3,6. However, under such conditions, it was not possible to reproduce the scale of the observed enhancement of even clone sizes in both the basal and suprabasal compartments. Furthermore, the model could not capture key features of the data, including the enrichment of even-sized jumps in clone size across consecutive time points (Extended Data Fig. 3a–l and Supplementary Note).

On the basis of these findings, we considered a two-proliferative compartment model, in which SC differentiation occurs through an intermediate TA-cell-like state (Fig. 1l and Supplementary Videos 1 and 2). Placing SC fate under the same mechanical feedback control, we proposed that SC division can also lead to three possible fate outcomes. However, in this case, TA cells are generated. TA cell division then leads to either terminal division or asymmetric fate outcome, which then gives rise to terminally differentiated basal cells. Finally, to capture even clone-size correlations (Extended Data Fig. 2e), we further proposed that the timing of sister SC and TA cell division and the timing of sister cell delaminations are correlated (Supplementary Note).

To fit the two-compartment model, we first took advantage of the temporal nature of the data to constrain model parameters. To that end, we segregated lineages into the following types of clones: clones comprising a maximum of a single cell (D-like), which represent cells induced in the differentiated basal cell compartment; clones that grow to a maximum of three cells (TA-like), which are enriched in clones derived from the putative TA cell compartment; and the remaining (SC-like) clones, which are derived from the putative SC compartment. Measurements of the average size of SC-like clones over the short term provided an estimate of the cell proliferation rate. By contrast, the long-term labelled cell densities provided an estimate of the SC fraction in the basal layer and the ratio of cell loss to proliferation rate. With these estimates in hand, a parameter scan showed that the two-proliferative compartment model could reproduce the clone-size dependencies of both male and female mice. These dependencies included the persistent even-odd clone size correlations and the transition probabilities of clone sizes between consecutive time points (Fig. 1i–n, Extended Data Fig. 4a–i and Supplementary Note).

We then questioned the origin of the differences observed between male and female animals. Notably, the average residency time of terminally differentiated basal cells was predicted to be shorter in male mice (3.0 ± 0.4 days) than in female mice (5.1 ± 0.3 days). By contrast, the frequency of terminal TA cell divisions was higher in female mice (>95%) than in male mice (70%). Moreover, the delamination of differentiated sister cells was highly coordinated for female mice, with a correlation of sister cell lifetimes of about 0.9 compared with about 0.6 in male mice (Supplementary Note).

These findings point to a spatially and temporally coordinated pattern of turnover in the ear skin IFE, in which most SC divisions result in an asymmetric fate outcome and the majority of TA cell divisions result in symmetric differentiation. The temporal coordination of sister cell divisions and delaminations leads predominantly to a local turnover of SC–TA cell pairs, which results in the observed long-term enhancement of even clone sizes.

To validate the model predictions, we closely monitored cell dynamics. To that end, we imaged clones every 12 h for 5 days and we recorded division and delamination events to reconstruct lineage trees and to assess correlations between sister cell events. Estimates of correlation coefficients revealed a strong temporal correlation between sister cell lifetimes. We calculated correlation coefficients between sister events, with a coefficient of 0 indicating no correlation, and 1 indicating events that occur in the same 12-h time window. In male mice, we observed coefficients of 0.89 ± 0.03 between delaminating sisters, 0.70 ± 0.15 between dividing sisters and 0.75 ± 0.18 when one sister divided and one sister delaminated. In female mice, the corresponding coefficients were 0.82 ± 0.14, 0.8 ± 0.12 and 0.3 ± 0.3, respectively (although the latter events were rare). Within a 12-h time window, 84% (42 out of 50) delaminations in male clones and 98% (45 out of 46) of delaminations in female clones involved sister cells. Moreover, the relative transition probabilities between clone sizes across consecutive time points were similar in short-term and long-term imaging data. Repeated daily imaging slightly increased cell division rates compared with the less invasive twice a month imaging schedule used for long-term clonal analyses and biophysical modelling. Overall, short-term imaging confirmed the tight temporal coordination predicted by the model to occur between sister cells (Extended Data Fig. 5a–d and Supplementary Note).

Sex differences in wound healing

We next investigated whether sex-specific differences in cell fate affect wound-healing kinetics. After a wound was made on the back skin, the kinetics of wound closure was significantly faster in female mice than in male mice. To determine whether this dimorphism was the consequence of different responses of keratinocytes, we assessed proliferation rates using EdU and clone fate mapping. As we and others have shown^36,37, keratinocyte responses to wound healing are divided into two areas. The first is a non-proliferative leading edge consisting of only migratory cells and the second is a proliferative ring, where cells actively divide to mediate wound closure. At D2 after wounding, a 4-h EdU pulse revealed a thicker proliferative ring in female than in male mice (800 ± 37 µm and 590 ± 60 µm, respectively) and a significantly higher density of proliferative cells (4,330 ± 320 and 2,860 ± 330 EdU cells per mm², respectively). By D4, this sex difference was markedly reduced (Fig. 2a–f and Supplementary Fig. 6a,b). Consistent with this difference in proliferation kinetics, lineage tracing showed that at D2, basal clone size was significantly larger in female than male mice (4.7 ± 0.2 and 3.0 ± 0.1 cells, respectively), a difference still present, albeit decreased, at D4. This difference was reflected in the enrichment of even clone sizes in female mice at D4 after wounding compared with males (Fig. 2g–h). These findings indicate that wound healing promotes a higher degree of cell coordination in female than in male mice. Together with more rapid SC and/or TA cell activation in the initial stage of wound healing, this process accelerates wound closure.

Sex hormones regulate basal cell fate

The hierarchical model of epidermal fate predicted that following commitment, basal cells delaminate faster in male than female mice. To test this prediction, we performed short-term EdU pulse–chase lineage tracing. At 96 h after injection, EdU doublets (2 EdU⁺ cells in contact, isolated from other EdU⁺ cells) confirmed that delamination occurred faster in male mice than in female mice (47 ± 6% and 21 ± 3%, respectively, for doublets with one delaminated cell). Single-housed males did not show significant differences from controls, excluding fighting-induced effects (Fig. 3a,b).

To test whether male sex hormones drive this sexual dimorphism, we castrated male mice and assessed delamination rates via short-term EdU pulse–chase. After castration, delamination rates were similar to female mice (21 ± 3% in female mice, 14 ± 2% in castrated male mice and 47 ± 6% in male control mice for doublets with one delaminated cell). The asymmetric division rate (estimated as EdU doublets with one KRT14⁺ and one KRT10⁺ cell) was higher in male (74 ± 4%) than in female (60 ± 3%) mice and decreased after castration (53 ± 2%). Ovariectomy in female mice did not lead to differences in basal cell delamination rates or cell fate outcomes (Fig. 3b,c).

To assess whether Y chromosome genes, beyond male sex hormones, affect epidermal cell differentiation and delamination, we used the four core genotype (FCG) mouse model³⁸. This model involves deletion of the testis-determining gene Sry from the Y chromosome and introducing an autosomal Sry transgene. Consequently, gonadal sex depends on the transgene rather than the Y chromosome. After 96 h of EdU pulse–chase, Sry-expressing mice, irrespective of their chromosomal status, showed a significantly increased delamination rate (17 ± 2% for Sry⁺ versus 8 ± 1% for Sry^–). XX and XY animals did not show significant delamination differences (13 ± 2% and 11 ± 3%, respectively). Similar data were observed for fate outcome of dividing cells (KRT14⁺KRT10⁺ EdU doublets of 60 ± 3% for Sry⁺ versus 42 ± 4% for Sry^–) (Fig. 3d,e). To further clarify the role of male sex hormones, we assessed the impact of hormone replacement using dihydrotestosterone (DHT), a biologically active form of testosterone, in castrated male mice. We confirmed that androgen receptors are expressed by IFE basal cells in these mice. DHT administration increased basal cell delamination at 96 h (10 ± 1% for castrated mice compared with 31 ± 4% for DHT-treated mice) and similarly affected dividing cell fate outcome (KRT14⁺KRT10⁺ EdU doublets of 41 ± 4% for castrated mice versus 58 ± 3% for DHT-treated mice) (Fig. 3f–h). Together, these data show that testosterone promotes epidermal differentiation and basal cell delamination.

To confirm that male sex hormones influence clonal dynamics of epidermal SC and TA cells, we performed long-term lineage tracing in castrated male mice. Mean basal clone size, clone persistence and labelled cell fractions were similar between control and castrated mice. However, clones from castrated mice showed an enrichment of even cell numbers compared with non-castrated male mice. Fitting of the two-proliferative compartment model to the clonal data independently confirmed that male sex hormones affect basal cell delamination and asymmetric TA cell division frequency (Fig. 3i and Extended Data Fig. 6a–e,k), which led to the restoration of female-like features. Together, these data demonstrate that male hormones influence epidermal SC and TA cell clonal dynamics to enhance the rate of cell delamination following terminal differentiation.

To better understand how keratinocyte proliferation is controlled in male and female mice, we performed 4-h EdU pulse experiments during homeostasis. In the ear skin IFE, female mice presented a higher degree of proliferation than male mice (7.2 ± 0.3% and 5.6 ± 0.4%, respectively). Modulation of male sex hormones (castration then with or without DHT treatment) did not influence proliferation. However, IFE basal cells expressed oestrogen receptors, and modulation of female sex hormones by ovariectomy decreased basal cell proliferation (7.2 ± 0.3% for control and 4.9 ± 0.3% after ovariectomy). Oestradiol administration after ovariectomy restored the proliferation rate in female mice (7.2 ± 0.3% for control and 7.6 ± 0.7% after oestradiol treatment; Fig. 3j–l). Finally, to assess how female sex hormones affect cell coordination, we performed long-term lineage tracing in ovariectomized mice. Absence of female sex hormones led to a near-total disappearance of even basal clone size enrichment (Fig. 3m and Extended Data Fig. 6f–k). This result suggests that regulation of proliferation by oestrogen affects the female-specific increase in basal cell coordination and the degree of asymmetric TA cell division (Supplementary Note).

Molecular heterogeneity of basal cells

To further characterize ear skin IFE heterogeneity and the molecular drivers of sexual dimorphism, we performed scRNA-seq of FACS-isolated basal cells from adult male and female mice. As previously reported^16,17,22, three G0/G1 clusters were identified (Methods, Extended Data Fig. 7a–j and Supplementary Fig. 7a–e). To investigate basal cell lineage specification, we visualized cell trajectories on a three-dimensional (3D) manifold based on a diffusion map. The 3D manifold preserved the continuity of cell cycle progression, which enabled us to interpret how different G0/G1 populations are connected to proliferative states while maintaining global organization of the differentiation manifold. Unsupervised clustering revealed three non-cell-cycle-related clusters that echoed the three G0/G1 clusters—SC, TA and differentiated cells—identified through clonal fate mapping. The SC cluster preferentially expressed SC-related genes^39,40,41, such as Lgr6, Trp73, Sox6 and Jag2, whereas the TA cluster was defined by high levels of basal keratins (Krt5 and Krt14), Wnt4, S100a6, Actb and Fth1. The differentiated cluster expressed differentiation markers such as Krt10, Grhl1, Notch1 and Krtdap (Fig. 4a–d and Extended Data Fig. 7k,n–p). Cells from both the SC and TA cell clusters were able to enter the cell cycle loop, whereas most of the differentiated cells were post-mitotic. RNA velocity analysis indicated a developmental progression from the SC to the TA cell cluster, followed by cell cycle exit and differentiation. SCENIC analysis of activated transcription factor (TF) regulons revealed enrichment of the Trp73 regulon in SCs, Klf4 in TA cells and Grhl1, Gata3 and Klf4 among differentiated cells (Fig. 4e–g and Extended Data Fig. 8a–f).

To link the clonal fate mapping, modelling and scRNA-seq results, we performed lineage tracing using different CreER mouse models based on the following putative SC and TA cell marker genes: Lgr6^creER for SC and Notch1^creER for TA and differentiated cells. To capture early fate decisions and clonal loss in TA and differentiated cells, we started imaging at D3 instead of D7. After combining clonal data from all animals, we observed a substantial difference in D28 clonal persistence between Lgr6^creER (44.6 ± 4.0%) and Notch1^creER (10.0 ± 0.3%). This finding indicated that Lgr6^creER (SC-enriched) and Notch1^creER (TA and differentiated-enriched) label functionally discrete basal populations. To assess whether these targeted cells were rooted in self-renewing populations, we looked at the change in average labelled cell number. Notably, an early decrease in labelled cell number in the Notch1^creER-targeted population indicated that initially labelled proliferative cells were not enriched in cells with self-renewal potential. This result is consistent with the identity of a putative TA cell population. Subsequently, the labelled cell density stabilized, which is most likely to be caused by residual SC labelling by Notch1^creER. By contrast, the Lgr6^creER-targeted population showed an early increase in average labelled cell number, followed by a convergence to a constant value, a result consistent with the labelling of SCs (Fig. 4h–m).

When following a single clone over time, 75.5 ± 1.0% of clones did not show division in Notch1^creER tracing, whereas in Lgr6^creER tracing, non-dividing clones represented only 18.7 ± 5.1% of the total. Of those clones that did divide, the frequency of clones that stayed within the maximum 3-cell threshold (putative TA cells) was 62.1 ± 3.7% in Notch1^creER-targeted clones and 38.6 ± 5.7% in Lgr6^creER-targeted clones (Fig. 4n–o). These findings were consistent with the fact that Notch1^creER-labelled cells mainly comprised terminally differentiated and TA cells compared with Lgr6^creER-labelled cells. By segregating clones into SC-like, TA-like and D-like populations, Notch1 and Lgr6 clones were consistent with our two-proliferative compartment model (Extended Data Figs. 9a–j and 10a–n and Supplementary Note).

The fraction of labelled basal cells expressing KRT10 at D3 was 7.7 ± 2.9% in Lgr6^creER-targeted cells from male mice and 9.2 ± 1.0% from female mice. By contrast, in Notch1^creER-targeted cells, the proportion was 38.2 ± 4.5% in male mice and 55.3 ± 3.8% in female mice. The proportion of Lgr6^creER-labelled cells that were suprabasal was 0% in male mice and 1.4 ± 1.4% in female mice, compared with 27.4 ± 4.4% of Notch1^creER-labelled cells in male mice and 21.3 ± 1.4% in female mice (Extended Data Fig. 9k–n). Altogether, these findings support the existence of two proliferative basal cell compartments, with SCs preferentially targeted by Lgr6^creER and TA and differentiated cells preferentially targeted by Notch1^creER.

Sex-specific transcriptomes

To define the molecular mechanisms that underlie sex-specific differences in epidermal fate, we assessed the most differentially expressed genes between male and female mice across the different clusters. Besides Y-expressed genes, the IFE SC cluster from male mice showed increased levels of Skint1 and differentiation-related genes such as Ehf and Krtdap. The TA cell cluster from male mice showed increased levels of differentiation markers such as Krtdap or Skint1. Similarly, differentiated basal cells had increased expression of the differentiation marker Krtdap. By contrast, female mice showed elevated expression of SC marker genes such as Postn and Smoc2, as well as specific expression of immune-related genes such as Il34 and Il20ra. In line with these results, SCENIC analysis identified sex-based enrichment of potential TF regulon activation. Among SC, TA and differentiated clusters, the well-established TF associated with epidermal differentiation, Klf4 (ref. ⁴²), was enriched in basal cells from male compared with female mice (Extended Data Figs. 7l,m and 8g–l and Supplementary Data 1–3).

To investigate in more detail the transcriptional differences mediated by sex hormones, we performed bulk RNA-seq on FACS-isolated basal cells from male and female mice and castrated male mice. As sex hormone regulation of epidermal cells could occur indirectly, we assessed whether castration influences other skin parameters. Dermal thickness was unchanged after castration. Hair follicle cycle interference was limited as—in the ear epidermis—hair follicles do not cycle after morphogenesis⁴³ and castration did not alter their resting stages (Extended Data Fig. 11a–d). Bulk RNA-seq is more sensitive than scRNA-seq in detecting differences in RNA expression levels, with a higher signal-to-noise ratio. However, it has the disadvantage of reflecting the average composition of the different cell states present in the basal layer at any given time. Nevertheless, similar to the scRNA-seq data, the most differentially expressed genes were linked to sex chromosomes. Some previously described human sex-specific genes were identified. Notably, osteoglycin (Ogn), a proteoglycan previously reported to be more highly expressed in male than in female mice⁴⁴, was differentially expressed in epidermal basal cells between male and female mice. Several extracellular matrix (ECM) genes such as proteoglycan glypican 2 (Gpc2), collagens and collagen-modifying proteins, basal lamina-specific proteins such as Lamb1 and actin regulators such as Ablim3 and Fhl1, were expressed at higher levels in epidermal cells in male than in female mice (Fig. 5a and Extended Data Fig. 12a,b,e–j).

By contrast, many differentiation markers were upregulated in basal cells in female compared with male mice. Genes associated with terminal differentiation of keratinocytes, including Krt1, involucrin (Ivl), filaggrin (Flg), the transglutaminase Tgm1, the TF Ovol1 and cornifins (Sprr2a1, Sprr2a2, Sprr1a, Sprr1b), were all expressed at higher levels in females. Gene set enrichment analysis (GSEA) of sex-specific differentially expressed genes, using commitment-associated genes from epidermal basal cells in vivo⁶, confirmed that basal cells from female mice are associated with more epidermal differentiation genes. Notably, many of these differentiation genes were upregulated in male mice after castration (Fig. 5b,c and Extended Data Fig. 12e–m). Note that the apparent difference between the scRNA-seq and bulk RNA-seq data reflects the differences in the relative proportion of SC, TA and differentiated cells in the basal layer between male and female mice owing to the faster cell delamination rate in males, which occurs in an androgen-regulated manner. To further test the enrichment of terminally differentiated basal cells in female mice, we assessed KRT14 and KRT10 expression in the basal layer. The frequency of epidermal cells expressing KRT10 but not KRT14 was increased in females, a finding consistent with an elevated proportion of differentiated basal cells (Fig. 5d–g). Consistent with results from the EdU pulse experiments (Fig. 3k), female mice showed increased expression of proliferation-associated genes, such as Aspm, Birc5 and Top2a (Fig. 5b), which indicated that sexual dimorphism affects both differentiation and proliferation. Moreover, some immune-related genes were expressed in a sexually dimorphic manner and modulated by castration (Fig. 5a,b).

To define how male sex hormones regulate sex-specific transcriptional programs, we assessed how genes differentially expressed between the sexes were regulated after castration. About 10% (45 out of 434) of genes preferentially expressed in male mice were downregulated following castration. Several sex dimorphic ECM components (Ogn, Gpc2, Clu, Egfl6 and Lamb1), actin and immunity-related genes (Ablim3, Fhl1, Cxcl12 and Ccl27b) were downregulated after castration (Extended Data Fig. 12d). This result reflects the fact that regulation of such genes is mediated by testosterone.

Almost 50% of the genes upregulated following castration were upregulated in female mice, which suggests that male sex hormones negatively regulate genes expressed at increased levels in female mice (Extended Data Fig. 12c,k–m). Most of the keratinization genes (Flg, Tgm1, Ivl, Dsc2 and Lce1g), which reflect the proportion of differentiated cells in the basal layer, were upregulated after castration, consistent with the decreased delamination rate observed following male castration (Fig. 3b).

Testosterone promotes differentiation

To assess the human relevance of our analysis, we studied the impact of male sex hormones on the differentiation and stratification of human epidermal cells using in vitro 3D skin equivalent (3DSE) keratinocytes⁴⁵. Primary human keratinocytes from male and female donors were seeded on 3T3-J2 feeder cells mixed with collagen gels. Human keratinocytes were cultured in the presence or absence of 100 nM DHT. At D1 after keratinocyte seeding, cells were moved to an air–liquid interface to facilitate proper stratification of the skin equivalent. At D10, once the skin equivalent was properly stratified, a pulse of EdU was administered and chased for 72 h. No significant difference in cell proliferation was observed between control and DHT-treated cultures. However, following a 72-h chase, DHT-treated human keratinocytes presented an increased fraction of EdU-labelled suprabasal cells (37.2 ± 2.0% for DHT-treated and 24.6 ± 2.0% for control). This result shows that testosterone promotes delamination of human keratinocytes (Fig. 5h–j and Extended Data Fig. 12n).

Finally, to assess whether DHT promotes human keratinocyte differentiation, we performed in vitro culture of human-derived skin epithelioids⁴⁶. In contrast to the 3DSE model, skin epithelioids allow keratinocytes to be cultured without a fibroblast layer (Extended Data Fig. 12o). Epithelioids presented a lower level of stratification than traditional skin equivalents, which precluded the quantification of stratification using EdU pulse–chase experiments. However, this technique enabled us to study the initial steps of differentiation. After 2 weeks of culture, epithelioids were starved in hormone-free medium for 72 h and subsequently treated for 72 h in the same medium with or without 100 nM DHT. Whole epithelioids, including basal and suprabasal layers, were collected and analysed by RT–PCR. Out of the three individuals (two male and one female), two samples presented an increase in differentiation (marked by upregulation of EVPL, IVL, KRTDAP, OVOL1 and TGM1) (Fig. 5k). Altogether, these results demonstrate that the male hormone DHT promotes epidermal differentiation and stratification in the absence of stroma or immune cells.

Discussion

Here using long-term time-lapse imaging, we uncovered sex-specific differences in the clonal dynamics of epidermal SC and TA cells during adult homeostasis of the mouse ear skin. Our results showed that SC differentiation occurs through an intermediate TA-cell-like compartment together with a high degree of coordination between sister cells, both in the timing of cell division and in delamination. These findings are resonant with previous studies. In particular, clonal tracing based on the Krt10 promoter provided evidence for basal progenitors with limited proliferative potential¹⁵. Furthermore, continuous live-imaging of ear and paw skin epidermis showed evidence of temporal correlations in the delamination of sister cells⁴. Studies of the back skin epidermis under perturbed or regenerative conditions also showed evidence for an enhancement in even clone sizes, which indicated that there is the same proliferative heterogeneity in other epithelial contexts⁹. Here through long-term live-imaging, we showed how the enrichment of even clone sizes is integrated into individual lineage dependences and, in doing so, revealed substantial functional differences between male and female mice.

To develop a minimal model of SC–TA cell fate, we imposed a mechanism of SC density regulation based on mechanical feedback and correlated fate decisions with local cell density^47,48. However, we note that mechanisms involving competition for niche factors⁴⁹ or a combination of both biochemical and mechanical cues can lead to the same qualitative dynamics. Moreover, provided the timing of sister cell events remains temporally coordinated, we cannot rule out the potential of more complex models of cell fate involving, for example, further heterogeneity of the basal cell compartment or a low frequency of reversion between TA and SC identity.

Previous studies have shown the positive effect of oestrogen^30,31 and the negative effect of androgens on skin wound-healing³². These phenomena have been attributed to the effect of sex hormones on the ECM and inflammatory response³⁰. However, some reports have indicated that keratinocytes could be involved in this sexual dimorphism, as conditional deletion of the gene encoding oestrogen receptor-β (ESR2) in the epidermis phenocopies the delay in wound healing observed in mice in which Esr2 is completely knocked out⁵⁰. Here we showed that wounds in female mice heal more rapidly than in male mice owing to more rapid basal cell activation in the initial stage of repair and a higher degree of cell coordination.

Alongside the development of reproductive systems, sex-specific differences have been reported in the number of cells that comprise tissues, neuronal connectivity and skin thickness both in animal models and in humans^25,26,29. Sex-specific differences in SC functions have also been reported in several tissues such as the mouse haematopoietic system, the fly intestine and in the regenerative capacity of skeletal SCs^51,52,53. In the fly gut, the difference is intrinsically regulated by a sex-specific transcriptional program⁵², whereas haematopoietic SCs are regulated by sex hormones, such as oestrogen produced by the ovaries⁵³. In the context of mouse skin, where the epidermis is thicker in females whereas the dermis is thicker in males, the proliferation rate in the female epidermis is mediated by oestrogen received by the canonical oestrogen receptor (oestrogen receptor-α, encoded by Esr1)²⁸. Whether this sex difference phenotype is the consequence of a different clonal dynamic of epidermal SC and TA cells is currently unknown. Here we observed a notable difference in the proliferative capacity of TA cells and the rate and coordination of basal cell delamination between male and female mice. This dynamic behaviour of TA cells depended on male hormones, as orchidectomy abolished the differences in clonal dynamics between male and female mice. Altogether, these data and previous studies indicate that there is sexual dimorphism in the homeostatic regulation of epidermal maintenance.

Differential gene expression between male and female mice uncovered sex dimorphic genes that were differentially regulated across all basal cell states, including X and Y chromosome genes, as well as sex dimorphic genes that were preferentially expressed in male or female mice in a particular cell state. Notably, this dimorphism was not limited to keratinization but extended to proliferation, ECM modification, actin network and inflammatory responses, with the inflammatory responses echoing the sex dimorphic presentation of highly prevalent inflammatory skin diseases such as atopic dermatitis⁵⁴. Finally, most of the sex dimorphic differentiation genes were affected by male sex hormone deprivation after castration. Consistent with the delamination rate and clonal dynamic of castrated mice, molecular profiling in vivo and in vitro indicated the ability of male sex hormones to modulate the transcriptome and basal layer composition of IFE basal cells. Altogether, these results demonstrate the direct impact of male sex hormones on basal cell fate. However, we cannot rule out that immune cells or ECM components contribute to the observed sex differences, as they have been demonstrated to be sensitive to androgens^29,55.

In conclusion, our long-term clonal analysis in vivo in living animals demonstrated the hierarchical organization of proliferative states in the IFE and uncovered a previously unknown hormonally regulated sexual dimorphism both in cell fate and temporal coordination of sister cells during homeostasis and wound healing. As epidermal SCs have a key role in skin cancer formation, it will be interesting to study the impact of sex-specific regulation of SCs during clonal expansion, field cancerization and tumour initiation.

Methods

Mouse

Krt14^creER transgenic mice were provided by E. Fuchs⁵⁶. iChr-control-mosaic transgenic mice were provided by R. Benedito⁵⁷. Notch1^creER (N1-CreERT2^SAT) mice were provided by S. Fre⁵⁸. Lgr6-GFP-creER mice were obtained from The Jackson Laboratory³⁹. Tissues from FCG mice were provided by the Shwartz Laboratory; the original FCG mice were obtained from A. Arnold under a material transfer agreement³⁸. CD1 mice were obtained from Janvier Labs.

All animals were raised on a CD1 background, with the exception of FCG mice, which were on a C57BL/6 background. Equal numbers of male and female mice were used for every experiment unless otherwise specified. No statistical methods were used to predetermine sample sizes. Randomization was not applicable for most analyses because of the inability to assign sex. For experimental procedures (castration, ovariectomy or hormonal replacement), mice were randomly selected from the same litter and allocated to experimental groups to minimize litter-specific confounding effects. Sex could not be blinded during experimental procedures owing to inherent biological differences. Downstream data collection and analyses were performed in a blinded manner. Mouse colonies (with the exception of FCG mice) were maintained in a certified animal facility in accordance with European guidelines. The housing conditions of all mice strictly followed ethical regulations. The room temperature ranged from 20 to 24 °C. The relative ambient humidity at the level of mouse cages was 55 ± 10%. Each cage was provided with ad libitum food, water and enrichment material types of nesting material. A semi-natural light cycle of 12–12 h (7:00–19:00) was used. The experiments were approved by the local ethics committee Commission d’Éthique du Bien-Être Animal under protocols 794N, 851N and 912N. The FCG mouse colonies were maintained in a specific pathogen-free facility at the University of Southern California (Department of Animal Resources, BCC facility) in accordance with institutional guidelines. All procedures were approved by the University of Southern California Institutional Animal Care and Use Committee protocol 21489. Mice were housed in individually ventilated cages with ad libitum access to food and water, provided with nesting enrichment, humidity at 55 ± 10% and maintained on a 12–12 h light–dark cycle.

Clonal induction and lineage tracing

For lineage tracing under experimental and control conditions, mice (aged P60) were treated with 0.2 mg tamoxifen per 30 g body weight by intraperitoneal injection. Clonal labelling was checked at the first imaging time point by intravital microscopy. In vivo lineage tracing in mouse ear epidermis was performed as previously described¹⁶. In brief, induction of anaesthesia was performed using 4% isoflurane in a chamber. For imaging, mice were under continuous low-dose of anaesthesia of around 2% isoflurane throughout the course of imaging sessions in a heated chamber at 32 °C to minimize hypothermia. Mouse ear skin was depilated using depilatory cream to clearly image the hair follicle and IFE. To image the same area over time, a tattoo on the ear skin was made. Ear skin was tattooed using a small needle with carbon ink. Images were acquired on a LSM 880 confocal laser (Zeiss) fitted on an AxioObserver 7 NLO inverted microscope (Zeiss). Multiphoton excitation from a tunable InSight X3 laser (Spectra Physics) was set at 880 nm to image SHG, Cerulean and GFP signals using emission filters (Chroma) SP 475 nm, BP 494/20 and BP 525/45 nm, respectively. mCherry was imaged using the fixed 1,045 nm laser with BP 580–650 nm emission filters, in front of GaAsP non-descanned detectors (Zeiss). Serial optical sections were acquired through a Plan-Apochromat ×40/0.95 dry objective (Zeiss) with a 2.5 µm step to image the entire thickness of the epidermis. Images were analysed using Zen 3.10 software (Zeiss). Quantification of clone size and composition of basal cells were defined as cells in close contact with the SHG signal emerging from the collagen fibres of the dermis and used as a proxy of the basal lamina (Extended Data Fig. 1a). Suprabasal cells were defined as cells resting over basal cells with an intact nucleus. When identification was unclear, orthogonal projection was used to precisely assess cell positions in the IFE (Extended Data Fig. 1b,c). For back and tail skin lineage tracing, mice were induced at P60 with 1 mg tamoxifen per 30 g body weight by intraperitoneal injection. Imaging was performed on epidermal whole-mounts as described below.

Mouse treatments

To perform castration, mice were anaesthetized via inhalation with 4% isoflurane. Analgesic treatment was administered preoperatively by subcutaneous injection of carprofen (5 mg kg^–1; maximum of 100 µl) and buprenorphine (0.05 mg kg^–1; maximum of 100 µl). The scrotum and perineal area were disinfected using Iso-Betadine. A 1 cm skin incision was made along the median raphe. The tunica, as well as the various fascia, were opened by means of an incision. The testicles were exposed and the spermatic cord cauterized to free the testicle. The skin was closed with three stitches using Ethilon 5.0 (Ethicon). To perform ovariectomy, mice were anaesthetized as described above and the back skin was shaved and disinfected using Iso-Betadine. A median incision was made to the skin, and the peritoneum was opened bilaterally. The ovaries were exposed and the vasculature was cauterized to free the ovaries. The peritoneum was closed with a single stitch using Vicryl 5.0 (Ethicon) and the skin was closed with three stitches using Ethilon 5.0 (Ethicon). Where indicated, hormonal replacement was conducted by inserting a hormonal implant under the neck skin via a 2 mm incision, which was closed with a single stitch using Ethilon 5.0 (Ethicon). DHT implants (DHT-M/60, 40 µg per 24 h, Belma Technologies) or 17β-oestradiol implants (E2-M/60, 1.6 µg per 24 h, Belma Technologies) were inserted at the time of surgery.

Wound healing

Mice were anaesthetized via inhalation with 4% isoflurane. Analgesic treatment was administered preoperatively by subcutaneous injection of carprofen (5 mg kg^–1; maximum 100 µl) and buprenorphine (0.05 mg kg^–1; maximum 100 µl). For macroscopic wound closure measurement, mice were wounded on shaved and disinfected back skin using a 6 mm biopsy punch (Kai). Images were acquired using an Axio Zoom.V16 microscope (Zeiss) with a PlanApo Z ×1.0/0.25 FWD 60 mm objective (Zeiss) at days 0, 2, 4, 7, 9 and 12 after wounding. Mice were anaesthetized for each image acquisition session. For lineage tracing, mice were induced with an intraperitoneal injection of tamoxifen (5 mg per 30 g body weight). At D7 after tamoxifen induction, mice were wounded on shaved and disinfected back skin using a 4 mm biopsy punch (Kai). On D2 or D4 after wounding, mice were injected with a single intraperitoneal injection of EdU (2.5 mg ml^–1 in PBS, 100 µl per 20 g) and euthanized after 4 h. Wounds were collected and processed as described above (whole mount). Staining was done for GFP clones as described below.

Proliferation and fate determination

To follow the fate of cells after division, mice were injected with a single intraperitoneal injection of EdU (2.5 mg ml^–1 in PBS, 100 µl per 20 g) at around 8:00. Mice were killed after a prespecified chase time (4 h, 96 h or 144 h). Tissue samples were collected and processed as described below (whole mount).

Epidermal whole-mount and immunostaining

Ear skin was dissected from experimental and control mice, cut into small pieces and incubated in 20 mM PBS–EDTA on a rocking plate at 37 °C for 1 h 15 min. The epidermis was separated from the dermis using forceps. The samples were washed twice with PBS and then fixed in 4% paraformaldehyde for 1 h at room temperature. They were then washed 3 times with PBS (5 min per wash) and then stored in PBS containing sodium azide (0.2%) at 4 °C. For immunofluorescence staining, pieces of epidermal sheet were used and incubated in blocking buffer (1% BSA, 5% horse serum and 0.8% Triton in PBS) for 3 h at room temperature on a rocking plate. For primary antibody staining, samples were incubated with the primary antibodies KRT14 (chicken, Thermo Fisher Scientific, custom batch, 1:2,000), KRT10 (rabbit, BioLegend, 905404, 1:1,000), GFP (goat, Abcam, ab6673, 1:3,000) and mCherry (rat, Invitrogen, M11217, 1:3,000) overnight at room temperature on a rocking plate (100 rpm). Samples were then washed 3 times with PBS containing 0.2% Tween for 20 min each and then incubated with secondary antibody and Hoechst (1:2,000) for nuclear staining for 2 h at room temperature on a rocking plate. The following secondary antibodies were used (dilution of 1:400): anti-rabbit and anti-chicken conjugated either to Alexa Fluor 488 (Molecular Probes) or Rhodamine Red-X (Jackson ImmunoResearch) or cyanine 5 (Cy5) (Jackson ImmunoResearch). Samples were washed 3 times with PBS 0.2% Tween for 20 min each and then mounted in DAKO mounting medium by keeping the basal layer against the coverslip. For EdU staining, permeabilization and EdU labelling were performed after secondary antibody staining following the protocol of the Click-iT EdU Cell Proliferation Kit for Imaging (Thermo Fisher Scientific). Hoechst (1:2,000) was added during the last wash of the Click-iT protocol.

Histology and immunohistochemistry

Ear skin samples were fixed overnight in 4% paraformaldehyde at 4 °C and routinely embedded in paraffin wax. Sections (4-μm thick) were placed on Superfrost glass microscope slides (Gerhard Menzel) and dried at 37 °C. Immunohistochemistry for the detection of androgen receptors was performed using a rabbit monoclonal antibody (rabbit, Abcam, ab108341, 1:500) and detection of oestrogen receptors was performed using a rabbit monoclonal antibody (rabbit, Cell Signaling, 13258, 1:200). In brief, deparaffinized and hydrated tissue sections were incubated at 97 °C in citrate buffer for 20 min for antigen retrieval and were treated with 3% H₂O₂ for 10 min to quench endogenous peroxidase. Tissue sections were incubated overnight with the primary antibody, washed and incubated for 30 min using a Vectastain Elite ABC-HRP kit (peroxidase, rabbit IgG) (Vector Laboratories). Immunolabelling was visualized by incubation with diaminobenzidine and hydrogen peroxide. Haematoxylin and eosin staining was routinely performed following similar deparaffinization steps. Slides were imaged using an Axioscan 7 microscope (Zeiss) with a Plan-Apochromat ×40/0.95 Korr M27 objective.

3DSE

Human healthy skin biopsy samples were obtained as anonymized surgical waste (abdominoplasty or mammoplasty) after informed consent (Ethics Committee of Area Vasta Emilia Nord, approval number 178/09 for healthy donors). In brief, skin biopsy samples were minced and treated with 0.05% trypsin and 0.01% EDTA for 4 h at 37 °C. Every 30 min, keratinocytes were collected, plated (2.5–3 × 10⁴ per cm²) on lethally irradiated 3T3-J2 cells (2.4 × 10⁴ per cm²) and grown at 37 °C, 5% CO₂ in a humidified atmosphere in keratinocyte growth medium (KGM)⁵⁹.

Primary human keratinocytes (3.5 × 10⁵ cells per scaffold) from male and female donors were seeded on 3T3-J2 cells (mouse fibroblast, 2 × 10⁵ cells per scaffold) mixed with collagen gels (final concentration of 4 mg ml^–1) in hormone-free keratinocyte growth medium (hfKGM) as previously described⁵⁹. hfKGM was composed of phenol-red-free DMEM (Gibco, 3105302) and phenol-red-free DMEM/F12 (Gibco, 21041025) (2:1 mixture) containing charcoal-stripped fetal bovine serum (10%) (A&E Scientific, FBS-CS-12A), 5 mg ml^–1 insulin (Lilly, HI0210), 0.18 mM adenine (Merck, 1152 CAS: 73-24-5), 0.4 mg ml^–1 hydrocortisone (Acros, AC352450010 CAS: 50-23-7), 0.1 nM cholera toxin (List Biological, 9100B), 2 nM triiodothyronine (liothyronine sodium) (Merck, T2877 Sigma), 10 ng ml^–1 EGF (Austral Biologicals, GF-010-8), 4 mM glutamine (Euroclone, ECB3000D) and 50 IU ml^–1 penicillin–streptomycin (Euroclone, ECB3001D). hfKGM was supplemented with 100 nM DHT (Merck, 31573) from initiation of the culture or with the same amount of ethanol for the control group. At D1 after keratinocyte seeding, samples were moved to an air–liquid interface to allow for proper stratification of the skin equivalent. At D9 after air lift, once the 3DSE was properly stratified, a pulse of 2 h of 10 μM EdU (in hfKGM and then wash with PBS) was performed and chased for a duration of 72 h. 3DSE samples were then fixed and embedded in paraffin, sectioned and stained for EdU for delamination quantifications as described above. Slides were imaged using an Axioscan 7 microscope (Zeiss) with an EC Plan-Neofluar ×20/0.50 objective.

Human epithelioids

Ethics approval for cultures of human skin samples was obtained from the Commissie voor Medische Ethiek U(Z) Gent (ONZ-2025-0531). Skin samples were collected during elective facial plastic surgery and preserved in cold PBS until processing. Samples were then cut into explants, placed on membrane inserts and cultured in cFAD medium (Dulbecco’s modified Eagle’s medium/Nutrient Mixture F12 at a ratio of 3.5:1.1 (Bio&Sell) supplemented with 5% of chelated (Bio-Rad, 1422832) fetal bovine serum (Sigma-Aldrich, S0615), 5% penicillin–streptomycin (Thermo Fisher Scientific, 15140122), 5 μg ml^–1 insulin (Sigma-Aldrich, I3536), 1.8 × 10⁻⁴ M adenine (Sigma-Aldrich, A2786), 1 × 10⁻¹⁰ M cholera toxin (Sigma-Aldrich C8052), 10 ng ml^–1 epidermal growth factor (Thermo Fisher Scientific, AF-100-15-500UG), 0.5 μg ml^–1 hydrocortisone (Sigma, H4001), 5 μg ml^–1 apo-transferrin (Merk, T2036-100MG), 40 µg ml^–1 gentamicin (Thermo Fisher Scientific, 15710049), 100 µg ml^–1 normocin (SAS InvivoGen, ant-nr-1), 100 µg ml^–1 plasmocin (SAS InvivoGen, ant-mpp) and 10 µM of the ROCK inhibitor Y-27632 (Stemcell Technologies, 72302)). Three to four skin explants (collected with a 2.5 mm biopsy needle) were placed on top of a transparent ThinCert 0.4 μm pore-size 12-well insert (Greiner Bio-One, 665641) with the epithelium facing upward. cFAD was added to the top and bottom compartments of the inserts, which were then incubated at 37 °C and 5% CO₂. Y-27632 was kept at 10 µM for the first 7 days, then decreased to 5 µM between D8 and D14, 2.5 µM between D15 and D21 and discontinued after confluency was reached. Once keratinocytes had migrated out of the explant and reached confluency, the explant was removed. The epithelioids were subsequently placed in an air–liquid interface (no medium in the insert) for around 14 days until further processing. To study the impact of DHT, confluent air–liquid interface epithelioids were cultured in hormone-free cFAD, in which FBS was replaced with charcoal-stripped FBS (A&E Scientific, FBS-CS-12A) and phenol-red-free DMEM (Gibco, 31053028) and DMEM/F12 (Gibco, 21041025) were used. Epithelioids were placed for 3 days in hormone-free medium followed by 3 days of hormone-free cFAD supplemented with 100 nM DHT (Merck, 31573) or hormone-free medium alone. Total RNA of epithelioids was extracted by adding 346.5 µl RLT (RNeasy micro kit, Qiagen) with 3.5 µl β-mercaptoethanol for 10 min on ice, collected, vortexed and snap frozen. All epithelioid experiments were performed in technical replicates using punch biopsy samples collected from slightly different surgical areas from the same patient.

Microscope image acquisition and analysis

All confocal images for whole-mount of epidermis were acquired at room temperature with a LSM780 confocal system fitted on an Axio Examiner.Z1 upright microscope equipped with C-Apochromat ×40/1.1 water-immersion objectives (Zeiss). Optical sections of 512 × 512 pixels with 1.5-µm step z stacks were collected sequentially for each fluorochrome. The datasets generated were merged and displayed using ZEN 3.10 software (Zeiss). Quantifications of 4-h EdU pulse experiments were performed using QuPath (v.0.6.0) with maximal intensity projection images and automated cell detection. Default parameters were used, except that cell expansion was modified to 2 µm and the intensity threshold was manually adjusted for correct segmentation of each image.

FACS analysis and cell sorting

To isolate skin epidermal cells, ears were dissected and opened in two along the cartilage with forceps. Fat tissue and cartilage were removed, and skin was floated on a solution of 0.25% trypsin and EDTA 1 mM (Gibco, Thermo Fisher Scientific) in DMEM (Gibco, Thermo Fisher Scientific) overnight at 4 °C. The day after, the epidermis was carefully separated from the dermis by using forceps and then floated on a rocking plate (100 rpm) at 37 °C for 10 min on the trypsin mix. Tissue was cut into pieces, and epidermal cells were mechanically dissociated by pipetting with a 5 ml pipette. Trypsin was neutralized by adding DMEM supplemented with 5% Chelex-treated fetal calf serum (CFCS), and samples were filtrated through a 70-μm filter then through a 40-μm filter and washed in PBS supplemented with 2% CFCS (FACS buffer). Single cells were incubated with a biotinylated anti-CD34 primary antibody (clone RAM34, eBiosciences, 13-0341-82, 1:50) diluted in FACS buffer, incubated for 30 min on ice and protected from light. Cells were washed with FACS buffer then incubated for 30 min on ice with streptavidin–APC secondary antibody (BD Biosciences, 554067, 1:100) and CD49f (integrin-α6) PE-Cyanine7 conjugated antibody (clone eBioGoH3 (GoH3), Thermo Fisher Scientific, 25-0495-82, 1:250). The cells were washed again and resuspended in a Hoechst (10 mg ml^–1) solution diluted at 1:4 000 in FACS buffer. Cell sorting was performed on a BD FACSAria III using FACSDiva software (BD Bioscience) at the ULB Flow Cytometry Platform. Live epidermal cells were gated by forward scatter, side scatter and negative for Hoechst. Bulge (hair follicle) cells were excluded by CD34 staining. Only basal cells (integrin-α6 positive) cells were sorted (Supplementary Fig. 7). A total of 100,000 cells were sorted for each sample (1 sample containing 2–3 different mice from the same litter).

scRNA-seq

Single-cell suspensions from two independent male mice and two independent female mice were loaded and multiplexed on a GEM-X OCM 3′ Chip kit (v.4 4-plex, PN-1000747, 10x Genomics) and barcoded with a Chromium X instrument according to the manufacturer’s recommendations (10x Genomics).

Libraries were prepared using a GEM-X Universal 3′ Kit (v4 4-plex, PN-1000779, 10x Genomics) and sequenced on an Illumina Novaseq 6000 (paired-end 100 bp reads).

Single-cell transcriptomic data analysis

Sequencing reads were aligned and annotated with the mm10-2020-A reference dataset as provided by 10x Genomics and demultiplexed using CellRanger (v.7.0.0)⁶⁰ with default parameters.

Quality control and downstream analyses were performed using the R package Seurat (v.5.1.0)⁶¹. For each sample, cells were filtered using the following criteria: more than 2,000 and fewer than 6,000 unique genes; less than 7.5% UMI counts belonging to mitochondrial sequences; and less than 5% UMI counts belonging to ribosomal sequences. Doublets were detected using DoubletFinder (v.2.0.6)⁶² and scds (v.1.26.1)⁶³ with default parameters. Cells with a doublet score on scds greater than or equal to 0.75 and labelled as doublet on DoubletFinder were filtered out.

To identify and retain IFE basal cells for downstream analyses, each sample was first processed and analysed independently. Read counts were normalized using the NormalizeData() function of Seurat, with the parameters ‘normalization.method = “LogNormalize” and scale.factor=10000’. A principal component analysis (PCA) for each sample was calculated using scaled expression data of the 2,000 most variable genes (identified as outliers on a mean/variability plot, implemented in the FindVariableFeatures() function). Uniform manifold approximation and projection (UMAP) calculation and graph-based clustering were performed for each sample using Seurat (default parameters) with the first 40 respective PCA components as input. The final clustering resolution for each sample was chosen to best reflect the underlying biological heterogeneity. On the basis of canonical marker gene expression, clusters corresponding to hair follicle SCs (Lgr5 and Lhx2) and infundibulum cells (Lrig1 and Sox9) were identified and excluded, whereas the remaining cells were retained as IFE basal cells for downstream analyses. Cell cycle phase scores (S phase and G2/M phase) were computed for each sample using the CellCycleScoring() function implemented in Seurat.

For comparison of transcriptomic profiles between male and female mice, the retained IFE basal cells from all samples were integrated using the standard CCA-MNN-based data integration workflow in Seurat. Feature selection was performed using the FindVariableFeatures() function from Seurat with default parameters, selecting the 2,000 most variable genes. Canonical correlation analysis (CCA) followed by integration anchor selection were then performed on the selected features using the FindIntegrationAnchors() function from Seurat, taking the first 40 dimensions from the CCA into account, as described in the standard workflow of Seurat. These anchors were then used to integrate data with the IntegrateData() function. This CCA-MNN workflow was used for two-dimensional visualization and cell-state annotation to maintain consistency with our previously established mouse IFE analysis framework^16,22, which enabled direct comparison with previous studies while preserving the established annotation of G0/G1 A, G0/G1 B, G0/G1 C, S, G2/M, committed and differentiated cell states (Extended Data Fig. 7d). The assessment of gene expression was based on known marker genes from previous scRNA-seq data^9,16,17,22 together with cluster-specific differentially expressed genes. In particular, the G0/G1-A cluster expressed stemness markers such as Lgr6, Trp73 or Smoc2, G0/G1-B was defined by S100a6 and high levels of Krt14, and G0/G1-C showed expression of Cdk8 or immune-related genes such as Il31ra. Committed and differentiated clusters expressed a gradient of differentiation genes with higher levels of Notch receptors (Notch1 and Notch2), Notch target genes (Hes1) and other TFs that regulate epidermal differentiation (Grhl1, Foxn1 and Gata3). Differentiated keratinocytes were characterized by upregulation of terminal differentiation makers, including Zfp750, Klf4, Irf6, Krt1, Krt10, Krtdap, Dmkn and Sbsn (Extended Data Fig. 7a–j).

Marker genes for each cluster were identified by Wilcoxon rank-sum test using the FindMarkerGenes() function in Seurat with Benjamini–Hochberg correction for multiple testing. Only genes expressed in at least 25% of cells with an average log₂ fold change of at least 0.35 and an adjusted P < 0.01 are reported. A small cell-dissociation-induced cluster marked by elevated levels of Jun, Fosb and Atf3 was subsequently removed from downstream analyses⁶⁴ (Extended Data Fig. 7a).

To investigate the global organization of epidermal cell state transitions, we constructed a 3D diffusion manifold using the standard analysis workflow implemented in Pegasus (v.1.9.0)⁶⁵. Highly variable genes were selected for per sample and used to calculate the first 50 principal components (PCs), which were batch-corrected using the run_harmony() function of Pegasus without cell cycle regression. Harmony was used in this context because it is natively integrated in the Pegasus workflow and provides efficient batch correction before construction of the diffusion map. The Harmony-corrected PCs were then used to compute 100 diffusion components via the diffusion map method (diffmap function of Pegasus). After excluding the first trivial diffusion component, the remaining components were embedded into three dimensions using the force-directed layout algorithm (fle function of Pegasus) to visualize the topology of the differentiation manifold while preserving cell cycle progression, which enabled interpretation of the relationships between cell identity, proliferation and lineage progression. The Python package scVelo (v.0.3.2)⁶⁶ was used to run RNA velocity analysis, which was visualized on a force-directed layout manifold. The Python package dynamo (v.1.4.1)⁶⁷ was used to estimate both the speed and acceleration of cellular transcriptomic dynamics.

To enable direct comparison with the three-state clonal modelling framework, Leiden clustering (resolution of 0.5), implemented in Scanpy (v.1.9.8)⁶⁸, was applied to Harmony-corrected PCs after cell cycle regression, partitioning the data into three major biological populations. This approach ensured that clustering is driven by transcriptional identity rather than cell cycle state or batch effects, which enabled the identification of biologically meaningful compartments independent of proliferation status. On the basis of canonical marker gene expression (for example, Trp73, Klf4 and Krt10) and transcriptomic profiles, these populations correspond to SC, TA and differentiated states, respectively. For downstream differential expression analyses, cells annotated as S and G2/M phases were retained but were assigned all cells to their corresponding biological populations (SC, TA and differentiated) based on the cell-cycle-regressed clustering, thereby minimizing confounding effects from cell-cycle-associated transcriptional programs.

Differentially expressed genes (DEGs) defining individual cell populations (SC, TA and differentiated) were first identified at the single-cell level using the Wilcoxon rank-sum test implemented in the “rank_genes_groups” function in Scanpy. Benjamini–Hochberg correction was applied to control the FDR, and genes expressed in at least 25% of cells (min.pct = 0.25), with a log₂ fold change > 0.5 (logfc.threshold = 0.5) and an adjusted P < 0.01, were considered differentially expressed. To leverage biological replicates and to improve the robustness of marker identification, we also performed pseudobulk differential expression analysis using both edgeR and DESeq2 with the same significance criteria (log₂ fold change > 0.5 and adjusted P < 0.01). Cell-type marker genes were defined on the basis of the consensus of these complementary analyses, thereby increasing the robustness of the identified marker genes. Given the substantially larger number of genes tested and the primary goal of defining robust cell-type markers, a more stringent significance threshold (adjusted P < 0.01) was applied.

DEGs between male and female mice were identified independently in each cell type (SC, TA and differentiated) using the Wilcoxon rank-sum test implemented in the “rank_genes_groups” function in Scanpy, followed by Benjamini–Hochberg correction for multiple testing. Given that sex-associated transcriptional differences were expected to be more subtle than differences between major cell populations, a less stringent adjusted P value threshold (0.05) was used to detect DEGs.

The R package SCENIC (v.1.3.1)⁶⁹ was used to infer regulons—sets of TFs and their putative target genes—based on co-expression and DNA motif enrichment, following the standard SCENIC workflow. Regulon activity scores were computed at the single-cell level using the AUCell module implemented in SCENIC. Differential regulon activity was assessed using the R function t.test, comparing both across cell types and between male and female mice in each cell type, followed by Benjamini–Hochberg correction. Because differential regulon activity was evaluated across a substantially smaller number of features (221 inferred regulons) than gene-level differential expression analyses, an adjusted P < 0.05 was used to define significant differences in regulon activity.

Bulk RNA-seq

FACS-isolated cells were collected into 350 µl lysis buffer (RLT buffer, RNeasy Microkit, Qiagen) supplemented with β-mercaptoethanol. RNA extraction was then carried out using a RNeasy Micro kit (Qiagen) according to the manufacturer’s recommendations. Before sequencing, the quality of RNA was evaluated by Fragment Analyzer (Agilent).

Indexed cDNA libraries were obtained using a NEBNext Single Cell/Low Input RNA Library Prep kit for Illumina (NEB E6420S) following the manufacturer’s recommendations. The multiplexed libraries were loaded on flow cells and sequences were produced using a NovaSeq 6000 S2 Reagent kit (200 cycles from Novaseq 6000 System, Illumina) on a NovaSeq 6000 system (Illumina).

The quality of each raw data was checked using FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Adaptor sequences and low-quality regions were trimmed using TrimmomaticPE (v.0.39)⁷⁰. Trimmed reads were mapped against the mouse reference genome (Grcm38/mm10) using STAR aligner (v.2.7.11b)⁷¹ to generate read alignments for each sample. Reference genome sequences and corresponding annotations for Grcm38.87 were obtained from ftp.Ensembl.org. Duplicate reads were removed using the Picard MarkDuplicates module (http://broadinstitute.github.io/picard/). After transcript assembly, gene-level counts were obtained using HTseq (v.0.11.1)⁷² and normalized to 20 million aligned reads per count. The average expression of each gene for the different cell populations was computed based on at least two biological replicates, and fold changes were calculated between the subpopulations.

For DEGs across the conditions, genes with a fold change of expression ≥ 2 were considered as upregulated and those with a fold change of expression ≤ 0.5 were considered downregulated. Genes were excluded from the differential expression analysis if they had more than a threefold difference in expression between replicates. Gene expression of less than 10 normalized read counts in at least one individual were not considered for differential expression.

Differential expression analysis was then performed using the R package DESeq2 (v.1.46.0)⁷³. Genes with very low expression were excluded by retaining only those with a total count greater than ten across all samples to reduce noise and to improve statistical power. Normalization was carried out using the median-of-ratios method implemented in DESeq2 to account for differences in sequencing depth. To stabilize variance across the dynamic range of expression values, regularized log (rlog) transformation was applied to the normalized counts. Quality assessment of the dataset was performed using multiple complementary approaches. Pairwise comparisons between biological replicates were examined to evaluate reproducibility, and sample relationships were further assessed using Euclidean distance-based clustering and correlation heatmaps. PCA was conducted on the transformed data to visualize global transcriptional differences between conditions, which confirmed appropriate clustering of samples according to experimental groups. Differential expression testing was performed using a negative binomial generalized linear model framework, as implemented in DESeq2. Dispersion parameters were estimated using empirical Bayes shrinkage, and significance was assessed using Wald tests. Pairwise contrasts were defined between male and female, male and castrated, and female and castrated conditions. Resulting P values were adjusted for multiple testing using the Benjamini–Hochberg procedure, and genes with an adjusted P < 0.1 were considered differentially expressed. To improve the interpretability and stability of effect-size estimates, log₂ fold changes were further shrunk using the normal shrinkage estimator.

In parallel, differential expression analysis was independently performed using the edgeR (v.4.4.2)⁷⁴ package to validate the robustness of the results. A DGEList object was constructed from the raw count matrix, and lowly expressed genes were filtered by retaining those with counts per million greater than or equal to 1 in at least two samples. Library size normalization was performed using the trimmed mean of M-values method. A design matrix was constructed based on experimental groups, and gene-wise dispersions were estimated using empirical Bayes methods before fitting a generalized linear model. Differential expression testing in edgeR was conducted using likelihood ratio tests for the same pairwise contrasts defined in the DESeq2 analysis. Resulting P values were adjusted using the FDR method. Differentially expressed genes were ranked on the basis of significance and effect size and annotated with gene identifiers before export for downstream analyses.

To assess the consistency between differential expression results obtained from DESeq2 and edgeR, overlap analyses were performed using Venn diagrams. The significance of the overlap between gene sets was evaluated using a hypergeometric test, which estimates the probability of observing a given number of overlapping genes by chance, given the total number of genes tested. Specifically, for each pairwise comparison, the total number of genes included in the analysis was defined as the background population, and the sizes of gene sets identified by each method, as well as their observed intersection, were used to compute one-sided P values corresponding to enrichment of overlap. For a given comparison, the P value was calculated using the cumulative hypergeometric distribution as implemented in R (function phyper), testing the probability of observing at least k overlapping genes between the two methods.

RT–qPCR

Total RNA from epithelioids was extracted using 350 µl lysis buffer (RLT buffer, RNeasy Microkit, Qiagen) supplemented with β-mercaptoethanol. RNA extraction was then carried out using a RNeasy Micro kit (Qiagen) according to the manufacturer’s recommendations. gDNA was eliminated via on-column DNase digestion. cDNA synthesis was performed using SuperScript II reverse transcriptase (Invitrogen, 18064022) according to the manufacturer’s recommendations. RT–qPCR was conducted using KAPA SYBR FAST (Roche, KK4600) on a Light Cycler 96 (Roche) with 2.5 ng cDNA per reaction. Gene expression levels were normalized to the housekeeping gene Tbp for relative quantification. The following probes from Eurogentec were used: TBP_Fw: GAGCTGTGATGTGAAGTTTCC, TBP_Rv: TCTGGGTTTGATCATTCTGTAG, EVPL_Fw: GCCCTGAAGCAGAAATTGGC, EVPL_Rv: AGCAGTCGGAAGTCGGAATC, FLG_Fw: TGAGGGCACTGAAAGGCAAA, FLG_Rv: TGGCCACATAAACCTGGGTC, IVL_Fw: CTGCCCACAAAGGGAGAAGT, IVL_Rv: AGCGGACCCGAAATAAGTGG, KRTDAP_Fw: CGACAAATTGCGATCTGCGTT, KRTDAP_Rv: CGCTCCTCAGTCCTTTCAGC, OVOL1_Fw: CAGGGAAGACGGCGACATT, OVOL1_Rv: CTCATGTTCAAAGCCAGCGG, TGM1_Fw: TCCCATCCATCCTGACCTGT and TGM1_Rv: GTCGTCATCTGCCGCATTTC.

Statistical analysis and reproducibility

All statistical analyses were performed using R (v.4.4.1) or Python (v3.12.3). Data are expressed as the mean ± s.e.m. Frequencies recorded in main text are shown as the mean ± s.e.m. unless otherwise indicated. Normality was tested using Shapiro–Wilk tests. For datasets that followed normality, P values were estimated with two-sided t-tests. For datasets that did not follow normality, P values were calculated using Wilcoxon rank-sum tests. For clonal analysis and transcriptomic data, P values calculated using different tests were corrected with the Benjamini–Hochberg procedure, unless otherwise specified.

To assess the bias to even clone sizes relative to a geometric (discrete exponential distribution), we defined the following metric: \({\Delta }_{{\rm{even}}}=\frac{{F}_{\text{e}}(t)-{F}_{\text{g}}(t)}{1-{F}_{\text{g}}(t)}\). Here F_e(t) denotes the fraction of even-sized clones at time point t and \({F}_{{\rm{g}}}(t)=\frac{\langle n(t)\rangle -1}{2\langle n(t)\rangle -1}\) is the fraction of even clones expected for a geometric (that is, a discrete exponential) distribution with the same mean \(\langle n(t)\rangle \) as given by the data. When there is no statistical bias towards even sizes, Δeven, whereas a bias towards even sizes relative to the exponential distribution is reflected by a positive value of Δeven. The normalization constant is chosen such that if all clones had even sizes, Δ_even = 1.

For all the figures, the number of mice, the number of clones and the number of areas analysed are indicated in the legend. All replicates shown are from distinct biological samples. For plots of averages and distributions for experimental and simulated data, uncertainty intervals were calculated by bootstrapping clonal trajectories. For all of the images, each experiment was repeated independently with similar results at least two times.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability

All the raw sequencing data have been deposited into the Gene Expression Omnibus (GEO) with the accession numbers GSE294042 (bulk RNA-seq) and GSE294426 (scRNA-seq). For bulk RNA-seq, the mouse reference genome assembly GRCm38/mm10 was used throughout this study. Genome sequence and gene annotation files were obtained from Ensembl release 87 (https://ftp.ensembl.org/pub/release-87/), including the mouse genome FASTA files (https://ftp.ensembl.org/pub/release-87/fasta/mus_musculus/dna/) and GTF annotation files (https://ftp.ensembl.org/pub/release-87/gtf/mus_musculus/). For scRNA-seq, reads were aligned and gene expression was quantified using the 10x Genomics pre-built Cell Ranger reference package refdata-gex-mm10-2020-A. The reference package and its build specifications are available from 10x Genomics (https://www.10xgenomics.com/support/software/cell-ranger/downloads/cr-ref-build-steps). Source data are provided with this paper.

Code availability

Codes for numerical simulations of the one-proliferative and two-proliferative compartment models have been deposited into Zenodo (https://doi.org/10.5281/zenodo.21997809)⁷⁵.

References

1. Clayton, E. et al. A single type of progenitor cell maintains normal epidermis. Nature 446, 185–189 (2007).

   Article  ADS  PubMed  CAS  Google Scholar 

2. Doupé, D. P., Klein, A. M., Simons, B. D. & Jones, P. H. The ordered architecture of murine ear epidermis is maintained by progenitor cells with random fate. Dev. Cell 18, 317–323 (2010).

   Article  PubMed  Google Scholar 

3. Doupé, D. P. et al. A single progenitor population switches behavior to maintain and repair esophageal epithelium. Science 337, 1091–1093 (2012).

   Article  ADS  PubMed  PubMed Central  Google Scholar 

4. Rompolas, P. et al. Spatiotemporal coordination of stem cell commitment during epidermal homeostasis. Science 352, 1471–1474 (2016).

   Article  ADS  PubMed  PubMed Central  CAS  Google Scholar 

5. Mesa, K. R. et al. Homeostatic epidermal stem cell self-renewal is driven by local differentiation. Cell Stem Cell 23, 677–686 (2018).

   Article  PubMed  PubMed Central  CAS  Google Scholar 

6. Mascré, G. et al. Distinct contribution of stem and progenitor cells to epidermal maintenance. Nature 489, 257–262 (2012).

   Article  ADS  PubMed  Google Scholar 

7. Sánchez-Danés, A. et al. Defining the clonal dynamics leading to mouse skin tumour initiation. Nature 536, 298–303 (2016).

   Article  ADS  PubMed  PubMed Central  Google Scholar 

8. Sada, A. et al. Defining the cellular lineage hierarchy in the interfollicular epidermis of adult skin. Nat. Cell Biol. 18, 619–631 (2016).

   Article  PubMed  PubMed Central  CAS  Google Scholar 

9. Aragona, M. et al. Mechanisms of stretch-mediated skin expansion at single-cell resolution. Nature 584, 268–273 (2020).

   Article  ADS  PubMed  PubMed Central  CAS  Google Scholar 

10. Ghuwalewala, S., Jiang, K., Ragi, S., Shalloway, D. & Tumbar, T. A transit-amplifying progenitor with biphasic behavior contributes to epidermal renewal. Development 151, dev202389 (2024).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

11. Blanpain, C. & Fuchs, E. Epidermal homeostasis: a balancing act of stem cells in the skin. Nat. Rev. Mol. Cell Biol. 10, 207–217 (2009).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

12. Alcolea, M. P. & Jones, P. H. Lineage analysis of epidermal stem cells. Cold Spring Harb. Perspect. Med. 4, a015206 (2014).

    Article  PubMed  PubMed Central  Google Scholar 

13. Marques-Pereira, J. P. & Leblond, C. P. Mitosis and differentiation in the stratified squamous epithelium of the rat esophagus. Am. J. Anat. 117, 73–87 (1965).

    Article  PubMed  CAS  Google Scholar 

14. Potten, C. S. The epidermal proliferative unit: the possible role of the central basal cell. Cell Prolif. 7, 77–88 (1974).

    Article  CAS  Google Scholar 

15. Cockburn, K. et al. Gradual differentiation uncoupled from cell cycle exit generates heterogeneity in the epidermal stem cell layer. Nat. Cell Biol. 24, 1692–1700 (2022).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

16. Bansaccal, N. et al. The extracellular matrix dictates regional competence for tumour initiation. Nature 623, 828–835 (2023).

    Article  ADS  PubMed  PubMed Central  CAS  Google Scholar 

17. Dekoninck, S. et al. Defining the design principles of skin epidermis postnatal growth. Cell 181, 604–620 (2020).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

18. Ghuwalewala, S. et al. Binary organization of epidermal basal domains highlights robustness to environmental exposure. EMBO J. 41, e110488 (2022).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

19. Joost, S. et al. The molecular anatomy of mouse skin during hair growth and rest. Cell Stem Cell 26, 441–457 (2020).

    Article  PubMed  CAS  Google Scholar 

20. Joost, S. et al. Single-cell transcriptomics reveals that differentiation and spatial signatures shape epidermal and hair follicle heterogeneity. Cell Syst. 3, 221–237 (2016).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

21. Lin, Z. et al. Murine interfollicular epidermal differentiation is gradualistic with GRHL3 controlling progression from stem to transition cell states. Nat. Commun. 11, 5434 (2020).

    Article  ADS  PubMed  PubMed Central  CAS  Google Scholar 

22. Sarate, R. M. et al. Dynamic regulation of tissue fluidity controls skin repair during wound healing. Cell 187, 5298–5315 https://doi.org/10.1016/j.cell.2024.07.031 (2024).

    Article  PubMed  CAS  Google Scholar 

23. Cheng, J. B. et al. Transcriptional programming of normal and inflamed human epidermis at single-cell resolution. Cell Rep. 25, 871–883 (2018).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

24. Negri, V. A. et al. Single-cell RNA sequencing of human epidermis identifies Lunatic fringe as a novel regulator of the stem cell compartment. Stem Cell Rep. 18, 2047–2055 (2023).

    Article  CAS  Google Scholar 

25. McCarthy, M. M. & Arnold, A. P. Reframing sexual differentiation of the brain. Nat. Neurosci. 14, 677–683 (2011).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

26. Ransick, A. et al. Single-cell profiling reveals sex, lineage, and regional diversity in the mouse kidney. Dev. Cell 51, 399–413 (2019).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

27. Stankiewicz, L. N. et al. Sex-biased human thymic architecture guides T cell development through spatially defined niches. Dev. Cell 60, 152–169 (2025).

    Article  PubMed  CAS  Google Scholar 

28. Movérare, S., Lindberg, M. K., Faergemann, J., Gustafsson, J.-A. & Ohlsson, C. Estrogen receptor α, but not estrogen receptor β, is involved in the regulation of the hair follicle cycling as well as the thickness of epidermis in male mice. J. Invest. Dermatol. 119, 1053–1058 (2002).

    Article  PubMed  Google Scholar 

29. Azzi, L., El-Alfy, M., Martel, C. & Labrie, F. Gender differences in mouse skin morphology and specific effects of sex steroids and dehydroepiandrosterone. J. Invest. Dermatol. 124, 22–27 (2005).

    Article  PubMed  CAS  Google Scholar 

30. Ashcroft, G. S. et al. Estrogen accelerates cutaneous wound healing associated with an increase in TGF-β1 levels. Nat. Med. 3, 1209–1215 (1997).

    Article  PubMed  CAS  Google Scholar 

31. Ashcroft, G. S., Greenwell-Wild, T., Horan, M. A., Wahl, S. M. & Ferguson, M. W. Topical estrogen accelerates cutaneous wound healing in aged humans associated with an altered inflammatory response. Am. J. Pathol. 155, 1137–1146 (1999).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

32. Ashcroft, G. S. & Mills, S. J. Androgen receptor-mediated inhibition of cutaneous wound healing. J. Clin. Invest. 110, 615–624 (2002).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

33. Blanpain, C. & Simons, B. D. Unravelling stem cell dynamics by lineage tracing. Nat. Rev. Mol. Cell Biol. 14, 489–502 (2013).

    Article  PubMed  CAS  Google Scholar 

34. Alt, S., Ganguly, P. & Salbreux, G. Vertex models: from cell mechanics to tissue morphogenesis. Philos. Trans. R. Soc. B Biol. Sci. 372, 20150520 (2017).

    Article  Google Scholar 

35. Bi, D., Yang, X., Marchetti, M. C. & Manning, M. L. Motility-driven glass and jamming transitions in biological tissues. Phys. Rev. X 6, 021011 (2016).

    PubMed  PubMed Central  Google Scholar 

36. Aragona, M. et al. Defining stem cell dynamics and migration during wound healing in mouse skin epidermis. Nat. Commun. 8, 14684 (2017).

    Article  ADS  PubMed  PubMed Central  Google Scholar 

37. Park, S. et al. Tissue-scale coordination of cellular behaviour promotes epidermal wound repair in live mice. Nat. Cell Biol. 19, 155–163 (2017).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

38. De Vries, G. J. et al. A model system for study of sex chromosome effects on sexually dimorphic neural and behavioral traits. J. Neurosci. 22, 9005–9014 (2002).

    Article  PubMed  PubMed Central  Google Scholar 

39. Snippert, H. J. et al. Lgr6 marks stem cells in the hair follicle that generate all cell lineages of the skin. Science 327, 1385–1389 (2010).

    Article  ADS  PubMed  CAS  Google Scholar 

40. Füllgrabe, A. et al. Dynamics of Lgr6⁺ progenitor cells in the hair follicle, sebaceous gland, and interfollicular epidermis. Stem Cell Rep. 5, 843–855 (2015).

    Article  Google Scholar 

41. Huang, S. et al. Lgr6 marks epidermal stem cells with a nerve-dependent role in wound re-epithelialization. Cell Stem Cell 28, 1582–1596 (2021).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

42. Segre, J. A., Bauer, C. & Fuchs, E. Klf4 is a transcription factor required for establishing the barrier function of the skin. Nat. Genet. 22, 356–360 (1999).

    Article  PubMed  CAS  Google Scholar 

43. Wang, Q. et al. A multi-scale model for hair follicles reveals heterogeneous domains driving rapid spatiotemporal hair growth patterning. eLife 6, e22772 (2017).

    Article  PubMed  PubMed Central  Google Scholar 

44. Woessner, M. N. et al. Osteoglycin across the adult lifespan. J. Clin. Endocrinol. Metab. 107, e1426–e1433 (2022).

    Article  PubMed  Google Scholar 

45. Tham, K.-C., Lim, S. S., Bonnard, C. & Common, J. E. A. Protocol to reconstruct an in vitro 3D full-thickness skin equivalent model with collagen I in 6- and 12-well inserts. STAR Protoc. 6, 103658 (2025).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

46. Herms, A. et al. Self-sustaining long-term 3D epithelioid cultures reveal drivers of clonal expansion in esophageal epithelium. Nat. Genet. 56, 2158–2173 (2024).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

47. Hannezo, E. & Heisenberg, C.-P. Mechanochemical feedback loops in development and disease. Cell 178, 12–25 (2019).

    Article  PubMed  CAS  Google Scholar 

48. Yamaguchi, H., Kawaguchi, K. & Sagawa, T. Dynamical crossover in a stochastic model of cell fate decision. Phys. Rev. E 96, 012401 (2017).

    Article  ADS  PubMed  Google Scholar 

49. Kitadate, Y. et al. Competition for mitogens regulates spermatogenic stem cell homeostasis in an open niche. Cell Stem Cell 24, 79–92 (2019).

    Article  PubMed  CAS  Google Scholar 

50. Campbell, L. et al. Estrogen promotes cutaneous wound healing via estrogen receptor β independent of its antiinflammatory activities. J. Exp. Med. 207, 1825–1833 (2010).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

51. Andrew, T. W. et al. Sexually dimorphic estrogen sensing in skeletal stem cells controls skeletal regeneration. Nat. Commun. 13, 6491 (2022).

    Article  ADS  PubMed  PubMed Central  CAS  Google Scholar 

52. Hudry, B., Khadayate, S. & Miguel-Aliaga, I. The sexual identity of adult intestinal stem cells controls organ size and plasticity. Nature 530, 344–348 (2016).

    Article  ADS  PubMed  PubMed Central  CAS  Google Scholar 

53. Nakada, D. et al. Oestrogen increases haematopoietic stem-cell self-renewal in females and during pregnancy. Nature 505, 555–558 (2014).

    Article  ADS  PubMed  PubMed Central  CAS  Google Scholar 

54. Kanda, N., Hoashi, T. & Saeki, H. The roles of sex hormones in the course of atopic dermatitis. Int. J. Mol. Sci. 20, 4660 (2019).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

55. Chi, L. et al. Sexual dimorphism in skin immunity is mediated by an androgen–ILC2–dendritic cell axis. Science 384, eadk6200 (2024).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

56. Vasioukhin, V., Degenstein, L., Wise, B. & Fuchs, E. The magical touch: genome targeting in epidermal stem cells induced by tamoxifen application to mouse skin. Proc. Natl Acad. Sci. USA 96, 8551–8556 (1999).

    Article  ADS  PubMed  PubMed Central  CAS  Google Scholar 

57. Pontes-Quero, S. et al. Dual ifgMosaic: a versatile method for multispectral and combinatorial mosaic gene-function analysis. Cell 170, 800–814 (2017).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

58. Fre, S. et al. Notch lineages and activity in intestinal stem cells determined by a new set of knock-in mice. PLoS ONE 6, e25785 (2011).

    Article  ADS  PubMed  PubMed Central  CAS  Google Scholar 

59. Hirsch, T. et al. Regeneration of the entire human epidermis using transgenic stem cells. Nature 551, 327–332 (2017).

    Article  ADS  PubMed  PubMed Central  CAS  Google Scholar 

60. Zheng, G. X. Y. et al. Massively parallel digital transcriptional profiling of single cells. Nat. Commun. 8, 14049 (2017).

    Article  ADS  PubMed  PubMed Central  CAS  Google Scholar 

61. Hao, Y. et al. Dictionary learning for integrative, multimodal and scalable single-cell analysis. Nat. Biotechnol. 42, 293–304 (2024).

    Article  ADS  PubMed  CAS  Google Scholar 

62. McGinnis, C. S., Murrow, L. M. & Gartner, Z. J. DoubletFinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Syst. 8, 329–337 (2019).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

63. Bais, A. S. & Kostka, D. scds: computational annotation of doublets in single-cell RNA sequencing data. Bioinformatics 36, 1150–1158 (2020).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

64. van den Brink, S. C. et al. Single-cell sequencing reveals dissociation-induced gene expression in tissue subpopulations. Nat. Methods 14, 935–936 (2017).

    Article  PubMed  Google Scholar 

65. Li, B. et al. Cumulus provides cloud-based data analysis for large-scale single-cell and single-nucleus RNA-seq. Nat. Methods 17, 793–798 (2020).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

66. Bergen, V., Lange, M., Peidli, S., Wolf, F. A. & Theis, F. J. Generalizing RNA velocity to transient cell states through dynamical modeling. Nat. Biotechnol. 38, 1408–1414 (2020).

    Article  ADS  PubMed  CAS  Google Scholar 

67. Qiu, X. et al. Mapping transcriptomic vector fields of single cells. Cell 185, 690–711 (2022).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

68. Wolf, F. A., Angerer, P. & Theis, F. J. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 19, 15 (2018).

    Article  PubMed  PubMed Central  Google Scholar 

69. Aibar, S. et al. SCENIC: single-cell regulatory network inference and clustering. Nat. Methods 14, 1083–1086 (2017).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

70. Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120 (2014).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

71. Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21 (2013).

    Article  PubMed  CAS  Google Scholar 

72. Anders, S., Pyl, P. T. & Huber, W. HTSeq—a Python framework to work with high-throughput sequencing data. Bioinformatics 31, 166–169 (2015).

    Article  PubMed  CAS  Google Scholar 

73. Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).

    Article  PubMed  PubMed Central  Google Scholar 

74. Chen, Y., Chen, L., Lun, A. T. L., Baldoni, P. L. & Smyth, G. K. edgeR v4: powerful differential analysis of sequencing data with expanded functionality and improved support for small counts and larger datasets. Nucleic Acids Res. 53, gkaf018 (2025).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

75. Hochstetter, J. & Simons, B. Code accompanying “Sexual dimorphism in epidermal stem cell dynamics”. Zenodo https://doi.org/10.5281/zenodo.21997809 (2026).

Download references

Acknowledgements

We thank staff at the ULB-Erasme Animal Facility (S. Delcambre and G. Lapouge), the ULB Light Microscopy Facility (J.-M. Vanderwinden and M. Martens), the ULB-Erasme flow cytometry platform (C. Dubois) and K. Bévant for their help with single-cell experiments; staff at the Brightcore Genomics Facility (F. Libert, A. Lefort and S. Paternot) and S. Vanuytven for their help with single-cell analyses; and M. T. Bejar and M. P. Alcolea for guidance on skin epithelioids. Access to human skin was provided by the surgeons P. Tonnard and A. Verpaele and their staff.

Funding

D.E. and A.B. are supported by a FRIA grant from the F.R.S.-FNRS. J.H. acknowledges support from the Cambridge Commonwealth Trust. B.D.S. acknowledges support from the Wellcome Trust (219478/Z/19/Z and 220379/Z/20/Z), the EPSRC (EP/V062654/1) and the Royal Society (EP Abraham Research Professorship, RP/R/231004). C.B. is supported by the WEL Research Institute, FNRS, TELEVIE, the Fond Erasme, the Fondation Contre le Cancer, the ULB Foundation, the Fondation Julie et Francoise Drion, Fonds Yvonne Boel, EOS FNRS/FWO and the European Research Council.

Author information

Author notes

1. These authors contributed equally: Dan Engelman, Joel Hochstetter, Anaïs Bauduin

Authors and Affiliations

1. Laboratory of Stem Cells and Cancer, Université Libre de Bruxelles (ULB), Brussels, Belgium

   Dan Engelman, Anaïs Bauduin, Yura Song, Nordin Bansaccal, Anne-Lise Delaunois, Mélanie Liagre & Cédric Blanpain

2. Department of Applied Mathematics and Theoretical Physics, Centre for Mathematical Sciences, University of Cambridge, Cambridge, UK

   Joel Hochstetter, Qiuyu Lian & Benjamin D. Simons

3. Gurdon Institute, University of Cambridge, Cambridge, UK

   Joel Hochstetter, Qiuyu Lian & Benjamin D. Simons

4. Eli and Edythe Broad CIRM Center for Regenerative Medicine, Los Angeles, CA, USA

   Sagar Patil & Yulia Shwartz

5. Centre for Regenerative Medicine “Stefano Ferrari”, Department of Life Sciences, University of Modena and Reggio Emilia, Modena, Italy

   Marta Valenti, Michele De Luca & Laura De Rosa

6. VIB Center for Inflammation Research, Ghent, Belgium

   Annagiada Toniolo, Mohammad Ghiasloo & Esther Hoste

7. Department of Biomedical Molecular Biology, Ghent University, Ghent, Belgium

   Annagiada Toniolo & Esther Hoste

8. Department of Human Structure and Repair, Ghent University, Ghent, Belgium

   Mohammad Ghiasloo

9. WEL Research Institute, Université Libre de Bruxelles (ULB), Brussels, Belgium

   Cédric Blanpain

Authors

1. Dan Engelman
2. Joel Hochstetter
3. Anaïs Bauduin
4. Sagar Patil
5. Marta Valenti
6. Yura Song
7. Qiuyu Lian
8. Annagiada Toniolo
9. Mohammad Ghiasloo
10. Nordin Bansaccal
11. Anne-Lise Delaunois
12. Mélanie Liagre
13. Michele De Luca
14. Esther Hoste
15. Laura De Rosa
16. Yulia Shwartz
17. Benjamin D. Simons
18. Cédric Blanpain

Contributions

D.E., J.H., B.D.S. and C.B. designed the experiments, performed data analyses and wrote the manuscript. D.E., J.H., A.B., B.D.S. and C.B. discussed the data and their interpretation. D.E. and A.B. performed most of the experiments. J.H. performed mathematical modelling and statistical analyses of the clonal data, under the supervision of B.D.S. S.P. performed the FCG experiments under the supervision of Y. Shwartz. M.V. performed the 3DSE experiments under the supervision of L.D.R. and M.D.L. Y. Song and Q.L. performed RNA-seq and scRNA-seq data analyses. A.T. performed the epithelioid experiments under the supervision and experimental guidance of E.H., with samples provided by M.G. who contributed to initial processing and culture maintenance. A.-L.D. and M.L. provided technical support. N.B. performed initial experiments. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Benjamin D. Simons or Cédric Blanpain.

Ethics declarations

Competing interests

The authors declare no competing interests.

Peer review

Peer review information

Nature thanks Salvador Benitah, G. Paolo Dotto, Qing Nie and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

Additional information

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Extended data figures and tables

Extended Data Fig. 1 Intravital imaging quantification of clone size.

a-c. Illustrative images of the method used to decipher the 3D composition of clones. Basal cells are touching the second harmonic generation (SHG), which indicates that the ECM is in contact with the basal layer and suprabasal cells are residing over the basal layer and presenting enlarged nuclei. Scale bar 10 µm.

Extended Data Fig. 2 Sex-specific differences in epidermal cell fate.

a-c, Average basal (black) and suprabasal (green) clone size (a) (n = 1171 clones), clone persistence (proportion of labelled clones at D7 maintaining at least one labelled cell; n = 1171 clones) (b) and labelled cell density (n = 1171 clones) (c) over time. The basal and suprabasal labelled cell densities correspond to the average number of basal and suprabasal cells per clone respectively including extinct clones. Points represent mean and SEM of experimental data, lines show the average results from 5 independent numerical simulations of the two-proliferative compartment model at the best fits parameters obtained for a tissue of approximately 1000 cells with the shading representing the SEM (for details of the model and simulations, see main text and Supplementary Note). d, Representative time-lapse live-imaging of the same clone over time in the IFE showing an enrichment in even numbers of basal cells. Dashed line represents hair follicles. Grey is SHG. Scale bar, 20 µm. e, Proportion of clones with given basal (dark) and suprabasal (light) cell numbers in male (red) and female (blue) mice presented as mean and SEM. f, Joint clone size distribution, representing 3D clone composition, at D42 in female and male mice. Panels (a-f) are based on n = 753 clones from N = 5 female mice and n = 418 clones from N = 5 male mice. g, Representative image of the basal layer (epidermis whole-mount) of the IFE labelled with 5-ethynyl-2-deoxyuridine (EdU) (yellow) in a pulse-chase experiment showing the progeny of a single cell (pseudoclone) at 144 h (D6) post-injection. Scale bar, 10 µm. h, Proportion of pseudoclones with different basal sizes in male and female mice at 144 h (D6) post-EdU injection. n = 100 pseudoclones per mouse from N = 3 female and N = 3 male mice, points represent means with SEM, grey lines indicate even values. i,k Representative clone in the IFE of the back skin (i) and the tail skin (k) of a female mice at D28 post tamoxifen showing an enrichment of even numbers of basal cells. Scale bar, 5 µm. j,l, Proportion of clones with given basal cell numbers in back skin (j) and tail skin (l) at D28. m,n, Sex-specific basal clone size in back skin (m) and tail skin (n) at D28, two-sided Wilcoxon rank-sum test. j,l,m,n, n = 179(tail) and 150(back) clones for female, n = 251(tail) and 160(back) clones for male from N = 3 female and 3 male mice.

Source data

Extended Data Fig. 3 Single-proliferative compartment model cannot capture clonal dynamics.

a, Schematic representation of the one-proliferative compartment model. Stem cells (red) can undergo symmetric duplication, asymmetric division (probability 1-r) or terminal division giving rise to basal differentiated cells (white), which in turn delaminate at a fixed rate. Symmetric stem cell fates occur with probability r with the balance between symmetric duplication and terminal division dependent on local cell density as indexed through the cell area (for details, see Supplementary Note). Differentiated cells (white) are lost from the basal layer through cell delamination. b-i, Simulations of clonal trajectories (b,f), average basal clone size (c,g: red and blue) and suprabasal clone size (c,g: pink and light blue) size, labelled cell density (d,h), and clone persistence (e,i), separated into SC-like clones (clones growing to 2 or more cells over time; shown as red for males and blue for females) and D-like clones (clones with maximum of one cell in total; shown as black). j-k, Basal (j) and suprabasal (k) clone size distribution in female (blue) and male (red) mice from data shown in Extended Data Fig. 2e. b-k, Points represent mean and SEM, lines show the average results from 10 independent numerical simulations of the one-proliferative compartment model at the best fits parameters obtained for a tissue of approximately 400 cells with the shading representing the SEM (for details of the model and simulations, see main text and Supplementary Note). l, Clone transition probabilities between configurations starting from 2, 3 or 4 basal cells to the next available time point in male and female mice, comparing results from the time-lapse imaging and modelling using the one-proliferative compartment model. Note that the one-proliferative compartment model is unable to capture the fine-grained features present in the clonal data for female mice.

Extended Data Fig. 4 Additional comparisons with hierarchical stem-TA model.

a-b, Predictions of the two-proliferative compartment model for clone persistence, breaking down clones into SC-like (dark; clones growing to a total of 4 cells or more over the time course), TA-like (light; clones growing to a maximum of 2 or 3 cells in total over the time course) and D-like (black; clones of one total cell maximum) clones for male (a) and female (b) mice. c-d, Model fits compared to the basal (c) and suprabasal (d) clone size distribution in female (blue) and male (red) mice for data shown in Extended Data Fig. 2e. e, Suprabasal labeled cell density of SC-like clones (dark) and TA-like (light) and D-like (black) clones in male (red) and female (blue) mice. a-e, Points represent mean and SEM, lines show the average results from 5 independent numerical simulations of the two-proliferative compartment model at the best fits parameters obtained for a tissue of approximately 1000 cells with the shading representing the SEM (for details of the model and fitting procedure, see main text and Supplementary Note). f,g, Snapshots showing typical clone dynamics of the basal layer of the model in male (f) and female (g) mice across time (see also Supplementary Videos 1 and 2). Red circles indicate SC identity, pink circles indicate TA cell identity, and unfilled circles indicate differentiated cells. Male IFE comprises a higher fraction of TA cells and a lower fraction of differentiated cells. Note that female clones are marginally bigger than male and comprised predominantly of even sizes, while male mice have a larger proportion of odd sizes. D42 is also presented of Fig. 1n. h, Clone transition probabilities between states starting from 2, 3 or 4 basal cells to the next available time point in male and female mice, comparing results from the time-lapse imaging and modelling using the two-proliferative compartment model. F: from 4b is also presented on Fig. 1m. i, Best fit parameters for the two-proliferative compartment model, with 95% uncertainty intervals. a-e,h,i, For experimental data n = 753 clones for F, n = 418 clones for M from N = 5 female and N = 5 male mice.

Extended Data Fig. 5 Short-term time-lapse live-imaging supports temporal coordination among divisions and delamination of sister cells.

a, Representative time-lapse live-imaging of the same clone over time in the IFE, on the basal and suprabasal layers, every 12 h for 72 h starting at D7 post tamoxifen, highlighting coordination in terms of divisions and delaminations. Scale bar, 5 µm. Grey is SHG. b, Schematic representation of cell divisions and putative fate outcomes based on the two-proliferative compartment model of the clone presented in a. c, Probabilities of cell divisions and delaminations starting from 2, 3 or 4 basal cells clones to the next available time point in male and female mice. Highlighting the increased probability of these events appearing in pairs when starting from even numbers of basal cells, or alone when starting from an odd number of basal cells. On n = 69 clones (188 pairs of unique events where both sisters’ fates were observed) in female mice and n = 60 clones (151 pairs of unique events) in male mice from N = 4 female and 4 male mice. d, Schematic representation of temporal coordination in basal cells division and basal cells delamination on computationally generated clonal dynamics. Simulations of the clone size distribution pattern using the absence of coordination between sister events; coordination coefficient of 1 for divisions; coordination coefficient of 1 for delaminations and coordination coefficient of 1 for both divisions and delaminations. Schematic in d created in BioRender; Blanpain, C. https://biorender.com/vjxkg4a (2026).

Source data

Extended Data Fig. 6 Additional lineage tracing on male and female sex hormones.

a-d,f-i, Average basal (magenta or yellow) and suprabasal (black or green) SC-like clone size (a,f), labelled cell density of SC-like (green or yellow), TA-Like (red or green) and D-like (black) basal (b,g) and suprabasal (c,h) clones as well as clone persistence (d,i) in castrated males (a-d) and ovariectomized female mice (f-i). e,j, Suprabasal clone size distribution of the ear skin IFE in castrated males (e) and in ovariectomized females (j). k, Best fit parameters for castrated males and ovariectomized females using the two-proliferative compartment model (for details see Supplementary Note). a-j, Points show mean values with SEM. Lines show average results of 5 independent simulations of the two-proliferative compartment model using best fit parameters based on a tissue of approximately 1000 cells with the shading representing the SEM (for details, see Supplementary Note). a-k, n = 326 clones from N = 3 male mice and n = 439 clones from N = 5 female mice.

Extended Data Fig. 7 Single-cell RNA-seq of IFE basal cells.

a,b, Uniform manifold approximation and projection (UMAP) of the unsupervised clustering analysis of integrated single-cell RNA-seq of ear skin IFE basal keratinocytes isolated from 2-month-old female and male mice (N = 3 mice per sample, 2 independent samples per sex, total of 3711 integrated cells). Gene expression is visualized as a colour gradient (Yellow: high, Purple: Low). Tissue dissociation-induced cell state, clustered away from IFE basal cells, was manually removed based on Egr1, Atf3, Fosb and Jun expression⁶⁴. Colours (left panel) represent cluster identities after dissociation-induced cluster removal. b,c, Representative proportion of all integrated IFE cells in male and female datasets based on the 2D (b) or 3D (c) plotting and clustering strategy described in Fig. 4 and Methods. d, Cluster identities based on gene expression (e-j). Diff, differentiated. e-j, UMAP plots showing the normalized level of gene expression for cluster specific genes (d). Gene expression is visualized as a colour gradient (Black: high, Yellow: Low). k, Visualization of cells on a 3D manifold based on diffusion map (see Methods). Projection of the clusters obtained in classical UMAP (d) on the 3D diffusion map. SC cluster colocalized with G0/G1 A, TA cluster colocalized with G0/G1 B and C and Diff cluster colocalize with committed and Diff. l-m, Illustrative differentially expressed genes between males and females shared among all clusters. n-p, Venn diagram of pseudobulk differential expression (DE) for each cell type using Wilcoxon test, DESeq2 and edgeR using the same thresholds (log2FC > 0.5 and adjusted pval <0.01). a-p, N = 3 mice per sample, 2 independent samples per sex, F1: 1144 cells, F2: 1063 cells, M1: 585 cells M2: 919 cells.

Extended Data Fig. 8 Single-cell RNA-seq analysis of cellular trajectories.

a-c, RNA velocity analysis (see Methods) of integrated single-cell RNA-seq of ear skin IFE basal keratinocytes isolated from 2-month-old females and males indicating a developmental progression (a) from the SC cluster to the TA cell cluster, followed by cell cycle exit and commitment to differentiation (arrows). Graphical representation of transcriptomic speed (b) and acceleration (c) in each cluster. Box plots: median (line), 25th–75th percentiles (box), min/max (whiskers). d-f, Differential regulon activity inferred using SCENIC, specific for SC (d), TA (e) and Diff (f) vs. other clusters using a two-sided t-test. g-i, Differentially expressed genes in SC (g), TA (h), differentiated clusters (i) between males and females Volcano plot of mRNA expression using significance margins of adjusted p-value < 0.005, using two-sided Wilcoxon rank sum test. j-l, Differential regulon activity inferred using SCENIC between males and females in SC (j), TA (k) and Diff (l) using a two-sided t-test. a-l, Total of 3711 integrated cells from N = 3 mice per sample, 2 independent samples per sex (n = 1504 cells for males and n = 2207 cells for females).

Extended Data Fig. 9 Lineage tracing of Lgr6CreER- and Notch1CreER-targeted basal cells.

a, Genetic strategy used to target Lgr6 and Notch1-expressing WT IFE basal cells with an inducible multicolor nuclear reporter that randomly labels nuclei with either H2B-mCherry, H2B-eGFP or H2B-cerulean fluorescent protein (Lgr6CreER/iChr-Control-Mosaic and Notch1CreER/iChr-Control-Mosaic). b, Protocol used for long-term time-lapse live imaging. c-d, Representative images of time-lapse analysis of the same Notch1CreER (c) and Lgr6CreER (d) traced clones over time. Dashed line around hair follicles. Scale bar, 20 µm. e-j, Sex-specific average basal (plain) and suprabasal (dashed) clone size (e-f), clone persistence (g-h) and labelled cell density (i-j) over time in female (blue) and male (red) mice in Lgr6CreER/iChr-Control-Mosaic (e,g,i) and Notch1CreER/iChr-Control-Mosaic (f,h,j) mice. For experimental data n = 114(F) and 152(M) clones for Lgr6CreER, n = 562(F) and 515(M) clones for Notch1CreER from N = 3 female and N = 3 male mice of each strain. Points represent mean values with SEM. k-n, Representative images (k) of the basal and suprabasal layers of the IFE and quantification (m) of Krt10 positive basal labelled cells or fraction of suprabasal labelled cells (n) using Lgr6CreER (k) and Notch1CreER (l) at D3 post-induction with tamoxifen in male and female mice. Scale bar = 10 µm. e-n, n = 200(F) and 177(M) clones for K14CreER, n = 85(F) and 86(M) clones for Lgr6CreER, n = 300(F) and 316(M) clones for Notch1CreER from N = 2 female and 2 male mice for K14CreER, N = 2 female and 2 male mice for Lgr6CreER and N = 3 F and 3 M mice for Notch1CreER.

Source data

Extended Data Fig. 10 Lgr6CreER and Notch1CreER with hierarchical stem-TA model.

a-l, Average basal (dark red and dark blue) and suprabasal (pink and light blue) clone size (a-d) and labelled cell density (e-h) and clone persistence (i-l) of SC-like clones (red and blue; clones growing to a total of 4 cells or more over the time course), TA-like (light red and light blue; clones growing to a maximum of 2 or 3 cells in total over the time course) and D-like (black; clones of one total cell maximum) clones in male (blue) and female (red) mice. m, Basal (left side) and suprabasal (right side) clone size distribution in female (blue, yellow/green) and male (red, purple) mice from Lgr6CreER and Notch1CreER mice. n, Best fit parameters for Lgr6CreER/iChr-Control-Mosaic and Notch1CreER/iChr-Control-Mosaic, male and female mice using the two-proliferative compartment model (for details on fitting procedure see Supplementary Note). a-m, Points represent mean and SEM of experimental data, lines show the average results from 5 independent numerical simulations of the two-proliferative compartment model at the best fits parameters obtained for a tissue of approximately 1000 cells with the shading representing the SEM (for details of the model and simulations, see main text and Supplementary Note). For experimental data n = 114(F) and 152(M) clones for Lgr6CreER, n = 562(F) and 515(M) clones for Notch1CreER from N = 3 female and N = 3 male mice of each strain. Points represent mean values with SEM mice.

Extended Data Fig. 11 Impact of sex hormones on dermal thickness and hair follicle cycle.

a-b, Representative images (a) and quantification (b) of dermal thickness in control males, castrated males, castrated males with dihydrotestosterone (DHT) implants, control females, ovariectomized females and ovariectomized females with estradiol (E2) implants. Mean of 5 measures per mouse from N = 3 mice per group. Scale bar, 20 µm, two-sided t-test c-d, Representative images (c) of hair follicles in back skin and ear skin and quantification (d) of hair follicle cycle stage in ear skin in control males and castrated males. Mean of 15 follicles per mouse from N = 3 mice per group. Scale bar, 100 µm.

Source data

Extended Data Fig. 12 Additional analysis of bulk RNA-seq of basal cells from female, male, and castrated male mice.

a-d, Relative mRNA expression of upregulated genes in females (a) and males (b) ear skin IFE basal cells isolated from 2-month-old females, males and castrated males, highlighting X- and Y-linked genes. (N = 2 mice per sample, 2 samples per condition.) c-d, Venn diagram of genes upregulated in females versus males and castrated males versus males (c) or upregulated in males versus females and upregulated in males versus castrated males (d). Genes were selected based on a FC > 2 or <0.5 (Methods). e-m, Venn diagram of differential expression (DE) of upregulated genes in males (e), females (f) and castrated males (k) using two-sided Wilcoxon test, DESeq2 (g,i,l) and edgeR (h,j,m) with the same thresholds (log2FC > 0.5 and adjusted pval <0.01). Limited to genes presented on Fig. 5 (g-j,l,m). For DESeq2, adjusted P values were calculated using two-sided Wald tests followed by Benjamini–Hochberg correction for multiple comparisons. For edgeR, differential expression was assessed using edgeR’s quasi-likelihood F-test based on a negative-binomial generalized linear model. The test was non-directional (two-sided with respect to log-fold change), and P values were adjusted for multiple comparisons using the Benjamini–Hochberg method to control the FDR. Statistically significant values are displayed in red. n-o, Schematic representation of 3D skin equivalents (3DSE) (n) and epithelioids (o) cultures, highlighting the presence of a fibroblast (fibro) gel in 3DSE and a higher level of stratification compared to epithelioids, both in air liquid interface (ALI). Schematics in BioRender: n, Blanpain, C. https://biorender.com/o4xdpm4 (2026); o, Engelman, D. https://biorender.com/ebs18vj (2026).

Supplementary information

Supplementary Information (download PDF )

Supplementary Information containing Supplementary Note (including Supplementary Figs. 1–4, Supplementary Tables 1–3 and Supplementary References) and Supplementary Figs. 5–7.

Reporting Summary (download PDF )

Supplementary Data 1 (download XLSX )

DEGs using edgeR in males and females from scRNA-seq data.

Supplementary Data 2 (download XLSX )

DEGs using DESeq2 in males and females from scRNA-seq data.

Supplementary Data 3 (download XLSX )

SCENIC AUCell scores for all regulons across all cells, differential regulon activity between cell types and sex-biased differences in regulon activity in each cell type from scRNA-seq data.

Peer Review File (download PDF )

Supplementary Video 1 (download MP4 )

Simulation of the clonal dynamics of the two-proliferative compartment STD model for male mice. SCs are marked red, TA cells are marked pink and differentiated cells with no marking. The coloured cells show a random selection of individual SC-derived clones. The simulations are shown based on the best-fit parameters of the model to the clonal data (see main text and Supplementary Note). The simulation shows a small degree of clone fragmentation, which occurs in the framework of the Voronoi model because of the rapid delamination time and local synchrony of events, which may lead to T1-type transitions following cell delamination. Note that clones show a small bias towards even sizes.

Supplementary Video 2 (download MP4 )

Simulation of the clonal dynamics of the two-proliferative compartment STD model for female mice. SCs are marked red, TA cells are marked pink and differentiated cells with no marking. The coloured cells show a random selection of individual SC-derived clones. The simulations are shown based on the best-fit parameters of the model to the clonal data (see main text and Supplementary Note). Note that with the majority of SC divisions resulting in asymmetric fate and the majority of TA cell divisions resulting in terminal differentiation, the clones show a strong bias towards even cell numbers.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions

About this article

Cite this article

Engelman, D., Hochstetter, J., Bauduin, A. et al. Sexual dimorphism in epidermal stem cell dynamics. Nature (2026). https://doi.org/10.1038/s41586-026-11159-2

Download citation

• Received: 14 April 2025

• Accepted: 11 September 2026

• Published: 14 October 2026

• Version of record: 14 October 2026

• DOI: https://doi.org/10.1038/s41586-026-11159-2