BACKGROUND
tp <dig> is an important transcription factor of the p53/p63/p <dig> family that dictates a wide range of cellular properties including but not limited to stem cell renewal, lineage choices and maintaining the balance between proliferation and differentiation
<cit> .
this diverse function of p <dig> is critical for morphogenesis during development, particularly for epithelial-enriched tissues such as the skin and its appendages such as the hair follicles and mammary glands.
indeed, p63-null mice die after birth and exhibit a dramatic agenesis of epithelial-rich structures and widespread developmental defects of the limb, orofacial region, and external genitalia
<cit> .
these p63-deficient structural defects are thought to be the result of a failed program of epithelial stratification and/or diminished capacity for stem cell renewal, both of which can jeopardize normal epithelial-stromal interactions needed during embryonic organ development
<cit> .
in agreement with the mouse phenotype, p <dig> mutations in humans lead to congenital abnormalities such as abnormal limb development and ectodermal dysplasia, which are associated with a spectrum of developmental disorders including aec or eec syndrome
<cit> .
the biological function of p <dig> is mediated by several isoforms derived from distinct transcripts
<cit> .
these include the longer tap <dig> isoforms and n-terminal deleted Δnp <dig> isoforms generated from an internal promoter located within intron  <dig>  furthermore alternative splicing can result in α, β, and γ isoforms, which differ in the c-terminus.
all p <dig> isoforms share the dna-binding and oligomerization domains, which are analogous to that of p <dig>  it is now well-established that ∆np <dig>  especially ∆np63α is the predominant isoform that is present in most epithelial cells such as the keratinocytes of the skin
<cit> .
importantly both gene complementation studies and isoform specific knockouts have conclusively affirmed that ∆np <dig> harbors most of the function and biological activity of p <dig>  particularly as it pertains to the epithelial tissues
.
the role of p <dig> in regulating transcription during development has been extensively studied in skin where ∆np <dig> is highly expressed and regulates the transition from simple ectodermal cells to stratified epithelium
<cit> .
given the master regulatory function of p <dig>  it is not surprising that the repertoire of p63-targets is vast and represents practically every crucial gene regulatory and signaling pathway.
this is evident from the ~ <dig>  binding sites for p <dig> in human keratinocytes as determined by chromatin immunoprecipitation with next-generation sequencing  studies
<cit> .
p <dig> controls expression of basal keratin genes k <dig> and k <dig> and regulates myc levels thereby controlling keratinocyte proliferation via the wnt/β-catenin and notch signaling pathways
<cit> .
the keratinocyte differentiation program is also regulated by p <dig>  in part via its effect on the znf750-klf <dig> regulatory axis
<cit> .
while the identification of p <dig> bound cis-regulatory elements in keratinocytes has received much attention, the mechanics of p63-dna interaction is still relatively unknown.
p <dig> binds a canonical motif, defined as closely spaced  <dig> decamers , although there is growing evidence that p <dig> can target sites that do not completely conform to this consensus sequence, including half sites
<cit> .
given the degenerate nature of the p <dig> binding motif, it is not surprising that by conservative estimates, there are more than  <dig> million such potential sites in the human genome.
however, as is the case with most other transcription factors , only a small subset of these sites are bound by p <dig> in vivo <cit> .
it is likely that the local chromatin architecture, among other factors plays an important deterministic role in dictating how and why p <dig> selects its target dna.
hence, this is an important area of future investigation; especially given the increasing evidence that p <dig> can play an important role in modulating the chromatin structure.
indeed, recent studies have demonstrated that p <dig> can functionally interact with several epigenetic factors in keratinocytes, which can in turn profoundly influence p63-dependent transcriptional activation and repression.
examples of such interactions include the reinforcement of p <dig> mediated repression of p <dig> by lsh, a member of the snf <dig> family of chromatin remodeling atpases
<cit> , direct recruitment of histone deacetylases, hdac <dig> and hdac <dig> by ∆np <dig> during repression of target genes in the embryonic epidermis
<cit>  and the crosstalk between p <dig> and chromatin organizer satb <dig> in regulating keratinocyte differentiation genes
<cit> .
p <dig> can also control higher-order chromatin structure in epidermal progenitor cells during skin development by regulating brg <dig>  a atp-dependent chromatin remodeler
<cit> .
given these emerging links between p <dig> and chromatin, it is important that any comprehensive studies on the mechanism of p63-genome interactions takes into account the underlying state of epigenetic modifications.
here we have utilized the p <dig> chip-seq dataset and available chromatin modification datasets for normal human epidermal keratinocytes  to investigate the rules that govern binding of p <dig> to its target dna.
we find that p <dig> binds to a canonical motif  at the majority  of its sites, whereas non-canonical motifs containing 1– <dig> spacer between decamers are present in only  <dig> % of the sites.
the chromatin at p <dig> binding sites is largely marked by active histone modifications .
moreover, chromatin accessibility with h3k4me <dig> can accurately predict bona-fide bound p <dig> sites without the need for any additional dna sequence information.
finally, using a comprehensive in silico approach, we identify several cooperating tfs that appear to define specific classes of p <dig> regulated genes.
RESULTS
underlying sequence patterns and chromatin architecture of p <dig> targets
several groups have determined global p <dig> binding locations in various primary and immortalized keratinocytes using chip-chip or chip-seq techniques
.
for our studies, we focused on the most comprehensive p <dig> chip-seq data
<cit>  available to date.
it had an added benefit of being generated from primary keratinocytes  and more importantly conforming to encode guidelines
<cit> .
to facilitate uniform comparisons across other encode datasets, we re-aligned the p <dig> chip-seq to the latest human genome build  with bowtie
<cit> .
in strong agreement with kouwenhoven et al., by using high stringency conditions , we identified a reliable and robust dataset of  <dig> p <dig> binding sites that were common among the three biological replicates.
on examining the underlying dna sequence of these p63-chiped elements, we found that  <dig> % of these sites have at least one p <dig> canonical motif.
among these, 32% show a close match to the p <dig> consensus  while the remaining  <dig> % are a weaker match .
both the strong and weak canonical motifs are significantly enriched at the p <dig> chiped regions compared to random genomic regions .
an additional  <dig> % of the binding sites show a close match to the non-canonical p <dig> motif, which has 1– <dig> base spacers separating the two half-sites .
interestingly, 10% of the p <dig> bound genomic sites do not have a recognizable p <dig> motif raising the possibility that p <dig> can perhaps also be recruited to target regulatory elements through indirect mechanisms such as protein-protein interactions.
our analysis also revealed that only a few of the p <dig> binding sites contained just a half site .
next we compared the chromatin architecture of the  <dig>  p <dig> bound sites and randomly selected  <dig>  unbound sites, which have at least  <dig> strong p <dig> canonical motif.
for this purpose we focused on  <dig> active  and  <dig> repressive  histone modifications profiles in nhek cells that have been generated by the encode consortium
<cit> .
as shown in figure
1a, the genomic segments of dna bound by p <dig> are characterized by a chromatin architecture consisting of high h3k27ac, h3k4me <dig>  h3k4me <dig>  and h3k4me <dig> and low h3k9me <dig> and h3k27me <dig> modifications.
in contrast, the randomly selected stretches of dna corresponding to the unbound p <dig> sites are completely lacking in these chromatin profiles .
hence chromatin architecture could contribute to selective targeting of p <dig> figure  <dig>
chromatin architecture differs between p <dig> bound and unbound sites.
average profiles of the  <dig> active  and  <dig> repressive  histone modifications are plotted for a 1 kb window centered at p <dig> motif.
average chromatin architecture at  <dig> p <dig> bound locations.
average chromatin architecture at  <dig>  unbound locations that have a strong p <dig> motif.
profiles were generated with archtex
<cit> .
p <dig> has been shown to target different types of regulatory regions  and involved in both activation and repression of gene expression
.
recent studies have demonstrated that specific states of chromatin modifications at the regulatory regions are strongly associated with the level of gene expression for the corresponding genes
<cit> .
to determine whether such differences exists for distinct classes  of p <dig> targets, we clustered the p <dig> bound regulatory regions by their underlying chromatin architecture.
we first divided p <dig> targets into two groups  based on the magnitude of signal  and then performed unsupervised clustering using the spatial arrangement of the histone modifications
<cit> .
using standard k-means clustering  we were able to capture all patterns within cluster a and b but this resulted in redundancy .
therefore, we compared all  <dig> sub-clusters to each other in both orientations, grouping together those with a pearson correlation above  <dig> .
this resulted in the classification of p <dig> targets into  <dig> distinct, non-redundant groups .
based on chromatin based segmentation
<cit> , the majority of p <dig> targets are predicted to represent enhancers   while only  <dig>  represent active promoters .
the distance of the regulatory regions to the nearest tss supports these predictions .
it is interesting to note that none of the clusters are characterized by repressive histone marks.
however this does not imply that p <dig> cannot act as a repressor.
in fact the largest sub-cluster b <dig>  has very low signal for active chromatin marks and predicted by chromatin based segmentation likely to represent repressed/inactive targets.
indeed, upon probing the rna-seq data from keratinocytes, it is evident that the regulatory regions of the b <dig> sub-cluster are associated with genes that are weakly expressed .
therefore this cluster might represent repressed/poised genes for which p <dig> may not play an active role in activation at least in the basal growth conditions of nhek cells.figure  <dig>
clustering of p <dig> targets by chromatin intensity and shape.
p <dig> bound locations are clustered by  <dig> active  and  <dig> repressive  histone modifications sequenced in nhek  cell-line.
they are first clustered by average signal intensity across a 1 kb window centered at p <dig> binding site and then clustered by spatial arrangement of the histone modifications.
annotating p <dig> targets clustered by chromatin profiles
p <dig> bound locations are clustered by  <dig> active  and  <dig> repressive  histone modifications sequenced in nhek cell-line.
the  <dig> groups are annotated by chromatin based segmentation and median distance to nearest tss.
great is used to find the top  <dig> go-bp terms and top  <dig> panther pathways that enrich for each cluster.
we next used genomic regions enrichment of annotations tool
<cit>  to determine whether the five clustered groups of p <dig> targets could be segregated into distinct classes of genes involved in specific biological pathways.
our analysis revealed that each cluster indeed was overrepresented by genes that were involved in closely related, yet disparate biological activity.
for example, cluster b <dig>  which primarily encodes for strong enhancers is enriched for gene-ontology biological process  of epidermis development, while cluster b <dig> representing repressed/inactive sites is enriched for go-bp of hair cycle.
this raises the intriguing possibility that p <dig> might play an important role in actively driving the epidermal developmental processes while keeping the hair cell fate repressed – a notion that is supported by data from prior transgenic mouse studies
<cit> .
the complete annotation for each group of p <dig> targets is provided in table
<dig>  thus, collectively the data obtained by chromatin architecture-dependent clustering of p <dig> targets allowed us to determine both the active and poised targets of p <dig> and annotate these to specific biological processes .
p <dig> binding can be accurately predicted from chromatin modifications
we next implemented two different statistical approaches to test whether epigenetic modifications can be exclusively utilized to identify functional p <dig> binding sites.
we generated a discriminant and a regression model of p63-binding events based on data from histone modifications as determined by chip-seq and chromatin accessibility as inferred from dnase-seq and faire-seq experiments.
dnase-seq and faire-seq assays generate both distinct and overlapping accessibility information about the genomic landscape
<cit> .
also the different histone modifications usually correspond to transcriptional regulation in a combinatorial fashion
<cit> .
as these various chromatin features are not mutually independent of each other, we created interaction datasets from each pair of chromatin datasets, such as dnase*faire .
these were also used as predictors for our statistical models.
for the discriminant approach we generated the model on a training dataset containing  <dig> p <dig> binding sites with  <dig>  unbound genomic locations and for the regression approach we generated the model on a training dataset containing  <dig> p <dig> binding sites with  <dig>  unbound genomic locations .
we then validated the models on a test dataset of  <dig>  genomic locations .
the discriminant model identifies boundaries between groups of objects , the boundaries being defined in terms of those variable characteristics that distinguish the objects into two groups.
this technique allows us to determine which chromatin modifications are strongly predictive of p <dig> binding status .
in addition this model can be used to predict if an unknown site is bound or unbound by p <dig>  we found that a three variable model was able to accurately classify p <dig> binding sites , specificity:  <dig> %) .
this model was based on df, h3k4me <dig>  and h3k27me <dig> datasets.
in this model both df and h3k4me <dig> had positive coefficients as anticipated by p <dig> binding to sites with high values for these chromatin datasets.
on the other hand, the repressive modification h3k27me <dig> had a negative coefficient suggesting that p <dig> is being excluded from genomic sites with this histone modification.figure  <dig>
p <dig> occupancy can be accurately predicted by chromatin marks.
<dig> active  and  <dig> repressive  histone modifications with chromatin accessibility as measured by dnase-seq and faire-seq were used to construct models of p <dig> binding.
discriminant model for p <dig> binding.
regression model for p <dig> occupancy.
the values for  <dig> datasets used for our models are shown with the data sorted in descending order by experimental p <dig> occupancy.
the standardized number of tags for the histone modifications and chromatin accessibility are displayed.
the dnase and faire interaction term  is calculated as product of the two and is displayed in standardized tag space.
the quality of the p <dig> sequence motif was determined by patser, good motif matches are purple and low quality motif matches are pink.
the results of our discriminant model is compared directly to the experimental determination of binding by macs.
our regression model is compared directly to experimental occupancy.
we also utilized the regression model to determine which chromatin characteristic best predicts the occupancy of p <dig> binding at any genomic location .
when applied to the test dataset, the regression model with only three variables has a mse  of  <dig>  and r <dig> of  <dig>   .
for the regression model h3k4me <dig>  dnase, and df variables were the most informative.
for both models, a chromatin accessibility dataset  with h3k4me <dig> were found to be very predictive .
furthermore, a regression model using only the df term has a r2 =  <dig> , which exemplifies the importance of chromatin accessibility in defining true in vivo binding sites.
however these models do not inform whether the accessibility is a prerequisite or a secondary effect of p <dig> binding.
we therefore integrated the nucleosome-dna interaction model in our analysis, which is based on the principle that certain dna sequences show a greater predisposition to wrap around the histone octamer and form nucleosomes
<cit> .
this model has previously been used to predict and experimentally validate that p <dig> preferentially binds nucleosome rich regions
<cit> .
our analysis showed that the p63-bound genomic segments in keratinocytes have higher sequence defined nucleosome occupancy than randomly chosen sequences .
this suggests that the increased chromatin accessibility at p <dig> bound sites is likely to be actively shaped by p63-dna interactions.figure  <dig>
p <dig> preferentially binds nucleosome rich regions.
predicted nucleosome occupancy as calculated by the nucleosome-dna interaction model
<cit> , was plotted for a 1 kb window centered at p <dig> binding site.
blue: average nucleosome occupancy at  <dig> p <dig> bound locations.
red: average nucleosome occupancy at  <dig>  random unbound locations.
for both the regression and discriminant models, p <dig> motif data was not included even though p <dig> motif is a statistically significant predictor .
indeed on its own the sequence information accounts for 14% of the variability in p <dig> binding .
we therefore performed additional analysis to ascertain the importance of p63′s motif in dictating dna binding, when taking into account epigenetic information.
we divided all accessible genomic locations as determined by dnase-seq into three groups based on the presence of a strong, weak, or the absence of a p <dig> motif.
our results indicate that the presence of a strong or weak p <dig> motif leads to stronger p <dig> occupancy .
this finding that p <dig> dna-sequence motif is an important determinant of in-vivo binding is difficult to reconcile with the fact that p <dig> motif was not considered a strong predictor in our statistical models.
one possible explanation for this discrepancy is that the chromatin datasets we use in our models have already embedded in it the dna sequence motif information.
indeed, using our regression model, we predicted p <dig> occupancy for each group and showed that predicted p <dig> occupancy is also higher for sites containing a strong or weak p <dig> motif .
this provides further support for the notion that h3k4me <dig>  dnase and df variables in our statistical models already account for the information represented by the p <dig> consensus motif.
thus accessible sites with a strong or weak motif are more accessible and have higher deposition of h3k4me <dig> than sites lacking a p <dig> motif .
taken together, these results lead us to postulate that p <dig> might exist as part of a chromatin remodeling complex, which creates a distinct epigenetic architecture at its binding sites .figure  <dig>
significance of p <dig> motif.
the genomic coordinates of all the accessible locations in human keratinocytes, as determined by dnase-seq in nhek cell-line, were obtained from encode via the ucsc genome browser.
the  <dig>  accessible locations are divided into  <dig> groups 1)  <dig> locations with strong p <dig> motif , 2)  <dig>  locations with weak motif  and 3)  <dig> locations with no motif.
box and whisker plots are made across 500 bp window for  experimental p <dig> occupancy,  predicted p <dig> occupancy by our best  <dig> regression model,  dnase tag density,  h3k4me <dig> tag density.
identification of p <dig> cooperating tfs
it is well established that many tfs often act in a combinatorial fashion to govern tissue-specific gene expression.
hence, we wanted to examine the repertoire of p <dig> associated tfs that might play such role in modulating p <dig> binding and possibly influencing p63-dependent gene expression in keratinocytes.
we therefore used a combination of in-silico analysis of the p <dig> binding sites and careful data mining of large-scale genomic datasets such as rna-seq and chip-seq from the encode project.
we posit that the tf that are likely to directly cooperate with p <dig> will have the following features.
first, their dna binding motifs will likely be enriched at p <dig> bound elements, keeping in mind that co-occurrence of such motifs adjacent to p <dig> sites is not a prerequisite for such interactions.
second, their in-vivo binding profiles will overlap with p <dig> binding profile.
here we reasoned that tf binding to some extent could be extrapolated from chip-seq in other cell types, as it has been shown that tfs share a large number of common binding sites across different cell lines
<cit> .
finally, we contend that the relevant tf should be expressed in keratinocytes.
applying this criterion, we examined  <dig> tfs for which chip and/or expression data was available from the encode project.
of these,  <dig> were enriched at p <dig> targets as determined by their in-vivo binding profile or in silico motif analysis.
however, we rejected  <dig> of these tfs as potential cooperating tfs because they were either not expressed in keratinocytes  or not enriched by both chip-seq and motif analysis .
an additional  <dig> tfs were also not considered from further evaluation due to ill-defined/unknown dna binding motifs or missing chip-seq data.
this shortlisted the potential p <dig> cooperating tfs to  <dig> tfs, which had both an enriched motif and in-vivo binding signal  at p <dig> bound locations in nhek cells .
these  <dig> tfs include: cebpb ; cfos, fosl <dig>  and jund ; bach <dig> ; tfap2c ; stat <dig> and stat <dig> ; max, c-myc, and usf <dig> , rfx <dig> ; elk <dig>  .
in this context, it is important to note that several of the predicted tf discovered by our analysis, such as ap <dig>  ap <dig>  myc, stat and c/ebp have been previously linked to p <dig> and keratinocyte biology
.
on the other hand, some of the other potential p <dig> cooperating tfs such as rfx <dig>  elk <dig> and bach <dig> are new members of this class.
they represent an interesting group for which much less is known about their expression and function in keratinocytes and any possible correlation with p <dig>  one particularly intriguing p <dig> cooperating tf is rfx <dig>  which is highly expressed in nhek cells as evident by the rna-seq data .
furthermore, a recent study has also shown that rfx <dig> mrnas and protein are preferentially expressed in the p63-rich basal cells of the human epidermis further raising the prospects of a functional interplay between these two tfs
<cit> .figure  <dig>
cooperating tfs of p <dig>   <dig> tfs corresponding to  <dig> motifs , bach <dig>  ap <dig>  stat <dig>  stat <dig>  mycmax , usf, rfx, elk1) were determined as potential cooperating tfs for p <dig>  the characteristics that determine these tfs as potential cooperating tfs are: 1) enrichment of tf binding at p <dig> targets  - as measured by chip-seq from any available cell-line, 2) high expression in keratinocytes  - as determined by rna-seq in nhek cell-line, 3) enrichment of the tf’s motif at p <dig> targets  and co-occurrence with p <dig> motif .
the tfs are sorted in descending order by an average percentage rank calculated across the different analysis.
the potential role of p <dig> in modulating different biological functions in coordination with various cooperating tfs
p <dig> is involved in myriad biological functions and it is possible that this diversity of p <dig> function is brought about by in part by different subtypes of p63-tf complexes.
we therefore constructed a correlation matrix of the  <dig> tfs based on their in-vivo binding signal and their dna-binding motif scores at the  <dig>  p <dig> bound locations .
the  <dig> tfs correspond to  <dig> unique dna-binding motifs .
there appears to be multiple complexes with unclear distinctions between them .
two distinct complexes containing overlapping tfs could be identified.
the first group comprised of c/ebpb, jund, cfos, and stat <dig>  whereas the second complex was primarily represented by elk <dig>  fosl <dig>  and myc-max.
we next asked if there were different sub-groups of genomic targets where p <dig> was bound in coordination with these distinct cooperating tfs.
upon clustering the  <dig> p <dig> targets,  <dig> clear groups emerged  based on the presence/absence of the cooperating tf motifs .
ap <dig>  stat <dig> and mycmax seem to be the driving force for the clustering, whereas cebpb, ap <dig> and elk <dig> motifs were quite ubiquitous in their presence.
further annotating the clusters using great allowed us to associate specific groups of p <dig> targets with distinct biological functions, based upon top go terms and significant pathways.
the complete annotation for each group of p <dig> targets and the type of complex it represents is shown in figure
7b.our clustering analysis revealed an interesting and complex correlation between p <dig>  its network of associated cooperating tfs and the nature and physiological role of the p <dig> target genes.
case in example is ap <dig>  it seems to be an important co-player of p <dig> in regulating genes involved in keratinocyte differentiation and epidermis development , modestly linked with hair follicle development and cell cycle arrest  but not at all involved in the activation of mapkk activity  or regulation of carbohydrate metabolic process .
similarly, while stat <dig> and stat <dig> seem to share similar co-regulatory function as evident by their enrichment in clusters c <dig> and c <dig>  there are also clear differences between these two closely related family members based on their enrichment patterns for clusters c <dig>  c <dig> and c <dig>  this point is further illustrated by the correlation matrices made in figure
7a.
in contrast, binding motifs for c/ebp, elk <dig> and ap <dig> did not seem to exhibit a preference for any of the particular clusters representing different go terms or biological pathways.
interestingly, none of the  <dig> clusters showed any difference in terms of p <dig> motif strength or the type of regulatory regions .figure  <dig>
p <dig> functions with multiple cooperating tfs.
a correlation matrix for the  <dig> probable p <dig> cooperating tfs at the  <dig>  p <dig> locations.
upper triangle shows pearson’s correlation coefficient calculated between chip-seq tag densities for each pair of tfs.
for bach <dig> and fosl <dig>  chip-seq signal are from k <dig> and a <dig> cell-lines respectively.
while for the rest of the  <dig> factors, chip-seq signal is from hela-s <dig> cell-line.
lower triangle shows correlation between motif scores generated by patser.
the  <dig> p <dig> bound locations are clustered by the presence  or absence  of the  <dig> cooperating tfs motifs , in a 1 kb window centered at p <dig> binding site.
one obvious limitation of our approach is the fact that for the purpose of the aforementioned studies, a specific p63-binding site was assigned to be regulatory element for the nearest gene, often without any supporting experimental data.
however it has been shown, that linking distal enhancer elements to the nearest gene, as is the norm, might not always be an accurate representation of the regulatory landscape
<cit> .
there is an absence of chromatin interactions data such as those obtained from chromosome conformation capture , in keratinocytes.
we therefore virtually linked the distal regulatory elements to their target genes, by exploiting the characteristic that regulatory regions become dnase hypersensitive  in synchrony with their target promoters.
to examine this further, we first retrieved the global map of distal dhs-to-promoter connections as generated by the encode consortium
<cit> .
overlapping this map with the p <dig> chip-seq data allowed us to link  <dig> p <dig> bound distal genomic regions to their putative target promoters/genes .
further analysis revealed interesting facets of long-range interactions between distal p63-bound elements and the promoter as evident by the krt <dig> gene, a known p <dig> target
<cit> .
while prior studies have focused primarily on the proximal regions 5′ of the krt <dig> gene, including a well-characterized p63-bound enhancer ~ <dig>  kb upstream
<cit> , we discovered  <dig> novel distal p63-bound elements that are located quite far away  as predicted by the dhs-to-promoter connections .
two of these p63-bound genomic segments, which are predicted to be enhancers for krt <dig> are shown in figure
8a.
these regulatory regions also contained potential binding sites for some of the p <dig> cooperating tfs such as cebpb, c-myc and ap <dig>  interestingly these two distal p63-bound elements are not exclusive to the krt <dig> gene since they are also putatively linked to related krt16/ krt <dig> genes, which are located in a relatively closely spaced genomic cluster .
this raises an interesting possibility that a specific p63-bound regulatory region  might be commonly utilized to coordinate the regulation of multiple closely related genes in keratinocytes.figure  <dig>
model of p <dig> binding at distal regulatory elements of keratin- <dig> and rfx <dig>
histograms in green show p <dig> tag density in a 1 kb window at the regulatory region.
loops links the distal regulatory regions to their target promoter as mapped by the encode project
<cit> .
<dig> p <dig> binding sites  <dig>  kb and  <dig>  kb upstream are linked to keratin- <dig> .
these regions also regulate other genes, among them keratin- <dig>  and keratin- <dig> .
motif analysis reveals other tfs binding within the 1 kb window.
p <dig> binds  <dig>  kb upstream to rfx <dig> at a distal regulatory element linked to rfx <dig> promoter.
this region is also linked to tuft <dig> and rorc promoters.
motif analysis shows p <dig> binding in complex with ap <dig>  e-box binding factors and bach <dig> on a weak canonical motif.
additional motifs are also found within the 1 kb window.
we also discovered possible new target genes under the control of p63-bound distal regulatory regions that otherwise would have been missed by the conventional strategy of assigning a p63-regulated element to its closest gene.
as shown in figure
8b, one such interesting candidate is rfx <dig>  a novel cooperating tf of p <dig> as discussed above.
p <dig> binds a regulatory element  <dig>  kb upstream of the rfx <dig> gene and motif analysis suggests that ap <dig>  e-box binding factors and bach <dig> are part of the transcriptional complex at this site .
hence, there might exist a transcriptional regulatory loop whereby p <dig> in cooperation with additional tfs activates rfx <dig>  which then in turn, modulates p <dig> binding to its target sites.
similar to the case with the krt <dig> gene, this distal regulatory element for rfx <dig> is predicted to be also linked to the tuft <dig> and rorc genes.
in the future in-depth studies such as 3c experiments will help to confirm these distal enhancer-promoter interactions and to firmly establish the true identity of the p63-driven gene network in keratinocytes.
discussion
this study aimed to decipher the mechanics of p <dig> binding by determining the minimal in-vivo motif required for binding, distinguishing between chromatin architecture of bound and unbound motifs and identifying cooperating tfs that modulate p <dig> biological activity.
p <dig> requires a full site for binding to target sites
the transcription factor p <dig> binds as a homotetramer to two decamers rrrcrwgyyy, rrrcwygyyy separated by a 0– <dig> base pair spacer region.
we found that the p <dig> sites containing the decamer pair with intervening spacers was much less prevalent  than those where the two half sites were juxtaposed to each other .
a very small subset  of p <dig> binding sites consisted of only a half site  .
this is in contrast to a similar study that found 3-4% of p <dig> binding sites having a half site
<cit> .
this observation can be explained by a slight difference in the consensus motifs of the two factors, which also results in  <dig> fold lower binding affinity of p <dig>  in comparison to p53
<cit> .
one possibility is that the dimer-dimer interactions are important for p <dig> dna binding specificity and therefore p <dig> requires a full site to bind dna efficiently.
indeed, such differences in the dna-protein interactions among p <dig> family members are quite evident from recent structural studies with p73
<cit> .
interestingly, not only the distance between two p <dig> half-sites influences the p <dig> quaternary structure, but tellingly transcriptional activity is also more affected by spacer length in p <dig> response element than in p <dig>  finally, it is worth noting that ~10% of p <dig> chiped sites in keratinocytes do not have a recognizable p <dig> binding sequence suggesting that the p <dig> binding at these sites is driven by indirect mechanisms that might involve other dna-binding tfs and/or non-canonical p <dig> motifs.
the p <dig> consensus motif is not required for predicting binding events in keratinocytes
we utilized two computational modeling approaches to uncover key characteristics defining p <dig> binding sites.
first, we utilized discriminant modeling that allowed us to predict p <dig> binding as a binary score .
second, in a parallel strategy, we used regression modeling that predicted the degree of p <dig> occupancy.
our approach was distinct from other published methods in that we trained our statistical models on random genomic sites which might or might not have p <dig> binding sequences
<cit> .
this allowed our models to include the 10% binding sites that do not have p <dig> dna-binding motif.
surprisingly, we found that our final models constituting only chromatin marks  predicted p <dig> binding with high accuracy  sensitivity, regression model –  <dig>  mse).
adding sequence information to the models, did not lead to any significant improvement even though p <dig> motif is a statistically significant predictor .
one possible explanation is that in keratinocytes, regions of the genome that have a functional p <dig> motif are on average more accessible and marked with active chromatin marks.
this result is not surprising if p <dig> is a key component of the regulatory complex that is involved in remodeling the chromatin at its binding sites.
indeed p <dig> target sequences dictate higher nucleosome occupancy than random genomic sequences according to the nucleosome-dna interaction model .
it can be then hypothesized that p <dig> binding shifts the nucleosomes creating an accessible active chromatin structure at its targets.
support for such a role for p <dig> comes from a recent study, which examined the selective loss versus gain of dhss targeted by lineage regulating tfs during lineage differentiation from escs
<cit> .
while the recognition landscape for p <dig> remains largely unchanged during development of all other lineages, there is a significant and selective gain of p <dig> binding elements in the dhss of the human skin keratinocytes, which represent the ectodermal lineage.
such deterministic function of p <dig> is further evident by data showing that p <dig> in combination with klf <dig> can efficiently convert human fibroblast into keratinocytes
<cit> .
these interesting correlative findings, together with our results presented here strongly suggest that p <dig> functions as part of a pioneering complex which can target and remodel chromatin at many of its sites.
binding of p <dig> in coordination with cooperating tfs
it is likely that the p63-depedent regulation of target genes in keratinocytes requires co-operation of other tfs.
we have used a multi-pronged approach to identify such p63-associated cooperating tfs by processing data from nhek rna-seq, available encode chip-seq and computational prediction methods based on tf motifs.
our analysis led to a few surprising observations about the identity of candidate tfs that were likely to be involved in p63-genomic interactions.
one striking result from our study is that many of the p63-associated factors belong to broadly expressed family of tfs such as ap <dig>  ap <dig>  myc and stat rather than highly tissue-specific factors.
although at first glance, this result may seem disappointing, we think that given the master regulatory role of p <dig>  such a finding makes biological sense.
indeed, given the fact that p <dig> is highly expressed in a lineage-restricted fashion and plays a crucial role in dictating keratinocyte cell fate, it is conceivable that some of the p63-associated cooperating tfs might just provide ancillary role in regulating gene expression.
another interesting possibility is that the keratinocyte-specific gene expression is mediated by a combinatorial interaction of multiple tfs as suggested by prior studies
<cit> .
however it is important to stress that many of the broadly expressed tfs such as ap <dig> and ap <dig> do have keratinocyte-specific roles that are often masked due to functional redundancy from expression of multiple family members
<cit> .
future functional studies on these tfs that are part of the p63-driven transcriptional network, including ones that are relatively under studied such as rfx <dig> will shed important insights into gene regulatory mechanisms in keratinocytes.
CONCLUSIONS
despite the wealth of information obtained from our data-mining studies, long term follow-up experimental studies are needed to better elucidate the p <dig> tf network and the role of chromatin in regulating myriad biological functions of p <dig>  unraveling the complex nature of the distal regulatory elements such as enhancers, which are by far the most common sites of p <dig> binding is a formidable challenge.
the new insight into the dynamic interplay between p <dig>  its many cooperating tfs and the local chromatin environment, as reported here is the first step towards tackling such challenges.

 tp is an important transcription factor of the p53/p63/p family that dictates a wide range of cellular properties including but not limited to stem cell renewal, lineage choices and maintaining the balance between proliferation and differentiation. this diverse function of p is critical for morphogenesis during development, particularly for epithelial-enriched tissues such as the skin and its appendages such as the hair follicles and mammary glands. indeed, p63-null mice die after birth and exhibit a dramatic agenesis of epithelial-rich structures and widespread developmental defects of the limb, orofacial region, and external genitalia. these p63-deficient structural defects are thought to be the result of a failed program of epithelial stratification and/or diminished capacity for stem cell renewal, both of which can jeopardize normal epithelial-stromal interactions needed during embryonic organ development. in agreement with the mouse phenotype, p mutations in humans lead to congenital abnormalities such as abnormal limb development and ectodermal dysplasia, which are associated with a spectrum of developmental disorders including aec or eec syndrome. the biological function of p is mediated by several isoforms derived from distinct transcripts. these include the longer tap isoforms and n-terminal deleted Δnp isoforms generated from an internal promoter located within intron furthermore alternative splicing can result in α, β, and γ isoforms, which differ in the c-terminus. all p isoforms share the dna-binding and oligomerization domains, which are analogous to that of p it is now well-established that ∆np especially ∆np63α is the predominant isoform that is present in most epithelial cells such as the keratinocytes of the skin. importantly both gene complementation studies and isoform specific knockouts have conclusively affirmed that ∆np harbors most of the function and biological activity of p particularly as it pertains to the epithelial tissues. the role of p in regulating transcription during development has been extensively studied in skin where ∆np is highly expressed and regulates the transition from simple ectodermal cells to stratified epithelium. given the master regulatory function of p it is not surprising that the repertoire of p63-targets is vast and represents practically every crucial gene regulatory and signaling pathway. this is evident from the ~ binding sites for p in human keratinocytes as determined by chromatin immunoprecipitation with next-generation sequencing studies. p controls expression of basal keratin genes k and k and regulates myc levels thereby controlling keratinocyte proliferation via the wnt/β-catenin and notch signaling pathways. the keratinocyte differentiation program is also regulated by p in part via its effect on the znf750-klf regulatory axis. while the identification of p bound cis-regulatory elements in keratinocytes has received much attention, the mechanics of p63-dna interaction is still relatively unknown. p binds a canonical motif, defined as closely spaced decamers, although there is growing evidence that p can target sites that do not completely conform to this consensus sequence, including half sites. given the degenerate nature of the p binding motif, it is not surprising that by conservative estimates, there are more than million such potential sites in the human genome. however, as is the case with most other transcription factors, only a small subset of these sites are bound by p in vivo. it is likely that the local chromatin architecture, among other factors plays an important deterministic role in dictating how and why p selects its target dna. hence, this is an important area of future investigation; especially given the increasing evidence that p can play an important role in modulating the chromatin structure. indeed, recent studies have demonstrated that p can functionally interact with several epigenetic factors in keratinocytes, which can in turn profoundly influence p63-dependent transcriptional activation and repression. examples of such interactions include the reinforcement of p mediated repression of p by lsh, a member of the snf family of chromatin remodeling atpases, direct recruitment of histone deacetylases, hdac and hdac by ∆np during repression of target genes in the embryonic epidermis and the crosstalk between p and chromatin organizer satb in regulating keratinocyte differentiation genes. p can also control higher-order chromatin structure in epidermal progenitor cells during skin development by regulating brg a atp-dependent chromatin remodeler. given these emerging links between p and chromatin, it is important that any comprehensive studies on the mechanism of p63-genome interactions takes into account the underlying state of epigenetic modifications. here we have utilized the p chip-seq dataset and available chromatin modification datasets for normal human epidermal keratinocytes to investigate the rules that govern binding of p to its target dna. we find that p binds to a canonical motif at the majority of its sites, whereas non-canonical motifs containing 1– spacer between decamers are present in only % of the sites. the chromatin at p binding sites is largely marked by active histone modifications. moreover, chromatin accessibility with h3k4me can accurately predict bona-fide bound p sites without the need for any additional dna sequence information. finally, using a comprehensive in silico approach, we identify several cooperating tfs that appear to define specific classes of p regulated genes.  underlying sequence patterns and chromatin architecture of p targets several groups have determined global p binding locations in various primary and immortalized keratinocytes using chip-chip or chip-seq techniques. for our studies, we focused on the most comprehensive p chip-seq data available to date. it had an added benefit of being generated from primary keratinocytes and more importantly conforming to encode guidelines. to facilitate uniform comparisons across other encode datasets, we re-aligned the p chip-seq to the latest human genome build with bowtie. in strong agreement with kouwenhoven et al., by using high stringency conditions, we identified a reliable and robust dataset of p binding sites that were common among the three biological replicates. on examining the underlying dna sequence of these p63-chiped elements, we found that % of these sites have at least one p canonical motif. among these, 32% show a close match to the p consensus while the remaining % are a weaker match. both the strong and weak canonical motifs are significantly enriched at the p chiped regions compared to random genomic regions. an additional % of the binding sites show a close match to the non-canonical p motif, which has 1– base spacers separating the two half-sites. interestingly, 10% of the p bound genomic sites do not have a recognizable p motif raising the possibility that p can perhaps also be recruited to target regulatory elements through indirect mechanisms such as protein-protein interactions. our analysis also revealed that only a few of the p binding sites contained just a half site. next we compared the chromatin architecture of the p bound sites and randomly selected unbound sites, which have at least strong p canonical motif. for this purpose we focused on active and repressive histone modifications profiles in nhek cells that have been generated by the encode consortium. as shown in figure  1a, the genomic segments of dna bound by p are characterized by a chromatin architecture consisting of high h3k27ac, h3k4me h3k4me and h3k4me and low h3k9me and h3k27me modifications. in contrast, the randomly selected stretches of dna corresponding to the unbound p sites are completely lacking in these chromatin profiles. hence chromatin architecture could contribute to selective targeting of p figure chromatin architecture differs between p bound and unbound sites. average profiles of the active and repressive histone modifications are plotted for a 1 kb window centered at p motif. average chromatin architecture at p bound locations. average chromatin architecture at unbound locations that have a strong p motif. profiles were generated with archtex. p has been shown to target different types of regulatory regions and involved in both activation and repression of gene expression. recent studies have demonstrated that specific states of chromatin modifications at the regulatory regions are strongly associated with the level of gene expression for the corresponding genes. to determine whether such differences exists for distinct classes of p targets, we clustered the p bound regulatory regions by their underlying chromatin architecture. we first divided p targets into two groups based on the magnitude of signal and then performed unsupervised clustering using the spatial arrangement of the histone modifications. using standard k-means clustering we were able to capture all patterns within cluster a and b but this resulted in redundancy. therefore, we compared all sub-clusters to each other in both orientations, grouping together those with a pearson correlation above. this resulted in the classification of p targets into distinct, non-redundant groups. based on chromatin based segmentation, the majority of p targets are predicted to represent enhancers while only represent active promoters. the distance of the regulatory regions to the nearest tss supports these predictions. it is interesting to note that none of the clusters are characterized by repressive histone marks. however this does not imply that p cannot act as a repressor. in fact the largest sub-cluster b has very low signal for active chromatin marks and predicted by chromatin based segmentation likely to represent repressed/inactive targets. indeed, upon probing the rna-seq data from keratinocytes, it is evident that the regulatory regions of the b sub-cluster are associated with genes that are weakly expressed. therefore this cluster might represent repressed/poised genes for which p may not play an active role in activation at least in the basal growth conditions of nhek cells.figure clustering of p targets by chromatin intensity and shape. p bound locations are clustered by active and repressive histone modifications sequenced in nhek cell-line. they are first clustered by average signal intensity across a 1 kb window centered at p binding site and then clustered by spatial arrangement of the histone modifications. annotating p targets clustered by chromatin profiles p bound locations are clustered by active and repressive histone modifications sequenced in nhek cell-line. the groups are annotated by chromatin based segmentation and median distance to nearest tss. great is used to find the top go-bp terms and top panther pathways that enrich for each cluster. we next used genomic regions enrichment of annotations tool to determine whether the five clustered groups of p targets could be segregated into distinct classes of genes involved in specific biological pathways. our analysis revealed that each cluster indeed was overrepresented by genes that were involved in closely related, yet disparate biological activity. for example, cluster b which primarily encodes for strong enhancers is enriched for gene-ontology biological process of epidermis development, while cluster b representing repressed/inactive sites is enriched for go-bp of hair cycle. this raises the intriguing possibility that p might play an important role in actively driving the epidermal developmental processes while keeping the hair cell fate repressed – a notion that is supported by data from prior transgenic mouse studies. the complete annotation for each group of p targets is provided in table  thus, collectively the data obtained by chromatin architecture-dependent clustering of p targets allowed us to determine both the active and poised targets of p and annotate these to specific biological processes. p binding can be accurately predicted from chromatin modifications we next implemented two different statistical approaches to test whether epigenetic modifications can be exclusively utilized to identify functional p binding sites. we generated a discriminant and a regression model of p63-binding events based on data from histone modifications as determined by chip-seq and chromatin accessibility as inferred from dnase-seq and faire-seq experiments. dnase-seq and faire-seq assays generate both distinct and overlapping accessibility information about the genomic landscape. also the different histone modifications usually correspond to transcriptional regulation in a combinatorial fashion. as these various chromatin features are not mutually independent of each other, we created interaction datasets from each pair of chromatin datasets, such as dnase*faire. these were also used as predictors for our statistical models. for the discriminant approach we generated the model on a training dataset containing p binding sites with unbound genomic locations and for the regression approach we generated the model on a training dataset containing p binding sites with unbound genomic locations. we then validated the models on a test dataset of genomic locations. the discriminant model identifies boundaries between groups of objects, the boundaries being defined in terms of those variable characteristics that distinguish the objects into two groups. this technique allows us to determine which chromatin modifications are strongly predictive of p binding status. in addition this model can be used to predict if an unknown site is bound or unbound by p we found that a three variable model was able to accurately classify p binding sites, specificity: %). this model was based on df, h3k4me and h3k27me datasets. in this model both df and h3k4me had positive coefficients as anticipated by p binding to sites with high values for these chromatin datasets. on the other hand, the repressive modification h3k27me had a negative coefficient suggesting that p is being excluded from genomic sites with this histone modification.figure p occupancy can be accurately predicted by chromatin marks. active and repressive histone modifications with chromatin accessibility as measured by dnase-seq and faire-seq were used to construct models of p binding. discriminant model for p binding. regression model for p occupancy. the values for datasets used for our models are shown with the data sorted in descending order by experimental p occupancy. the standardized number of tags for the histone modifications and chromatin accessibility are displayed. the dnase and faire interaction term is calculated as product of the two and is displayed in standardized tag space. the quality of the p sequence motif was determined by patser, good motif matches are purple and low quality motif matches are pink. the results of our discriminant model is compared directly to the experimental determination of binding by macs. our regression model is compared directly to experimental occupancy. we also utilized the regression model to determine which chromatin characteristic best predicts the occupancy of p binding at any genomic location. when applied to the test dataset, the regression model with only three variables has a mse of and r of. for the regression model h3k4me dnase, and df variables were the most informative. for both models, a chromatin accessibility dataset with h3k4me were found to be very predictive. furthermore, a regression model using only the df term has a r2 = , which exemplifies the importance of chromatin accessibility in defining true in vivo binding sites. however these models do not inform whether the accessibility is a prerequisite or a secondary effect of p binding. we therefore integrated the nucleosome-dna interaction model in our analysis, which is based on the principle that certain dna sequences show a greater predisposition to wrap around the histone octamer and form nucleosomes. this model has previously been used to predict and experimentally validate that p preferentially binds nucleosome rich regions. our analysis showed that the p63-bound genomic segments in keratinocytes have higher sequence defined nucleosome occupancy than randomly chosen sequences. this suggests that the increased chromatin accessibility at p bound sites is likely to be actively shaped by p63-dna interactions.figure p preferentially binds nucleosome rich regions. predicted nucleosome occupancy as calculated by the nucleosome-dna interaction model, was plotted for a 1 kb window centered at p binding site. blue: average nucleosome occupancy at p bound locations. red: average nucleosome occupancy at random unbound locations. for both the regression and discriminant models, p motif data was not included even though p motif is a statistically significant predictor. indeed on its own the sequence information accounts for 14% of the variability in p binding. we therefore performed additional analysis to ascertain the importance of p63′s motif in dictating dna binding, when taking into account epigenetic information. we divided all accessible genomic locations as determined by dnase-seq into three groups based on the presence of a strong, weak, or the absence of a p motif. our results indicate that the presence of a strong or weak p motif leads to stronger p occupancy. this finding that p dna-sequence motif is an important determinant of in-vivo binding is difficult to reconcile with the fact that p motif was not considered a strong predictor in our statistical models. one possible explanation for this discrepancy is that the chromatin datasets we use in our models have already embedded in it the dna sequence motif information. indeed, using our regression model, we predicted p occupancy for each group and showed that predicted p occupancy is also higher for sites containing a strong or weak p motif. this provides further support for the notion that h3k4me dnase and df variables in our statistical models already account for the information represented by the p consensus motif. thus accessible sites with a strong or weak motif are more accessible and have higher deposition of h3k4me than sites lacking a p motif. taken together, these results lead us to postulate that p might exist as part of a chromatin remodeling complex, which creates a distinct epigenetic architecture at its binding sites.figure significance of p motif. the genomic coordinates of all the accessible locations in human keratinocytes, as determined by dnase-seq in nhek cell-line, were obtained from encode via the ucsc genome browser. the accessible locations are divided into groups 1) locations with strong p motif, 2) locations with weak motif and 3) locations with no motif. box and whisker plots are made across 500 bp window for experimental p occupancy, predicted p occupancy by our best regression model, dnase tag density, h3k4me tag density. identification of p cooperating tfs it is well established that many tfs often act in a combinatorial fashion to govern tissue-specific gene expression. hence, we wanted to examine the repertoire of p associated tfs that might play such role in modulating p binding and possibly influencing p63-dependent gene expression in keratinocytes. we therefore used a combination of in-silico analysis of the p binding sites and careful data mining of large-scale genomic datasets such as rna-seq and chip-seq from the encode project. we posit that the tf that are likely to directly cooperate with p will have the following features. first, their dna binding motifs will likely be enriched at p bound elements, keeping in mind that co-occurrence of such motifs adjacent to p sites is not a prerequisite for such interactions. second, their in-vivo binding profiles will overlap with p binding profile. here we reasoned that tf binding to some extent could be extrapolated from chip-seq in other cell types, as it has been shown that tfs share a large number of common binding sites across different cell lines. finally, we contend that the relevant tf should be expressed in keratinocytes. applying this criterion, we examined tfs for which chip and/or expression data was available from the encode project. of these, were enriched at p targets as determined by their in-vivo binding profile or in silico motif analysis. however, we rejected of these tfs as potential cooperating tfs because they were either not expressed in keratinocytes or not enriched by both chip-seq and motif analysis. an additional tfs were also not considered from further evaluation due to ill-defined/unknown dna binding motifs or missing chip-seq data. this shortlisted the potential p cooperating tfs to tfs, which had both an enriched motif and in-vivo binding signal at p bound locations in nhek cells. these tfs include: cebpb; cfos, fosl and jund; bach; tfap2c; stat and stat; max, c-myc, and usf, rfx; elk. in this context, it is important to note that several of the predicted tf discovered by our analysis, such as ap ap myc, stat and c/ebp have been previously linked to p and keratinocyte biology. on the other hand, some of the other potential p cooperating tfs such as rfx elk and bach are new members of this class. they represent an interesting group for which much less is known about their expression and function in keratinocytes and any possible correlation with p one particularly intriguing p cooperating tf is rfx which is highly expressed in nhek cells as evident by the rna-seq data. furthermore, a recent study has also shown that rfx mrnas and protein are preferentially expressed in the p63-rich basal cells of the human epidermis further raising the prospects of a functional interplay between these two tfs.figure cooperating tfs of p tfs corresponding to motifs, bach ap stat stat mycmax, usf, rfx, elk1) were determined as potential cooperating tfs for p the characteristics that determine these tfs as potential cooperating tfs are: 1) enrichment of tf binding at p targets - as measured by chip-seq from any available cell-line, 2) high expression in keratinocytes - as determined by rna-seq in nhek cell-line, 3) enrichment of the tf’s motif at p targets and co-occurrence with p motif. the tfs are sorted in descending order by an average percentage rank calculated across the different analysis. the potential role of p in modulating different biological functions in coordination with various cooperating tfs p is involved in myriad biological functions and it is possible that this diversity of p function is brought about by in part by different subtypes of p63-tf complexes. we therefore constructed a correlation matrix of the tfs based on their in-vivo binding signal and their dna-binding motif scores at the p bound locations. the tfs correspond to unique dna-binding motifs. there appears to be multiple complexes with unclear distinctions between them. two distinct complexes containing overlapping tfs could be identified. the first group comprised of c/ebpb, jund, cfos, and stat whereas the second complex was primarily represented by elk fosl and myc-max. we next asked if there were different sub-groups of genomic targets where p was bound in coordination with these distinct cooperating tfs. upon clustering the p targets, clear groups emerged based on the presence/absence of the cooperating tf motifs. ap stat and mycmax seem to be the driving force for the clustering, whereas cebpb, ap and elk motifs were quite ubiquitous in their presence. further annotating the clusters using great allowed us to associate specific groups of p targets with distinct biological functions, based upon top go terms and significant pathways. the complete annotation for each group of p targets and the type of complex it represents is shown in figure  7b.our clustering analysis revealed an interesting and complex correlation between p its network of associated cooperating tfs and the nature and physiological role of the p target genes. case in example is ap it seems to be an important co-player of p in regulating genes involved in keratinocyte differentiation and epidermis development, modestly linked with hair follicle development and cell cycle arrest but not at all involved in the activation of mapkk activity or regulation of carbohydrate metabolic process. similarly, while stat and stat seem to share similar co-regulatory function as evident by their enrichment in clusters c and c there are also clear differences between these two closely related family members based on their enrichment patterns for clusters c c and c this point is further illustrated by the correlation matrices made in figure  7a. in contrast, binding motifs for c/ebp, elk and ap did not seem to exhibit a preference for any of the particular clusters representing different go terms or biological pathways. interestingly, none of the clusters showed any difference in terms of p motif strength or the type of regulatory regions.figure p functions with multiple cooperating tfs. a correlation matrix for the probable p cooperating tfs at the p locations. upper triangle shows pearson’s correlation coefficient calculated between chip-seq tag densities for each pair of tfs. for bach and fosl chip-seq signal are from k and a cell-lines respectively. while for the rest of the factors, chip-seq signal is from hela-s cell-line. lower triangle shows correlation between motif scores generated by patser. the p bound locations are clustered by the presence or absence of the cooperating tfs motifs, in a 1 kb window centered at p binding site. one obvious limitation of our approach is the fact that for the purpose of the aforementioned studies, a specific p63-binding site was assigned to be regulatory element for the nearest gene, often without any supporting experimental data. however it has been shown, that linking distal enhancer elements to the nearest gene, as is the norm, might not always be an accurate representation of the regulatory landscape. there is an absence of chromatin interactions data such as those obtained from chromosome conformation capture, in keratinocytes. we therefore virtually linked the distal regulatory elements to their target genes, by exploiting the characteristic that regulatory regions become dnase hypersensitive in synchrony with their target promoters. to examine this further, we first retrieved the global map of distal dhs-to-promoter connections as generated by the encode consortium. overlapping this map with the p chip-seq data allowed us to link p bound distal genomic regions to their putative target promoters/genes. further analysis revealed interesting facets of long-range interactions between distal p63-bound elements and the promoter as evident by the krt gene, a known p target. while prior studies have focused primarily on the proximal regions 5′ of the krt gene, including a well-characterized p63-bound enhancer ~  kb upstream, we discovered novel distal p63-bound elements that are located quite far away as predicted by the dhs-to-promoter connections. two of these p63-bound genomic segments, which are predicted to be enhancers for krt are shown in figure  8a. these regulatory regions also contained potential binding sites for some of the p cooperating tfs such as cebpb, c-myc and ap interestingly these two distal p63-bound elements are not exclusive to the krt gene since they are also putatively linked to related krt16/ krt genes, which are located in a relatively closely spaced genomic cluster. this raises an interesting possibility that a specific p63-bound regulatory region might be commonly utilized to coordinate the regulation of multiple closely related genes in keratinocytes.figure model of p binding at distal regulatory elements of keratin- and rfx histograms in green show p tag density in a 1 kb window at the regulatory region. loops links the distal regulatory regions to their target promoter as mapped by the encode project. p binding sites  kb and  kb upstream are linked to keratin-. these regions also regulate other genes, among them keratin- and keratin-. motif analysis reveals other tfs binding within the 1 kb window. p binds  kb upstream to rfx at a distal regulatory element linked to rfx promoter. this region is also linked to tuft and rorc promoters. motif analysis shows p binding in complex with ap e-box binding factors and bach on a weak canonical motif. additional motifs are also found within the 1 kb window. we also discovered possible new target genes under the control of p63-bound distal regulatory regions that otherwise would have been missed by the conventional strategy of assigning a p63-regulated element to its closest gene. as shown in figure  8b, one such interesting candidate is rfx a novel cooperating tf of p as discussed above. p binds a regulatory element  kb upstream of the rfx gene and motif analysis suggests that ap e-box binding factors and bach are part of the transcriptional complex at this site. hence, there might exist a transcriptional regulatory loop whereby p in cooperation with additional tfs activates rfx which then in turn, modulates p binding to its target sites. similar to the case with the krt gene, this distal regulatory element for rfx is predicted to be also linked to the tuft and rorc genes. in the future in-depth studies such as 3c experiments will help to confirm these distal enhancer-promoter interactions and to firmly establish the true identity of the p63-driven gene network in keratinocytes. discussion this study aimed to decipher the mechanics of p binding by determining the minimal in-vivo motif required for binding, distinguishing between chromatin architecture of bound and unbound motifs and identifying cooperating tfs that modulate p biological activity. p requires a full site for binding to target sites the transcription factor p binds as a homotetramer to two decamers rrrcrwgyyy, rrrcwygyyy separated by a 0– base pair spacer region. we found that the p sites containing the decamer pair with intervening spacers was much less prevalent than those where the two half sites were juxtaposed to each other. a very small subset of p binding sites consisted of only a half site. this is in contrast to a similar study that found 3-4% of p binding sites having a half site. this observation can be explained by a slight difference in the consensus motifs of the two factors, which also results in fold lower binding affinity of p in comparison to p53. one possibility is that the dimer-dimer interactions are important for p dna binding specificity and therefore p requires a full site to bind dna efficiently. indeed, such differences in the dna-protein interactions among p family members are quite evident from recent structural studies with p73. interestingly, not only the distance between two p half-sites influences the p quaternary structure, but tellingly transcriptional activity is also more affected by spacer length in p response element than in p finally, it is worth noting that ~10% of p chiped sites in keratinocytes do not have a recognizable p binding sequence suggesting that the p binding at these sites is driven by indirect mechanisms that might involve other dna-binding tfs and/or non-canonical p motifs. the p consensus motif is not required for predicting binding events in keratinocytes we utilized two computational modeling approaches to uncover key characteristics defining p binding sites. first, we utilized discriminant modeling that allowed us to predict p binding as a binary score. second, in a parallel strategy, we used regression modeling that predicted the degree of p occupancy. our approach was distinct from other published methods in that we trained our statistical models on random genomic sites which might or might not have p binding sequences. this allowed our models to include the 10% binding sites that do not have p dna-binding motif. surprisingly, we found that our final models constituting only chromatin marks predicted p binding with high accuracy sensitivity, regression model – mse). adding sequence information to the models, did not lead to any significant improvement even though p motif is a statistically significant predictor. one possible explanation is that in keratinocytes, regions of the genome that have a functional p motif are on average more accessible and marked with active chromatin marks. this result is not surprising if p is a key component of the regulatory complex that is involved in remodeling the chromatin at its binding sites. indeed p target sequences dictate higher nucleosome occupancy than random genomic sequences according to the nucleosome-dna interaction model. it can be then hypothesized that p binding shifts the nucleosomes creating an accessible active chromatin structure at its targets. support for such a role for p comes from a recent study, which examined the selective loss versus gain of dhss targeted by lineage regulating tfs during lineage differentiation from escs. while the recognition landscape for p remains largely unchanged during development of all other lineages, there is a significant and selective gain of p binding elements in the dhss of the human skin keratinocytes, which represent the ectodermal lineage. such deterministic function of p is further evident by data showing that p in combination with klf can efficiently convert human fibroblast into keratinocytes. these interesting correlative findings, together with our results presented here strongly suggest that p functions as part of a pioneering complex which can target and remodel chromatin at many of its sites. binding of p in coordination with cooperating tfs it is likely that the p63-depedent regulation of target genes in keratinocytes requires co-operation of other tfs. we have used a multi-pronged approach to identify such p63-associated cooperating tfs by processing data from nhek rna-seq, available encode chip-seq and computational prediction methods based on tf motifs. our analysis led to a few surprising observations about the identity of candidate tfs that were likely to be involved in p63-genomic interactions. one striking result from our study is that many of the p63-associated factors belong to broadly expressed family of tfs such as ap ap myc and stat rather than highly tissue-specific factors. although at first glance, this result may seem disappointing, we think that given the master regulatory role of p such a finding makes biological sense. indeed, given the fact that p is highly expressed in a lineage-restricted fashion and plays a crucial role in dictating keratinocyte cell fate, it is conceivable that some of the p63-associated cooperating tfs might just provide ancillary role in regulating gene expression. another interesting possibility is that the keratinocyte-specific gene expression is mediated by a combinatorial interaction of multiple tfs as suggested by prior studies. however it is important to stress that many of the broadly expressed tfs such as ap and ap do have keratinocyte-specific roles that are often masked due to functional redundancy from expression of multiple family members. future functional studies on these tfs that are part of the p63-driven transcriptional network, including ones that are relatively under studied such as rfx will shed important insights into gene regulatory mechanisms in keratinocytes.  despite the wealth of information obtained from our data-mining studies, long term follow-up experimental studies are needed to better elucidate the p tf network and the role of chromatin in regulating myriad biological functions of p unraveling the complex nature of the distal regulatory elements such as enhancers, which are by far the most common sites of p binding is a formidable challenge. the new insight into the dynamic interplay between p its many cooperating tfs and the local chromatin environment, as reported here is the first step towards tackling such challenges. 
furthermore, we use a comprehensive in silico approach that leverages encode data to identify several known tfs such as ap ap and novel tfs that can potentially cooperate with p to modulate its myriad biological functions in keratinocytes. collectively, our results suggest that p might actively remodel and/or influence chromatin dynamics at its target sites and in the process dictate its own dna binding and possibly that of adjacent tfs. furthermore, ectopic expression of p alone is sufficient to convert simple epithelium into stratified epithelial tissues in vivo and p is required for efficient transdifferentiation of fibroblasts into keratinocytes. but, in human keratinocytes less than of these sites are bound in vivo suggesting that underlying chromatin architecture and cooperating tfs mediate p63-genome interactions. we find that the chromatin architecture at p63-bound targets possess distinctive features and can be used to categorize p targets into proximal promoters, enhancers and repressed or inactive regulatory elements. p which acts both as an activator and repressor of transcription, recognizes a canonical binding motif that occurs over million times in the human genome. our analysis shows that p bound genomic locations in keratinocytes are accessible, marked by active histone modifications, and co-targeted by other developmentally important transcriptional regulators. however, little is known about the molecular mechanisms of p function, in particular how it selects its target sites in the genome.
 the transcription factor p is a master regulator of epidermal development and differentiation as evident from the remarkable skin phenotype of p mouse knockouts. furthermore, ectopic expression of p alone is sufficient to convert simple epithelium into stratified epithelial tissues in vivo and p is required for efficient transdifferentiation of fibroblasts into keratinocytes. however, little is known about the molecular mechanisms of p function, in particular how it selects its target sites in the genome. p which acts both as an activator and repressor of transcription, recognizes a canonical binding motif that occurs over million times in the human genome. but, in human keratinocytes less than of these sites are bound in vivo suggesting that underlying chromatin architecture and cooperating tfs mediate p63-genome interactions.  we find that the chromatin architecture at p63-bound targets possess distinctive features and can be used to categorize p targets into proximal promoters, enhancers and repressed or inactive regulatory elements. our analysis shows that the chromatin modifications h3k4me h3k27me along with overall chromatin accessibility status can accurately predict bonafide p63-bound sites without a priori dna sequence information. interestingly, however there exists a qualitative correlation between the p binding motif and accessibility and h3k4me levels. furthermore, we use a comprehensive in silico approach that leverages encode data to identify several known tfs such as ap ap and novel tfs that can potentially cooperate with p to modulate its myriad biological functions in keratinocytes.  our analysis shows that p bound genomic locations in keratinocytes are accessible, marked by active histone modifications, and co-targeted by other developmentally important transcriptional regulators. collectively, our results suggest that p might actively remodel and/or influence chromatin dynamics at its target sites and in the process dictate its own dna binding and possibly that of adjacent tfs. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords p63chromatinchip-seqtranscriptionkeratinocyteencodeissue-copyright-statement© the author 2014 