Characterization of chicken p53 transcriptional function via parallel genome-wide chromatin occupancy and gene expression analysis

The tumor suppressor p53, which acts primarily as a transcription factor, can regulate infections from various viruses in chickens. However, the underlying mechanisms of the antiviral functions of chicken p53 (chp53) remain unclear due to the lack of detailed information on its transcriptional regulation. Here, to gain comprehensive insights into chp53 transcriptional regulatory function in a global and unbiased manner, we determined the genome-wide chromatin occupancy of chp53 by chromatin immunoprecipitation, which was followed by sequencing and chp53-mediated gene expression profile by RNA sequencing using chemically immortalized leghorn male hepatoma (LMH) cells with ectopic expression of chp53 as the model. The integrated parallel genome-wide chromatin occupancy and gene expression analysis characterized chp53 chromatin occupancy and identified 754 direct target genes of chp53. Furthermore, functional annotation and cross-species comparative biological analyses revealed the conserved key biological functions and DNA binding motifs of p53 between chickens and humans, which may be due to the consensus amino acid sequence and structure of p53 DNA-binding domains. The present study, to our knowledge, provides the first comprehensive characterization of the chp53 transcriptional regulatory network, and can possibly help to improve our understanding of p53 transcriptional regulatory mechanisms and their antiviral functions in chickens.


INTRODUCTION
The tumor suppressor p53, is well known as 'the guardian of the genome' and is a transcription factor that controls a complex signal transduction network referred to as the p53 pathway (Vogelstein et al., 2000;Harris and Levine, 2005;Levine and Oren, 2009;Lawrence et al., 2014). Among the various biological processes of the p53 network, p53 has been shown to control many aspects of host immune responses, such as innate immunity, nitric oxide-induced apoptosis, and histocompatibility complex I expression (Albina et al., 1993;Messmer and Brune, 1996;Zhu et al., 1999;Thiery et al., 2005;Xue et al., 2007;Lujambio et al., 2013;Schwitalla et al., 2013;Wang et al., 2013), all of which are important defence mechanisms of the host against viral infection. To date, the antiviral function of p53 has been confirmed in many viruses, such as vesicular stomatitis virus (Munoz-Fontela et al., 2005;Garijo et al., 2014), poliovirus (Pampin et al., 2006), hepatitis C virus (Dharel et al., 2008), and influenza A virus (Yan et al., 2015). For example, p53 plays an essential role in enhancing the type I IFN-mediated immune response against Influenza A virus (IAV) infection. Knockdown of p53 expression by RNAi enhanced IAV replication, which is associated with reduced expression of antiviral IFN-stimulated genes (ISGs), such as IRF7, IRF9, ISG15, ISG20, GBP1, RIG-I, and OAS1. In addition, pretreatment of p53-knockdown cells with IFN-a failed to inhibit IAV replication, showing impaired antiviral activity (Zhu et al., 2014). In chickens, Marek's disease virus must inhibit chicken p53 (chp53)-mediated transcriptional activity and apoptosis by interacting directly with chp53 via its Meq oncoprotein to replicate successfully (Deng et al., 2010). Efficient inhibition of the replication of J subgroup avian leukosis virus (ALV-J)  and infectious bursal disease virus (Ouyang et al., 2017) by chp53 has been reported, which provides a direct link between chp53 function and the antiviral immune response in chickens. Our previous study also provided experimental evidence for the antiviral effect of chp53, which identified chp53 as the key determinant of avian infectious laryngotracheitis virus (ILTV) infection by suppressing ILTV replication in infected cells and protecting uninfected standby cells from ILTV-induced paracrine-regulated apoptosis, which reduces pathological damage and delays viral transmission (Li et al., 2018). However, the underlying molecular mechanisms by which chp53 represses infections of these avian viruses remain unclear.
At the molecular level, p53 acts primarily as a transcription factor, inducing a network of genes that contributes to its biological responses, including immune responses (Bieging and Attardi, 2012). Currently, p53 has been identified as a transcription factor for several important antiviral immune genes, such as genes encoding retinoic acid inducible gene I (Hsu et al., 2012), a ligand for NK-cell activating receptor NKG2D, Toll-like receptor 3, and NKG2D ligands ULBP1 and ULBP2 (Taura et al., 2008;Textor et al., 2011;Shatz et al., 2012). The transcriptional regulatory networks of human and mouse p53 proteins and their target genes have been extensively studied Nikulenkov et al., 2012;Li et al., 2014;Jiang et al., 2015;Zhang et al., 2017), which demonstrates that the pleiotropic character of the p53 network makes it difficult to predict the consequence of p53 activation. Thus, discovery of the mechanisms governing p53mediated biological outcomes is extremely important for application of p53-based therapies . In chickens, the global characterization of the p53 transcriptional regulatory network is still lacking, which limits the in-depth illustration of the antiviral function of chp53.
In this study, to gain comprehensive insights into chicken p53 transcriptional regulatory function, we performed chromatin immunoprecipitation sequencing (ChIP-seq) and RNA sequencing (RNA-seq) to characterize p53 transcriptional responses in a global and unbiased manner. We characterized the important features of p53 chromatin occupancy and identified novel target genes.

Cell Culture
The leghorn male hepatoma (LMH) cells (ATCC CRL-2117) were maintained in Dulbecco's modified Eagle's medium supplemented with 10% fetal bovine serum, 100 units/mL penicillin, 100 mg/mL streptomycin, and 2 mM L-glutamine. Cell cultures were incubated at 37°C in 5% CO 2 . After the cells are recovered, they are passed to the fourth passage and then plated. For RNA-seq and ChIP-seq, LMH cells were seeded in 100 mm dishes 12 h before transfection. Cells were harvested 24 h after transfection.

Plasmids and Transfection
Chicken p53 cDNA was subcloned into the expression vectors p3xFLAG-CMV-7.1 to generate the recombinant plasmid Flag-chp53. These 2 plasmids were kindly provided by Zhiyong Ma and Yafeng Qiu (Shanghai Veterinary Research Institute, Chinese Academy of Agricultural Science, Shanghai, China). The LMH cells were plated onto tissue culture plates 12 h before transfection. Transfection was performed using the Turbofect transfection reagent (R0531, Thermo Scientific, Rockford, IL) according to the manufacturer's instructions.The transfection dose of both plasmids was 1 mg/1.0 £ 10 6 cells.

Flow Cytometry and Immunofluorescence
Fluorescence-activated cell sorting analyses were conducted using a BD FACScan cell sorter and CellQuest software version 4.0.2 (Becton Dickinson-Pharmingen, San Diego, CA). Cell death was assayed by examining cells in the sub-G1 phase of the cell cycle with propidium iodide-staining of permeabilized cells, as described previously . For immunofluorescence staining, the samples were washed with PBS and fixed with 4% paraformaldehyde for 10 min. After quenching excess aldehyde, the samples were permeabilized with 0.1% Triton X-100. Nonspecific antibody binding was blocked with 2% bovine serum albumin for 1 h, and then the samples were incubated with a mouse anti-DYKDDDDK antibody (flag epitope tag, 1 mg/mL; A00187, GenScript Corporation, Piscataway, NJ) against overexpressed chicken p53 or a mouse IgG2b isotype control antibody (1 mg/mL; 53484, Cell Signaling Technology, Boston, MA), followed by a secondary FITC-conjugated goat anti-mouse antibody (1:64 dilution; F9006, Sigma-Aldrich, St. Louis, MO). All cell nuclei were stained with DAPI (2 mg/mL; C1002, Beyotime Biotech, Jiangsu, China). Fluorescent signals were detected with an EVOS FL fluorescence microscope (AMG, Bothell, WA).

Chromatin Immunoprecipitation Assays
Chromatin immunoprecipitation (ChIP) experiments were conducted as previously described (Wang et al., 2021) with some minor modifications. LMH cells were fixed in 1% formaldehyde for 10 min and then with 0.125 M glycine. DNA was sheared using a sonicator (6 mm; Cole Parmer, Chicago, IL) to an optimal DNA fragment size of 200 to 500 bp. Samples were kept ice cold throughout the sonication. Thirty cycles (30 s pulses/ 30 s rest) were set. Sonication amplitude was 30%. Each ChIP experiment was performed with sheared chromatin samples from LMH cells (5 £ 10 6 cells) using 5 mg of one of the following antibodies: anti-flag or isotype control IgG2b. Protein A/G PLUS−agarose beads were used for pull-down according to the manufacturer's instructions (sc-2003, Santa Cruz Biotechnology, Santa Cruz, CA). The immunoprecipitated DNA was purified using a QIAquick PCR Purification Kit (28106, QIA-GEN, Valencia, CA). Chromatin immunoprecipitation followed by quantitative PCR (ChIP-qPCR) was performed using Luna Universal qPCR Master Mix (M3003L, NEB, Ipswich, MA) following the manufacturer's protocol with a Bio-Rad CFX96 instrument. The primer sequences are shown in Table 1. All samples were measured in triplicate.

RNA Sequencing and ChIP Sequencing
Genome-wide gene expression profiling and genomewide chromatin occupancy of LMH cells was performed via RNA or DNA deep sequencing by Annoroad Gene Technology Co., Ltd. (Beijing, China), according to the manufacturer's instructions. Four biological replicates were performed for both ChIP-seq and RNA-seq experiments. For RNA-seq, RNA was isolated from cells using RNeasy Plus mini kit (74314, QIAGEN, Hilden, Germany) according to the manufacturer's instructions. For ChIP-seq, 10 ng ChIP DNA was prepared for HiSeq 2500 sequencing with TruSeq ChIP Sample Prep Kit (IP-202-1012, Illumina, San Diego, CA) following manufacturer instructions. Libraries for RNA-Seq were prepared with the TruSeq RNA Sample Prep Kits v2 (Illumina) following manufacturer instructions starting from the RNA fragmentation step. For sequencing, paired-end reads were obtained with the length of 150 bp. RNA-seq and ChIP-seq raw data were uploaded to the National Center for Biotechnology Information database under the accession number GSE193188 and GSE200405. Related statistics of all sequenced libraries (RNA-seq and ChIP-seq) have been added as supplementary data (Supplementary Table S3).

RNA-seq Data Analysis
RNA sequencing (RNA-seq) data was analyzed with the Galaxy web-based tool (Blankenberg et al., 2010). Gene ontology and pathway analysis was performed with DAVID (gene enrichment analysis using the EASE score with a modified Fisher's exact test P value of < 0.05 as the threshold (Huang et al., 2009).

ChIP-seq Data Analysis
Chromatin immunoprecipitation sequencing (ChIPseq) data analysis was performed with a web-based tool, Galaxy (Blankenberg et al., 2010), and R (http://www.rproject.org/). Venn diagrams were plotted with TBtools . The annotation and visualization of the ChIP peaks were performed using the ChIPseeker package (Yu et al., 2015). For de novo motif discovery analysis, we used the MEME program (Bailey et al., 2015). As input data, we used sequences at a distance § 50 bp from the position of the peak maximum height.
The transcription factors (TFs) of the putative chp53 direct target genes were predicated using oPOSSUM v 3.0 (Ho Sui et al., 2005) based on the DNA sequences in the interval of § 100 bp of the ChIP-seq peaks we identified. The parameters were as follows: conservation cutoff = 40%; matrix score threshold = 85%; results sorted by Fisher score; amount of upstream/downstream = 5000/ 5000. Protein interaction information was acquired using the STRING database (Szklarczyk et al., 2019). Cytoscape software (Shannon et al., 2003) was applied to visualize the protein interaction networks.

RNA Extraction and Reverse Transcription-Quantitative PCR
Total RNA was extracted from cells using TRIzol reagent (Takara, Dalian, China). The quality and concentration of total RNA were measured using a microvolume spectrophotometer (Implen GmbH, Munich, Germany). A One-Step SYBR PrimeScript RT-PCR Kit (RR086A, Takara, Tokyo, Japan) was used for reverse transcription-quantitative PCR (RT-qPCR) according to the manufacturer's instructions. Primers sequences were designed using Oligo 7 software (version 7.6.0, Molecular Biology Insights, Inc., Colorado Springs, CO). The primer sequences are shown in Table 2. All samples were measured in triplicate. The relative quantification of target gene expression was calculated by 2 ÀDDCt method.

Western Blot Analysis
Western blot was performed strictly according to previously described procedures (Wang et al., 2019). Briefly, cells were washed with ice-cold PBS and soluble proteins were extracted with cell lysis buffer (100 mM Tris-HCl [pH = 8], 150 mM NaCl, 1% NP-40, phosphatase and protease inhibitor cocktail tablets [S8820, Sigma-Aldrich]) according to the manufacturer's protocol. The protein concentration of each sample was determined by a BCA Kit (P0012, Beyotime Biotech). An equal amount of protein was separated by SDS-PAGE and transferred onto a nitrocellulose membrane (Millipore, Billerica, MA). The membrane was then blocked with 5% non-fat milk for 2 h at room temperature and incubated with primary antibodies overnight at 4°C. Antibodies were obtained from GenScript Corporation (A00187) and Proteintech (1:3000 dilution; b-actin, 60008-1-Ig, Wuhan, China). Signals were visualized using infrared imaging systems (Odyssey CLX, LiCor Biosciences, Lincoln, NE).

p53 Amino Acid Sequence and Structure Analysis
The p53 amino acid sequences were aligned using MegAlign program (DNAstar 7.1, DNASTAR Inc., Madison, WI). Next, Conserved Domain Databasesearch tool in NCBI (https://www.ncbi.nlm.nih.gov/ cdd/) and MEME (Bailey et al., 2015) was used to predict the conserved domains of proteins. Finally, we use the online website Phyre2 (Kelley et al., 2015) to predict the three-dimensional structure of the protein.

Statistical Analysis
The SPSS software package (SPSS for Windows version 13.0, SPSS Inc., Chicago, IL) was used for all statistical analyses. Data obtained from several experiments are reported as the mean § standard deviation (SD).
The significance of differences between 2 groups was determined with two-tailed unpaired Student's t test. For all analyses, a probability (P) value of < 0.05 was considered statistically significant.

Genome-Wide Map of Chicken p53 Binding to Chromatin Upon Its Ectopic Overexpression
Due to the lack of an ideal monoclonal antibody for the ChIP of chicken p53 (chp53), the ectopic overexpression of chp53 was conducted in LMH cells which have been used as an in vitro model to study the viral-host interactions of ILTV in our previous studies (Li et al., , 2018Wang et al., 2019;Qiao et al., 2020;Wang et al., 2021).Chp53 was overexpressed in LMH cells by transfecting a Flag-chp53 plasmid (Deng et al., 2010). The efficiency of overexpression was determined by indirect immunofluorescence ( Figure 1A) and flow cytometry ( Figure 1B) analyses using an antibody targeting Flag, which revealed the highest transfection efficiency at 24 h after transfection. Furthermore, Western blot assays revealed ectopic chp53 expression starting 12 h after transfection, with the highest level of chp53 expression occurring at 24 h ( Figure 1C). In this in vitro model, LMH cells were transiently transfected with Flag-chp53 plasmid or negative control empty vector for 24 h. Chp53 overexpression induced significant levels cell death but not cell cycle arrest, as assayed by fluorescence-activated cell sorting of propidium iodide-stained cells ( Figure 1D).
CDKN1A and MDM2 have been experimentally identified as the direct target genes of chicken p53 (Deng et al., 2010). The amount of overexpressed p53 bound to the transcription regulating regions of these two classic chp53 target genes were quantified by ChIP-qPCR in LMH cells upon ectopic expression of chp53. Input DNA was used as a positive control, and ChIP with an isotype IgG2b antibody was used to determine the level of nonspecific binding of the antibody to DNA. Increased binding of overexpressed chp53 to the transcriptionregulating regions of these 2 genes was observed upon ectopic expression of p53 by ChIP-qPCR ( Figure 1E). The transcriptional activity of chp53 bound to these 2 genes was validated by the increased transcription levels of these 2 genes after chp53 overexpression as assayed by RT-qPCR ( Figure 1F). These results demonstrate that in our in vitro model, overexpressed chp53 could successfully interact with its binding sites on chromatin with an activated transcriptional regulatory function. Using this model, ChIP-seq and RNA-seq were conducted to obtain the genome-wide chromatin occupancies of chp53 in LMH cells and related gene expression profiles ( Figure 1G).
We performed ChIP-seq of chp53 in LMH cells upon chp53 overexpression. Thresholds of P < 0.01, q < 0.01, and fold change > 2 were applied, with 65,803 peaks corresponding to chp53-bound DNA fragments identified upon chp53 overexpression in LMH cells. The locations of chp53-bound DNA fragments around genes are summarized in Figure 2A. The promoter region accounted for 39.89%, 5 0 UTR accounted for 0.59%, 3 0 UTR accounted for 2.07%, exon region accounted for 5.74%, intron region accounted for 28.19%, immediate downstream region, which is the proximal intergenic region (≤ 300 bp) accounted for 1.78%, and the distal intergenic region accounted for 21.73%. Although chp53 binding sites were located in varying parts of the genes, approximately 60% of chp53-bound fragments were located within § 10 kb of the transcription start site (TSS; Figure 2B).
We generated a consensus p53-binding site based on p53-bound fragments by applying the de novo motif discovery algorithm, MEME, for the 65,803 peaks (Evalue = 2.1e-056; Figure 2C), which fits well with the original human p53 consensus site as well as the one refined later (Wei et al., 2006;Riley et al., 2008;Smeenk et al., 2008;Menendez et al., 2009). Figure 2D summarizes the amounts of p53-binding sites in different pileup (pileup height at peak summit) regions. As shown in the Venn diagram, approximately 28.69% of these p53bound genes have multiple p53 binding sites.
For validation purposes, 31 putative chp53 binding sites were selected randomly from each pileup region and the occupancy of chp53 at each site in LMH cells upon ectopic expression of p53 was determined by ChIP-qPCR. As shown in Figure 2E, all these selected putative chp53 binding sites were successfully validated with no significant differences in validation rates among different pileup regions.

Identification of Chicken p53 Target Genes
Next, we explored the effects of the ectopic expression of chp53 on the gene expression profiles in LMH cells via RNA-seq. Bioinformatics analysis identified 923 genes that were differentially expressed upon ectopic expression of chp53 based on the following criteria: 1) P value of < 0.01, 2) q value of < 0.01, and 3) fold change of > 1.5. Hierarchical clustering analysis using these differentially expressed genes (DEGs) indicated efficient clustering of biological replicates for each group of cells ( Figure 3A). Compared with the control cells, 114 genes were downregulated and 809 were upregulated in LMH cells upon ectopic expression of p53 (Supplementary Table S1). For validation, the transcript levels of 13 genes randomly selected from the 923 DEGs were examined by RT-qPCR analysis. The directions of change determined by RT-qPCR detection and transcriptome sequencing analysis were generally identical ( Figure 3B), suggesting a correspondence between these 2 methods in our study. Additional Gene Ontology analysis revealed that the upregulated DEGs were significantly enriched in biological processes such as regulation of the canonical Wnt signaling pathway, apoptosis, gene expression, and cell growth ( Figure 3C). The downregulated DEGs were significantly enriched in biological processes such as signal transduction, cytosolic calcium ion concentration, and the Notch signaling pathway ( Figure 3D).
Analysis of chp53 occupancy showed that chp53 bound to 40,243 sites corresponding to 15,108 genes ( Figures 2D and 4A). Among these genes, the transcription of only 675 genes was induced (83.44% of upregulated DEGs) and 79 genes were repressed (69.30% of downregulated DEGs) upon ectopic expression of p53 ( Figure 4A), suggesting that these 754 DEGs are putative direct chp53 target genes (Supplementary Table  S2). Among these putative direct chp53 target genes, 8 induced and 5 repressed genes were selected randomly for validation. The binding of chp53 to these genes was confirmed by ChIP-qPCR and the transcriptional regulation of chp53 to these genes was validated by RT-qPCR in LMH cells upon ectopic chp53 expression ( Figure 4B), suggesting that these 754 genes are bona fide chp53 target genes. A total of 169 genes that changed their expressions upon ectopic expression of chp53 were not bound by chp53 and most likely represent secondary effects of p53. The transcriptions of the majority of the genes bound by chp53 remained unaffected, suggesting that the binding of chp53 per se is not sufficient to regulate gene expression. Consistent with the finding that the majority of DEGs induced by ectopic chp53 expression were bound by chp53 ( Figure 4A), further functional annotation of the putative chp53 target genes that were induced or repressed by chp53 showed similar results to those of induced or repressed DEGs (Figures 4C and D), indicating direct transcriptional control of these biological processes by chp53. We generated two chp53-binding motifs based on the p53-bound fragments of the chp53-target genes we identified (754 genes; E-value = 5.9e-008) and the only-binding p53-bound fragments (14,354 genes; E-value = 4.2e-004), respectively, using the de novo discovery algorithm MEME (Bailey et al., 2015). Two identical motifs were obtained, both of which are highly similar to the classic human p53 consensus motif ( Figure 4E), indicating high conservation of the p53 binding motif among species.

Cross-Species Analysis of p53 Transcriptional Regulation Networks
To determine the central modulators of the molecular network in LMH cells upon ectopic expression of p53, we performed TFs prediction using oPOSSUM v 3.0 (Ho Sui et al., 2005) based on the DNA sequences in § 100 bp intervals of the ChIP-seq peaks we identified. Twenty-nine TFs, including p53, were predicted to be the key regulators of the transcriptional shift induced by chp53 overexpression in LMH cells (Z score ≥10 and Validation of chromatin occupancies and mRNA expressions in a selected panel of newly identified p53 direct target genes. ChIP using the antibody Flag, followed by qPCR, was used to determine the relative p53 occupancies in the promoters of 13 selected genes in LMH cells upon ectopic expression of chp53 (upper panel). The mRNA levels of 13 selected genes were determined by RT-qPCR in LMH cells upon ectopic expression of p53 (bottom panel). Data are presented as the means § SDs (n = 3). (C and D) Gene ontology analysis was performed with DAVID (P < 0.05). (C) Biological process results from upregulation of putative chp53 direct target genes. (D) Biological process results from the downregulation of putative chp53 direct target genes. (E) The p53 consensus binding site was analyzed using the only-binding p53-bound fragments (14354-motif) or p53-bound fragments of chp53-target genes (675+79-motif), respectively, using the MEME algorithm. Abbreviation: LMH, leghorn male hepatoma.  Figures S1 and S2) reveals p53 as the most promising central modulator of the molecular events. The TF nodes were arranged in order of their functional connectivity to other proteins in the TF network (right). (A) TFs predicted by the p53-bound fragment sequences in chicken p53 direct target genes. (B) Human TFs were derived from a published human p53 dataset (Nikulenkov et al., 2012). (C) Venn diagram shows that seven of the human and chicken TFs are identical (left). The correlations analysis of the seven common TFs (middle and right). The Z score of chicken and human TFs in the TF network or TF together with p53 direct target gene network were calculated by normalizing the degree (middle and right). (D) Venn diagram showing the intersections of chicken p53 direct target genes and human p53 direct target genes (left). Gene ontology analysis of chicken p53 direct target genes (middle) and human p53 direct target genes (right) (P < 0.05, count > 5).
Fisher score ≥7) and were submitted for protein interaction analysis using the STRING database (Szklarczyk et al., 2019). The protein-protein interaction networks of these TFs or these TFs together with the chp53 target genes we identified were visualized by Cytoscape program (Shannon et al., 2003), with the layout of group attributes determined according to the connectivity degrees of the nodes ( Figure 5A and Supplementary Figure S1). A similar analysis of a published human p53 dataset was also conducted (Nikulenkov et al., 2012), and 45 TFs were predicted to be the key regulators of the transcription network of human p53 ( Figure 5B and Supplementary Figure S2). As expected, both chp53 and human p53 were predicted to be some of the most central transcriptional regulators with high number of connections with other key TFs and their target genes in their own transcriptional networks. Among these TFs, only 6 were common between chickens and humans, except p53 itself ( Figure 5C), which is consistent with the great difference in p53 target genes between the 2 species ( Figure 5D). Similar findings were observed when comparing human and mouse data (Tonelli et al., 2017;Supplementary Figure S3).
However, many well-known p53-mediated biological processes related to cell proliferation, apoptosis, and transcription are commonly enriched by both chicken and human p53 target genes ( Figure 5D). Consistently, among the 6 common TFs, SP1, EGR1, PAX5, and NFKB have been identified as the essential coregulators of p53 in these p53-mediated biological processes ( Figure 5C; Liu et al., 2011;Zwang et al., 2012;Li et al., 2014;Huang et al., 2016), suggesting the conservation of p53 key functions between species.
To gain an in-depth understanding of the above findings, protein alignments of chicken, human, and mouse p53 proteins were performed using ClustalW (MegAlign, DNAStar 7.1.0) with their full amino acid sequences. The sequence distances indicated the divergence and percent identity values of each sequence pair in the current alignment. Divergence is calculated by comparing sequence pairs in relation to the phylogeny reconstructed by MegAlign. Percent identity compares sequences directly, without accounting for phylogenetic relationships ( Figure 6A, lower left panel). These analyses revealed obvious differences in the sequences among chickens and the other species ( Figure 6A, upper panel). In addition, a phylogenetic tree was constructed using the neighbor-joining method to view the evolutionary relationships predicted from the above multiple sequence alignments. The length of each pair of branches represents the distance between sequence pairs, while the units at the bottom of the tree indicate the number of substitution events. A dotted line on a phenogram indicates a negative branch length, a common result of averaging ( Figure 6A, lower right panel). The results show that the evolutionary relationship between chicken and human or between chicken and mouse p53 proteins is distant. This finding was further proven by 3D in silico models generated using the amino acid sequences of p53 (e.g., chicken NP_990595.1, human NP_000537.3 and mouse NP_035770.2; Figure 6B). However, analysis of the conserved domains of p53 proteins from different species indicates the structural and functional conservation of p53 proteins from different species, suggesting that these p53 proteins from different species all perform regulation through similar signaling pathways ( Figure 6C). Furthermore, prediction of the secondary structure suggests that these p53 proteins in different species can maintain relatively conservative and stable functional topological folding structures, which is consistent with the higher homology and domain conservation at the level of their primary structure ( Figure 6D), indicating that p53 proteins from different species may have similar functions.
Thus, alignments of the amino acid sequences of the DNA binding domain of p53 proteins were also conducted using ClustalW, which revealed better identities and lower divergences among chicken p53 and p53 of other species compared with those of the full sequences of p53 proteins ( Figure 7A). The phylogenetic tree analysis also indicated closer evolutionary relationships between the chicken and human or between chicken and mouse p53 DNA-binding domains. Meanwhile, the conserved motifs of the DNA binding domains of p53 proteins were further analyzed using MEME (Bailey et al., 2015) and 2 conserved motifs in p53 proteins were identified among species ( Figure 7B), which further confirmed the conservation of the DNA binding domains of the p53 protein sequence. The 3D in silico models were generated using the amino acid sequences of these p53 DNA-binding domains ( Figure 7C). The predictions of the secondary structures suggests that these p53 DNAbinding domains from different species can maintain highly conserved and stable functional topological folding structures, which is consistent with the higher homology and domain conservation at the level of their primary structures. This proves that p53 proteins from different species have similar DNA binding domains, leading to their identical DNA binding motifs and may have conserved transcriptional regulating functions ( Figure 7D).

DISCUSSION
The effects of p53 on viral infections have been demonstrated in chickens, including our previous studies, which identified chp53 as a key determinant of ILTV infection (Li et al., 2018). However, due to the lack of fundamental information on the chp53 transcription regulatory network, the underlying molecular mechanisms of the antiviral function of chp53 remain unclear. Here, to gain deeper insights into chicken p53 functions, we integrated ChIP-seq and RNA-seq analyses to describe the genome-wide binding and direct regulation of genes by chp53 in chickens using LMH cells with ectopic expression of chp53 as the model and identified 754 direct chp53 target genes. Additional functional analysis and comparative biological analysis among species revealed conserved key biological functions and DNA binding motifs of p53, which may be due to the consensus amino acid sequence and structure of p53 DNA-binding domains among species.
P53 has a wide range of biological functions, many target genes, and very complex transcriptional regulation mechanisms. Therefore, it is necessary to clarify its target genes and cofactors that regulate the corresponding functions. To date, the antiviral functions of chp53 have been confirmed in many viruses, such as Marek's disease virus (Deng et al., 2010), avian leukosis virus subgroup J , and avian infectious laryngotracheitis virus (Li et al., 2018). However, only a few genes such as CDKN1A and MDM2, have been identified as direct target genes of chp53. Other target genes and cofactors are unknown, which hinders further mechanistic research on chicken p53. Parallel analysis of chicken p53 genomic occupancy and gene expression profiles allowed us to identify 752 novel p53 target genes, which greatly expanded our knowledge of chicken p53 transcriptional regulation and provided data support for further research in this field.
Our cross-species analysis found that many classic p53-mediated biological processes were enriched by the direct target genes of both chicken p53 and human p53 ( Figure 5D), indicating that the key functions of p53 may be conserved among species. However, the direct target genes and cofactors of regulation are quite different ( Figures 5C and D), suggesting differences among species in the p53 target genes that perform the same function. The transactivation domain or transactivating domain (TAD) is a transcription factor scaffold domain that contains binding sites for other proteins such as transcription coregulators. These binding sites are frequently referred to as activation functions (AFs) (Anette et al., 2003). The absence of the TAD domain in chicken p53 compared to human and mouse p53 ( Figure 6C) may partially explain the difference in genes targeted to p53 among species. In addition, among the 2.5% of the human sequences aligned with chickens (International Chicken Genome Sequencing Consortium, 2004), 44% are located in protein-coding regions, 25% are intronic, and 31% are intergenic, indicating that promoter/genic regions tend to be more conserved than intronic and intergenic regions. Our analysis of the locations of chp53 genomic binding sites revealed that 39.89% are located in promoter regions, 28.19% are intronic, and 23.51% are intergenic (Figure 2A), which are generally identical to the proportions of these regions in human-chicken conserved sequences. Thus, our data indicate that chp53 genomic binding sites tend to occur in conserved regions. However, the core cofactors and the targeted motif of chp53 are similar to those of humans ( Figure 4E and Figure 5C). This is because, despite the great differences in the full amino acid sequences of p53 among chickens, humans, and mice ( Figure 6), both the amino acid sequences and the simulated three-dimensional structures of the DNA binding domain of p53 are highly conserved (Figure 7). Therefore, from the perspective of comparative biology, the sequences and structures of functional domains, the sequences of p53 target sites on the genome, and at least part of its key biological functions are relatively conserved across species.
In the present study, genome-wide transcriptional regulation by chicken p53 was conducted using the chicken Figure 7. Comparative biological analysis of the predicted p53 DNA-binding domains. (A) The amino acid sequence analysis of p53 DNA-binding domains was performed using MegAlign 7.0.1, which included alignment reports, sequence distances, and phylogenetic trees. (B) The conserved motifs were determined using MEME with the p53 DNA-binding domain protein sequences. (C) In silico 3D prediction analyses for the three p53 DNA-binding domains were performed by the Phyre2 server. (D) Secondary structure and disorder prediction were also conducted by the Phyre2 server. hepatoma cell line, LMH, which is a permissive cell line for ILTV infection and has been used in vitro to study the viral-host interactions of ILTV in our previous studies. In this model, we discovered an important regulatory role of chp53 in ILTV replication by p53 activation and inhibition Li et al., 2018;Wang et al., 2019;Qiao et al., 2020;Wang et al., 2021;Xu et al., 2022). In the chicken fibroblast cell line, DF-1, which is nonpermissive for ILTV, remarkable activation of chp53 upon IBDV infection was reported, and overexpression of chp53 inhibited IBDV replication and upregulated the expressions of multiple chicken antiviral innate immunity genes (e.g., IPS-1, IRF3, PKR, OAS, and MX), whereas suppression of chp53 led to the opposite effect (Ouyang et al., 2017). Our present study provides a comprehensive view of the antiviral functions of chp53. While, as the cell line we used here, the ChIP-seq experiments for human p53 were mainly applied in cancer cells too (Koeppel et al., 2011;Nikulenkov et al., 2012;Hafner et al., 2017). Given the high tumor-suppressor potency of p53 and the induction of cell death of LMH cells by ectopically expressing p53 that we observed ( Figure 1D), our data may also contribute to the exploration of chp53 targets and networks related to responding to the cancerogenesis in chickens.
In conclusion, our study identified 754 chicken p53 target genes and conducted a preliminary comparative biological analysis with humans and mice. These results are expected to help harness p53 activity by targeting p53 modulators and achieve a desired p53 transcriptional response in chickens in the future.