Transcription factor p63 bookmarks and regulates dynamic enhancers during epidermal differentiation

Gene expression dynamics during epidermal differentiation

Using an in vitro differentiation model of human epidermal keratinocytes ^16,20, we generated a complete map of dynamically regulated genes during epidermal differentiation at a genomewide scale using RNA polymerase II (RNAPII) ChIP-seq that measures the occupancy of RNAPII in mainly gene-associated genomic regions, as well as RNA-seq that measures the steady-state mRNA level. With both methods, the detected expression pattern of epidermal marker genes such as KRT5, KRT1, and involucrin (IVL) (Fig1A), and KRT14, KRT10 and loricrin (LOR) (Supplementary Fig S1A) was consistent with our RT–qPCR and Western blotting analyses (Supplementary Fig S1B and C). Expression of these marker genes was also in agreement with their expression dynamics in vivo ^20,21. At the global level, our RNA-seq analysis showed a good correlation (Spearman rho > 0.8 for all comparisons) with a previous RNA-seq study that used a slightly different in vitro keratinocyte differentiation method (Supplementary Fig S2A) ²². Thus, our results show that this in vitro differentiation model of epidermal keratinocytes recapitulates the epidermal stratification process. Of note, RNAPII ChIP-seq and RNA-seq profiles at the TP63 locus showed that ΔNp63 is the main isoform expressed in keratinocytes throughout the epidermal differentiation process, and the expression of the TAp63 isoform was not detected with both RNAPII ChIP-seq and RNA-seq methods (Supplementary Fig S1D). Western blot analysis confirmed the decreased protein level of ΔNp63 during differentiation and that the expression of TAp63 was undetectable (Supplementary Fig S1B). The decreased ΔNp63 expression at the mRNA and protein levels is in agreement with previous reports ^16,23.

To capture the gene expression dynamics during differentiation, we clustered the genes according to their expression levels quantified as reads per kilobase per million mapped reads (RPKM) in both the RNAPII ChIP-seq and RNA-seq profiles. For RNAPII ChIP-seq analyses, 7,223 genes with an RPKM value (expression level) of at least 1.0 were clustered using K-means clustering into nine groups with distinct expression patterns (Fig1B; Supplementary Table S1A). To assess the function of these gene clusters and their association with human diseases, we performed functional annotation analyses using the web tools DAVID Gene Ontology (GO) ²⁴ (Fig1B; Supplementary Table S1B) and the Human Phenotype Ontology (HPO) database ²⁵ (Supplementary Table S1C), respectively. For examples, genes that are expressed highly in proliferating keratinocytes (cluster 4) are associated with skin blistering and hyperkeratosis (e.g., KRT5). Genes that are highly expressed at late differentiation stages (cluster 3, e.g., LCE genes and cluster 6, e.g., TGM1) are specifically enriched for late differentiation genes and those involved in the regulation of apoptosis and cell death. These genes are associated with dry and scaly skin and non-blistering phenotypes. Cluster 1 contains genes with high expression at early differentiation stages (day 2 and day 4) and reduced expression at the late stage (day 7), for example KRT1 and KRT10, and is enriched with genes involved in epidermis development, blistering, and hyperkeratosis. For the RNA-seq analyses, we identified 3,813 genes that were significantly differentially expressed between keratinocyte differentiation stages (P-value < 0.05, cuffdiff2 software) and with an RPKM > 10 in at least one stage. Six distinct gene clusters were defined based on the expression dynamics of these genes using K-means clustering (Supplementary Fig S2B; Supplementary Table S2A). Similar to the RNAPII ChIP-seq data, the RNA-seq data also showed that genes with distinct expression patterns during epidermal differentiation are associated with distinct roles during epidermal differentiation and with different skin disease phenotypes (Supplementary Fig S2B; Supplementary Table S2B and C).

We further investigated and compared gene expression detected by RNAPII ChIP-seq and by RNA-seq. The global analysis showed that the range of expression levels (RPKM values) is narrower in the RNAPII ChIP-seq profile than in the RNA-seq profiles (Supplementary Fig S2C). Nevertheless, the relative expression levels of genes (gene ranks based on their RPKM) are generally similar in both RNAPII ChIP-seq and RNA-seq profiles (Supplementary Fig S2D, Supplementary Table S3A), as the relationship between RNAPII occupancy and steady-state mRNA level remains largely unchanged during epithelial differentiation for most genes. However, a small proportion of genes showed substantially higher RNAPII occupancy than mRNA levels (Supplementary Fig S2D, the bottom cluster), and these genes are strongly enriched for non-coding RNAs (mean-rank gene set enrichment, P-value < 10⁻⁹⁵ at all stages) (Supplementary Table S3C and D). Despite the global similarity, closer inspection reveals different dynamics of gene expression patterns between these two methods for many genes (Fig1A and C). For example, the highest expression level of KRT10 was detected at differentiation day 2 in RNAPII ChIP-seq, while its expression in the RNA-seq profile only reached its highest level at day 4 (Supplementary Fig S1A). The expression pattern detected in the RNA-seq analysis correlates with the protein levels (Supplementary Fig S1B). Consistent with the expression of individual marker genes, expression patterns of genes in the RNAPII ChIP-seq clusters were often more dynamic than their corresponding RNA-seq expression patterns, with the expression changes in the RNAPII ChIP-seq profiles preceding those in the RNA-seq profiles (Fig1C). These observations are consistent with RNAPII-mediated transcription preceding mRNA buildup and highlight differences between transcription rate control and steady-state mRNA level maintenance. Interestingly, Spearman correlation analysis showed a slightly better correlation of H3K27ac levels with the RNAPII occupancy at gene bodies than with the steady-state mRNA abundance (Supplementary Fig S2E). This may reflect the influence of post-transcriptional processes such as RNA stability and processing of mRNA by miRNAs on mRNA steady-state levels. These data suggest that RNAPII ChIP-seq is a preferred method to study the direct effect of transcription factors on gene regulation as compared to RNA-seq. Based on this rationale, we used the RNAPII occupancy at the gene body as the measure of gene expression to study the p63-regulated transcription in subsequent analyses.

MicroRNAs (miRNAs) have been shown to play an important role in epidermal differentiation ^26,27,28, and therefore we analyzed our RNA-seq and RNAPII ChIP-seq for expression of miRNAs. Examining the RNAPII occupancy at annotated miRNA genes using K-means clustering analysis showed dynamically regulated expression of 103 primary miRNA (pri-miRNA) transcripts during differentiation (Supplementary Fig S3A; Supplementary Table S4A) ²⁹. Some of these miRNAs have been reported to play a role in epidermal differentiation such as miR-205 and miRNA-203 (Supplementary Fig S3B) ^27,30,31. For most of these miRNAs, RNA-seq signals at miRNA regions are low and the correlation between RNAPII ChIP-seq and RNA-seq signals is poor (Supplementary Table S4C). Thus, RNAPII occupancy appears to give a better readout of pri-miRNA transcription.

p63-bound enhancers during epidermal differentiation

Next, we assessed the role of p63 in regulating dynamic gene expression during epidermal differentiation. We performed ChIP-seq of p63 from keratinocytes at different stages of differentiation. Among a total number of 38,980 p63 binding sites detected from all stages (Supplementary Table S5E), more than 90% were found to contain at least one p63 binding consensus motif using our previously established p63scan algorithm ¹⁸, and there is a significant enrichment of p63 binding in promoter regions and within 25 kb distance of genes compared to genomic distributions of all binding sites (hypergeometric test, P-value < 1E-323 and P-value = 6.97E-174, respectively) (Supplementary Fig S4A). Surprisingly, in contrast to the dynamic gene expression pattern observed in RNAPII ChIP-seq and RNA-seq analyses (Fig1C), the number of p63 binding sites and the p63 binding signals only showed a minor overall decrease during differentiation (Fig2A; Supplementary Fig S4B). This is probably due to the decreasing p63 protein levels during epidermal differentiation (Supplementary Fig S1B) ^16,23.

To resolve the seeming inconsistency between the rather static p63 binding and the dynamic gene expression during differentiation, we examined the epigenetic state of p63 binding sites using ENCODE data of primary human keratinocytes (NHEK) ³². Of all p63 binding sites, 50.5% are co-localized with the H3K27ac histone mark associated with active enhancers (Supplementary Table S6A). This suggests that approximately half of the identified p63-bound regions are active enhancer elements in human keratinocytes. H3K27me3 and H3K9me3, repressive histone modification marks associated with poised enhancers in embryonic stem cells and heterochromatin ^32,33, respectively, were not detected at any of the identified p63 binding sites (Fig2A), even though these marks are present at promoters of inactive genes (examples in Supplementary Fig S5A). As ENCODE data of NHEKs were collected at only one condition and provide limited information for the transcription program during epidermal differentiation, we further investigated the H3K27ac occupancy at all four stages during epidermal differentiation to gain insight into the activity dynamics of the p63-bound regions. Our data showed that 56.9% of all identified p63 binding sites were co-occupied by H3K27ac at any of the four stages (Supplementary Table S6B–F), consistent with the ENCODE data analysis (Fig2A; Supplementary Table S6A). The p63 binding sites with a detectable level of H3K27ac were considered as active p63 binding sites. When active p63 binding sites were mapped to adjacent genes (Supplementary Table S6G), a strong correlation was observed between the RNAPII occupancy at most genes and H3K27ac levels at the nearby p63 binding sites across the four stages of differentiation (median Pearson's R = 0.77) (Fig2B). In contrast, the p63 levels at these binding sites showed no tendency to correlate with nearby gene body RNAPII levels during differentiation (median Pearson's R = 0.09). These results indicate that H3K27ac occupancy at p63 binding sites, rather than p63 occupancy itself, reflects the regulatory activity of p63 binding sites in the transcription of nearby genes. Additionally, H3K27ac occupancy also correlates well with gene expression in general apart from at p63 binding sites (median Pearson's R = 0.72) (Fig2B), which is consistent with previous reports of the correlation of H3K27ac with gene transcription ³⁴. These findings suggest that p63-regulated gene expression can be predicted by dynamic H3K27ac occupancy at p63 binding sites.

To further investigate the dynamics of H3K27ac occupancy at p63 binding sites, we clustered p63 binding sites into seven groups based on the co-localized H3K27ac signals using K-means clustering (Fig3A; Supplementary Table S9A). Subsequently, we compared genes associated with these H3K27ac clusters (H3K27ac occupancy clusters in Fig3B) with those in the nine clusters based on gene body RNAPII occupancy (Fig1B, RNAPII occupancy clusters in Fig3B), to assess the correlation between the H3K27ac occupancy at the p63 binding sites and transcription of the nearest genes. We observed that RNAPII and H3K27ac clusters with similar dynamics shared more genes than expected, while those with divergent dynamics shared fewer genes than expected (Fig3B). For example, genes associated with H3K27ac cluster 2 with increased H3K27ac signal are enriched for genes present in RNAPII cluster 3 that show increased expression during late differentiation (e.g., KRT23) (Fig3B and C; Supplementary Tables S1A and S9A and B). In contrast, genes associated with H3K27ac cluster 2 are depleted for genes present in RNAPII clusters 1, 4, 5, and 7, all of which show decreased expression during differentiation (e.g., DSG3) (Fig3B and C; Supplementary Tables S1A and S9A and B). Genes associated with H3K27ac cluster 3 with first increased and then decreased H3K27ac signal are enriched for genes from RNAPII cluster 1 that show similar dynamics (e.g., SFN) (Fig3B and C; Supplementary Tables S1A and S9A and B).

Taken together, both general and per-cluster correlation analyses of RNAPII and H3K27ac occupancy indicate that H3K27ac signal at p63 binding sites rather than p63 binding signal alone is informative for the prediction of the transcriptional program regulated by p63 and that assessing H3K27ac occupancy at p63 binding sites could facilitate the identification of p63 target genes in keratinocytes.

Potential co-activators of p63

As p63 binding alone does not explain epidermal gene expression, we searched for potential p63 co-regulators by transcription factor motif enrichment analysis using position weight matrices (PWM) of known transcription factor family motifs (threshold, P-value = 0.01, 1.5-fold enrichment versus background) ³⁵. An enrichment of motifs from the ZFX, RUNX, and EBF1 families was detected in active p63 binding sites with the H3K27ac mark (Fig4A; Supplementary Tables S7 and S8A). Among transcription factors that can bind to these motifs, RUNX1 is expressed in proliferating and differentiating keratinocytes (Supplementary Fig S5B; Supplementary Table S7). ChIP-qPCR analysis showed that RUNX1 was indeed bound to a number of tested p63 binding sites with the RUNX motif (Fig4B; Supplementary Table S8B), indicating that RUNX1 can cooperate with p63 in epidermal cells. For “inactive” p63 binding sites that were not co-localized with the H3K27ac mark, an enrichment of the NR2 nuclear receptor family motif was found (Fig4A; Supplementary Table S7), and among transcription factors that can bind to this motif, RXRA is expressed in keratinocytes. Some motifs were enriched in both active and inactive p63 binding sites. For example, a strong enrichment of the motif sequence bound by the p53 family proteins including p63 was detected (Supplementary Table S7). This is consistent with the high percentage of p63 motif-containing p63 binding sites detected by p63scan (90%). Other motifs that were enriched in both active and inactive binding sites include the motif of TFAP2α that has been reported to be a co-regulator of p63 (Supplementary Table S7). TFAP2α binding to p63-bound regions in keratinocytes was confirmed by ChIP-qPCR (Fig4B; Supplementary Table S8C).

To search for potential p63 co-regulators during specific stages of epidermal stratification, we determined motif enrichment in the seven clusters of p63 binding sites that show differential H3K27ac occupancy patterns during differentiation. Several transcription factor family motifs were enriched in specific clusters, for example, FOXQ1 in H3K27ac cluster 1, MAFB and SP100 in cluster 2, E2F in cluster 5, and VDR in cluster 6 (Fig4C; Supplementary Table S9C). Although the precise function of these transcription factors in epidermal stratification is not yet clear, it is likely that they cooperate with p63 in keratinocyte differentiation as they have been shown to regulate cell proliferation and differentiation ^{36,37,38,39,40}.

Clustered p63 binding sites (cl-p63BS) and clustered epidermal enhancers

Clustered enhancers, also referred to as super-enhancers, are defined as stretches of active enhancers that are bound by high levels of lineage-specific master transcription factors and co-activators that drive gene expression to define cell lineage identity ^41,42. The essential role of p63 for the commitment to the epithelial lineage has been well established, and it is likely that p63 binds to clustered enhancers. Therefore, we investigated whether p63 can bind to clustered enhancers for its function in epithelial lineage determination. We observed an enrichment of clusters of five or more p63 binding sites (hypergeometric test, P-value = 9.06E-13) within 100 kb of the transcription start site of transcribed genes (RNAPII occupancy, RPKM > 1) relative to all genes (Supplementary Table S10). We subsequently performed analyses using an established algorithm for detecting clustered enhancers ⁴² on two of our datasets: on p63 ChIP-seq data to identify regions with clustered p63 binding sites with high p63 binding signals, referred to as clustered p63 binding site regions (cl-p63BS) (Supplementary Table S14A), and on H3K27ac ChIP-seq data to identify regions with clustered H3K27ac-marked enhancers. Regions with clustered H3K27ac-marked enhancers have been reported to drive cell-type-specific gene expression ^42,43 and are therefore referred to as clustered epidermal-specific enhancers (cl-epidermal enhancers) hereafter (Supplementary Table S11A). Our analyses identified 1,293 cl-epidermal enhancers (Fig5D and E) and 1,991 regions as cl-p63BS (Fig5E). Of all 1,293 identified cl-epidermal enhancers, 1,182 (91.4%) sites overlapped with at least one p63 binding site (Fig5A and D). In addition, p63 preferentially binds to cl-epidermal enhancers (clustered H3K27ac-marked enhancers) rather than to single active enhancers (single H3K27ac peaks) and random controls (hypergeometric test and normal distribution probability test, respectively, P-value < 1E-323 in all cases) (Fig5A). Three p63 binding sites on average were identified within each cl-epidermal enhancer, whereas randomized control sets contained zero p63 binding sites on average (Fig5B). These observations show that p63 binding sites are enriched in cl-epidermal enhancers in epidermal keratinocytes, and underscore the importance of p63 as a key transcription factor in epithelial lineage commitment. Out of a total number of 1,991 cl-p63BS, 617 share genomic regions with 589 cl-epidermal enhancers out of a total number of 1,293 (Fig5E) and are hereafter considered as active epidermal cl-p63BS (Supplementary Table S15A and B). Interestingly, approximately two-thirds of the cl-p63BS do not overlap with active cl-epidermal enhancers. Next, we used K-means clustering to identify six gene clusters associated with cl-epidermal enhancers showing different patterns of H3K27ac occupancy during differentiation (Fig5C; Supplementary Table S11A). These clusters were enriched for genes with concordant RNAPII occupancy patterns during differentiation and depleted for genes with discordant RNAPII differentiation dynamics (Fig5C). Compared to the association of single p63-bound enhancers with gene expression (Fig3B), these p63 binding site-containing cl-epidermal enhancers showed higher fold enrichment of gene clusters with similar expression patterns (Fig5C). In summary, our data showed that p63 binds to clustered enhancers, consistent with the role of p63 in epithelial commitment, and that these clustered enhancers play roles not only in epithelial commitment but also during epidermal differentiation.

p63-regulated genes and their relevance to disease

Knockdown or knockout of p63 in epithelial cells to identify de-regulated genes ^14,17 has been the major approach in previous p63 target gene studies. The main drawback of this approach is that loss of p63 results in epithelial differentiation defects, and genes that are expressed at a later stage of differentiation and regulated by p63 may not be identified. In this study, we took a different approach and used differentiating normal keratinocytes. To identify p63 target genes, we first associated all identified p63 binding sites (38,980) with their nearest gene (11,890) (Supplementary Table S12A) and analyzed gene function using DAVID Gene Ontology annotation ²⁴. This showed an enrichment of genes involved in intracellular signaling cascade, phosphorylation, vascular development, cell migration, and cell death (Supplementary Table S12B). Altogether, 10,020 p63 binding sites were identified in the proximity of differentially regulated genes (RNAPII ChIP-seq clusters, RPKM > 1) during epidermal differentiation and 74.2% of these binding sites are co-localized with H3K27ac (Supplementary Table S1D). These differentially regulated genes with active p63 binding sites nearby were considered potential p63 target genes (1,961 based on RNAseq and 2,760 based on RNAPII) (Supplementary Tables S1D and S2D). For example, p63 binding sites with increasing levels of the H3K27ac mark during differentiation were observed near the IVL, LCE1B, and LCE1C genes that are highly expressed in late differentiating keratinocytes (Supplementary Fig S5C). Several of these potential p63 target genes are known to be causative genes involved in epidermal diseases, such as CDH3, DSG1, KRT5, KRT10, and FLG (Supplementary Fig S5C; Supplementary Table S13A and B) ^{10,11,44,45,46}. Of note, many of these potential p63 target genes such as KRT1 and PRDM1 were not previously considered as direct p63 target genes (Supplementary Fig S5D; Supplementary Table S13A–D). Furthermore, we compared the potential p63 targets identified in this study with previous p63 RNAi knockdown experiments (Supplementary Table S13C) ^14,17. Among the 120 p63-regulated genes previously reported in two p63 RNAi knockdown experiments ^14,17, we found that 70 of them (57.5%) are among potential direct p63 targets in this study (Supplementary Table S13D). Interestingly, many co-regulators of p63 including RUNX1 ⁴⁷, VDR ⁴⁸, RXRA ⁴⁹, and FOXQ1 in epidermal keratinocytes are either known p63 targets genes or identified as potential p63 targets in this study (Fig4A and C). Taken together, our data show that p63 targets are expressed in different layers of the epidermis and associated with distinct disease phenotypes (Supplementary Fig S2B).

Intriguingly, p63 binding sites were also identified in the vicinity of 80 pri-miRNAs of the 130 pri-miRNAs that are differentially regulated during keratinocyte differentiation (Supplementary Fig S3A; Supplementary Table S4B). Among these pri-miRNAs, 74 of them contain at least one p63 binding site with H3K27ac occupancy, including pri-miR-203 and pri-miR-205 (Supplementary Fig S3B; Supplementary Table S4B), which suggests that these miRNAs are targets of p63. Among these p63-regulated miRNA genes, miR-205 is a known p63 target ^50,51. Another potential p63-regulated miRNA is miR-203 that is known to regulate p63 at the mRNA level during keratinocyte differentiation ³⁰, and thus, our data suggest a regulatory feedback loop of p63 and miRNA-203 during differentiation.

In comparison with genes mapped to single p63 binding sites, genes mapped to the 1,991 cl-p63BS (clustered p63 binding sites) (Fig5E) showed a stronger enrichment of genes involved in ectoderm and epithelium development (Supplementary Table S14A–D). This may be due to genes in the vicinity of cl-p63BS being more likely to be true p63 target genes compared with those near single p63 binding sites. Genes associated with nearby cl-epidermal enhancers (clustered epidermal enhancers marked by H3K27ac) (Fig5D and E) are enriched for genes with epidermal-related functions (Supplementary Table S11C). Genes mapped to these active epidermal cl-p63BS (overlap of cl-p63BS and cl-epidermal enhancers) likely play a role in epidermal development (Supplementary Table S15C), similar to genes mapped to cl-epidermal enhancers (Supplementary Table S11C). Human Phenotype Ontology (HPO) analysis showed that this group of genes including CLDN1, p63, and KRT17 are associated with skin-related abnormalities such as epidermal thickening, aplasia/hypoplasia of the nails, and hyperkeratosis (Fig5E; Supplementary Table S15D). Genes in the vicinity of cl-epidermal enhancers that do not overlap cl-p63BS were enriched with genes that play a role in epidermal development (Supplementary Table S16C), similar to the genes that were mapped by all cl-epidermal enhancers (Supplementary Table S11C). In contrast, the 69% of cl-p63BS that do not overlap cl-epidermal enhancers mapped to genes with broader functions including regulation of cell proliferation, tube development, vascular development, and regulation of apoptosis (Supplementary Table S17D). HPO analysis revealed that these genes are associated with epithelial-related non-skin phenotypes such as limb and tooth defects (Fig5E; Supplementary Table S17D), suggesting that inactive enhancers regulate gene expression during the development of other epithelial-related structures.

Ethics statement

All procedures regarding establishing human primary keratinocytes were approved by the ethics committee of the Radboud University Nijmegen Medical Centre (“Commissie Mensgebonden Onderzoek Arnhem-Nijmegen”). Informed consent was obtained.

Human primary keratinocyte culture

Skin biopsies were taken from the abdomen of two healthy volunteers to set up the primary keratinocyte culture, referred to as HKC1 and HKC2 ⁶³. Keratinocyte cultures in keratinocyte growth medium (KGM) under undifferentiated condition were previously described ⁶⁴, and differentiation of keratinocyte cultures was induced by cell contact inhibition and excluding several growth factor supplements: bovine pituitary extract (Bio Whittaker), EGF (Sigma), insulin (Sigma), and hydrocortisone (Calbiochem) from the medium. The medium was changed every second day, and before harvesting of the RNA and chromatin. Cells were collected on day 0, and on days 2, 4 and 7 after induction of differentiation.

ChIP and ChIP-seq

Human primary keratinocytes under proliferating, early, mid-, and late differentiation condition generated from HKC1 were used for ChIP and ChIP-seq analysis. ChIP-sequencing datasets were generated with an antibody recognizing the alpha isoforms of p63 (H129, Santa Cruz), an antibody that detects all forms of the large subunit of RNA polymerase II (8WG16, Santa Cruz), and a H3K27ac antibody recognizing the acetylation of the lysine 27 residue of histone H3 (H3K27ac, C15410174, pAb-174-050, Diagenode). ChIP-qPCR experiments were performed using 2 μg RUNX1 antibody (ab23980, Abcam) and 1 μg TFAP2 antibody (C18X, Santa cruz). RNAPII, p63, H3K27ac, RUNX1, and TFAP2 ChIP experiments were performed as previously described ⁶⁵, with a minor change of using magnetic beads (Novex by Life Technologies, 10008D and 10009D).

RNA extraction and RNA-seq

Human primary keratinocytes under proliferating, early, mid-, and late differentiation condition were used for RNA-seq analysis. RNA extraction was carried out using the NucleoSpin RNA II kit (740,955.50; MACHEREY-NAGEL, Düren, Germany). cDNA synthesis from 1 μg total RNA as control for RNA-seq was carried out using the iScript cDNA synthesis kit (Bio-Rad; 170-8891, CA, Hercules, USA). Starting amount of RNA-seq sample was 1 μg, and rRNA was depleted using the Ribo-Zero rRNA Removal Kit (Human/Mouse/Rat; MRZH11124, Epicentre, Madison, WI, USA) according to the manufacturer's instructions. RNA fragmentation reactions were performed using fragmentation buffer (5×; 200 mM Tris-Ac, 500 mM potassium-Ac, 150 mM magnesium-Ac) in a final concentration of 1× per reaction. Each 50-μl fragmentation reaction was incubated at 95°C for 1.5 min on a thermal cycler and placed on ice for 10 min. Ethanol precipitation was used to purify the reactions. The fragmented rRNA-depleted RNA was combined with 5 μg of random hexamers (Roche) in a final volume of 11 μl, incubated at 65°C for 5 min, and afterwards immediately placed on ice. The remaining reagents were added to the reaction: 4 μl 5× first-strand buffer (Invitrogen), 2 μl DTT (0.1 M Invitrogen), 1 μl dNTP mix (10 mM Invitrogen), 2 μl actinomycin D (0.1 μg/μl Sigma), 1 μl of RNase H (40 U/ml Ambion), and 1 μl SuperScript III (200 U Invitrogen). The first-strand reaction was incubated at 25°C for 10 min, 50°C for 90 min, followed by deactivation at 70°C for 15 min, followed by MinElute Reaction Cleanup Kit (Qiagen), according to the manufacturer's protocol. The second strand was synthesized by adding 20 μl 5× second-strand buffer (Invitrogen), 4 μl 5× first-strand buffer (Invitrogen), 2 μl DTT (0.1 M Invitrogen), 1 μl random hexamers (5 mg/ml Roche), dUTP mix (12.5 mM Invitrogen), 1 μl RNase H (8 U/ml Ambion), 1 μl E. coli DNA polymerase I (10 U/μl Invitrogen), and 1 μl E. coli DNA ligase (10 U/μl NEB) to the purified sample (100 μl total volume). After 2 h at 16°C, 1 μl of T4 DNA polymerase (10 U/μl Promega) was added and the reaction was incubated at 16°C for 10 min. The ds-cDNA was purified using the MinElute Reaction Cleanup Kit (Qiagen), according to the manufacturer's protocol.

Illumina library preparation

DNA samples were prepared for sequencing by end repair of 6 ng total DNA (H129 ChIP-seq and RNAPII ChIP-seq samples), 10 ng total DNA (H3K27ac ChIP-seq samples), and 5 ng total cDNA (RNA-seq samples) as measured by Qubit (Invitrogen). Adaptors were ligated to the DNA fragments, followed by a pre-PCR of 4 cycles, size selection (∼300 bp), and subsequently 11 cycles of PCR amplification. Cluster generation and sequencing (50 bp) was performed with the hiSeq Solexa Genome analyzer according to standard Illumina protocols.

Quantitative real-time (qRT–PCR)

Quantitative PCR primers were designed using Primer3 (http://frodo.wi.mit.edu), and qPCRs were performed in the CFX96 Real-Time system (C1000 Touch-Thermal Cycler, Bio-Rad) by using GoTaq qPCR mix (Promega) according to the manufacturer's protocol. For qPCR of cDNA analysis, human acidic ribosomal protein (hARP) was used as a housekeeping gene to normalize the amount of cDNA. Differences in the expression of each gene during differentiation (relative expression) were calculated by 2^ΔΔCt method ⁶⁶. For qPCR of ChIP analysis, one primer set was used for each tested binding region (Supplementary Table S18) and ChIP efficiency of each binding site was calculated by comparison of the percentage of ChIPped DNA against input chromatin to this percentage at a negative control region. Validation experiments were performed by RT–qPCR with primers (Biolegio BV, Nijmegen, the Netherlands) as shown in Supplementary Table S18.

ChIP-seq data analysis

All 50-bp sequence reads were uniquely mapped to the human genome NCBI build 37 (hg19) using bwa 0.6.1 with standard parameters ⁶⁷ (Supplementary Table S19). All p63, H3K27ac, and RNAPII ChIP-seq datasets were normalized to the same sequencing depth by randomly removing aligned reads. Duplicated reads were removed before normalization. Peak recognition for p63 ChIP-seq and RNAPII ChIP-seq datasets was performed using MACS2, an updated version of MACS ⁶⁸ that is specifically designed to process mixed signal types (https://github.com/taoliu/MACS) with default settings and a P-value threshold of 1E-9 using a genomic DNA as background control, and results are summarized in Supplementary Table S19. Peaks were mapped to RefSeq genes, downloaded from UCSC Genome Browser (hg19), to determine genomic location of the p63 binding sites at promoter regions (defined as 5 kb upstream of the transcription start site and the end of the first intron of genes), intragenic, intergenic, or 25 kb away from the transcription start site of the gene. Peak recognition for the H3K27ac ChIP-seq dataset was performed using MACS2 using default settings for broad peak calling and a P-value threshold of 1E-9 using a genomic DNA as background control, and results are summarized in Supplementary Table S19. The identified p63 peaks (38,980) that overlap with the identified H3K27ac peaks by MACS2 are considered active p63 binding sites (22,190), while those that do not overlap the MACS2-detected H3K27ac peaks are considered inactive p63 binding sites (16,790).

RNA-seq data processing

All 50-bp single-end RNA-seq reads from two keratinocyte cell lines (biological replicates) were mapped to the human genome NCBI build 37 (hg19) using gsnap version 2013-03-31 ⁶⁹. The generated data per condition are summarized in Supplementary Table S19. Spliced reads were mapped using the Ensembl 68 reference transcriptome ⁷⁰. Transcript quantification and differential expression analysis were done with cufflinks and cuffdiff v2.1.1 ⁷¹, using the RefSeq transcriptome ⁷², downloaded from the UCSC genome browser website ⁷³ in June 2013.

RNA-seq clustering

RNA-seq RefSeq genes were filtered for those determined to be differentially expressed in at least one comparison by cufflinks/cuffdiff (P-value < 0.05, beta negative binomial model). These were filtered for genes with RPKM expression level of at least 10 in at least one condition, resulting in a set of 3,639 distinct genes, and 3,813 when considering different isoforms. RPKM of gene body from RNA-seq was clustered into 6 clusters using the PAM clustering algorithm in R (R core team (2013) R: A language and environment for statistical computing. R foundation for Statistical Computing. Vienna, Austria, http://www.R-project.org/, Bioconductor ⁷⁴) using the z-scores calculated from the log10 (RPKM) values. GO analysis was performed with DAVID (http://david.abcc.ncifcrf.gov/).

RNAPII and H3K27ac gene body occupancy quantification

The occupancy of RNAPII and H3K27ac of the ChIP-seq datasets from HKC1 at Refseq gene transcripts (hg19) was quantified by the detection of the sequenced reads, and the reads per kilobase per million mapped reads (RPKM) was calculated for each gene.

RNAPII ChIP-seq clustering

RNAPII ChIP-seq RefSeq genes were filtered for genes with RPKM expression level of at least 1 in at least one condition, resulting in a set of 7,223 distinct genes, and 8,741 when considering different isoforms. These genes were clustered into 9 clusters using the PAM clustering algorithm in R ⁷⁴ using the z-scores calculated from the log₂ (RPKM) values.

RNAPII occupancy versus RNA steady-state comparison

Rank-based comparisons of expression levels were used to compare RNAPII occupancy and RNA steady state due to the different distributions of RNAPII occupancy and RNA-seq expression levels. Genes were ranked separately according to RNAPII ChIP-seq occupancy and RNA-seq RPKM. For each gene, the difference between the RNA-seq and ChIP-seq ranks was taken, and the genes were ranked according to this difference. Subsequently, the enrichment of non-coding RNAs with high RNAPII occupancy relative to RNA-seq levels was determined using mean-rank gene set enrichment analysis ⁷⁵ as implemented in the “limma” R package. Overall correlations between the RNA-seq and RNAPII levels of different samples were calculated using the Spearman rank correlation coefficient, due to the fact that the genomewide gene expression profiles per sample do not follow a normal distribution.

Correlation between enhancer H3K27ac, p63 signal, and gene body RNAPII ChIP-seq signal

For all expressed genes mapped to nearby enhancers, the Pearson correlation coefficient was calculated between the H3K27ac and p63 ChIP-seq read counts at the associated enhancers (normalized for per-sample total read count) and the gene body RNAPII RPKM score across the four differentiation stages. H3K27ac enhancer signal was defined as a 4-kb window around the identified p63 ChIP-seq peak center to accommodate its broader binding pattern. When multiple enhancers mapped to a single gene, their read counts were summed.

Clustered enhancers calling

Clustered epidermal enhancers were called using the ROSE tool, version 0.1 ⁴². The identified H3K27ac peaks were used to call the cl-epidermal enhancers using the same settings as described in Hnisz et al, ⁴², with an enhancer stitching distance of 12-kb and a 4-kb exclusion window around promoters. H3K27ac (cl-epidermal) enhancer differentiation profiles were clustered into 6 clusters using PAM clustering, applied to the z-scores of the normalized read counts.

p63-bound enriched regions were determined using the same approach and settings as the detection of cl-epidermal enhancers, except that the identified p63 peaks from p63 ChIP-seq files were used. Additionally, the 4-kb exclusion window around promoters was not used. Enriched regions consisting of only one single p63 peak were removed from the set.

Association between H3K27ac signal dynamics at single enhancers and gene body RNAPII dynamics

The H3K27ac signal at co-localized single p63 binding sites was defined by the calculation of RPKM using a 4-kb window around the identified p63 ChIP-seq peak center from the H3K27ac ChIP-seq datasets of all four differentiation stages of HKC1. To analyze H3K27ac dynamics, z-scores of log₂ (RPKM) were determined and clustered into 7 clusters by K-means clustering. These enhancers were mapped to the closest gene (human genome NCBI build 37, hg19). The detected RNAPII gene body occupancy (RPKM > 1) and clustering were defined as previously described from the RNAPII ChIP-seq datasets of all four differentiation stages of HKC1. Over- and underrepresentation of shared genes associated with p63 binding sites with H3K27ac occupancy dynamics (for all 7 clusters) were compared with the gene clusters defined by RNAPII occupancy dynamics (for all 9 clusters). Significant fold enrichment or depletion was defined using hypergeometric distribution analysis (P-value < 0.001) combined with a fold change higher than 1.5 or smaller than 0.5, respectively.

Association between H3K27ac signal dynamics at cl-epidermal enhancers and gene body RNAPII dynamics

The H3K27ac signal at cl-epidermal enhancers was defined by summing the reads falling within each cl-epidermal enhancer region from the count-normalized H3K27ac ChIP-seq datasets of all four differentiation stages of HKC1. The BAM files containing the H3K27ac ChIP-seq reads were normalized to the same number of reads per dataset by randomly removing excess reads from the larger BAM files. To analyze H3K27ac dynamics, the z-scores of the normalized counts across the four differentiation stages were calculated and clustered into 6 clusters using K-means clustering. These cl-epidermal enhancers were mapped to the closest gene or any overlapping genes (human genome NCBI build 37, hg19). The detected RNAPII gene body occupancy (RPKM > 1) and clustering were defined as previously described from the RNAPII ChIP-seq datasets of all four differentiation stages of HKC1. Over- and underrepresentation of shared genes associated with cl-epidermal enhancer H3K27ac occupancy dynamics (for all 6 clusters) were compared with the gene clusters defined by RNAPII occupancy dynamics (for all 9 clusters). Significant fold enrichment or depletion was defined using hypergeometric distribution analysis (P-value < 0.001) combined with a fold change higher than 1.5 or smaller than 0.5, respectively.

Enrichment of p63 binding sites in clustered epidermal enhancers

Enrichment of p63 binding sites in clustered epidermal enhancers was determined by randomly shuffling the clustered enhancers across the genome using the “bedtools shuffle” command from the BEDTools suite ⁷⁶, excluding unsequenced genomic gaps larger than 10 kb. Overlap with p63 binding sites was determined using the “bedtools intersect” command. The P-value of p63 binding site enrichment in clustered enhancers was determined by fitting a normal distribution to the overlap percentages for the shuffled regions (which were normally distributed, Shapiro–Wilks P-value = 0.52) and determining the area under the cumulative probability curve to the right of the actual overlap percentage.

Detected transcribed miRNAs by RNAPII occupancy

The occupancy of RNAPII ChIP-seq datasets at the transcribed miRNAs was quantified by RPKM calculation of the detected sequenced reads for each identified RNAPII peak. A cutoff of RPKM > 1 was set, and only miRNAs, smaller than 200 bp, intersecting these RNAPII peak regions were considered.

Motif enrichment analysis

Motif enrichment was determined using Gimme scan of the GimmeMotifs package ⁷⁷, which was used to scan the p63BS datasets for motif enrichment with a motif score of at least 99% of known family motifs ¹⁶, including the RUNX family motifs. The p63BS datasets were centered on the peak summit of the identified p63 peaks (Supplementary Table S5) and extended 300 bp to both directions. A background file of 10,000 genomically matched random peaks was generated using Gimme background, of which a validation set (sequences not used for motif prediction) was generated using Gimme threshold of all known family motifs with a false discovery rate P-value of 0.01 (considering a motif score of at least 99%). Motif enrichment in active versus inactive p63 binding site and motif enrichment for the different p63 binding site clusters with H3K27ac occupancy was determined. For each set of p63BS sequences, a background set was generated, and the validation set was used as a cutoff file to determine enrichment in the family motifs of all p63BS datasets. We used 1.5-fold as a cutoff to determine enrichment of motifs present in the datasets versus matched genomic background datasets, and between the different datasets.

ENCODE data

H3K27ac ChIP-seq (GSM733674), H3K4me1 ChIP-seq (GSM733698), H3K27me3 ChIP-seq (GSM733701), H3K9me1 ChIP-seq (GSM733655), H3K9me3 ChIP-seq (GSM1003528), and Open chromatin by DNaseI HS (GSM816635) data from NHEK was downloaded from the ENCODE UCSC browser (http://genome.ucsc.edu/ENCODE/).

Accession codes

All RNA-seq and ChIP-seq data have been deposited in the Gene Expression Omnibus ⁷⁸ (GEO GSE59827).