Mouse
Krt14creER transgenic mice were provided by E. Fuchs56. iChr-control-mosaic transgenic mice were provided by R. Benedito57. Notch1creER (N1-CreERT2SAT) mice were provided by S. Fre58. Lgr6-GFP-creER mice were obtained from The Jackson Laboratory39. Tissues from FCG mice were provided by the Shwartz Laboratory; the original FCG mice were obtained from A. Arnold under a material transfer agreement38. 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 described16. 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% H2O2 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 × 104 per cm2) on lethally irradiated 3T3-J2 cells (2.4 × 104 per cm2) and grown at 37 °C, 5% CO2 in a humidified atmosphere in keratinocyte growth medium (KGM)59.
Primary human keratinocytes (3.5 × 105 cells per scaffold) from male and female donors were seeded on 3T3-J2 cells (mouse fibroblast, 2 × 105 cells per scaffold) mixed with collagen gels (final concentration of 4 mg ml–1) in hormone-free keratinocyte growth medium (hfKGM) as previously described59. 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−4 M adenine (Sigma-Aldrich, A2786), 1 × 10−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% CO2. 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)60 with default parameters.
Quality control and downstream analyses were performed using the R package Seurat (v.5.1.0)61. 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)62 and scds (v.1.26.1)63 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 framework16,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 data9,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 log2 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 analyses64 (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)65. 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)66 was used to run RNA velocity analysis, which was visualized on a force-directed layout manifold. The Python package dynamo (v.1.4.1)67 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)68, 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 log2 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 (log2 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)69 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)70. Trimmed reads were mapped against the mouse reference genome (Grcm38/mm10) using STAR aligner (v.2.7.11b)71 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)72 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)73. 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, log2 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)74 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 Fe(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.
