- Open Access
p63 cooperates with CTCF to modulate chromatin architecture in skin keratinocytes
Epigenetics & Chromatinvolume 12, Article number: 31 (2019)
The transcription factor p63 regulates epidermal genes and the enhancer landscape in skin keratinocytes. Its molecular function in controlling the chromatin structure is, however, not yet completely understood. Here, we integrated multi-omics profiles, including the transcriptome, transcription factor DNA-binding and chromatin accessibility, in skin keratinocytes isolated from EEC syndrome patients carrying p63 mutations, to examine the role of p63 in shaping the chromatin architecture. We found decreased chromatin accessibility in p63- and CTCF-bound open chromatin regions that potentially contributed to gene deregulation in mutant keratinocytes. Cooperation of p63 and CTCF seemed to assist chromatin interactions between p63-bound enhancers and gene promoters in skin keratinocytes. Our study suggests an intriguing model where cell type-specific transcription factors such as p63 cooperate with the genome organizer CTCF in the three-dimensional chromatin space to regulate the transcription program important for the proper cell identity.
Skin development and homeostasis requires tightly regulated epidermal keratinocyte proliferation, differentiation and apoptosis, and these processes are governed by cooperation of cis-regulatory elements, transcription factors (TFs), chromatin accessibility as well as higher-order chromatin organization [1,2,3,4,5]. The TF p63 is a key regulator of epidermal development. At the molecular level, p63 regulates a large number of genes, for example, p21 (CDKN1A) in cell cycle arrest , Fras1 in maintaining basement membrane integrity  and key genes such as keratins, filaggrin and loricrin required for epidermal morphogenesis and differentiation [1, 8]. Furthermore, p63 directly regulates chromatin factors: Satb1, Lsh and Brg1, to control the chromatin remodeling during epidermis development [9,10,11]. Recent studies showed that p63 exerts a crucial role in establishing the enhancer landscape [4, 5, 12, 13]. Through active enhancers, p63 cooperates with its co-regulating TFs to modulate transcriptional program required for epidermal homeostasis .
In human, heterozygous mutations of TP63 encoding p63 cause a spectrum of ectoderm-related disorders . For example, ectrodactyly–ectodermal dysplasia–cleft lip/palate (EEC) syndrome is caused by point mutations located in the p63 DNA-binding domain and manifests ectodermal dysplasia with defects in the epidermis and epidermal-related appendages, limb malformation and cleft lip/palate. Five hotspot mutations affecting amino acids, R204, R227, R279, R280 and R304, cover approximately 90% of the all EEC syndrome cases . Our previous study showed that mutant p63 resulted in a genome-wide redistribution of enhancers in keratinocytes established from EEC patients . Consistently, the gene network analysis identified a significant co-expression gene module of ‘nucleosome assembly’, implying a less-organized chromatin structure in EEC syndrome keratinocytes . How mutant p63 affects the chromatin structure is, however, not yet clear.
In this study, we characterized the chromatin accessibility using ATAC-seq in keratinocytes established from EEC patients carrying p63 mutations, in comparison with control keratinocytes. A clear difference in chromatin accessibility that correlated with the transcriptional dynamics was detected. Unexpectedly, strong enrichment of CTCF binding sites was observed at control-specific open chromatin regions. By combining published promoter capture Hi-C seq data, we found that CTCF and p63 were cooperatively involved in DNA loops to regulate epidermal genes. Our findings provide new insights into the coordinated regulatory role of CTCF and p63 in chromatin interactions in epidermal keratinocytes.
Differential chromatin accessibility between control and p63 mutant keratinocytes
We performed assay for transposase-accessible chromatin with high-throughput sequencing (ATAC-seq) to characterize the accessible genome as well as the nucleosome position in both control and p63 mutant keratinocytes at the proliferation stage. Two replicas of ATAC-seq analyses showed high correlation (Fig. 1a), indicating high reproducibility. The principle component analysis (PCA) plot based on the chromatin accessibility displayed a clear separation between control and p63 mutant keratinocytes (Fig. 1b), which is highly consistent with results from gene expression profiles (Fig. 1b). The two p63 mutant lines showed high similarity in both genome accessibility and transcription level, as they are close on the PC1 axis that represents major variations. As the goal of this work is to examine the difference between control and EEC mutant lines, these two mutant lines were considered as one group, termed as p63 mutant keratinocytes, in the following analyses.
Subsequently, we identified the differential open chromatin regions (OCRs) marked by ATAC-seq signals between control and p63 mutant keratinocytes. In total, there were 2492 open chromatin regions that showed higher signal in control keratinocytes, termed as control-specific open chromatin regions (Ctr-OCRs); in parallel, there were 3716 regions that showed higher signal in mutant keratinocytes, termed as mutant-specific open chromatin regions (Mt-OCRs) (Fig. 1c). As expected, we found that differential OCRs were positively correlated with gene expression, when they were assigned to the nearest differentially expressed (DE) genes (Fig. 1c). Genes associated with Ctr-OCRs were mainly involved in the regulation of ‘cell cycle’ (Fig. 1d), e.g., CDC7 (Fig. 1e), while genes nearby Mt-OCRs were mainly enriched in ‘actin filament-based process’ (Fig. 1d), e.g., FAP (Fig. 1e).
p63 and CTCF occupancy at Ctr-OCRs
To identify the potential TFs involved in these differential OCRs, we first performed a comparative motif analysis between control and p63 mutant keratinocytes using HOMER. Among the Ctr-OCRs, we found that p63 motif is most enriched (Fig. 2a), which is consistent with our previous finding that loss of p63 binding resulted in decreased enhancers in p63 mutant keratinocytes . Unexpectedly, the motif of CTCF that is a well-known genome organizer  was identified as the second highly enriched motif in the Ctr-OCRs. Among the Mt-OCRs, AP-1 motif family was most enriched (Fig. 2b), consistent with findings on gained enhancers in p63 mutant keratinocytes .
To validate our motif analysis result and explore the role of p63 in controlling Ctr-OCRs, we used published p63 ChIP-seq data from our previous work . In control keratinocytes, a total of 1223 p63 binding sites (BSs) were covered by the 2492 Ctr-OCRs, showing a significant overlap (p < 0.001). When ATAC-seq and p63 ChIP-seq reads were mapped to the 1223 p63-bound Ctr-OCRs, we found a decrease in both chromatin accessibility and p63 binding signals at the p63-bound Ctr-OCRs in p63 mutant keratinocytes (Fig. 2c).
To confirm the role of CTCF in Ctr-OCRs, we performed CTCF ChIP-seq in both control and p63 mutant keratinocytes. In total, we found 12,116 CTCF BSs, among which only 363 were differential in CTCF binding between control and p63 mutant keratinocytes (p < 0.05), indicating that CTCF binding is rather stable. In addition, among all putative CTCF BSs, we found 376 overlapped with the 2492 Ctr-OCRs, which is statistically significant (p < 0.001). We then mapped ATAC-seq and CTCF ChIP-seq reads to the 376 CTCF-bound Ctr-OCRs. As expected, decrease in ATAC-seq signals was observed at the CTCF-bound Ctr-OCRs in p63 mutant keratinocytes. However, unexpectedly, there was no significant change in CTCF binding signals at the CTCF-bound Ctr-OCRs in p63 mutant keratinocytes, as compared to the control keratinocytes (Fig. 2d). Another unexpected finding was that there is few overlap (n = 435) out of 45,350 p63 BSs and 12,116 CTCF BSs (Additional file 1: Fig. S1A), although both motifs were identified in Ctr-OCRs. Consistently, there was no visible CTCF binding signal at all p63 BSs, and vice versa, no visible p63 binding at the CTCF BSs (Additional file 1: Fig. S1B). We also observed decreased p63 binding and largely unchanged CTCF binding signal in p63 mutant keratinocytes at the co-localized binding sites (Additional file 1: Fig. S1C).
Subsequently, we analyzed the nucleosome positioning at p63-bound and CTCF-bound Ctr-OCRs by mapping the distribution of ATAC-seq fragments centered on either p63 motif or CTCF motif (Fig. 2e). In control keratinocytes, we found an enrichment of short fragments (< 150 bp) surrounding both p63 motif and CTCF motif. The enrichment of such short reads clearly decreased in p63 mutant keratinocytes (Fig. 2e and Additional file 2: Fig. S2). We also observed an enrichment of fragments that had lengths slightly shorter than 200 bp which represent the stable flanking nucleosomes in control keratinocytes. However, in p63 mutant keratinocytes, the approximate 200 bp length fragments seemed to be less present (Additional file 2: Fig. S2), indicating an altered nucleosome organization.
Taken together, we observed decreased accessibility at p63-bound and CTCF-bound Ctr-OCRs in p63 mutant keratinocytes when compared to control keratinocytes. Notably, there was a clear decrease in p63 binding at the p63-bound Ctr-OCRs, whereas CTCF binding did not seem to be affected at the CTCF-bound Ctr-OCRs in p63 mutant keratinocytes.
Characterization of p63-bound OCRs and CTCF-bound OCRs
We reasoned that the chromatin states and the genomic locations of these Ctr-OCRs might be informative to understand the difference between p63-bound and CTCF-bound Ctr-OCRs. To this end, we firstly analyzed the chromatin states (CSs) of all CTCF and p63 BSs with the combination of six histone marks from Roadmap project  (Fig. 3a). We found that p63 bound to both promoters and enhancers marked by H3K27ac and H3K4me1, whereas CTCF was more enriched at promoters marked by H3K4me3 (Fig. 3a and Additional file 3: Fig. S3). We then performed analyses on the p63- and CTCF-bound Ctr-OCRs (Fig. 3b). When comparing to all Ctr-OCRs, p63-bound Ctr-OCRs were overrepresented at enhancers (CSs 8, 9 and 11), while the majority of CTCF-bound Ctr-OCRs showed overrepresentation at promoter-proximal regions (CSs 2 and 4, Fig. 3b), consistent with all p63 and CTCF BSs.
Next, to analyze the localization of p63-bound and CTCF-bound Ctr-OCRs in regulatory chromatin interactions, we utilized the promoter capture Hi-C (PCHiC) data and TAD (topologically associating domain) regions from published Hi-C data in keratinocytes . In the PCHiC, promoters were used as baits to capture the interacting genomic regions. In total, this PCHiC dataset captured 119,648 interactions, including 103,497 (86.5%) promoter–enhancer loops and 16,151 (13.5%) promoter–promoter loops (Fig. 3c). As p63 plays a more prominent role at enhancers (Fig. 3a) , we focused on promoter–enhancer loops. These loops had a median loop size of 250.9 kb (Fig. 3d). Among the previously defined 1223 p63-bound Ctr-OCRs (Fig. 2c), 388 of them were located at the anchors of 1674 loops, at either the promoter or the enhancer, and on average one p63-bound Ctr-OCR was connected by ~ 4.3 loops. Similarly, among the 376 CTCF-bound Ctr-OCRs (Fig. 2d), 191 of them were located at the anchors of 1093 loops; on average, one CTCF-bound Ctr-OCR was connected by ~ 5.7 loops. Therefore, CTCF-bound Ctr-OCRs at anchors were connected by more loops.
Furthermore, we found that CTCF-bound Ctr-OCRs located at loop anchors were more enriched for promoters (fourfold) (Fig. 3e), consistent with all CTCF BSs and all CTCF Ctr-OCRs. In contrast, p63-bound Ctr-OCRs associated with loops were enriched in both promoters and enhancers to a similar extent (Fig. 3e), which is different from all p63 BSs and all p63-bound Ctr-OCRs that were only enriched at enhancers (Fig. 3b). We also examined whether p63- and CTCF-bound Ctr-OCRs were localized at the TAD boundaries. As expected, CTCF-bound Ctr-OCRs showed ~ 2.5-fold enrichment for TAD boundaries, whereas p63-bound Ctr-OCRs were depleted of TAD boundaries (Fig. 3e). Randomization by reshuffling genomic regions with same sizes did not show any enrichment in TAD boundaries, promoters or enhancers. As CTCF-bound Ctr-OCRs were enriched at the promoters and p63-bound Ctr-OCRs were enriched at both promoters and enhancers, it is plausible that p63-bound enhancers interact with CTCF-bound promoters to regulate genes through DNA looping and a subset of these loops may be affected by p63 mutations.
p63 and CTCF mediate a subset of loops to regulate transcription
To assess how p63 and CTCF cooperates via DNA looping to regulate gene expression, we performed further analyses using previously published gene expression data  to examine the relation between gene expression and DNA looping. In control keratinocytes, we defined highly expressed genes by quartiles to four groups using their transcriptional levels (Q1-4, FPKM ≥ 1), and another group with lowly expressed genes (Q0, FPKM < 1). We detected an increasing number of associated loops with elevated gene expression levels (Fig. 4a). This suggests that highly expressed genes may be more indispensable of regulation through loops, which is line with previous study showing multi-loop activation hubs at key regulatory genes . Moreover, the differentially expressed genes between control and p63 mutant keratinocytes (Fig. 1b) were associated with significantly more loops (p < 0.001) than all expressed genes (Fig. 4b). We also performed the same analysis among loops with either p63-bound or CTCF-bound Ctr-OCRs. Similar results were found, although in these comparisons the difference was not significant, probably due to the low number of associated loops (Fig. 4c, d).
We then mapped 388 p63-bound and 191 CTCF-bound Ctr-OCRs which are located at loop anchors to their target genes by identifying their interacting gene promoters located at anchors of the other end of the loops. In some cases, Ctr-OCRs were located in promoters, and the associated genes were considered as target genes. In total, 388 p63-bound Ctr-OCRs interacted with 304 gene promoters, while 191 CTCF-bound Ctr-OCRs interacted with 500 gene promoters through looping interactions, indicating that on average a CTCF-bound Ctr-OCR is involved in the regulation of more genes than a p63-bound Ctr-OCR. There was a significant (p < 0.001) overlap (n = 96) of genes that were potentially regulated by p63- and CTCF-bound Ctr-OCRs (Fig. 4e). These 96 genes represent genes connected by the promoter–enhancer loops that were anchored by both p63-bound and CTCF-bound Ctr-OCRs and that may be affected in p63 mutant keratinocytes. As the three-dimensional chromatin landscape is relatively stable, we reasoned that the difference of the chromatin landscape in control and p63 mutant keratinocytes may influence gene expression during keratinocyte differentiation . Therefore, we examined the expression of these 96 genes using previously reported differential gene expression in control and p63 mutant keratinocytes during differentiation . Indeed, we found that 39 of these 96 genes showed deregulation (p < 0.001). Notably, many of the 96 genes were involved in ‘epidermal cell differentiation,’ e.g., KRT5, KRT75 and B9D1 (Fig. 4f).
Next, we dissected how p63-bound and CTCF-bound Ctr-OCRs at loop anchors co-regulated these 96 genes: whether p63-bound and CTCF-bound Ctr-OCRs were involved in the same or different promoter–enhancer loops. We found that 60 genes were regulated by promoter–enhancer loops where p63-bound and CTCF-bound Ctr-OCRs were both located at the anchors in the same loops, although they may co-localize in different manner (Fig. 4g). Among these 60 genes, 26 showed differential expression between control and p63 mutant keratinocytes during differentiation, with PIGV and KRT5 shown as examples (Fig. 4i, j). Among the loops associated with the 60 genes, the case where p63-bound Ctr-OCR was located at one anchor, while CTCF-bound Ctr-OCR was located at the other anchor of the loop counted for 12.6% (Fig. 4g, (1)). In the majority of the cases, both p63-bound and CTCF-bound Ctr-OCRs were present either at the promoter (70.6%) (Fig. 4g, (2)) or at the enhancer (16.8%) (Fig. 4g, (3)). Interestingly, for the rest 36 genes gene promoters were connected with both p63-bound and CTCF-bound Ctr-OCRs from different loops (Fig. 4h), suggesting that enhancers occupied by p63-bound and CTCF-bound Ctr-OCRs can co-regulate the same gene through different loops. Among these 36 genes, 13 were deregulated in p63 mutant keratinocytes during differentiation.
Over the years, molecular understanding of p63 function in keratinocytes has led to the identification of numerous target genes including chromatin factors [20,21,22] as well as of a solid role of p63 in orchestrating the enhancer landscape [4, 5, 23]. However, whether p63 is involved in the direct regulation of chromatin architecture is not known. In this study, we characterized the chromatin accessibility with ATAC-seq in both control and EEC syndrome patient keratinocytes carrying p63 mutations (R204W and R304W). Interestingly, we found decreased chromatin accessibility in p63- and CTCF-bound open chromatin regions in mutant p63 keratinocytes. These less accessible chromatin regions interact with epidermal gene promoters and potentially regulate gene expression. Collectively, we proposed an unreported role of p63 cooperating with CTCF in chromatin interactions.
It is envisaged that p63 is a key factor in skin keratinocytes [4, 5, 24]. Ample studies have shown that the proper regulation of gene expression by p63 through the precise control of enhancers is essential for maintaining the epidermal cell identity [4, 13, 25, 26]. It is equally important that p63 requires additional co-regulating TFs to regulate transcription. Several studies reported p63 co-regulating TFs, e.g., AP1, AP2, STAT5 and RFX5 [20, 22, 23, 27], which are essential for the epidermal gene expression in a temporal and spatial manner [4, 28]. Besides, p63 interacts with the chromatin remodeling factor BAF1 to maintain the open chromatin regions  and with Dnmt3a to locate enhancers . In this study, we found that CTCF, a well-known TF maintaining the chromatin architecture, is also a p63 co-regulating TF modulating a subset of chromatin loops. However, this cooperation is probably not through direction protein–protein interactions but mainly through bridging enhancers and promoters, given the few overlap observed between p63 and CTCF BSs. It is known that the majority of p63-bound regions are located in intergenic or intragenic regions (enhancers) rather than promoters , and therefore, how p63-bound enhancers interact with gene promoters is an interesting question. Our current study provides the first line of evidence that the cooperation with CTCF assists the looping of p63-bound enhancers to gene promoters.
The observation of the enriched CTCF motif in Ctr-OCRs detected in our ATAC-seq analyses was somewhat unexpected. Our previous enhancer-centered analysis that focused on H3K27ac enriched genomic regions did not capture CTCF motifs [4, 13], probably due to the fact that CTCF is more enriched in open chromatin regions depleted of H3K27ac but enriched for H3K4me3 signal that are often promoters (Fig. 3a and Additional file 3: Fig. S3). Interestingly, detecting the CTCF motif in Ctr-OCRs in keratinocytes is in line with the previous finding that the motif of YY1 was enriched in p63-bound sequences . YY1 has been shown to co-localize with CTCF and stabilize chromatin looping [29, 30].
CTCF maintains the genome organization by occupying the boundaries of megabase-scale TADs and functions as an insulator to block enhancer–promoter interactions between different TADs . Our data showed that, in addition to demarcating TADs, CTCF mediates promoter–enhancer loops, often located in promoter-proximal regions (Fig. 3e), to facilitate the promoter–enhancer interactions within one TAD. This is in line with the concept that a subpopulation of CTCF associates with the RNA polymerase II (Pol II) protein complex to activate transcription . It is likely that CTCF helps to bridge the p63-bound enhancers to transcription start site-proximal regulatory elements and to initiate transcription by interacting with Pol II, thus supporting a role of CTCF in facilitating contacts between transcription regulatory sequences. This model has been demonstrated by the previous work on the beta-globin locus . To get a clearer picture, technologies that are similar to ChIA-PET , such as HiChIP , would be useful to understand long-range contacts associated with CTCF or p63 in skin keratinocytes.
The three-dimensional chromatin landscape is relatively stable once established in a specific cell type, probably during cell commitment, and cell type-specific looping structures mainly control the accessibility of enhancers to their specific targets . Recent studies showed that cell type-specific TFs are involved in regulation of DNA looping in macrophage development, neural development as well as during cell reprogramming, which shed light on the role of cell type-specific TFs in regulating chromatin interactions [19, 36, 37]. Here, we proposed that a number of loci nearby epidermal genes were organized into a ‘regulatory chromatin hub’ within the chromatin interactions mediated by CTCF in epidermal keratinocytes (Fig. 5a). Such hubs contain multiple connecting DNA loops that require not only CTCF binding that is rather static but also binding of cell type-specific TFs (e.g., p63) for the transcriptional activity. In this model, cell type-specific TFs may be essential to make the DNA loops active in transcription. This hypothesis is consistent with our observation that unchanged CTCF binding signal but decreased accessibility in the CTCF-bound Ctr-OCRs in p63 mutant keratinocytes (Fig. 2d).
Moreover, our observations shed light on the chromatin-based mechanisms underlying EEC syndrome caused p63 mutations. These EEC mutations, such as R204W and R304W, were shown to disrupt p63 DNA-binding, resulting in impaired transactivation activity and loss of epidermal cell identity [13, 38, 39]. Our data suggest that deregulated function of DNA loops mediated by p63 and CTCF represents an additional layer to the disease mechanism (Fig. 5b). Given that CTCF binding did not change, it is likely that the looping interactions are static, but the looping is no longer active due to the absence of p63 binding (Fig. 5b (1)). Therefore, the epidermal genes were down-regulated and p63 mutant keratinocytes showed loss of epidermal cell identity . On the other hand, it is also possible that the loop ablation as a result of lost p63 binding leads to the decreased accessibility at CTCF-bound promoters and deregulated gene expression in p63 mutant keratinocytes (Fig. 5b (2)). The latter scenario suggests that p63 is the key factor determining the regulatory activity of such promoter–enhancer interactions. Nevertheless, how p63 mutations affect the cooperation between p63 and CTCF and gene expression is still unclear. Among the 96 genes with loops connected by p63- and CTCF-bound Ctr-OCRs, only 39 of them were deregulated in p63 mutant keratinocytes. It is possible that ‘shadow enhancers’ that are composed of clusters of enhancer and contribute to robust gene transcription [40, 41], are involved in regulation of the 57 genes that are not deregulated in mutant keratinocytes. Further experiments are required to test these hypotheses.
Taken together, our results highlight the cooperation of p63 and CTCF in mediating DNA loops as well as gene regulation. These observations provide new insights into the general regulatory role of chromatin looping which depends on the coordination of CTCF and cell type-specific TFs such as p63 in skin keratinocytes.
Materials and methods
All procedures for establishing and maintaining human primary keratinocytes were approved by the ethical committee of the Radboud university medical center (“Commissie Mensgebonden Onderzoek Arnhem-Nijmegen”). Informed consent was obtained from all donors of a skin biopsy.
Human primary keratinocyte culture
Primary keratinocytes were established previously from skin biopsies of three EEC syndrome patients carrying heterozygous mutations in the p63 DNA-binding domain, R204W  and R304W , as well as of non-EEC volunteers (Dombi23, referred to as control) . R204W and R304W are treated as a p63 mutant group for further analyses to minimize the effect from individual difference. As previously described , primary keratinocytes were cultured in Keratinocyte Basal Medium (KBM, Lonza #CC-4131) supplemented with 100 U/mL penicillin/streptomycin (Gibco Life Technology #15140122), 0.1 mM ethanolamine (Sigma-Aldrich #141-43-5), 0.1 mM O-phosphoethanolamine (Sigma-Aldrich #1071-23-4), 0.4% (vol/vol) bovine pituitary extract, 0.5 μg/mL hydrocortisone, 5 μg/mL insulin and 10 ng/mL epidermal growth factor (Lonza #CC-4131). Medium was refreshed every other day. When cells were more than 90% confluent, cells were collected at the proliferation stage. No mycoplasma contamination is found during cell culture.
ATAC libraries of the control and p63 mutant keratinocytes were prepared by a documented protocol . In brief, keratinocytes were treated with Accutase® solution on plates and well re-suspended into single cells. Pellet cells for 5 min at 500 g. After twice wash with ice-old PBS, pellet the cells again. Re-suspend cells and take out 100,000 cells for lysis with ice-old freshly made lysis buffer (containing 10 mM Tris–HCl pH 7.5, 10 mM NaCl, 3 mM MgCl2 and 0.1% IGEPAL CA-630 detergent). Then perform tagmentation using 2 μL of Tn5 transposase and 12.5 uL 2 × TD buffer (Illumina #FC-121-1031) at 37 °C for 1 h with 650 rpm shaking. The resulted DNA fragments underwent two sequential seven-cycle PCR amplifications, and in between the libraries were selected for < 500 bp fragments using SPRI beads. The final PCR products were purified with QIAquick PCR Purification Kit (QIAGEN #28106) and quantified with the KAPA Library Quantification Kit (Kapa Biosystems #KK4844), and then sequenced in a paired-ended manner using the NextSeq 500 (Illumina) according to standard Illumina protocols.
Chromatin for ChIP was prepared as previously described . ChIP assays were performed following a standard protocol  with minor modifications. On average, 0.5 M keratinocytes were used in each ChIP. 4x ChIP reactions are pooled to prepare one ChIP-seq sample. Antibodies against CTCF (Millipore #07-729, 5uL) were used in each ChIP assay. Resulted DNA fragments from four independent ChIP assays were purified with QIAquick PCR Purification Kit (QIAGEN #28106). Afterward, 5 ng DNA fragments were pooled and proceeded on with library construction using KAPA Hyper Prep Kit (Kapa Biosystems #KK8504) according to the standard protocol. The prepared libraries were then sequenced using the NextSeq 500 (Illumina) according to standard Illumina protocols.
RNA-seq data analysis
Paired-end RNA-seq reads were mapped to the human genome hg19 using STAR aligner  in two-pass mode with default parameters and enumerate stranded gene-level read counts at the same time. The generated count matrix was used as input for DESeq 2 package  to distinguish differential expressed genes between control and mutant keratinocytes. These genes greater than 1.5-fold changed at adjusted p value < 0.05 were considered significantly deregulated.
ChIP-seq and ATAC-seq data analysis
Sequenced reads were aligned against the UCSC hg19 human reference genome with Burrows–Wheeler Aligner (BWA) program with default parameters . The potential PCR and optical duplicates were removed using Picard MarkDuplicates option. The filtered BAM files were inputted to MACS2  for peak calling. The p63 ChIP-seq data were published and re-analyzed in this study . The p63, CTCF and H3K27ac were called using the narrow setting (default) with a q value of 0.01. For ATAC-seq, the open chromatin regions were predicted using default parameters except for using -f BAMPE option. Peaks overlapping with the consensus excludable ENCODE blacklist were dropped to avoid confounding by repetitive regions. All alignment files were extended to 200-bp and scaled to RPKM-normalized read coverage files using deepTools  for visualization. To compare binding profiles between different samples unbiasedly, we applied library size factors estimated from DESeq2  on RPKM values.
Differentially accessible regions were detected using DESeq 2 package  with fold change less than 2.0 and p value below 0.05. Differential motif analysis in the differential DHSs was employed by the findMotifs function in HOMER tool (http://homer.salk.edu/homer/motif/) with other default parameters, which can normalize the background sequences to remove GC-bias. The BEDtools suite (https://bedtools.readthedocs.io/en/latest/content/bedtools-suite.html) was used to test overlap and enrichment between different intervals.
Dimensionality reduction and functional annotation
For visualization, top 1000 variable genes were first selected based on interquartile range (IQR) of normalized expression values, and further used to reduce dimensionality (principal component analysis) of the dataset by pca function in R. Functional enrichment was evaluated by Metascape online tool , to gain insight into the biological functions for deregulated genes. Only these functional terms with Benjamini-adjusted p value < 0.05 were considered significantly overrepresented.
Capture Hi-C data processing
Raw promoter capture Hi-C sequencing reads from GSE84662 were processed using the HiCUP pipeline , like quality control, alignment to hg19 and reads filtering. Technical replicates were merged and de-duplicated, and then, pooled biological replicates were then jointly used for interaction identification with CHiCAGO  and the associated chicagoTools suite. All significant interactions were defined based on a strict interaction threshold (CHiCAGO score ≥ 5).
Briefly, CHiCAGO predicts interactions based on a convolution background model reflecting both ‘Brownian’ (real, but expected interactions) and ‘technical’ (assay and sequencing artifacts) components. The putative p values are corrected using a weighted false discovery control procedure that specifically accommodates the fact that increasingly larger numbers of tests are performed at regions where progressively smaller numbers of interactions are expected. The weights were learned based on the decrease in the reproducibility of interaction calls between the individual replicates of macrophage samples with distance. Interaction scores were then computed for each fragment pair as log-transformed, soft-thresholded, weighted p values.
Truong AB, Kretz M, Ridky TW, Kimmel R, Khavari PA. p63 regulates proliferation and differentiation of developmentally mature keratinocytes. Genes Dev. 2006;20(22):3185–97.
Lippens S, Denecker G, Ovaere P, Vandenabeele P, Declercq W. Death penalty for keratinocytes: apoptosis versus cornification. Cell Death Differ. 2005;12(S2):1497.
Botchkarev VA, Flores ER. p53/p63/p73 in the epidermis in health and disease. Cold Spring Harbor Perspect Med. 2014;4(8):a015248.
Kouwenhoven EN, Oti M, Niehues H, van Heeringen SJ, Schalkwijk J, Stunnenberg HG, van Bokhoven H, Zhou H. Transcription factor p63 bookmarks and regulates dynamic enhancers during epidermal differentiation. EMBO Rep. 2015;16(7):863–78.
Bao X, Rubin AJ, Qu K, Zhang J, Giresi PG, Chang HY, Khavari PA. A novel ATAC-seq approach reveals lineage-specific reinforcement of the open chromatin landscape via cooperation between BAF and p63. Genome Biol. 2015;16(1):284.
LeBoeuf M, Terrell A, Trivedi S, Sinha S, Epstein JA, Olson EN, Morrisey EE, Millar SE. Hdac1 and Hdac2 act redundantly to control p63 and p53 functions in epidermal progenitor cells. Dev Cell. 2010;19(6):807–18.
Koster MI, Dai D, Marinari B, Sano Y, Costanzo A, Karin M, Roop DR. p63 induces key target genes required for epidermal morphogenesis. Proc Natl Acad Sci. 2007;104(9):3255–60.
Romano R-A, Ortt K, Birkaya B, Smalley K, Sinha S. An active role of the ΔN isoform of p63 in regulating basal keratin genes K5 and K14 and directing epidermal cell fate. PLoS ONE. 2009;4(5):e5623.
Fessing MY, Mardaryev AN, Gdula MR, Sharov AA, Sharova TY, Rapisarda V, Gordon KB, Smorodchenko AD, Poterlowicz K, Ferone G, et al. p63 regulates Satb1 to control tissue-specific chromatin remodeling during development of the epidermis. J Cell Biol. 2011;194(6):825–39.
Mardaryev AN, Gdula MR, Yarker JL, Emelianov VU, Poterlowicz K, Sharov AA, Sharova TY, Scarpa JA, Joffe B, Solovei I, et al. p63 and Brg1 control developmentally regulated higher-order chromatin remodelling at the epidermal differentiation complex locus in epidermal progenitor cells. Development. 2014;141(1):101–11.
Keyes WM, Pecoraro M, Aranda V, Vernersson-Lindahl E, Li W, Vogel H, Guo X, Garcia EL, Michurina TV, Enikolopov G. ΔNp63α is an oncogene that targets chromatin remodeler Lsh to drive skin stem cell proliferation and tumorigenesis. Cell Stem Cell. 2011;8(2):164–76.
Rinaldi L, Datta D, Serrat J, Morey L, Solanas G, Avgustinova A, Blanco E, Pons JI, Matallanas D, Von Kriegsheim A, et al. Dnmt3a and Dnmt3b associate with enhancers to regulate human epidermal stem cell homeostasis. Cell Stem Cell. 2016;19(4):491–501.
Qu J, Tanis S, Smits J, Kouwenhoven EN, Oti M, Logie C, Stunnenberg H, Mulder K, Zhou H. Mutant p63 affects epidermal cell identity through rewiring the enhancer landscape. Cell Rep. 2018;25:3490–503.
Rinne T, Brunner HG, van Bokhoven H. p63-associated disorders. Cell Cycle. 2007;6(3):262–8.
Rinne T, Hamel B, Bokhoven HV, Brunner HG. Pattern of p63 mutations and their phenotypes—update. Am J Med Genet A. 2006;140(13):1396–406.
Ong C-T, Corces VG. CTCF: an architectural protein bridging genome topology and function. Nat Rev Genet. 2014;15(4):234.
Kundaje A, Meuleman W, Ernst J, Bilenky M, Yen A, Heravi-Moussavi A, Kheradpour P, Zhang Z, Wang J, Ziller MJ. Integrative analysis of 111 reference human epigenomes. Nature. 2015;518(7539):317.
Rubin AJ, Barajas BC, Furlan-Magaril M, Lopez-Pajares V, Mumbach MR, Howard I, Kim DS, Boxer LD, Cairns J, Spivakov M. Lineage-specific dynamic and pre-established enhancer–promoter contacts cooperate in terminal differentiation. Nat Genet. 2017;49(10):1522.
Phanstiel DH, Van Bortle K, Spacek D, Hess GT, Shamim MS, Machol I, Love MI, Aiden EL, Bassik MC, Snyder MP. Static and dynamic DNA loops form AP-1-bound activation hubs during macrophage development. Mol Cell. 2017;67(6):1037–48.
Kouwenhoven EN, van Heeringen SJ, Tena JJ, Oti M, Dutilh BE, Alonso ME, de la Calle-Mustienes E, Smeenk L, Rinne T, Parsaulian L, et al. Genome-wide profiling of p63 DNA-binding sites identifies an element that regulates gene expression during limb development in the 7q21 SHFM1 locus. PLoS Genet. 2010;6(8):e1001065.
Vigano MA, Lamartine J, Testoni B, Merico D, Alotto D, Castagnoli C, Robert A, Candi E, Melino G, Gidrol X, et al. New p63 targets in keratinocytes identified by a genome-wide approach. EMBO J. 2006;25(21):5105–16.
McDade SS, Henry AE, Pivato GP, Kozarewa I, Mitsopoulos C, Fenwick K, Assiotis I, Hakas J, Zvelebil M, Orr N, et al. Genome-wide analysis of p63 binding sites identifies AP-2 factors as co-regulators of epidermal differentiation. Nucleic Acids Res. 2012;40(15):7190–206.
Sethi I, Sinha S, Buck MJ. Role of chromatin and transcriptional co-regulators in mediating p63-genome interactions in keratinocytes. BMC Genom. 2014;15:1042.
Cavazza A, Miccio A, Romano O, Petiti L, Malagoli Tagliazucchi G, Peano C, Severgnini M, Rizzi E, De Bellis G, Bicciato S, et al. Dynamic transcriptional and epigenetic regulation of human epidermal keratinocyte differentiation. Stem Cell Reports. 2016;6(4):618–32.
Fan X, Wang D, Burgmaier JE, Teng Y, Romano RA, Sinha S, Yi R. Single cell and open chromatin analysis reveals molecular origin of epidermal cells of the skin. Dev Cell. 2018;47:21–37.
Soares E, Xu Q, Li Q, Qu J, Zheng Y, Raeven HHM, Brandao K, van den Akker WMR, Tang F, Zhou H. Single-cell RNA-seq identifies a reversible epithelial-mesenchymal transition in abnormally specified epithelia of p63 EEC syndrome. bioRxiv 2019;437632.
Yang A, Zhu Z, Kapranov P, McKeon F, Church GM, Gingeras TR, Struhl K. Relationships between p63 binding, DNA sequence, transcription activity, and biological function in human cells. Mol Cell. 2006;24(4):593–602.
Kouwenhoven EN, van Bokhoven H, Zhou H. Gene regulatory mechanisms orchestrated by p63 in epithelial development and related disorders. Biochim Biophys Acta. 2015;1849(6):590–600.
Schwalie PC, Ward MC, Cain CE, Faure AJ, Gilad Y, Odom DT, Flicek P. Co-binding by YY1 identifies the transcriptionally active, highly conserved set of CTCF-bound regions in primate genomes. Genome Biol. 2013;14(12):R148.
Beagan JA, Duong MT, Titus KR, Zhou L, Cao Z, Ma J, Lachanski CV, Gillis DR, Phillips-Cremins JE. YY1 and CTCF orchestrate a 3D chromatin looping switch during early neural lineage commitment. Genome Res. 2017;27:1139–52.
Chernukhin I, Shamsuddin S, Kang SY, Bergström R, Kwon Y-W, Yu W, Whitehead J, Mukhopadhyay R, Docquier F, Farrar D. CTCF interacts with and recruits the largest subunit of RNA polymerase II to CTCF target sites genome-wide. Mol Cell Biol. 2007;27(5):1631–48.
Wendt KS, Grosveld FG. Transcription in the context of the 3D nucleus. Curr Opin Genet Dev. 2014;25:62–7.
Fullwood MJ, Wei C-L, Liu ET, Ruan Y. Next-generation DNA sequencing of paired-end tags (PET) for transcriptome and genome analyses. Genome Res. 2009;19(4):521–32.
Mumbach MR, Rubin AJ, Flynn RA, Dai C, Khavari PA, Greenleaf WJ, Chang HY. HiChIP: efficient and sensitive analysis of protein-directed genome architecture. Nat Methods. 2016;13(11):919.
Jin F, Li Y, Dixon JR, Selvaraj S, Ye Z, Lee AY, Yen C-A, Schmitt AD, Espinoza CA, Ren B. A high-resolution map of the three-dimensional chromatin interactome in human cells. Nature. 2013;503(7475):290.
Bonev B, Cohen NM, Szabo Q, Fritsch L, Papadopoulos GL, Lubling Y, Xu X, Lv X, Hugnot, JP, Tanay A, Cavalli G. Multiscale 3D genome rewiring during mouse neural development. Cell 2017;171(3):557–72.
Stadhouders R, Vidal E, Serra F, Di Stefano B, Le Dily F, Quilez J, Gomez A, Collombet S, Berenguer C, Cuartero Y, Hecht J. Transcription factors orchestrate dynamic interplay between genome topology and gene regulation during cell reprogramming. Nat Genet 2018;50(2):238.
Browne G, Cipollone R, Lena AM, Serra V, Zhou H, van Bokhoven H, Dotsch V, Merico D, Mantovani R, Terrinoni A, et al. Differential altered stability and transcriptional activity of DeltaNp63 mutants in distinct ectodermal dysplasias. J Cell Sci. 2011;124(Pt 13):2200–7.
Celli J, Duijf P, Hamel BC, Bamshad M, Kramer B, Smits AP, Newbury-Ecob R, Hennekam RC, Van Buggenhout G, van Haeringen A, et al. Heterozygous germline mutations in the p53 homolog p63 are the cause of EEC syndrome. Cell. 1999;99(2):143–53.
Hong J-W, Hendrix DA, Levine MS. Shadow enhancers as a source of evolutionary novelty. Science. 2008;321(5894):1314.
Barolo S. Shadow enhancers: frequently asked questions about distributed cis-regulatory information and enhancer redundancy. BioEssays. 2012;34(2):135–41.
van Bokhoven H, Hamel BC, Bamshad M, Sangiorgi E, Gurrieri F, Duijf PH, Vanmolkot KR, van Beusekom E, van Beersum SE, Celli J, et al. p63 Gene mutations in eec syndrome, limb-mammary syndrome, and isolated split hand-split foot malformation suggest a genotype-phenotype correlation. Am J Hum Genet. 2001;69(3):481–92.
Rheinwald JG, Green H. Epidermal growth factor and the multiplication of cultured human epidermal keratinocytes. Nature. 1977;265(5593):421–4.
Rinne T, Clements SE, Lamme E, Duijf PH, Bolat E, Meijer R, Scheffer H, Rosser E, Tan TY, McGrath JA, et al. A novel translation re-initiation mechanism for the p63 gene revealed by amino-terminal truncating mutations in Rapp-Hodgkin/Hay-Wells-like syndromes. Hum Mol Genet. 2008;17(13):1968–77.
Liu NQ, Ter Huurne M, Nguyen LN, Peng T, Wang S-Y, Studd JB, Joshi O, Ongen H, Bramsen JB, Yan J. The non-coding variant rs1800734 enhances DCLK3 expression through long-range interaction and promotes colorectal cancer progression. Nat Commun. 2017;8:14418.
Novakovic B, Habibi E, Wang S-Y, Arts RJ, Davar R, Megchelenbrink W, Kim B, Kuznetsova T, Kox M, Zwaag J. β-glucan reverses the epigenetic state of LPS-induced immunological tolerance. Cell. 2016;167(5):1354–68.
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15–21.
Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106.
Li H, Durbin R. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics. 2009;25(14):1754–60.
Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, Nussbaum C, Myers RM, Brown M, Li W, et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 2008;9(9):R137.
Ramírez F, Dündar F, Diehl S, Grüning BA, Manke T. deepTools: a flexible platform for exploring deep-sequencing data. Nucleic Acids Res. 2014;42(W1):W187–91.
Tripathi S, Pohl MO, Zhou Y, Rodriguez-Frandsen A, Wang G, Stein DA, Moulton HM, DeJesus P, Che J, Mulder LC. Meta-and orthogonal integration of influenza “OMICs” data defines a role for UBR4 in virus budding. Cell Host Microbe. 2015;18(6):723–35.
Wingett S, Ewels P, Furlan-Magaril M, Nagano T, Schoenfelder S, Fraser P et al. HiCUP: pipeline for mapping and processing Hi-C data. F1000Res 2015; 4:1310.
Cairns J, Freire-Pritchett P, Wingett SW, Várnai C, Dimond A, Plagnol V, Zerbino D, Schoenfelder S, Javierre B-M, Osborne C. CHiCAGO: robust detection of DNA looping interactions in Capture Hi-C data. Genome Biol. 2016;17(1):127.
We thank Eva Janssen-Megens, Siebe van Genesen and Rita Bylsma for operating the Illumina analyzer. We thank the ENCODE Consortium for sharing their data.
JQ, GY and HZ conceived and designed the experiments. JQ, GY, HZ wrote and revised the manuscript. JQ performed the experiments. JQ, GY and HZ analyzed the data. All authors read and approved the final manuscript.
This research was supported by Netherlands Organization for Scientific Research (NWO/ALW/MEERVOUD/836.12.010, HZ), Radboud University fellowship (HZ) and Chinese Scholarship Council grant 201406330059 (JQ).
To review our complete dataset GEO accession GSE123711 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE123711). All data supporting the findings of the study and in-house codes are available on request.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.