ZFP42 helps keep human skin stem cells active and controls their daily gene activity through CRY1

To verify whether the keratinocytes we isolated from fetal and adult skin tissues were epidermal stem and progenitor cells (EPSCs), we analyzed the expression of several well-known surface markers for human EPSCs, including CD49f²⁶, CD34²⁷, CD71²⁸, and keratin 14 (a basal cell marker), using flow cytometry (Supplementary Fig. 1A, see Supplementary Data 4 for gating strategy). Our results showed that both groups of cells exhibited high expression levels of these markers. We also confirmed minimal expression of keratin 10 (an early differentiation marker) in both young and adult EPSCs (Supplementary Fig. 1B). Additionally, immunofluorescence staining for Mki67, along with the clonogenic assay, demonstrated that EPSCs from both fetal and adult skin tissues stained positively for Mki67 and retained the ability to form colonies (Supplementary Fig. 1C, D). These findings confirm that the isolated cells exhibit key characteristics of human EPSCs. Furthermore, our results align with previous studies that fetal EPSCs display enhanced stem cell-like properties, including higher expression of CD71 and Mki67, as well as larger colony sizes²¹.

To characterize the rhythmic genes in fetal and adult EPSCs, we initiated the process by synchronizing the cells with a brief serum pulse and subsequently harvested RNAs at intervals of 3 h within a 24-h period, as outlined in a previous study by Janich et al. ²⁹. Following this, RNA-Seq was conducted to profile the transcriptome of these samples (Fig. 1A). Our results indicated that all samples of fetal and adult EPSCs expressed genes at similar levels and exhibited greater similarity to themselves by Principal Component Analysis (PCA) analysis (Supplementary Fig. 1E, F).

Subsequently, we employed the Jonckheere–Terpstra–Kendall (JTK) algorithm to identify rhythmic genes in both fetal and adult EPSCs, resulting in the detection of 1291 rhythmic genes in fetal EPSCs and 905 rhythmic genes in adult EPSCs, with an overlap of 208 genes³⁰ (Supplementary Data 1A, B). These rhythmic genes also show significant overlap with findings from other studies in primary keratinocytes and human epidermis^29,31,32 (Supplementary Data 2A, B). We also assessed circadian rhythm parameters including Phase, Amplitude, and Mesor using a comprehensive and interactive tool, DiffCircaPipeline (Fig. 1B and Supplementary Data 1A, B)^29,31–33. Notably, several core clock genes, including BMAL1, PER1, NR1D2, and CRY1, showed significant rhythmicity in both sample sets, whereas NR1D1 in adult EPSCs (p = 0.062) and PER2 in fetal EPSCs (p = 0.065) exhibited only near-significant oscillatory patterns (Fig. 1C and Supplementary Fig. 1G). These findings support overall synchronization of the cultures, with gene-specific and context-dependent differences in the robustness of rhythmic oscillations. In both fetal and adult EPSCs, BMAL1 and CLOCK displayed similar peak times at different phases, with a phase delay of ~12 h. Conversely, the peak times of PER1 and PER2 were antiphasic to BMAL1, in accordance with BMAL1’s known antiphasic relationship with its core clock targets²⁹. Furthermore, a sequential pattern of peak times emerged, with BMAL1 initially peaking around CT5 post-synchronization, followed by NR1D1/2 and CLOCK peaking around CT10-CT15. Subsequently, PER1/2 reached their peaks around CT15-CT20, and CRY1/2 at CT20-CT24. Notably, BMAL1 and PER1/2 shared a similar period time of 24 h, while CLOCK and NR1D1/2 exhibited longer period times exceeding 24 h. These findings indicate that the successive oscillation expression patterns of core clock components were generally conserved throughout human epidermal development.

Interestingly, we observed a notable reduction in the overall amplitude of rhythmic genes, in adult EPSCs when compared to fetal EPSCs. This suggests a dampening effect on circadian rhythm due to aging (Fig. 1D). In addition to this, the overlapped 208 rhythmic genes in adult EPSCs exhibited an average delay of ~3 h in peak times, and significant phase disparities were detected within the CT8–CT14 time frame (Fig. 1E, F). These findings underscore the existence of substantial disparities in rhythmic gene expression patterns between fetal and adult EPSCs, warranting further investigation.

To comprehend the function of the conserved rhythmic genes between fetal and adult EPSCs, we generated a heatmap for these 208 overlapped genes, organized based on the order of fetal EPSCs (Supplementary Fig.). Subsequently, we conducted Gene Ontology (GO) and KEGG analyses on the set of 208 overlapping genes. Both the GO terms and KEGG pathways revealed significant enrichment in terms such as DNA repair, DNA replication, cell cycle, and DNA synthesis, which are biological processes closely related to circadian regulation (Supplementary Fig.). 2A 2B

To delve deeper into the differences between fetal and adult EPSCs regarding the Phase shift, Amplitude, and Mesor of these 208 rhythmic genes, we employed a comprehensive and interactive tool, DiffCircaPipeline³³. Intriguingly, we identified 128 genes with significant phase shift differences between fetal and adult EPSCs, and these genes were enriched in GO terms and KEGG pathways related to antiviral mechanism by IFN-stimulated genes, defense response to virus, DNA damage response, and cell cycle (Supplementary Fig. 2C). The amplitude of these 208 genes exhibited significant overall differences (Supplementary Fig. 2D). Among them, 112 genes with altered amplitudes were enriched in GO terms associated with the cell cycle, DNA damage, and regulation of cellular senescence, as well as KEGG pathways such as the cell cycle, transcriptional regulation by TP53, and apoptosis (Supplementary Fig. 2E). Furthermore, 155 genes displayed differences in Mesor and were also in terms related to the cell cycle, DNA replication, and associated GO terms and KEGG pathways. Notably, GO terms like “S-adenosylmethionine-dependent methyltransferase”, “tricarboxylic acid (TCA) cycle”, “deubiquitination”, and “TP53 regulation of transcription of genes involved in G2 cell cycle arrest” were also observed (Supplementary Fig. 2F). Several rhythmic genes (ARG2, BARD1, CLSPN, NEIL3) common in fetal and adult EPSCs associated with cell senescence, apoptosis, DNA replication, and DNA repair are profiled here (Supplementary Fig. 2G). These findings suggest that essential biological processes are under the regulation of rhythmic genes during development, which may be conserved during adult stem cell aging.

To elucidate the functional characteristics of the distinctive rhythmic genes with successive peaks in fetal and adult EPSCs, we employed hierarchical clustering for the time series genes and analyzed the enriched GO terms within each group using the R package ClusterGVi, respectively. As depicted in Fig. 2A, the earliest peak genes in fetal EPSCs were associated with GO terms related to wound healing and acute inflammatory responses. These were followed by genes linked to ribonucleoprotein complex biogenesis, the regulation of metabolic processes, and mRNA processing. Subsequently, a substantial cluster of genes (n = 478) related to DNA replication, mitotic nuclear division, and nuclear chromosome segregation exhibited peaks around CT6 to CT15. Afterward, genes associated with the regulation of mitochondrial membrane potential and cytoplasmic translation displayed their peaks.

In contrast, the initial group of genes with the earliest peaks in adult EPSCs was connected to epidermis development and keratinization at around CT3 (Fig. 2B). This was followed by genes related to DNA replication, nuclear division, and double-strand break repair (n = 338) that peaked from CT9 to CT15. Then, a cluster of 161 genes associated with protein localization to the nucleus and cell projection assembly showed peak expression at a similar time as another group of 175 genes linked to the non-canonical Wnt signaling pathway. Finally, GO terms encompassing carboxylic acid biosynthetic processes and unsaturated fatty acid biosynthetic processes, along with genes related to mitochondrial ATP synthesis coupled with electron transport, exhibited peaks around CT18 and CT24, respectively. Notably, the pattern where keratinocyte differentiation-related genes peaked earlier than genes associated with DNA replication and cell division (around late-night to early-morning hours vs. afternoon and evening) was consistent with earlier findings by Janich et al.²⁹. The profiles of several representative rhythmic genes unique in fetal EPSCs and adult EPSCs were displayed. Cell division-related gene PRC1 and mitosis-related gene BUB1 are rhythmic only in fetal EPSCs. While epidermal differentiation-associated markers, including ZNF750 and GRHL3, are rhythmic only in adult EPSCs (Fig. 2C, D). In contrast to the emphasis on biological functions in the DNA replication and cell cycle in fetal EPSCs, the functional segregation of pro-proliferative and differentiation cues enables these adult epidermal stem cells to maintain daily homeostasis and prepare for daylight UV radiation or external stressors.

To investigate the underlying molecular mechanisms that contribute to the differences in circadian rhythmicity between fetal and adult EPSCs, we performed the ChEA analysis for the unique rhythmic genes in each group. The analysis showed that transcription factors such as YY1 and GABP are potential regulators of rhythmic gene expression in fetal EPSCs, while AF4, FOXA1, and MITF may exert a more prominent role in regulating rhythmic gene expression in adult EPSCs (Fig. 2E, F). It’s worth noting that MYC is an oncoprotein known to influence the expression of molecular clock genes, potentially disrupting the molecular clock oscillation in various cancer types³⁴. Additionally, GABP proteins, characterized as 12-hour pacemakers in hepatic cells are significant transcription factors in this context³⁵. Moreover, FOXA1 has recently been associated with its role in prostate cancer therapy as a downstream factor of ARNTL³⁶. MITF has also been observed to exhibit rhythmic expression with a 24-hour periodicity in the presence of BMAL1³⁷. The promoter motif analysis of the rhythmic genes in fetal EPSCs also confirmed the binding motif of YY1 family protein (ZFP42) as one of the top motifs enriched (Fig. 2G). In contrast, GRHL2 and NF-κB were among the top promoter motifs enriched in rhythmic genes of adult EPSCs (Fig. 2H). These findings indicate the existence of an intricate regulatory network involving transcription factors and their potential impact on circadian rhythmicity in EPSCs at different developmental stages in the skin.

Next, we examined the significantly differentially expressed genes (DEGs) between fetal and adult EPSCs by comparing the paired time-series data from these two groups. Our analysis revealed a total of 691 genes exhibited high expression in fetal EPSCs, while 929 genes were highly expressed in adult EPSCs with a cutoff at a fold change ≥2 and a p-value < 0.05. The most enriched GO terms in fetal EPSCs included extracellular matrix organization, cell junction assembly, and collagen metabolic processes. In contrast, highly expressed genes in adult EPSCs demonstrated enrichment in GO terms like epidermal cell differentiation, epidermis development, and negative regulation of cell population proliferation (Supplementary Fig. 3A, B). To validate these findings, we performed qPCR on representative genes associated with extracellular matrix organization and epidermis development in multiple fetal, young and older adult EPSC donors (Supplementary Fig. 3C, D). These results suggest that young and older adult EPSCs are more differentiated with less proliferative potential than fetal EPSCs. Furthermore, upon examining the overlap between highly expressed genes and their respective rhythmic genes, we observed that 64 rhythmic and highly expressed genes in fetal EPSCs were enriched for cellular processes, developmental processes, and the regulation of biological processes (Supplementary Fig. 3E, F). In contrast, the 38 rhythmic and highly expressed genes in adult EPSCs were enriched for skin development, establishment of skin barrier, and immune response (Supplementary Fig. 3G, H). However, it’s essential to highlight that a substantial portion of the rhythmically expressed genes are not necessarily highly expressed. These findings further underscore the critical role of genes with circadian rhythm may be essential in the control of human epidermis development and epidermal function.

Among the top 100 DEGs between fetal and adult EPSCs, we observed the presence of ZFP42, a YY1 subfamily protein and well-known pluripotent transcription factor³⁸ that was highly enriched in fetal EPSCs. In contrast, cell senescence gene CDKN2A and epidermal differentiation genes (SPRR1B and IVL) were highly expressed in adult EPSCs (Supplementary Fig. 3I). We further validated the ZFP42 expression in fetal, adult and old EPSCs from different donors using western blot and flow cytometry (Supplementary Fig. 4A–C, see Supplementary Data 4 for gating strategy). ZFP42, also known as Rex1, has been extensively studied as a pluripotent stem cell marker that is regulated by Oct3/4. More recently, it has been studied as a member of the YY1 subfamily. It was discovered to originate from a retroposition duplication of YY1 in eutherian mammals, and they share some similarities at their core regions³⁹. Our ChEA analysis and promoter motif analysis of the fetal rhythmic genes showed significant enrichment for ZFP42 binding sites, which suggests that it may regulate the unique rhythmic genes in fetal EPSCs (Fig. 2E, G). Thus, we suspect that ZFP42 could potentially be involved in the circadian regulation of fetal EPSCs.

Due to limited availability of fetal EPSCs, we proceeded to perform ZFP42 knockdown experiments in young EPSCs (derived from young children) using siRNAs or shRNAs. Depletion of ZFP42 significantly inhibited cell growth and decreased cell viability as indicated by cell proliferation and cell viability measured by CCK8 assay (Fig. 3A, B and Supplementary Fig. 4D). Cell cycle analysis also showed that the knockdown of ZFP42 resulted in an accumulation of cells in the G1 phase and a decrease of cells in the S phase, suggesting DNA synthesis in the ZFP42 knockdown cells was significantly reduced (Fig. 3C, see Supplementary Data 4 for gating strategy). Knockdown of ZFP42 also increased cell apoptosis (Fig. 3D, see Supplementary Data 4 for gating strategy). Furthermore, a decrease in the expression of proliferation-related genes, such as FGFBP1, KRT15, MYC, PRC1, and BUB1, was observed (Fig. 3E). Meanwhile, the knockdown of ZFP42 promoted the expression of differentiation-related genes, including SPRR1A, SPRR1B, and SPINK5 (Fig. 3F). Immunostaining also confirmed the increase of differentiation protein, Keratin 1, and the decrease of MKI67 positive cells in ZFP42 knockdown cells (Fig. 3G and Supplementary Fig. 4E).

To gain a deeper understanding of the molecular changes resulting from ZFP42 knockdown in fetal EPSCs, we conducted RNA-Seq analysis to profile the transcriptome alterations. As expected, we found that 583 upregulated genes were enriched for processes related to tissue development, epidermal cell differentiation, and cornification, while 764 downregulated genes were associated with multicellular organism signaling and chemical synaptic transmission (Fig. 3H–J and Supplementary Data 3A). Remarkably, approximately 20% of these DEGs following ZFP42 knockdown overlapped with DEGs observed when comparing fetal and adult EPSCs (Fig. 3K). This set of 261 overlapping genes was enriched in functions related to epidermal cell differentiation, epidermis development, and positive regulation of nuclear division (Fig. 3L).

More importantly, overexpression of ZFP42 rejuvenated the colony formation abilities in EPSCs from older adults (harvested from subjects ≥ 60 years old) (Fig. 4A, B), demonstrating ZFP42’s pivotal role as a stemness regulator in human EPSCs. Restoration of ZFP42 exprtession also rescued the expression of cell cycle-associated genes (Fig. 4C). Notably, however, ZFP42 itself is not rhythmic in either fetal (p = 0.080) or adult EPSCs (p = 0.509) (Supplementary Fig. 4F). In fetal EPSCs, ZFP42 showed modest positive correlations KLF4 (r = 0.31, p = 0.404) and NANOG (r = 0.31, p = 0.415), whereas no meaningful correlation was detected in adult EPSCs (KLF4: r = −0.06, p = 0.869; NANOG: r = −0.02, p = 0.961) (Supplementary Fig. 4G). Because these p-values exceed the conventional significance threshold, we interpret these trends cautiously and do not claim statistical significance. Together, our result suggests that ZFP42 may influence fetal EPSCs’ stemness genes by regulating clock genes rather than possessing inherent rhythmicity.

We then closely examined the rhythmic genes that were regulated by ZFP42, and interestingly, we found that 5.27% of genes were rhythmic in adult EPSCs, while 6.60% were rhythmic in fetal EPSCs. Additionally, 5.04% of genes were rhythmic in both sample sets, including the core circadian gene CRY1 (Fig. 4D). Furthermore, we assessed the expression of core clock genes (BMAL1, CLOCK, PER1/2, CRY1/2, and NR1D1/2) following ZFP42 knockdown and discovered that only CRY1 exhibited significant downregulation at both the mRNA and protein levels (Fig. 4E, F). Also, the decrease of CRY1 upon ZFP42 knockdown was not time-specific (Supplementary Fig. 4I). Consistent with our hypothesis that ZFP42 transcriptionally regulates CRY1, we observed a significant positive linear correlation between ZFP42 and CRY1 expression in fetal EPSCs (r = 0.49, p = 0.168), whereas no correlation was detected in adult EPSCs (r = −0.14, p = 0.717) (Supplementary Fig. 4G). Although the fetal correlations do not meet the p < 0.05 criterion, we report it as a positive biological trend that requires further validation with an increased sample size. For completeness, Spearman’s rank correlation is also provided in the panel (Supplementary Fig. 4H). Knockdown of ZFP42 also disrupts the rhythmicity of CRY1 and CLSPN in EPSCs, which is similar to changes (increase in amplitude and a phase shift) from fetal to adult EPSCs (Fig. 4G–I, and Supplementary Fig. 4J). We also confirmed the binding of ZFP42 to the promoter region (–420 to –172 bp upstream of transcription start site) of CRY1 by ChIP-QPCR (Fig. 4J). A promoter activity assay using a Dual-Luciferase reporter system was conducted to analyze CRY1 promoter activity in both ZFP42 knockdown and overexpression cells, as well as their respective controls (Fig. 4K, L). The results showed a significant reduction in CRY1 promoter activity following ZFP42 knockdown, whereas its activity increased with ZFP42 overexpression. As a positive control, we used CCNB1, a previously reported ZFP42 target gene⁴⁰, to validate the assay (Supplementary Fig. 4K, L). These findings suggest that CRY1 may act downstream of ZFP42, functioning as a potent stemness regulator in human EPSCs.

To investigate this, we knocked down CRY1 using three independent siRNAs. Similar to the effects of ZFP42 knockdown, the loss of CRY1 significantly inhibited cell growth and promoted epidermal differentiation (Fig. 5A and Supplementary Fig. 5A–C). Notably, cell cycle analysis indicated that the knockdown of CRY1 led to an accumulation of DNA synthesis in the G1 phase and a decrease in the S and G2 phases (Fig. 5B, see Supplementary Data 4 for gating strategy). Immunofluorescence staining of CRY1 knockdown cells confirmed an increased expression of epidermal differentiation marker Keratin 1 and a decreased expression of Mki67, mirroring the effects observed in ZFP42 knockdown cells (Fig. 5C and Supplementary Fig. 5D). Furthermore, our RNA-Seq analysis of CRY1 knockdown cells identified a total of 689 upregulated genes and 752 downregulated genes, with a cutoff at a fold change >=2 and a p-value < 0.05. The GO analysis of the upregulated genes indicated enrichment in terms related to epidermal cell differentiation, keratinocyte differentiation, and skin development. Conversely, the downregulated genes were enriched in GO terms like mitotic sister chromatid segregation, mitotic spindle organization, and positive regulation of the cell cycle process (Fig. 5D–F and Supplementary Data 3B). A significant portion (24.3%) of the DEGs resulting from CRY1 knockdown, also overlapped with DEGs observed when comparing fetal and adult EPSCs (Supplementary Fig. 5E). Once again, these 350 overlapped genes were enriched in processes related to epidermal development, extracellular matrix organization, and the regulation of cell population proliferation (Supplementary Fig. 5F). Collectively, these findings indicate the significant role of CRY1, a core circadian gene, in maintaining the self-renewal and differentiation of epidermal stem cells.

To validate, QPCR was employed to confirm the decreased expression of the cell cycle-related genes, such as MKI67, PRC1, and BUB1, and the increase of the differentiation-inducing genes and markers, including GRHL3, KRTDAP, SCEL, and FLG upon CRY1 knockdown (Fig. 5G, H). By comparing the DEGs between ZFP42i and CRY1i RNA-Seq, we found 171 (~13%) genes overlapped (Supplementary Fig. 5G). Moreover, overexpression of CRY1 did not alter the expression of ZFP42 but was able to rescue the decrease in the expression of proliferation gene expressions resulting from ZFP42 knockdown, including MKI67, BUB1, CDK1, CDK2, and YAP1 (Fig. 5I and Supplementary Fig. 5H). To confirm the rescue of ZFP42 knockdown phenotype by CRY1 overexpression, we performed colony formation, cell cycle analysis, and immunostaining for Mki67. The results showed that CRY1 overexpression partially restored the EPSCs stemness, as evidenced by the increased number of colonies, a higher percentage of Mki67 positive cells and an increase of S phase, consistent with a decrease of CDKN1B expression (Fig. 5I–K, Supplementary Fig. 5I–K, see Supplementary Data 4 for gating strategy). However, it did not prevent the increase in differentiation genes upon ZFP42 knockdown (Supplementary Fig. 5L). These findings suggest that CRY1 plays a more prominent role in regulating cell proliferation rather than differentiation in EPSCs.

Furthermore, we analyzed the rhythmic genes that were influenced by CRY1 manipulation and found that 14.22% of the DEGs were rhythmic in fetal EPSCs, compared to 10.51% in adult EPSCs, with 13.8% displaying rhythmicity in both fetal and adult EPSCs (Supplementary Fig. 5M). In total, ~40% of the DEGs upon CRY1 knockdown exhibited rhythmic patterns, underscoring the crucial role that core clock genes play in the biological processes related to human epidermis homeostasis. Interestingly, restoring CRY1 expression in ZFP42 knockdown cells rescued the rhythmicity of proliferation genes (BUB1 and PRC1), but did not restore rhythmicity for differentiation-inducing genes (ZNF750 and GRHL3) (Fig. 5L, Supplementary Fig. 6A). Our data suggest that ZFP42 regulates circadian output partly through CRY1 and partly via CRY1-independent pathways.

Recent studies have characterized the chromatin binding of CRY1 in both human and mouse genomes. These studies also demonstrated non-traditional role of CRY1 in terms of regulating and reprogramming induced pluripotent stem cells and DNA repair in cancer^41–43. We were able to identify CRY1 binding to the promoter regions of proliferation genes such as CDK1, CDK2, and YAP1 in both human and mouse datasets (Supplementary Fig. 6B, C). We further validated these binding sites by ChIP-QPCR in both control and CRY1 knockdown cells (Fig. 5M). In summary, our data suggested that CRY1 is regulated by ZFP42 and also functions in maintaining proliferating genes and inhibiting differentiation in human EPSCs.

Primary human epidermal keratinocytes were isolated from the discarded foreskin of young children (aged between 5–10 years old, 5 donors) who underwent circumcision or adult normal skin after plastic surgery (aged between 24 and 69 years old, 5 donors). Fetal epidermal keratinocytes were isolated from fetal skin of natural aborted fetus (~20 weeks, 2 donors). Skin from both genders were included. Collection was performed with patient and/or parental consent under an Institutional Review Board protocol approved by The First Affiliated Hospital of Sun Yat-sen University (reference no. 2022-350). In all instances, samples were obtained in de-identified manner and constituted non-human subjects research. The study was conducted in accordance with the Declaration of Helsinki principles. All ethical regulations relevant to human research participants were followed.

The skin’s subcutaneous tissue was excised using scissors and washed twice with Ca²⁺ and Mg²⁺-free PBS. The excised tissue was then incubated overnight at 4 °C in a 3 mg/mL Dispase solution (Sigma: D4693-1g). The following day, the epidermis was carefully separated from the dermis using forceps. The epidermal sheets were rinsed with PBS and then treated with 0.25% trypsin solution at 37 °C for 5–7 min before being neutralized with 10% fetal bovine serum (FBS). The resulting epidermal cells were filtered through a 70 µm cell strainer, centrifuged at 350 × g for 6 min, and the supernatant was discarded.

Primary fetal and adult keratinocytes were propagated in low calcium Epilife medium (Thermo Fisher Scientific: MEPI500CA) with recommended supplements to maintain an undifferentiated, stem and progenitor-enriched cell population. The same passage (passage 3 or 4) of fetal and adult keratinocytes were cultured in Epilife medium and characterized for epidermal stem cell markers, including CD49f, CD34, CD71, KRT14, KRT10, and ZFP42 before being used for all experiments.

Passage 3 keratinocytes were resuspended using 100 µL staining buffer (PBS containing 3% BSA) at a concentration of 1 ×10⁵ cells/mL. The cells were labeled separately with FITC-conjugated CD49f (BioLegend: 313605), APC-conjugated CD71 (BioLegend: 334107), APC-conjugated CD34 (BioLegend: 128611), keratin 14 (Abcam: ab181595), keratin 10 (Abcam: ab76318) and ZFP42 (Invitrogen: 710190) primary antibodies or corresponding isotype controls. For keratin 14, keratin 10 and ZFP42 staining, goat anti-rabbit Alexa Fluor 488 (Invitrogen: A11034) was used as the secondary antibody. 5 μL of each antibody was added to the cell suspensions, followed by incubation at 4 °C for 1 h in the dark. After incubation, cells were centrifuged at 300 × g for 5 min, the supernatant was discarded, and the cells were washed twice with PBS containing 3% BSA. Finally, the cells were resuspended in PBS and analyzed using a CytoFlex flow cytometer (Beckman Coulter, USA). Data were analyzed using FlowJo v10.4 software (BD Biosciences).

Primary fetal and adult keratinocytes were propagated in low calcium Epilife medium (Thermo Fisher Scientific: MEPI500CA) with recommended supplements to maintain an undifferentiated, stem and progenitor enriched cell population. Same passage (passage 3 or 4) of fetal and adult keratinocytes were cultured in Epilife medium and used for all experiments. For synchronization, cells were treated with 20% FBS for 3 h. After that, cells were rinsed with PBS twice before changed back to the Epilife medium. Cells were harvested every 3 h within 24 h for rhythmic gene expression detection.

Control or knockdown cells were seeded in 96-well plates (2000 cells per well) and measured using Cell Counting Kit-8 (CCK-8, Dojindo) according to the manufacturer’s instructions at different timepoints.

For clonogenic assays with primary human keratinocytes, 500 cells were plated in a 6-well plate and were grown for 10–14 days in Epilife medium and then fixed in 4% paraformaldehyde and stained with 0.1% crystal violet/0.1% rhodanile blue (Beyotime: C0121).

Control and ZFP42 or CRY1 knockdown cells were stained with Annexin V conjugated to Alexa Fluor 488 (Life Technologies: A13201) and analyzed using the Beckman flow cytometer (Beckman Coulter, USA) according to manufacturers’ instructions.

Control and ZFP42 or CRY1 knockdown cells were detached and washed with PBS containing 1% FBS, and then fixed with 80% ethanol at 4 °C for 15 min. Fixed cells were stained with 50 µg/mL propidium iodide in a PBS buffer containing 1% FBS, 0.05% Triton X-100, and 10 µg/mL RNaseA. After a 30-min incubation at room temperature, a flow cytometer was used to determine the percentage of cells in different phases of the cell cycle (G1, S, and G2/M phase) with CytoFLEX software (Beckman Coulter, CA).

Total RNA from cells was extracted using Trizol reagent (Thermo Fisher Scientific: 15596026) according to the manufacturer’s instructions. RNA quality and quantity were analyzed using Nanophotometer (Implen). One μg of total RNA was reversed transcribed using the PrimeScript™ RT Reagent Kit (Takara: RR047). cDNA was synthesized using the ABI Veriti 96-Well Thermal Cycler (Applied Biosystems). Quantitative PCR was carried out with TB Green® Premix Ex Taq™ II (Takara: RR820) on an Applied Biosystems™ QuantStudio™ 5 system. Gene expression was calculated using the ^ΔΔCt method, with L32 used as internal control for normalization. The sequences of primers used in this study are provided in Supplementary Table 1.

Fetal or adult EPSCs were fixed in 4% paraformaldehyde for 11 min followed by blocking in PBS with 2.5% normal goat serum, 0.3% triton X-100, and 2% bovine serum albumin for 30 min. Primary antibodies used were Keratin 1 (Biolegend: 905204) at 1:200, MKI67 (Abcam: AB16667) at 1:200. The secondary antibodies used were Alexa 488 Goat anti-rabbit IgG (Thermo Fisher Scientific: A11034) at 1:500. Nuclear dye was stained with DAPI Fluoromount-G (Southern Biotech: 0100-20). Images were acquired using an Olympus/BX63 Upright Fluorescence Microscope.

Total protein was extracted using RIPA lysis buffer with protease inhibitors and quantified using BCA™ Protein Assay Kit. Fifty micrograms of protein from each sample were denatured and separated using BioRad Mini-PROTEAN Tetra System. The resolved proteins were transferred onto PVDF membranes (Milipore: IPVH00010) at 4 °C for two hours. Primary antibodies used include beta-actin (Santa Cruz: Sc-47778) at 1:30000, ZFP42 (Invitrogen: 710190) at 1:2000, CRY1 (Proteintech: 13474-1-AP) at 1:3000 overnight at 4 °C. Secondary antibodies including HRP conjugated anti-mouse IgG or anti-rabbit IgG (Cell Signaling Technology), were used at 1:2000.

Protein detection was performed using the enhanced chemiluminescence method (EpiZyme: SQ202) and analyzed with an Amersham ImageQuant 800 (Cytiva: 29399481).

To construct the ZFP42 overexpression vector, the ZFP42 cDNA fragment was cloned into the BamHI (Thermo Fisher Scientific: FD0054) and EcoRI (Thermo Fisher Scientific: FD0274) sites of the pCDNA3.1(+) vector to generate the recombinant pCDNA3.1-ZFP42 plasmid. For the luciferase reporter plasmids containing the CRY1 or CCNB1 promoters, promoter regions of CRY1 (-2000 to -80 bp relative to the translation start site) and CCNB1 (−2000 to +7 bp relative to the translation start site) were amplified from genomic DNA of primary human EPSCs. The amplified sequences were inserted into the XhoI (Thermo Fisher Scientific: FD0694) and HindIII (Thermo Fisher Scientific: FD0504) cloning sites of the pGL4.10[luc2] basic luciferase vector (Promega: E6651).

For the promoter activity assay, primary fetal EPSCs were seeded at a density of 50,000 cells per well in a 24-well plate. Cells were co-transfected with 1 μg of the pGL4.10[luc2] plasmid containing either the CRY1 or CCNB1 promoter, along with 1 μg of the pCDNA3.1-ZFP42 overexpression plasmid (or the pCDNA3.1-empty vector as a control), and 20 ng of pRL-SV40P plasmid (Addgene: 27163) as a normalization control. After 24 h of transfection, the culture medium was replaced with fresh medium. Forty-eight hours post-transfection, cells were harvested and lysed in luciferase lysis buffer and the luciferase activity was measured using the Dual-Luciferase Reporter Assay Kit (Vazyme: DL101-01) following the manufacturer’s instructions. Firefly and Renilla luciferase activities were quantified using the Dual-Luciferase Reporter Assay System (Promega). Firefly luciferase activity obtained from each sample was normalized to the Renilla luciferase activity from the same sample. Relative luciferase activity was calculated by comparing the normalized luciferase activity of knockdown or overexpressed cells to the control cells. Promoter sequences were summarized in Supplementary Table 1.

Two biological replicates of both fetal and adult EPSCs were pooled to generate RNA-Seq data across circadian time points. RNA-Seq was performed using the DNBSEQ platform (BGI, China) with PE150. About 40 million pair-ended reads per sample were obtained. Reads were aligned to the GENCODE v19 transcriptome hg19 using Bowtie2 (v2.3.4.3) with default settings⁶⁰. RSEM (v1.3.1) was used to generate gene and isoform expression levels and pheatmap (v1.0.8) was used to plot heatmap. Differential expression among samples was calculated using DESeq2⁶¹. Analysis of the read count distribution indicated that a threshold of ten reads per gene generally separated expressed from unexpressed genes, so all genes with fewer than ten reads were excluded from DESeq2 analysis. Gene lists for significantly upregulated or downregulated genes were created using p < 0.05 and 2-fold change. Enriched GO terms and KEGG analysis for RNA-seq differentially expressed gene sets were identified using Enrichr and Metascape^62–65. For de novo promoter motif analysis, transcription factor motif finding was p erformed with motif finder program HOMER with the command findMotifs.pl on the interested genes. The motif with the most significant p-value predicted by HOMER was presented. To conduct motif analysis on non-rhythmic genes, we performed a matched case-control enrichment analysis: for each rhythmic gene in fetal EPSCs, we selected a matched set of non-rhythmic genes, ensuring they were similar in factors such as promoter length and GC content. We repeated this analysis several times to ensure robustness.

To identify genes with a rhythmic expression from the fetal and adult EPSCs, the R package with Jonckheere–Terpstra–Kendall (JTK)_cycle (version:3.0) was used³⁰. A p < 0.05 threshold is applied to identify rhythmic genes. The expression levels for all genes were normalized as TPM values and used for JTK analysis. DiffCircaPipeline (version:0.99.0) R package was used to detect different circadian characteristics³³. Kronos/Cosinor (version:1.0.0) R package was used to produce the plots⁶⁶. R package ClusterGVis (version: 0.0.6) was used for better clustering and visualizing the time series gene expression data⁶⁷. Amplitude and phase estimations of oscillating genes were extracted from the DiffCircaPipeline algorithm.

Five million control and CRY1 or ZFP42 knockdown cells and 3 μg of antibody were used for each antibody pulldown experiment for ChIP^14,19,68. Cells for the CRY1 or ZFP42 pulldown ChIP-QPCR were fixed in both formaldehyde (1% final concentration, Thermo Fisher Scientific: 28908) and disuccinimidyl glutarate (2 mM final concentration, Thermo Fisher Scientific: 20593). As described before^18,19,69, the cells were treated with Farnham lysis buffer and sheared with a syringe to assist lysis. After shearing, the cells were centrifuged and resuspended in SDS lysis buffer and sonicated using a Bioruptor® Plus system (Diagenode). Once the appropriate fragment size was achieved, the lysate was incubated with CRY1 (Proteintech: 13474-1-AP) or ZFP42 (Invitrogen: 710190) and Rabbit IgG (Millipore: 12-370). 50 μl of dynabeads™ Protein G (Invitrogen: 10004D) were added to each sample and rotated for 4 h at 4 °C. Wahses were then carried out to reduce non-specific binding. The beads were then placed in an elution buffer and the supernatant was isolated. Finally, the DNA was de-crosslinked and used for downstream ChIP-qPCR analysis. qPCR results are represented as enrichment of pulldown over IGG. Primers for ChIP-qPCR were summarized in Supplementary Table 1.

Results are expressed as mean±s.d. unless stated otherwise. Statistical comparisons between two groups were evaluated by unpaired Student’s t test, two-tailed. Statistical comparisons between three groups were evaluated by one-way ANOVA followed by Tukey’s multiple comparision. A probability (p) value of <0.05 was considered to indicate statistical significance.

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

Communications Biology thanks Leo Kurian and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editors: Jun Wei Pek and Kaliya Georgieva.

The raw RNA-Seq datasets generated from this study have been deposited in the Genome Sequence Archive at the National Genomics Data Center (https://ngdc.cncb.ac.cn/gsa-human/↗) with accession code HRA007097. CRY1 ChIP-Seq datasets were downloaded from Gene Expression Omnibus Repository (GEO) with the accession codes GSE130602↗⁴¹ and GSE230321↗⁴². The source data for the graphs in the study are found in Supplementary Data 4. Uncropped and unedited blot images are presented in Supplementary Fig. 7. Any other data that support the findings of this study, are available within the article, its Supplementary Information and from the corresponding author upon reasonable request.