Introduction
Breast cancer is the most frequent form of carcinoma in females. Approximately 230,480 new cases of invasive breast cancer and 39,520 breast cancer mortalities are expected to occur among females in the USA in 2011 (1). The etiology of breast cancer is multifactorial, the period of development can span decades and the clinical course is highly variable. Even though novel strategies in the treatment of breast cancer are emerging (2–4), overall results remain unfavorable. In the past few years, numerous studies have been performed on the identification of cytokines as prognostic factors. Innate and acquired arms of the immune system were discovered to play crucial roles in the anti-tumor response (5,6). Evidence in recent years has proven that cell-specific enhancers also play an important role in the development of cancer (7,8).
p300, a histone acetyltransferase and transcription coactivator, is a near-ubiquitously expressed component of enhancer-associated protein assemblies and is critically required for embryonic development (9). Recently, p300 has been studied for its roles in the development of human disease. Investigators have demonstrated that a high expression of p300 in human breast cancer correlates with tumor recurrence and predicts an adverse prognosis (10). Furthermore, simultaneous modulation of p300/nuclear factor-κB signaling can inhibit proliferation and induce apoptosis in breast cancer cells (11).
Chromatin immunoprecipitation followed by sequencing (ChIP-seq) is a technique for genome-wide profiling of DNA-binding proteins, histone modifications or nucleosomes (12). This technique accurately predicts tissue-specific activity of enhancers (8). In this study, we present the results of chromatin immunoprecipitation coupled with massively parallel sequencing with the enhancer-associated protein p300, and several thousand binding sites of p300 were mapped in the estradiol (E2) stimulation and control groups in the MCF7 human breast cancer cell line. The present results indicate that p300 binding is a highly accurate process for identifying enhancers and their associated activities, and suggests that these data sets are useful for investigating the role of specific enhancers in human biology and disease on a genome-wide scale.
Materials and methods
ChIP-seq data of p300 in human breast cancer cell line MCF-7
ChIP-seq is a technique for genome-wide profiling of DNA-binding proteins, histone modifications or nucleosomes (12). The ChIP-seq data of GSE39623 (13), including E2 stimulation data of GSM986085 and control data of GSM986087 of p300 in MCF-7 cell, was obtained from National Center of Biotechnology Information Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo), which is based on the Illumine GAIIx to sequence the binding fragment of extracted p300.
Alignment ChIP-seq reads with the reference genome
Comparison of the reads was performed using Bowtie (14) (version 0.12.9) on the hg19 Refseq RNA sequences downloaded from the University of California, Santa Cruz Genome Browser (http://genome.ucsc.edu). Bowtie was run with i) the uniquely aligned reads and ii) the mismatch numbers were <2. Only reads matching the above requirements were used for subsequent analysis. The detailed statistical information is shown in Table I.
Identification of p300 binding sites among the whole genome
SAMtools is a library and software package for parsing and manipulating alignments in the SAM/BAM format (15). Firstly, SAMtools was utilized to remove PCR duplicates. For multiple reads on a single chromosome, a flattened operation was adopted to ensure that the number of reads in a single comparison site was <1. Secondly, the two-cross-correlation software of SPP and MaSC was employed to confirm the average length of ChIP-seq fragments, half of which were selected as translational distance of reads 5′ and 3′ direction. Thirdly, MACS1.4.0 software was used to identify enriched ChIP regions of p300 and peak calling was employed to infer the actual binding loci from the positional distribution of sequenced DNA fragments. The default setting of the significant level was 0.00001. Only when all the reads of the enriched region were P>0.00001 were they regarded as p300 binding sites.
Potential target gene annotations of p300 binding sites and gene ontology enrichment analysis
ChIPpeakAnno, a bioconductor package within the statistical programming environment R to facilitate batch annotation of the enriched peaks identified from ChIP-seq and ChIP-chip data (16), was used to detect target genes of the p300 binding sites. Gene Ontology (GO) analysis has become a commonly used approach for functional studies of large-scale genomic or transcriptional data (17). Based on the GO database, functional enrichment analysis of these target genes in the biological process was conducted (P<0.00001).
Differential analysis of dynamic variation of p300 binding sites
To eliminate impact bias resulting from the length of peak region and difference numbers of the available reads, the random particle-mesh method was employed to standardize the read signal of the peak region in any sample. NOIseq, a non-parametric approach and an RNA-seq differential expression method robust for sequencing depth biases (18), was used to identify the differential expression of the peak of reads enrichment in two samples. The significant Q-value of NOISeq was <0.8.
Heat map of the ChIP-seq signal of p300 binding sites
The cluster heat map is an original display that simultaneously reveals row and column hierarchical cluster structure in a data matrix (19). SeqMINER-1.3.3e, an integrated user-friendly platform that allows the central questions in the ChIP-seq analysis workflow to be addressed, was used to generate the heat map of the ChIP-seq signal of the p300 binding sites and the mean value of K for all the grid points was utilized to conduct clustering analysis.
Motif finding of transcription factors (TFs) of p300 binding region
Online Seqpos was applied to conduct motif finding analysis. The search scope of a motif is 300 base pairs upstream and downstream of the hub of the peak. The P-value was <0.00001.
Results
Screening out the p300 binding sites in MCF7 cells in the E2-treated (p300_E2) and control groups (p300_Veh) and identifying the distribution characteristics
MACS peak calling was employed to identify the p300 binding sites, including 32,249 in p300_E2 and 43,156 in p300_E2. Subsequently, the distribution characteristics of the p300 binding sites on the chromosome were analyzed (Fig. 1). From the result, the p300 binding sites were found to focus on two flanks of 5′untranslated region (UTR), 3′UTR, transcriptional start sites (TSS) and transcriptional termination sites, while the distribution density in the open reading frame internal sites was low. In addition, in p300_E2, the distribution density in the two flanks of TSS slowly increased, whilst the other regions were stable.
Identification of p300-regulated target genes
Based on the gene detection around the p300 binding sites, 12,883 target genes in p300_E2 and 15,511 target genes in p300_Veh were found. There are 2–3 binding sites around the average target gene.
GO enrichment analysis of p300 target genes
In combination with the GO database, GO functional enrichment analysis was conducted for p300-regulated target genes (Figs. 2 and 3). In p300_Veh, p300-regulated target genes mainly participated in the neural cell differentiation-related biology process, including neuron differentiation, neuron projection morphogenesis and forebrain development. Furthermore, the functions of the transmembrane receptor protein tyrosine kinase signaling pathway, axonogenesis and circulatory system development were significantly enriched in MCF7 cells (P<0.00001). Whereas in p300_E2, partial functions of p300 regulatory genes were changed, including cell-cell junction organization and regulation of cell adhesion, which were significantly increased. Furthermore, the GO-based annotation similarity of genes in neurogenesis and circulatory system development remained unchanged and is referred to as the p300-regulated target gene.
Differential signal analysis of p300 binding sites in E2-treated and control groups in MCF7 cells
Based on the NOISeq results, 24,899 differential binding sites were identified in total. A ChIP-seq signal of 39% (9,812) binding sites was significantly increased in p300_E2, while the remaining 61% (158077) binding sites maintained a high level in p300_Veh. In addition, the other 28,090 binding sites had no statistically significant change in the two groups (Fig. 4).
Identifying target genes of the differential p300-bound sites and GO enrichment analysis
In total, 5,642 and 6,952 sample-specific target genes of differential binding sites in p300_E2 and p300_Veh were identified, respectively. Subsequently, GO enrichment analysis was performed and found that all these target genes were associated with neural differentiation. However, cell locomotion emerged only in p300_E2 (Figs. 5 and 6).
Exploring p300 concordant potential TFs
To analyze the effect of E2 in MCF-7 cells, motif finding analysis was utilized to identify the sample-specific p300 binding region in the two groups, involving 225 and 176 TFs of significantly enriched motif clustering around the binding sites (P<0.00001). Removing 98 identical TFs, 127 TFs were searched and colocalized with p300 in p300_E2, while there were 78 TFs in p300_Veh. GO enrichment analysis of the newly-presented TFs in p300_E2 showed that these concordant TFs mainly participated in the biology processes of cell proliferation regulation, negative regulation of cell differentiation and cell migration (Table I). Furthermore, ESR1 and BRCA1 are involved in the steroid hormone receptor-signaling pathway and response to the estrogen stimulus.
Discussion
Transcriptional coactivator p300 has been shown to play a critical role in the transcription process (20). Several studies revealed that p300 is a positive regulator of cancer progression and is associated with tumorigenesis of various types of human cancer (21,22). However, the expression dynamics of p300 in breast cancer and its effect on the prognosis of breast cancer are poorly understood. In the present study, the genome-wide distribution of p300 binding sites was identified and the distribution characteristics in p300_E2 and p300_Veh were analyzed using ChIP-seq directly from the human breast cancer cell line, MCF7 (23). Despite two flanks of TSS, p300 binding sites of other regions were not affected by E2 stimulation.
To develop an objective analysis of the main processes altered in p300_E2 and p300_Veh, p300-regulated target genes were identified and GO enrichment analysis was conducted. The results showed that in p300_Veh, the largest and most noteworthy is the neural cell differentiation-related gene ontology terms. They are primarily associated with neurogenesis system development, axonogenesis, neuron projection development and neuron differentiation. Differentiation programs are abnormal in cancer cells, which allows them to express differentiation markers in addition to their tissue of origin (24). This indicates that terminal differentiation to the anticipated cellular type is compromised in the cancerous state, but aberrant multi-lineage trans-differentiation or lineage infidelity may be wide spread in cancer phenomenon. Metastasis to distant tissues is the main driver of breast cancer-related mortality. Investigators have demonstrated that activation of the sympathetic nervous system is a novel neural regulator of breast cancer metastasis (25). The role of axonogenesis and neurogenesis in cancer remains relatively unexplored (26).
In p300_E2, the biological processes of cell locomotion and cell adhesion regulation have significantly enriched through studies. E2 enhances cocaine-stimulated locomotion in mice predominantly through estrogen receptor-α (ERα) (27) and in a direct manner that cancer cell locomotion is in fact important in invasion of biological tissues (28). A recent study indicates that random cell locomotion plays an important role in C5α-triggered cancer cell invasion, which may provide a useful therapeutic option for cancer treatment (29). Hein et al (30) indicate that overexpression of the activated leukocyte cell adhesion molecule in breast cancer may be relevant for the outcome in ductal carcinomas.
To study the effect of E2 in MCF7 cells, motif-finding analysis to the p300 specific binding region was conducted. The results suggested that E2 stimulation affected binding affinity of DNA and TFs in MCF7 cells in a large scale. These TFs can regulate transcription of relevant target genes in association with p300. The GO functional analysis of p300 concordant TFs discovered that these TFs, such as ESRRA, ESR1 and BRCA1 participated in the biological process of cell proliferation regulation, negative regulation of cell differentiation, cell migration, steroid hormone receptor signaling pathway and response to estrogen stimulus.
ESRRA has been most studied in the context of breast cancer. ESRRA is a negative prognostic marker in ER(−) tumors (31,32), and it induces VEGF mRNA expression and contributes to the malignant phenotype of a breast cancer cell line (33,34). Aecent study has proved that ESRRA-C11orf20 is a recurrent gene fusion in serous ovarian carcinoma (35).
Approximately 70% of breast cancers are positive for the nuclear receptor ERα and require estrogen for cell proliferation. Gene expression studies and TF-binding maps have revealed that a number of estrogen-regulated genes are subsequently used by the ER complex for efficient gene transcription (36,37). Estrogen-ER can transcribe genes that ultimately contribute to ER function.
BRCA1 has been indicated in numerous important cellular functions, such as DNA damage repair, cell cycle checkpoint control, apoptosis and transcriptional regulation (38). BRCA1 was initially identified as one of the genes conferring genetic predisposition to breast and ovarian cancer. One of the noteworthy aspects of BRCA1-associated cancers is the observed specificity for estrogen-responsive tissues, including the breast and ovary. Recent advances in the understanding of BRCA1-linked breast cancers have revealed a complex association between BRCA1 and ERα signaling. Estrogen stimulation increases BRCA1 expression at the mRNA and protein level and conversely BRCA1 functions to induce ERα mRNA expression and act as a negative regulator of ERα signaling (39). Estrogen stimulation increases BRCA1 mRNA expression levels through a mechanism potentially involving p300 and aromatic hydrocarbon receptor. In addition, BRCA1 competes with p300 and cyclin D for binding to ERα and negatively regulates ERα-mediated transactivation of its downstream target genes. A more recent study showed an alternative model of regulation, in which ERα and its cofactor p300 are recruited to an activator protein site on the BRCA1 promoter following E2 stimulation (40).
Overall, the present study confirms an integral role of p300 and concordant TFs. The organization of mammalian enhancers into correlated networks is possibly mediated by the joint action of TFs through shared motifs. As the predominant hormone in breast cancer biology is estrogen, it is important that the interactions between ER and p300 be determined under these specific conditions. Therefore, p300 concordant TFs should be considered to further confirm other carcinogenic factors. This may be useful in the study of the biology of human disease and provide theoretical basis for cancer therapy.