Lineage commitment of dermal fibroblast progenitors is controlled by Kdm6b‐mediated chromatin demethylation

Abstract Dermal Fibroblast Progenitors (DFPs) differentiate into distinct fibroblast lineages during skin development. However, the epigenetic mechanisms that regulate DFP differentiation are not known. Our objective was to use multimodal single‐cell approaches, epigenetic assays, and allografting techniques to define a DFP state and the mechanism that governs its differentiation potential. Our initial results indicated that the overall transcription profile of DFPs is repressed by H3K27me3 and has inaccessible chromatin at lineage‐specific genes. Surprisingly, the repressive chromatin profile of DFPs renders them unable to reform the skin in allograft assays despite their multipotent potential. We hypothesized that chromatin derepression was modulated by the H3K27me3 demethylase, Kdm6b/Jmjd3. Dermal fibroblast–specific deletion of Kdm6b/Jmjd3 in mice resulted in adipocyte compartment ablation and inhibition of mature dermal papilla functions, confirmed by additional single‐cell RNA‐seq, ChIP‐seq, and allografting assays. We conclude that DFPs are functionally derepressed during murine skin development by Kdm6b/Jmjd3. Our studies therefore reveal a multimodal understanding of how DFPs differentiate into distinct fibroblast lineages and provide a novel publicly available multiomics search tool.


Introduction
Fibroblasts are mesenchymal cells that create and maintain the extracellular matrix of the connective tissues supporting the functions of various organs (Watt & Fujiwara, 2011;Driskell & Watt, 2015;Plikus et al, 2021).They are distinct across different organs but are also diverse within the same tissue (Dhouailly & Sengel, 1975;Fernandes et al, 2004;Muhl et al, 2020;Buechler et al, 2021).Mammalian skin serves as a critical model to study development, homeostasis, and the role of fibroblast heterogeneity in tissue-resident stem cell niche, tissue regeneration, fibrosisassociated diseases, and the effects of chronic aging (Singer & Clark, 1999;Gurtner et al, 2008;Chen et al, 2013;Cole et al, 2018;Griffin et al, 2020;Gruber et al, 2020).Recent studies have demonstrated the distinct roles of dermal fibroblast subpopulations during embryonic development and skin wound healing.These studies have revealed that manipulations of specific fibroblast subpopulations lead to hair follicles and adipocytes regeneration in different wounding scenarios (Rinkevich et al, 2015;Plikus et al, 2017;Jiang et al, 2018;Correa-Gallegos et al, 2019;Guerrero-Juarez et al, 2019;Abbasi et al, 2020;Phan et al, 2020;Shook et al, 2020;Mascharak et al, 2021;Sinha et al, 2022).However, the establishment of fibroblast heterogeneity from a common progenitor during embryonic development remains poorly characterized.
During early embryonic development, mesenchymal stem cells migrate from different embryonic origins to subepidermal regions of the skin, where the dorsal dermis arises from the somite, and the facial dermal fibroblasts derive from the neural crest (Dhouailly, 1973;Fernandes et al, 2004;Atit et al, 2006;Wong et al, 2006;Tran et al, 2010;Thulabandu et al, 2018).Approximately E13.5-14.5 in murine skin, reciprocal signaling between the overlying epidermis and the Dermal Fibroblast Progenitors (DFPs) leads to condensation of upper fibroblasts to form the placodes (Kishimoto et al, 2000;Millar, 2002;Chen et al, 2012;Fu & Hsu, 2013).Between E14.5 and E16.5, the DFPs then differentiate and commit to distinct lineages of upper papillary and lower reticular fibroblasts (Driskell et al, 2013).Upper papillary fibroblasts differentiate into the dermal papilla, neonatal papillary fibroblast, and the arrector pili muscles postnatal.Lower reticular fibroblasts differentiate into collagen-secreting reticular and the dermal white adipocyte (DWAT; Festa et al, 2011;Wojciechowicz et al, 2013).Moreover, classic chimeric skin grafting studies between avian and murine skin have shown that DFPs need to differentiate before responding to epidermal signals, and that fibroblast identity is also region-specific (Dhouailly, 1973).Altogether, during embryonic development, we propose that DFPs must undergo an intrinsic reprogramming process to commit to fibroblast lineages and differentiate toward the terminal fates.
Epigenetic regulators have been extensively studied in the context of cellular differentiation (Boyer et al, 2006;Lee et al, 2006;Kang et al, 2019;Flora & Ezhkova, 2020;Plikus et al, 2021).A proposed mechanism of chromatin regulation is through the addition of methyl groups to Histone 3 Lysine 27 (H3K27), which is associated with condensed and silenced chromatin while the removal of methyl group is associated with accessible and active chromatin (Cao et al, 2002;Agger et al, 2007;Xiang et al, 2007;Kang et al, 2019).The most widely known regulators of methylation on H3K27 are part of the Polycomb repressive complex 2 (PRC2), EZH1 and EZH2.In the epidermis, the addition and removal of the trimethyl group (me3) on H3K27me3 were shown to be critical in controlling the differentiation process of epidermal stem cells during development (Sen et al, 2008;Ezhkova et al, 2009).Demethylases, such as KDM6A/UTX and KDM6B/JMJD3, specifically remove methyl group and derepress chromatin (Agger et al, 2007;Xiang et al, 2007).In dermal fibroblasts, EZH2 has been shown to act as an epigenetic rheostat to control the Wnt and Retinoic Acid signaling in DFPs (Thulabandu et al, 2021).However, a full understanding of how the methylation status of H3K27 regulates fibroblast differentiation has not been shown.In addition, other studies have also investigated the roles of different epigenetic marks in regulating and maintaining fibroblast lineages in developed the skin.Acetylation of H3K27 and HDAC2 activities are required for the maintenance of Papillary Fibroblast lineages postnatally (Kim et al, 2022).H3K4me3, H3K9me3, and H3K27me3 all display reduced methylation levels in hair follicle stem cells prior to hair growth, and inhibition of hypomethylation delays skin wound healing (Lee et al, 2016;Kang et al, 2020).Nonetheless, an epigenetic regulation that programs dermal fibroblast lineage commitment and differentiation during embryonic development remains unknown.
The recent emergence of single-cell technology has advanced the study of fibroblast heterogeneity drastically (Plikus et al, 2017).The transcriptomic changes in restricted fibroblast lineages, specifically in the hair follicle-associated fibroblasts, have been characterized at different timepoints that regulate Dermal Papilla establishment (Biggs et al, 2018;Gupta et al, 2019;Mok et al, 2019).In our recent study of transcriptomic and chromatin landscapes at the single-cell level, we showed that the transcriptomic profiles of fibroblast lineages could only reflect the cellular states, while the accessible chromatin profiles can predict the lineage restriction and terminal fate commitment (Thompson et al, 2022).Overall, we hypothesize that the differentiation of dermal fibroblast progenitors could be defined by incorporating chromatin accessibility profiles and epigenetic regulation.
To establish the multimodal reference point for dermal fibroblast progenitors, we compared the molecular states of dermal fibroblast progenitors (DFP) at E14.5 to late embryonic fibroblasts at E18.5 using single-cell transcriptomics, single-cell ATAC, ChIP-seq, genetic ablation of epigenetic regulator Kdm6b/Jmjd3, and functional chamber allografting assays.We found that DFPs at E14.5 are predicted to differentiate into hair follicle-associated fibroblast, but they are not intrinsically programmed as they failed to support hair follicle formation ex vivo.E14.5 DFPs possess a closed-off chromatin profile, marked by a high level of H3K27me3, and they require epigenetic programming mediated by Kdm6b/Jmjd3 to commit to specific fibroblast lineages.Finally, we share our processed data for all our experiments on our website skinregeneration.org.

Results
DFPs differentiate into distinct fibroblast lineages in vivo but fail to regenerate functional skin in grafting assays Using lineage tracing assays, we have previously shown that DFPs are present throughout the dermis at E14.5 and that a DFP has the potential to contribute to the formation of all layers and fibroblast subtypes (Driskell et al, 2013;Rognoni et al, 2016).We hypothesized that E14.5 fibroblasts are multipotent DFP progenitors in developing murine skin.We performed scRNA-seq from whole skin samples of E14.5, E18.5, and P5 mice (Fig 1A -E).Altogether, we captured and sequenced a total of 43,740 cells, of which 17,356 cells were classified as dermal fibroblasts (Fig EV1A and B).To evaluate DFP differentiation and developmental potential, we integrated all three time points together and embedded the single-cell clusters as a PAGA map to improve interpretation but also to preserve the global topology of cells transition (Fig 1D -E; Materials and Methods) (Hie et al, 2019;Wolf et al, 2019).DFPs at E14.5 were classified into two distinct subsets based on their horizontal anatomical positions as previously published (Gupta et al, 2019;preprint: Jacob et al, 2022).Upper DFPs (clusters 0, 2, 11) expressing Crabp1, Cav1, Nkd1, and Lef1 resided closer to the epidermis, while lower DFPs expressed Thbs1, Ptn, and Mfap5 (clusters 1, 5, 6) and EV1A and B).Furthermore, fibroblast subsets in E18.5 and P5 populations were identified based on gene expression profiles that have been previously defined (Fig EV1A and B) (Guerrero-Juarez et al, 2019;Phan et al, 2020;Thompson et al, 2022).We identified clusters of Dermal Papilla, Papillary Fibroblasts, Reticular Fibroblasts, Fascia, and Preadipocytes in both E18.5 and P5 skin.In addition, clusters of mature papillary and ECM-producing fibroblasts were detected only in P5 skin, reflecting the maturation of postnatal dermal fibroblasts.Our PAGA integrated analysis revealed that most E18.5 fibroblasts were closely associated with E14.5 DFPs while P5 fibroblasts clustered distinctly based on their lineage functions, except in the lineage-committed specialized fibroblasts such as Dermal Papilla and Preadipocytes, where E18.5 and P5 clustered together (Fig 1D -E).
To further establish the continuous transition between DFPs and fibroblast lineages, we performed RNA velocity analysis on the PAGA-embedded fibroblast subsets, which accounted for the transcriptional dynamics of splicing kinetics within each cell (Bergen et al, 2020).Upper DFPs are shown to differentiate into papillary fibroblasts (clusters 7, 8, 16, 21), dermal papillae (clusters 4, 12, 20), dermal sheath (clusters 10), and mature fibroblasts (clusters 3, 17), while lower DFPs are projected to become upper DFPs and reticular (clusters 9,13,15,18; (Driskell et al, 2014;Correa-Gallegos et al, 2019;Guerrero-Juarez et al, 2019;Joost et al, 2020;Shook et al, 2020).Unexpectedly, fascia and preadipocytes showed inverted RNA velocity trajectories toward replicating reticular fibroblasts.This observation could be due to the multiple rate kinetics genes (MURK) (Barile et al, 2021), which are associated with assumptions and caveats for velocity analysis in scRNA-seq.For example, dynamical MURK genes can violate the assumptions that velocity analysis is based on (Barile et al, 2021).In addition, it has been proposed that current batch correction implementations may not be accurate (Bergen et al, 2021).Consequently, it is important to test in-silico transcriptome-based results through robust independent systems.Moreover, we have noticed that a large portion of fascia and preadipocytes is actively replicating and clustered closely with replicating lower reticular (cluster 15).The newly generated unspliced mRNA is higher in the replicating population, contributing to the inversion of RNA velocity trajectory.RNA velocity analysis also provided us with a list of putative driver genes that need to be activated for cells to differentiate toward a specific cluster.Using these genes, we generated a heatmap depicting the pseudotime transitions of upper and lower DFPs into two distinct lineages: Dermal Papilla and Preadipocytes (Fig EV1C).We found that multiple velocity driver genes of E18.5 and P5 Dermal Papilla were expressed as early as E14.5 Dermal Condensate (Fig EV1C).Overall, our results suggested that at E14.5, DPFs have the potential to differentiate into specialized cell types to support and maintain the different dermal niches in postnatal skin.
Grafting assays, in which cells are transplanted into host tissues and allowed to interact with the surrounding microenvironment, provide a powerful tool to assess the differentiation potential and lineage restrictions of stem cells ex vivo (Jensen et al, 2010;Driskell et al, 2013).To functionally challenge the potential of E14.5 DFPs to support the formation of fully functional skin, we performed a skin reconstitution experiment utilizing the chamber grafting technique.We collected whole, unsorted DFPs from E14.5 skin, and compared their ability to reform hair follicles and skin from fibroblasts from E18.5 and P5 ROSA-CMV-tdTomato-positive mice to trace their contribution within grafted tissue.A combination of dermal cells and epidermal cells was injected into a silicone chamber grafted on top of an immunocompromised Foxn1 À/À mouse (Fig 1G).Skin from areas of the graft was collected at 28 days postsurgery to assess the reformation and development of the skin.The analysis and quantification of the grafting assay revealed that E18.5 and P5 differentiated fibroblast populations were able to support the reformation of fully functional hair follicles and hair fibers while E14.5 DFPs failed to regenerate and develop hair follicles despite their contribution to skin reformation (Fig 1H -J).This result suggests that, despite the potential trajectory and presence of early dermal papilla cells, DFPs from E14.5 skin are not intrinsically programmed to differentiate ex vivo to support fully functional skin reformation.
DFP programming occurs by specifically increasing the accessibility of chromatin architecture in distinct lineages The chromatin accessibility landscape regulates cellular gene expression and is highly informative of its states and fates of cells (Buenrostro et al, 2013(Buenrostro et al, , 2015;;Cusanovich et al, 2015Cusanovich et al, , 2018;;Trapnell, 2015;Cao et al, 2017;Thompson et al, 2022).We hypothesized that the intrinsic differentiation program within E14.5 DFPs could be reflected in its chromatin accessibility profiles.We performed single-cell ATACseq assays on murine DFPs at E14.5 and differentiating fibroblasts at E18.5.After quality control, filtering, and preprocessing, 10,913 fibroblasts from both time points were used for downstream analyses including dimension reduction and unsupervised clustering (4,187 from E14.5 and 6,006 from E18.To investigate the relationship between chromatin depression and the establishment of fibroblast heterogeneity from DFPs, unsupervised clustering of scATAC-seq data was performed resulting in 15 distinct clusters from both time points (Figs 2E and EV2A).Cross-modality integration and label transfer were used to integrate scATAC-seq with scRNA-seq data for cluster identification (Materials and Methods; Fig EV2A) (Stuart et al, 2021).Consistent with our transcriptomic analysis, DFPs at E14 were split into Upper DFPs and Lower DFPs, while E18 fibroblasts contained Dermal Papilla, Papillary, Reticular, and Preadipocytes clusters (Fig 2E).Interestingly, even though fascia fibroblasts were detected in our scRNA analysis, we did not detect them in our scATAC-seq data (Fig 2E).Differential peak analysis revealed distinct chromatin landscapes among fibroblasts at E18.5, while the E14.5 Upper and Lower DFPs possessed less unique accessible peaks profiles, which was consistent with the opening of chromatin during DFPs differentiation processes (Fig EV2C).Similar to P0 fibroblasts (Thompson et al, 2022), E18.5 fibroblasts also shared a gradient of significant peaks within their respective lineages, with the differentiated Dermal Papilla and Preadipocytes fates possessing the most distinct chromatin profiles (Figs 2F and G and EV2B and C).Considering the predicted RNA velocity trajectories of DFPs to differentiate into Dermal Papilla and Preadipocyte, we hypothesized that the chromatin accessibility profiles of the equivalent cell types will also reflect a similar differentiation trajectory.We constructed a pseudotime heatmap based on our scATAC-seq data but defined by RNA velocity driver genes for each

Chromatin programming of DFPs is associated by the removal of the repressive histone modification H3K27me3
Posttranslational modifications of histones can regulate gene accessibility and expression during the differentiation process in various tissue models (Boyer et al, 2006;Lee et al, 2006;Heintzman et al, 2009;Karli c et al, 2010).Specifically, the Histone 3 Lysine 27 trimethylation mark (H3K27me3), which is associated with the repression of gene transcription, plays a key role in the differentiation process of epidermal keratinocytes (Sen et al, 2008;Ezhkova et al, 2009;Kang et al, 2019).Our scATAC-seq data suggested that the opening of chromatin regions at lineage-specific genes was critical for fibroblast differentiation.To study the role of H3K27me3 in regulating open chromatin regions of DFPs during embryonic differentiation, we performed H3K27me3 ChIP-seq comparing E14.5 and E18.5 Fibroblasts (Fig 3A).In conformity with the increase in open chromatin regions between E14.5 and E18.5 fibroblasts, we found a sharp decrease of 52% in the number of H3K27me3 peaks (from 19,060 peaks at E14.5 to 8,768 peaks at E18.5, n = 2) by visualizing the heatmaps and coverage plots of peaks called at each time point (Figs 3B-D and EV2D).In addition, we also confirmed that the majority of the genes associated with H3K27me3 peaks at E18.5 were present at E14.5 (95%), indicating that there was a significant reduction in H3K27me3-associated genes and not a dynamic change where undifferentiated genes were repressed and differentiated genes were derepressed (Fig 3D).We next utilized DiffBind to identify the significant changes between E14.5 and E18.5 H3K27me3 peaks in DFs, which resulted in 4699 significant peaks that decreased from E14.5 to E18.5 and only 29 peaks that increased (Fig 3E).To investigate their functions, Gene Ontology analysis of genes associated with the significantly derepressed peaks from E14.5 revealed the biological processes associated with skin differentiation and development such as multicellular organism development, cell differentiation, neuron differentiation, Canonical Wnt signaling pathway, hair follicle development, and glucose homeostasis (Fig 3F The EMBO Journal Quan M Phan et al Overall, our data revealed that the changes in chromatin landscapes in DFPs between E14 and E18 coincided with the demethylase of H3K27me3 markers within genes required for DFP differentiation.Furthermore, we also found an active, high level of expression for the H3K27me3 specific demethylase-Kdm6b/Jmjd3-in both E14 and E18 fibroblasts (Fig EV2E ).

Genetic ablation of histone demethylase Kdm6b/Jmjd3 inhibits DFP differentiation
Previous studies have identified the roles of Ezh2 and PRC2 complex in controlling epidermal differentiation by adding the H3K27me3 marker at specific genes to maintain the stemness of the progenitor cells.In contrast, the H3K27me3 demethylase Kdm6b/Jmjd3 has also been shown to be critical in epidermal stem cell differentiation, where ablation of Jmjd3 leads to inhibited differentiation of this cell type (Sen et al, 2008;Ezhkova et al, 2009).Based on the significant decrease in H3K27me3 markers of DFPs from E14 to E18, we hypothesize that Kdm6b/Jmjd3 demethylase activity is required for the differentiation of dermal fibroblasts.To test the role of Kdm6b/ Jmdj3 in dermal differentiation, we generated a transgenic dermalspecific knockout of Jmjd3 mouse model driven by activation of Dermo1 in dermal fibroblasts (Fig 4A).The Dermo1Cre-Kdm6b fl/fl (Kdm6b-cKO) mice are neonatally lethal.Consequently, we focused our studies on the changes that occur between E14.5 and E18.5.Histological analysis of E14.5 skin showed similar stages of skin development between Kdm6b-cKO skin and the WT control where a thin layer of the dermis was found in conjunction with the formation of dermal condensates (Fig 4B).However, at E18.5, we observed drastic changes as Kdm6b-cKO skin appeared to be underdeveloped with smaller hair follicles and significantly thinner dermal layers (Fig 4C and D).We also found that the  EV3A-C).Interestingly, we did not observe a significant change in the number of Dermal Papilla between WT (Kdm6b fl/fl ) and Kdm6b-cKO.Instead, we found that there was a dramatic reduction in number of preadipocytes detected in Kdm6b-cKO skin compared to the WT control (1,081 in WT, 238 in Kdm6b-cKO), along with a large decrease in the number of lymphatic vasculatures (Figs 5C and D,and EV3A and B).This observation coincides with previous studies reporting the relationship between dermal fibroblasts, dermal preadipocytes, and endothelial vasculature niche (Franc ¸ois et al, 2008;Tang et al, 2008;Gupta et al, 2012;Jiang et al, 2014;Martinez-Corral et al, 2015;Pichol-Thievend et al, 2018;Stone & Stainier, 2019;Oliver et al, 2020).We next performed a differential gene expression analysis comparing the two conditions using Diffxpy, a more optimized package for statistical testing of single-cell     To assess the role of epigenetic regulation in Dermal Papilla functions, we hypothesized that the ablation of Kdm6b/Jmjd3 in DFPs will lead to inaccessibility of transcription factor binding sites that regulate fibroblast lineage functions.To address this hypothesis, we utilized the motif finder Homer to screen for binding motif enrichment in the significant peaks regulated by Kdm6b (Fig 6F;Heinz et al, 2010).The known Homer motifs revealed an enrichment of Ebf1, Tcf4/Tcf7l2binding motifs in peaks that were increased in Kdm6b-cKO fibroblasts, while the de novo Homer motifs detected the significant enrichment of dermal adipocyte differentiation transcription factors Med1 and Zfp423 (Fig 6F;Gupta et al, 2010Gupta et al, , 2012;;Harms et al, 2015;Shao et al, 2016;Biferali et al, 2021).Interestingly, the expression of these transcription factors was not specific to just Dermal Papilla or Preadipocyte, but these were expressed throughout fibroblasts populations (Fig 6G).This result suggested that they might have distinct roles in both Adipogenesis and Dermal Papilla functions.
Since Kdm6b-cKO mice are neonatally lethal, we performed a chamber grafting assay to evaluate the postnatal potential of the mutant dermis to forming functional dermal papilla that support hair reformation.In order to confirm the contribution of mutant dermis within the grafting assay, we generated a triple transgenic mouse line that harbors the dermal-specific Kdm6b deletion with the mTmG mouse line (Fig 7A).Wild-type cell dermal populations were tracked utilizing the mTmG mouse line.Our analysis of the chamber grafting experiment after 28 days post surgery revealed that hair follicle reformation was robust when wild-type E18.5 dermis was used to seed the dermal component of the graft (Fig 7B and D).However, when E18.5 Kdm6b-cKO dermis was used to seed the chamber grafts (Fig 7C and D) hair follicle reformation was almost ablated despite significant contribution to the graft (Fig 7E and F).Overall, we concluded that DFPs differentiation into distinct cell fates and lineages was programmed by the establishment of unique open chromatin profiles, which is mediated by Kdm6b demethylase activity removing H3K27me3 markers as specific, functional genes loci.

Defining the epigenetic state of a dermal fibroblast
Recent studies have demonstrated that stem cells require extensive chromatin accessibility priming prior to lineage commitment (Buenrostro et al, 2018;Ranzoni et al, 2021).The accessible chromatin profiles of dermal fibroblasts and other cell types could reflect their potential to differentiate into a terminal fate (Li et al, 2017;Adam et al, 2018;Ma et al, 2020;Kim et al, 2022;Thompson et al, 2022).Here, we discovered that DFPs possess inaccessible chromatin profiles at differentiation specific genes required for lineage commitment.The inaccessible and unprogrammed chromatin profiles of E14.5 DFPs compared to the differentiating fibroblasts helped explain the inability of E14.5 dermis to support HF reformation in the chamber grafting assay (Figs 1 and 2D).Moreover, our scATAC-seq data confirmed the bifurcation of distinct fibroblast lineages during embryonic development (Sorrell & Caplan, 2004;Driskell et al, 2013).We found that Upper DFPs possessed the most stem-like chromatin profiles with the least number of accessible peaks.Papillary and Dermal Papilla are closely associated with each other, where significant differential peak analysis and integration with RNA velocity driver genes showed the large overlapping regions between them (Figs 2F and EV1D).In addition, our integrative RNA velocity and scATAC-seq analysis revealed a stepwise arrangement of accessible chromatin for DFP differentiation, where Dermal Papilla possess the most unique chromatin architecture, but genes specific for Papillary fibroblasts also remain accessible in the Dermal Papilla (Fig 2F).This result confirmed our previous finding of fibroblasts states and fates, with the Dermal Papilla and Preadipocytes being two defined terminal fates for DFP differentiation (Thompson et al, 2022).
Our investigation of the regulation of chromatin architecture by epigenetic histone modifications revealed that DFP differentiation occurs through specific derepression of chromatin to support fibroblast lineage function.We found that DFPs required a significant decrease in H3K27me3 levels in genes associated with differentiating fibroblasts functions such as Sox18, Alpl, Cebpa, and Pparg.Notably, during the differentiation process, we only observed a small increase in H3K27me3 markers (24 peaks) (Fig 3B-F).Our result suggested that the opening of chromatin accessibility profiles is correlated with the demethylation of H3K27me3 from E14.5 and during fibroblast lineage commitment.This finding confirms the epigenetic regulation of stem cell differentiation where the PRC2 complex is needed to trimethylate H3K27 and repress genes required for cellular differentiation (Boyer et al, 2006;Lee et al, 2006;Ezhkova et al, 2009;Wang et al, 2010;Thulabandu et al, 2018).Our model suggests that after Ezh2 establishes an epigenetic profile of early dermal fibroblast progenitors that requires Kdm6b/Jmjd3 activity to initiate chromatin derepression at differentiation genes.

Mechanisms regulating adipogenesis and dermal papilla lineage commitment
Lineage tracing and single-cell studies have demonstrated the importance of fibroblast lineages commitment of DFPs in skin development ◀ Figure 7. Kdm6b regulates DFP differentiation into fully functional Dermal Papilla.

A
Experimental strategy for skin reconstitution assay (chamber grafting assay) to analyze the hair forming potential of E18.5 Fibroblast of WT and Kdm6ab-cKO dermis.B, C Images of hair follicle formations within grafted area comparing WT and Kdm6b-cKO Fibroblasts.D Quantification of hair follicles reformed in chamber grafting assays.E, F Immunofluorescence analysis of WT and Kdm6b-cKO grafted area to verify contribution to skin reconstitution.Source data are available online for this figure.

Ó 2023 The Authors
The EMBO Journal 42: e113880 | 2023 and homeostasis (Driskell & Watt, 2015;Griffin et al, 2020;Plikus et al, 2021).Genetic ablation of Jmdj3/Kmd6b in dermal fibroblasts revealed an inhibition of DFPs to differentiate into terminal fates such as functional Dermal Papilla and Adipocytes.Interestingly, the fibroblast states of papillary and reticular lineages appeared to be unaffected.We found that even though hair follicle morphogenesis was initiated and that Dermal Papilla were observed, Kdm6b-cKO hair follicles appeared stunted and underdeveloped (  Festa et al, 2011;Lien et al, 2014).After the integration of scATAC-seq data with our ChIPseq data, we discovered that these transcription factors might regulate different sets of downstream genes in different cell types (Fig EV5B).Overall, our data suggested that the chromatin accessibility profiles can dictate the potential of fibroblast commitment to the terminal fates, and that Kdm6b/Jmdj3 demethylation regulates this process.The differential activity of Kdm6b/Jmjd3 in diverse fibroblast lineages acts as a mechanism to regulate the dual-function of non-unique transcription factors.Finally, we have also reconstructed a high-resolution molecular blueprint of dermal fibroblast differentiation processes across multiple modalities including transcriptomic, accessible chromatin profiles, and epigenetic markers to define DFPs and distinct fibroblast fates, freely accessible at our webtool.

scATAC-seq provides insights into the regulation of progenitor cell differentiation and adult tissue homeostasis
Recent studies have characterized the single-cell chromatin profiles of fibroblasts in murine skin wound healing.In regular scarring wounds, Foster et al (2021) found different populations of fibroblasts with distinct putative roles (Foster et al, 2021).While comparing large WIHN and small scarring wounds, another study defined the distinct chromatin landscapes of regenerative and scarring fibroblasts (Abbasi et al, 2020), although the competent regenerative fibroblasts might arise from the periphery regenerative fibroblasts (Phan et al, 2021).Biernaskie and colleagues have also discovered the inaccessible inflammatory regulome in reindeer fibroblasts which could lead to the regenerative states found in reindeer's velvet (Sinha et al, 2022).Altogether, these studies indicated that the dynamic cellular state of dermal fibroblasts is defined by accessible chromatin profiles, and that manipulation of distinct regulatory mechanisms to achieve a regenerative fibroblast state might hold the key to promoting scarless wound healing.Alternatively, there may not be an advantage for achieving a 'stem-like' fibroblast state but to utilize and achieve optimal populations of different committed fibroblast populations that arise from their distinct lineages.Our current study provides an in-depth analysis looking at the changes in chromatin architecture happening during embryonic development while establishing the ground state of accessible chromatin, transcriptomic, and epigenetic profiles for DFPs and differentiating fibroblasts.
In the context of aging, dermal fibroblasts have been shown to have lost their cellular identity at transcriptomic and epigenetic levels (L opez-Ot ın et al, 2013;Salzer et al, 2018;Yang et al, 2023).Future studies will require the identification of the homeostatic chromatin states in different fibroblast populations to activate fibroblasts to resist the effect of chronic aging.Our molecular blueprint will provide the reference points of differentiating fibroblasts in embryogenesis and postnatal dermis to help determine the chromatin profiles of how fibroblast heterogeneity was established.

Single-cell suspension digestion
Tissues were collected from embryonic E14.5, E18.5, and P5 mice skins.Cells isolation procedure described previously was followed to obtain single-cell suspension from tissues (Jensen et al, 2010).

Histological analysis
Tissues were collected from mice back skins.For formalin-fixed paraffin-embedded, tissues were fixed with 4% paraformaldehyde in 1× PBS overnight.Hematoxylin and Eosin staining were done with 5 lm-thick sections.Colored images were taken with Nikon E600 and Nikon DS-Fi3 camera.Horizontal wholemount preparation was done as previously published (Salz & Driskell, 2017).Immunofluorescent images were taken with Leica SP8 confocal microscope at 20× objective.Analysis and quantification of H&E images were given unique ID and scored by a blinded researcher.

Chamber allografting assay
This procedure was performed as previously described (Jensen et al, 2010).A 8 × 10 6 WT neonatal epidermal cells were combined with 10 6 dermal cells from different timepoints or from different transgenic lines.Combination of cells were injected into silicone chambers grafted on 6 mm circular dorsal wounds in Foxn1 À/À mice.Collection of grafted areas were done 28 days post-surgery.
Low-quality and doublet cells were filtered by larger than 5,000 genes or less than 700 genes, less than 10 percent mitochondrial genes.Data were normalized to 10,000 reads per cells, log transformed, then scaled to maximum standard deviation 10.Integration and batch-correction by Time were done using Scanorama (Hie et al, 2019).Differential gene expression was calculated using Diffxpy (https://github.com/theislab/diffxpy).
Partition-based graph abstraction (PAGA) was used to infer the trajectory of continuous cell transitions while preserving the global topology of data (Wolf et al, 2019).RNA velocity analysis was done using dynamical model to retrieve latent time (Bergen et al, 2020).
For the time points scRNA-seq in Fig 1, we aimed to investigate the differentiation trajectories of dermal progenitor fibroblasts across three different time points that were sequenced separately.To establish the continuous transition states of fibroblasts among these development time points, Scanorama was used to integrate and correct for batch effect.We next utilized PAGA to provide.
All codes used to analyze this data is available on our github page.
Single-cell assay for transposase-accessible chromatin Dermal cells were collected and digested into the single cell suspension from 3 individual mouse embryos pooled together for each sample.Nuclei isolation from single cell suspension were prepared targeting 10,000 nuclei per library.Isolated nuclei were processed using 10× Genomics Single Cell ATAC v2 (CG000496).Libraries were sequenced on Illumina NovaSeq 6000 (100 bp Paired-end).Fastq files were aligned to mm10 reference genome using Cellranger-atac version 2.1.0.Cellranger outputs including Peaks called, Fragments, and metadata were used for downstream analysis with Signac in R (Stuart et al, 2021).Basic filtering and preprocessing were done screening for Peak region fragment < 100,000 and Percent reads in peaks > 40%.Clustering analysis was done using latent semantic indexing (LSI) (Cusanovich et al, 2015), dimension reduction was carried out by UMAP (McInnes et al, 2018).scRNAseq integration with scATAC-seq data for cell type classification was done by Signac Label Transferred methods using scRNA-seq fibroblasts object and scATAC-seq fibroblast object (Stuart et al, 2021).To improve the peaks calling sensitivity and accuracy, we performed Peak Call within Signac using MACS2 function (Zhang et al, 2008).Differential accessible peaks between clusters were calculated using FindMarkers() function within Signac.Motif binding sites analysis of specific cell type was done using Motif-MatchR (https://greenleaflab.github.io/motifmatchr/).
To generate the peaks heatmap in Fig 2, we performed the differential peaks analysis to isolate the list of peaks associated with distinct fibroblast clusters in scATAC-seq data in Signac.Within each cluster, we filtered and retained the list of peaks that are associated with genes found in the velocity driver genes of the same cell type from scRNA-seq analysis.
Codes used to analyze this data is available on our github page.

Flow cytometry
E18.5 single-cell suspension from WT and Kdm6bcKO mice were stained with primary antibody and secondary antibody for 1 h each with 3× PBS washes in between and after.Flow sorting of cells were done on Sony SH800 Flow Sorter.Flow analysis were done on Ther-moFisher Attune Flow Cytometer.

Statistical analysis
Histology data comparing two groups were tested for normality using Shapiro-Wilk, equal distribution using Levene, and t-test for independent samples if assumptions are met, otherwise we used Wilcoxon Rank Sum test.Single-cell RNA and ATAC data statistical tests were done using nonparametric Wilcoxon Rank Sum test.

Figure 1 .
Figure 1.DFPs are not intrinsically programmed to support skin reformation.A-C H&E histology of murine skin sections at Embryonic Day E14.5 (A), E18.5 (B), and Post-natal day 5 (P5) (C).D-E PAGA-initialized single-cell embedding of fibroblasts of 3 timepoints integrated.(F) Expression of Fibroblast Marker Genes.G-J (G) Ex vivo chamber grafting assay testing the ability of fibroblast populations at E14.5, E18.5, and P5 to support hair follicle formation.(H) Photographs of area of reconstituted skin.(I) Histology of reconstituted skin in chamber grafting assay.(J) Quantification of the number of hair formations in grafted area.Source data are available online for this figure.
5) (Fig 2A and B).Standard quantification of the ATAC peak number per cell comparing E14.5 and E18.5 fibroblasts revealed a large increase in open chromatin regions during the differentiation process (Fig 2C).Differential peak analysis was performed using Wilcoxon testing between time identified 6,200 peaks in E14.5 and 13,566 peaks in E18.5 fibroblasts.Coverage plots depicting the significant changes in open chromatin regions between E14.5 and E18.5 highlighted the increase of accessible peaks during the differentiation process of DFPs (Fig 2D).Our results indicate that chromatin reprogramming of DFPs during differentiation is associated with a selective, derepression of the genome.
fibroblast cell type from DFPs to Papillary to Dermal Papilla and DFPs to Reticular to Preadipocyte (Materials and Methods; Fig 2F).We found that as Upper DFPs differentiated toward Papillary and Dermal Papilla, the velocity driver genes of Papillary and Dermal Papilla became more accessible (Fig 2F).Notably, the driver genes of Upper DFPs were accessible across all fibroblast subtypes.Dermal Papilla had accessible chromatin regions for Papillary driver genes, but DP-specific genes were not accessible in Papillary, indicating the differentiated fate of Dermal Papilla as previously reported (Fig 2F; REF).This finding confirmed the continuous transition of Upper DFPs to Papillary to Dermal Papilla, where Upper DFPs require making their chromatin landscape specifically accessible during their differentiation process.The closed chromatin profiles of DFPs at E14.5 without accessible regions for Dermal Papilla driver genes might help explain the inability of DFPs to reform hair in chamber grafting assays.Similarly, the Lower DFPs differentiation trajectory was also observed using Lower DFPs, Reticular, and Preadipocytes RNA velocity driver genes (Fig 2F).Overall, our analysis revealed that during differentiation, DFPs specifically make their chromatin accessible in regions associated with the specific driver genes at E18.5 in a stepwise manner (Fig 2F).We also concluded that the differentiated fibroblast fates of Dermal Papilla and Preadipocytes required unique reprogramming of dermal fibroblasts with specifically accessible regions such as Sox18, Stmn2 for Dermal Papilla and Fabp4, AdipoQ for Preadipocytes (Fig 2G).

Figure 2 .
Figure 2. DFP differentiation is defined by inaccessible chromatin architecture.A Experimental design of scATAC-seq experiments comparing E14.5 and E18.5 dermal fibroblasts.B UMAP plot showing the distinct clustering of E14.5 and E18.5 fibroblasts.C Coverage Plots highlighting the significant increase in open accessible chromatin regions from E14.5 to E18.5 fibroblasts.D Quantification of Total ATAC Peaks per cell.E UMAP plot showing the classification of fibroblast cell types in scATAC-seq data.F Pseudotime-heatmap of scATAC-seq peaks associated with RNA velocity genes required for Dermal Papilla and Preadipocyte differentiation.Dermal Papilla and Preadipocyte populations are flanking the E14.5 DFPs which are located in the centrally.G Tracks-plot arranged in pseudotime showing specific accessibility regions of Sox18 and Stmn2 in Dermal Papilla, Fabp4 and AdipoQ in Preadipocytes.Source data are available online for this figure.
development of one differentiated cell type-Dermal Papilla-was also inhibited, with fewer Dermal Papilla per hair follicle (Fig 4E).Flow cytometry staining for CD133/ Prom1, a surface marker of Dermal Papilla, confirmed a significant 20% decrease in their number from whole dermal preparation (Fig 4F).To extensively characterize the changes induced by Jmjd3 ablation, we performed a scRNA-seq experiment comparing the E18.5 dermis of Kdm6b-cKO mice to the littermate genetic control (Fig 5A).After preprocessing and filtering steps, 12,947 fibroblasts were retrieved for downstream analysis (Fig 5B).PAGA embedding of single-cell clusters was performed to visualize the trajectory and continuous transition among different fibroblast cell types.In addition, we performed cluster classification as in Fig 1D-F (Figs 5B and

Figure 3 .
Figure 3. Differentiation of DFPs during embryonic development is coupled with demethylation of H3K27me3.A Experimental design of H3K27me3 ChIP-seq comparing E14.5 and E18.5 Fibroblasts.B Heatmaps of H3K27me23 Peaks comparing E14.5 and E18.5 Fibroblasts.C Quantification showing the decrease in H3K27me3 peaks (19,060 peaks at E14.5 and 8,768 peaks at E18.5, n = 2).D Venn diagram revealing the distribution of genes associated with H3K27me3 peaks changes between fibroblasts from E14.5 to E18.5.E Quantification of significantly differentiate H3K27me3 binding between E14.5 and E18.5 fibroblasts (4,699 peaks at E14.5 and 24 peaks at E18.5).F GO Analysis of Significantly De-repressed Genes revealed important Fibroblast differentiation processes.G Track plot of important de-repressed genes associated with fibroblast differentiation (Sox18 and Alpl for Dermal Papilla, Cebpa and Pparg for Preadipocyte).Source data are available online for this figure.

Figure 5 .
Figure 5. scRNA-seq analysis of Wild Type and Kdmb6b-cKO mouse fibroblasts.A Experimental strategy to generate single-cell RNA-seq data from E18.5 WT and Kdm6b-cKO mice.B, C PAGA embedded plots showing cell clusters of integrated E18.5 WT and Kdm6b-cKO fibroblasts.D Bar plots showing the contribution of each condition (WT and Kdm6b-cKO) to each fibroblast sub-cell types.E GO analysis of Genes significantly downregulated in Kdm6b-cKO fibroblasts.F, G PAGA plots showing the expression levels of Adipogenic genes AdipoQ and Fabp4 between WT and Kdm6b-cKO.Source data are available online for this figure.
transduction, Cell differentiation, Multicellular organism development, Cell-cell signaling, and Neuron differentiation (Fig EV4C).We did not observe processes directly related to the skin, hair follicle, and dermal white adipocyte differentiation being enriched, even in the present of important genes for Dermal Papilla and Preadipocyte such as Cebpa, Hes5, Rspo1, and Eya2 (Fig 6E; Sennett et al, 2015; Mok

Figure 6 .
Figure 6.Kdm6b regulates epigenetic de-repression of dual-function Adipogenic and Dermal Papilla transcription factors.A Experimental design for H3K27me3 ChIP-seq comparing E18.5 WT and Kdm6b-cKO.B Heatmaps depicting the changes in H3K27me3 peaks in Kdm6b-cKO fibroblasts.C Bar plots showing the number of H3K27me3 peaks called in E18.5 WT and Kdm6b-cKO Fibroblasts.D Track plots of relevant transcription factors loci for Dermal Papilla and Preadipocyte lineage commitment (Cebpa, Rspo1, Eya2, and Hes5).E Homer motif analysis of significantly repressed peaks in Kdm6b-cKO fibroblasts.F scRNA-seq expression levels of Motifs regulated by Kdm6b activity.Source data are available online for this figure.
Fig 4C).Our single-cell RNA-seq experiment comparing WT and Kdm6b-cKO skin revealed a major decrease of Preadipocytes in Kdm6b-cKObut not Dermal Papilla (Fig 4I).It is possible that the existence of other epigenetic regulators and histone marks such as bivalency with H3K4me3, or the acetylation of H3K27me3, may regulate the formation but not the function of Dermal Papilla.Our scRNA-seq comparing WT and Kdm6b-cKO skin also detected the significant upregulation of Chromatin Modification GO process, indicating a compensation mechanism to remodel and program the chromatin profile bypassing H3K27me3 demethylation (Fig EV4D and E).Our H3K27me3 ChIP-seq evaluating E18.5 WT and Kdm6b-cKOfibroblasts validated the ablation of Kdm6b/Jmdj3 which led to the increase in H3K27me3 levels across the genome.Motif analysis of significantly increased H3K27me3 peaks in Kdm6b-cKOfibroblasts revealed an array of transcription factors that might possess dual-functions.These transcription factors such as Ebf1, Zfp423, Tcf7l2 are not distinctly expressed across fibroblast lineages, even though they have been shown to play important roles in adipogenesis and Dermal Papilla function (Fig EV5A-D; Gupta et al, 2010; Figure EV1.scRNA-seq analysis.

Figure
Figure EV3.scRNA-seq analysis of all cells in WT and Dermo1Cre-Kdm6b fl/fl skin and differential expression analysis of genes from GO Analysis of WT and Kdm6b-cKO.A UMAP projections of WT and Dermo1Cre-Kdm6b fl/fl skin.B Quantification of the number of cells within each cluster type based on genotype represented as a percentage in the cluster.C Dotplot of gene utilized to identify cell clusters.D GO analysis of differentially regulated genes comparing WT and Dermo1Cre-Kdm6b fl/fl skin.E Gene expression of Chromatin Organization genes from GO Analysis.F-I Differential expression of GO-Terms genes from (F) Lipid Metabolic Process (G) Fatty Acid Metabolic Process (H) Brown Fat Cell Differentiation and (I) White Fat Cell Differentiation between WT and Kdm6b-cKO fibroblasts.

Figure EV5 .
Figure EV5.Kdm6b regulates dual function transcription factor bindings sites.A Model for Dual-Function Transcription Factor regulation of Dermal Papilla activity regulated by Kdm6b/Jmjd3 de-repression.B Dual-Function Transcription Factor, Ebf1, binding activities represented by a GO-analysis of Ebf1 motifs in scATAC-seq peaks from Dermal Papilla and Preadipocytes fibroblast populations.C scATAC-seq track-plot of Dermal Papilla genes, Sox18, Lamc3, and Fgf10.Predicted Ebf1 motifs are highlighted in a gray box on top of the peak that is predicted to have the motif.D Motif analysis of scATAC peaks for GO analysis.