Global Egr1-miRNAs binding analysis in PMA-induced K562 cells using ChIP-Seq.

Although much is known about microRNAs' regulation in gene expression and their contributions in cell fate, to date, globally lineage-(cell-) specific identification of the binding events between a transcription factor and its targeting microRNA genes is still waiting for elucidation. In this paper, we performed a ChIP-Seq experiment to find the targeting microRNA genes of a transcription factor, Egr1, in human erythroleukemia cell line K562. We found Egr1 binding sites near the promoters of 124 distinct microRNA genes, accounting for about 42% of the miRNAs which have high-confidence predicted promoters (294). We also found EGR1 bind to another 63 pre-miRNAs. We chose 12 of the 187 microRNAs with Egr1 binding sites to perform ChIP-PCR assays and the positive binding signal from ChIP-PCR confirmed the ChIP-Seq results. Our experiments provide the first global binding profile between Egr1 and its targeting microRNA genes in PMA-treated K562 cells, which may facilitate the understanding of pathways controlling microRNA biology in this specific cell line.

1. Introduction

MicroRNAs (miRNAs) are a family of ~22-nucleotide small noncoding RNAs in eukaryotes and mainly involve in regulation at posttranscriptional level by translational repression or degradation their target mRNAs [1, 2]. More than 700 human miRNAs have been identified up till now [3] and they are estimated to control about one-third of human known genes [4]. miRNAs have been reported to regulate hematopoietic lineage differentiation, angiogenesis, cell adhesion and so on [1, 5]. K562 is a cell line deriving from chronic myeloid leukemia, which is a common progenitor of megakaryocytic and erythroid lineages of the hematopoietic stem cell differentiation. It can be induced to differentiate into erythrocytes or megakaryocytes (MK) by hemin and Phorbol-12-myristate 13-acetate (PMA), respectively, [6-9]. Recently, miRNAs have been found to play a key role in K562 differentiation. miR-27a, miR-34a, miR-223 were up-regulated when K562 was induced to MK status while miR-27a, miR-223, miR-103, miR-130a, miR-210, and miR-18b were downregulated when K562 was induced to erythroid differentiation [10-13]. The changes of miRNAs expression level gave clues for their functions in hematopoietic lineage differentiation but a more detailed regulation pathway is anticipated to be understood: by which targets miRNAs realize their functions in hematopoietic differentiation and miRNAs are subjected to which factors controlling their transcription? Garzon et al. confirmed that miR-130 targets the transcription factor MAFB and participates in MK differentiation by up-regulating its expression level in CD34+ hematopoietic progenitors [14]. Navarro et al. found that independently of p53, miR-34a directly regulates expression of MYB facilitating megakaryocytic differentiation of K562 cells and of CDK4 and CDK6, to inhibit the G1/S transition [11]. Lu et al. found that miR-150 regulates megakaryocyte-erythrocyte progenitors (MEPs) differentiation and is preferentially expressed in megakaryocytic lineage. Besides, Lu et al. identified that transcription factor MYB is also a critical target of miR-150 in this regulation [15]. A feedback loop between Runx1 and miR-27a was found by Ben-Ami et al. that miR-27a plays a regulatory role in megakaryocytic differentiation by attenuating Runx1 expression, and during megakaryopoiesis, Runx1 exerts positive regulation of miR-27a expression [10]. While the researches, which deepen our understanding of the underlying mechanism of cellular differentiation involving miRNAs, are based on individual specific miRNA, currently the interactions between a TF and its target miRNAs in this cellular process on a large scale have been sparsely investigated.

Egr1 is an immediate-early response protein which is rapidly and transiently induced by various stimuli, such as different growth factors, cytokines, mechanical injury, shear stress [16]. As a transcription factor, many of its known transcription target genes are protein coding genes. Among the known transcriptional targets of Egr1, a part of genes were implicated in the pathogenesis of vascular disease, including PDGF-A, PDGF-B, FGF-2, SOD1, p53, CD44 (see commentary [17]). Besides, Egr1 may regulate cell interaction with the extracellular matrix by coordinated induction of TGF-β1, FN, and PAI-1 in human glioblastoma cells [18]. Moreover, the Egr1 gene is functionally implicated in cell proliferation and in the regulation of apoptosis and is considered as a potential target for prostate cancer therapy [19]. Few studies in the past, on the other hand, has been focused on the Egr1's regulation role on noncoding target genes except a finding of Egr1's role in hsa-miR-106a transcription. Through hsa-miR-106a, Egr1 indirectly regulates the IL-10 expression [20]. Previous studies showed that PMA-induced activation of Egr1 expression is involved in megakaryocytic differentiation of K562 cells but the regulation pathway has not been investigated. As an early response factor, Egr1 may regulate both downstream protein-coding and noncoding genes. ChIP-Seq, a technique based on next-generation sequencing, is able to capture the genome-wide snapshot of protein-DNA interactions in vivo [21]. The interactions between a TF and its target genes can be identified by ChIP-Seq and miRNA genes are also the targets of TFs. Here we identified the transcription factor Egr1 binding sites near miRNA promoters and premiRNAs in PMA-induced K562 cells using ChIP-Seq. Combining with the Navarro's miRNA expression profile [11] before and after PMA or hemin treatment in K562 cells, we analyzed the potential role of Egr1 in differentiation fate of K562 cells through targeting miRNA genes.

2. Materials and Methods

2.1. Cell Culture

K562 cell was cultured in RMPI 1640 medium containing 10% fetal bovine serum (GIBCO) at 37°C with 5% CO2 incubator. K562 cells were treated by PMA (10 ng/ml) for 2 hours.

2.2. Chromatin Immunoprecipitation (ChIP)

ChIP experiment was carried out as described previously [22]. Briefly, PMA treated K562 cells were crosslinked with 1% formaldehyde and incubated at room temperature (RT) for 12 minutes, prior to the addition of glycine to 0.125 M followed by incubation at RT for 5 minutes. Chromatin were sonicated and was capture with the beads coated with EGR1 antibody (588) (from Santa Cruz Biotechnology) or no IgG. The beads were then precipitated and washed. Bead complexes were eluted with elution buffer. Samples were incubated overnight at 65°C for reversing crosslink. Samples were then incubated with Proteinase K and RNase A for 2 hours at 50°C. DNA was purified by two rounds of phenol-chloroform extraction and ethanol precipitation and resuspended in 50 μl dH2O.

2.3. Sequencing

Briefly, EGR1 ChIP DNA was sheared into 90–150 bp and then was repaired by DNA end-repair kit (Epicentre). SOLiD adaptors were then added and DNA was subjected to 19 cycles of PCR (2 × enhancer (invitrogen) was added to increase the amplification for GC-rich sequences). PCR product of 180–220 bp size was incised from gel and then was subjected to emulsion PCR. Template positive beads were collected and extended in the presence of terminal transferase and Bead Linker (SOLiD reagents, Applied Biosystems). Extended beads were then deposited and covalently attached onto slides at >20,000 beads per panel (SOLiD reagent, Applied Biosystems). Template bead slides were then loaded onto a SOLiD Analyzer, and cycled ligation sequencing using the SOLiD Sequencing System was performed.

2.4. Mapping of the SOLiD Sequencing Data

Sequence reads from ChIP fragments were aligned to human reference genome using standard SOLiD mapping pipeline tool (Corona Lite, [23]), allowing up to 3 mismatches out of 35 bp. All mapping was performed in the SOLiD color space corresponding to dinucleotide encoding of the sequenced DNA.

2.5. Peak Detection and Motif Analysis

The uniquely mapped reads were retained for subsequent analysis. An integrated software CisGenome was used to detect peaks which are enriched from background reads [24]. The algorithm applied the sliding window method to count the reads and the window size was set 500 bp, a tradeoff size which had best coverage for the known Egr1 targets and was consistent with length of DNA fragments from ChIP experiments. The position with highest density in a window should contain at least nine reads otherwise the regions will be discarded. We used the strand filtering function to make sure each strand meet the minimum three reads requirement.

All the motif analysis, including known motif mapping and new motif discovery, in this study were completed by CisGenome. The known motif sequence (degenerate pattern: 5′-G[T]C[A]GG[T]GGGCG[A,T]-3′) for EGR1 was curated from literature [17, 25].

2.6. miRNA Promoters and PremiRNAs Annotation of Egr1 Binding Sites

We directly applied 294 predicted high-confidence miRNA promoters [26] with minor revisions to annotate our Egr1 ChIP-Seq dataset. Firstly, the coordinates of the predicted promoter regions were converted to the current version of the human genome (hg18) from hg17 using the liftover utility provided by UCSC Genome Browser [27]. Then, the miRNA gene names were revised based on miRBase v13 [3] because they are supposed to be criterion to link the miRNA expression profile in the subsequent analysis.

For sequence-specific factors such as Egr1, previous studies used a more relaxed region of 8 kb surrounding the promoters [26, 28]. As the length of those predicted promoters varies, here we define that the enriched regions which are located within 5 kb regions of predicted promoters and premiRNAs are candidate Egr1 binding sites.

2.7. Gene Ontology Analysis of the Targets of miRNA

The experimentally supported miRNA targets were downloaded from TarBase v5.0 [29]. We chose the Batch-Genes tools provided by GOEAST to do Gene Ontology(GO) enrichment analysis [30]. Given a set of genes, Batch-genes can find statistically significantly enriched GO terms among them by using a species-specific random background. Hypergeometric test was chose and the P-value was set to  .01 (Yekutieli FDR adjusted) to call an enriched GO term.

2.8. ChIP-PCR Analysis

We randomly selected seven miRNAs together with another five up-regulated miRNAs after PMA treatment (hsa-miR-135b, hsa-miR-141, hsa-miR-152, hsa-miR-199b, hsa-let-7g) [11] among Egr1 target miRNAs to perform ChIP-PCR analysis and their corresponding ChIP DNA fragments resulted from two different cell conditions: before PMA treatment and after PMA treatment of K562 cells. The designed PCR primers are listed in supplementary Table 1at supplementary material available online at doi:10.1155/2010/867517.

3. Results

We performed a ChIP-Seq experiment on next-generation sequencer SOLiD system for the transcription factor Egr1 in PMA-induced K562 cells. Relying on the hypothesis that chromatin signatures (H3K4me3) can be used to locate transcription start sites of most genes in genome, Marson et al. [26] predicted 294 distinct high-confidence miRNA promoters utilizing the H3K4me3 dataset from ChIP-chip and ChIP-Seq experiment. We adopted these promoters to annotate our Egr1 binding sites and anticipated to find Egr1-miRNA gene interactions. We also surveyed the premiRNA with Egr1 binding signal. ChIP-PCR was conducted to assess the reliability of part of those binding events. Besides, a recent published miRNAs expression profile in the similar cell condition facilitated to analyze the expression changes of those miRNAs which are bound by Egr1 in PMA-induced MK differentiation of K562 cells [11].

3.1. Bioinformatics Evaluation of the Quality of ChIP-Seq Experiment

By setting the cutoff 9 during the peak detection, we got 14,636 enriched regions. Although the negative binominal model adopted by CisGenome is able to estimate a stringent false discovery date (FDR) and the FDR was less than 10% when the peak height is equal to and higher than 17 at the window size of 500 bp, we chose to combine the FDR estimation with further experimentally validations to determine the threshold and the cutoff 9 is supposed to compromise the coverage and specificity of the results.

Here we evaluate our ChIP-Seq experiments by motif analysis. The enriched regions were scanned by the canonical EGR1 motif (degenerate pattern: 5′-G[T]C[A]GG[T]GGGCG[A,T]-3′ [17, 25]) within different peak height rank. Among the peaks with height ≥17, 60.8% possessed the canonical EGR1 motifs. We observed that the percentages of the canonical EGR1 binding motifs decreased with peak rank. 43.2% of the full list of peaks can be found canonical EGR1 binding motifs. Besides, the sequences extracted from 2042 peaks with high quality (peak height more than 20) were used to do de novo motif finding analysis. We found two novel motifs, one (5′-GCGG[T]GGGC[T]GG-3′) resembling the known consensus of EGR1 binding site [31], and one is the other (5′-GGGGC[T]GGGG-3′) resembling the consensus of EGR1 binding sites which has been identified by Kubosaki using ChIP-chip [32]. These information are presumed to give an overview of the technical quality of this ChIP-Seq experiments in a bioinformatics way.

We next examined the EGR1 binding sites distribution relative to protein-gene structure. We found the majority of the EGR1 binding loci (40.7% of the 14,636 sites) are within the range from 10 kb upstream of the transcription start site (TSS) to 10 kb downstream of the transcription end site. About 12.2% of the peaks are located within 1 kb upstream of the TSS of the known genes, suggesting the transcriptional relevance of EGR1. A more detailed discussion of protein-coding gene targets of EGR1 using this ChIP-Seq experiment can be found in [33]. Due to the unbiased identification of sequencing method, it is unsurprising to find that many EGR1 binding sites are located in exons, introns, and intergenic regions.

3.2. Identification of the Occupancy of miRNA Promoters and PremiRNAs by Egr1 in K562 Cells

That 37.3% and 50.2% of the 14,636 enriched regions are located in intron regions and intergenic regions, respectively, elicited our curiosity to check if EGR1 is able to bind to noncoding genes like miRNAs, as 5,091 (91%, 5,604 miRNAs' context is known in miRBase v13 [3]) miRNAs' context is intron. With the emergence of the predicted miRNA promoters, we firstly defined those regions, which are located within 5 kb regions upstream and downstream of the predicted promoters, are candidate binding sites by Egr1. The promoters of 124 different miRNAs were found with ChIP-Seq signals and they accounts for 42% of those miRNAs with known predicted promoters. hsa-mir-10a is one of the miRNA genes and its Egr1 binding signal is visualized in Figure 1. hsa-mir-10a has two Egr1 ChIP-Seq enriched binding signals (chr17:44008884-44011068, reads count:51; chr17:44013523-44015182, reads count:9) near its predicted promoters (chr17:44014310-44014709).

Figure 1

Visualization of Egr1 ChIP-Seq peaks enriched around the promoter of miRNA hsa-mir-10a. Two peaks (chr17:44008884-44011068, reads count:51; chr17:44013523-44015182, reads count:9) around the promoter of hsa-mir-10a can be detected. Raw signal represents the aligned tags. Total IP, 5′ IP and 3′ IP are the signal profile built by sliding windows. 5′ IP and 3′ IP means immunoprecipitationed DNA signal which are clustered by the forward and reverse aligned tags, respectively. The figure is drawed by CisGenome browser [24].

Although most miRNA genes (469) are still lacking high-confidence predicted promoters, to get a comprehensive view of the target miRNA genes, we also surveyed the EGR1 binding signal around all precursor miRNAs (miRBase v13 [3]). 125 precursor miRNAs were found EGR1 binding signal by the same distance standard for the miRNAs with promoters (also listed in supplementary Table 1). Among the 125 target miRNAs identified by precursor coordinates and 124 target miRNAs identified by promoter coordinates, 62 miRNA genes were found present in both lists (see Figure 2), thus, 187 miRNAs were found EGR1 binding signal in all. The reason for those 62 miRNA genes which can be identified by both coordinates is that the relative distance to mature miRNAs is one of elements to score the predicted miRNA promoters. Among the 12 miRNA genes whose ChIP-Seq EGR1 binding were confirmed by ChIP-PCR (see Figure 3), 9 (mir-10a, mir-130b, mir-95, mir-23b, mir-17, mir-152, mir-141, mir-135b, and mir-199b) were among the 62 miRNAs whose EGR1 binding signal can be identified by both the promoters and precursors of miRNAs and the left 3 ones were only found near the predicted promoters of miRNAs.

Figure 2

The overlap of miRNA list with EGR1 binding. 125 precursor miRNAs were found EGR1 binding signal in addition that the promoters of 124 different miRNAs were found with ChIP-Seq signals. 62 miRNA genes were found present in both lists.

Figure 3

ChIP-PCR results of 12 miRNAs before (−) and after (+) PMA treatment. All the 12 miRNAs showed positive Egr1-binding signals after PMA treatment, which verified the ChIP-Seq results (+Ab).

A region containing canonical consensus of a transcription factor has great potential to be the targets of this factor. To determine if the extra 63 miRNAs should be accepted as the targets of Egr1, the peaks associated with miRNAs were scanned by Egr1 canonical consensus sequence (degenerate pattern: 5′-G[T]C[A]GG[T]GGGCG[A,T]-3′ [17, 25]). 61.2% of the peaks associated with miRNA genes (220 peaks in all) were found at least one Egr1 canonical motif (see Table 1), which is significantly more than the ratio of the full set peaks (Fisher's exact test, P-value = 9.8 × 10−8). However, the difference of the motif frequency between the peaks associated with the promoters of miRNAs and the peaks linking to premiRNAs is not significant (Fisher's exact test, P-value = .1), so we decided to view the extra 63 miRNAs as the putative targets of Egr1. All miRNAs with Egr1 binding sites enriched from ChIP-Seq signals are listed in supplementary Table 1. It is noteworthy that a miRNA gene may have two or more predicted promoters and a promoter may have more than one Egr1 binding sites.

Table 1

Genomic features (especially intergenic region and intron region) and motif properties of the ChIP-Seq enriched regions associated with miRNA genes.

	No. of peaks located in intergenic region (ratio)	No. of peaks located in intron region (ratio)	total peak no.	No. of peaks with canonical motif (ratio)
All peaks enriched	7340 (50.2%)	5460 (37.3%)	14,636	6318 (43.2%)
Peaks associated with miRNA genes (187)	95 (43.2%)	90 (40.1%)	220	135 (61.2%)
Peaks associated with the promoter of miRNA genes (124)	64 (46.0%)	67 (48.2%)	139	94 (67.6%)
Peaks associated with pre miRNAs (125)	60 (46.9%)	46 (35.9%)	128	74 (57.8%)

Among the 14,636 peaks, most are located in intergenic and intron region as mentioned above, so we asked the genomic features of those peaks relating to miRNA. 220 unique peaks are associated with the 187 miRNA genes. Among them, 40.1% of peaks are located in intron region (see Table 1), which is a little more than the ratio of the same genomic peaks (37.3%) among the whole set of peaks (14,636). Higher ratio of peaks (48.2%) located in intron region was observed for peaks near the predicted promoter of miRNA genes. However, lower ratio of peaks (35.9%) located in intron region was observed for peaks near the premiRNAs. This is because most premiRNA genes are already located in intron region and the average length of human intron regions is 6,206 bp. In this case, peaks around the −5 kb ~ +5 kb of the precursor of miRNA genes can easily be outside of the intron region. In general, no significant increase of the amount of peaks as putative Egr1-miRNA binding are located in introns (Fisher's exact test, P-value = .29).

3.3. Gene Ontology Enrichment Analysis of Egr1-Bound miRNAs' Targets

In order to elucidate the functions of the 187 Egr1-bound miRNAs, we examined gene ontologies (GO) using the web-based analysis toolkit GOEAST [30] for the targets of those miRNAs as miRNAs regulate cellular activities via pairing to mRNAs. Since our experimental results are nonmicroarray-based high throughput, we chose the Batch-Genes tools provided by GOEAST to do GO enrichment analysis. Given a set of genes, Batch-genes can find statistically significantly enriched GO terms among them by using a species-specific random background.

Although most mammalian mRNAs are conserved targets of miRNAs, a significant proportion of false positives are thought to be included in the predicted miRNA targets. Here, we chose to analyze those experimentally supported miRNA targets. Among the 187 Egr1-bound miRNAs, 38 miRNAs were found experimentally validated targets (247) according to TarBase v5.0 [29]. In the GOEAST analysis, the 247 genes were compared to 18,596 genes as background and the P-value of a GO term less than.01 (Yekutieli FDR adjusted) was thought to be enriched using hypergeometric test. We present the first 25 enriched biological process terms in Table 2 and the full enriched molecular function terms in Table 3. Interestingly, the statistically significantly overrepresented GO biological process terms are highly enriched for hemopoiesis, immune system development, cell differentiation and so on, which are the cases supposed to be occurring in this cell condition. Moreover, among enriched GO molecular function terms, the targets of miRNA genes with EGR1 binding signal include binding, protein binding, similar to the previous GO enrichment finding of gene targets of EGR1 using ChIP-chip [32]. Besides, molecular function GO terms are significant for transcription factor transcription, transcription factor activity, transcription activator binding, and transcription activator activity. This results is consistent with the discovery that transcription factors also prevail among human miRNA targets [34]. It also supports the notion that EGR1 acts as an early response protein in cell events and may trigger a series of transcription factor activity. Three graphical outputs showing the hierarchical relationships of enriched GO terms (P ≤ .01, Yekutieli FDR adjusted) in each GO category (can be seen in supplementary Figures 1–3).

Table 2

Enrichment of Gene Ontology biological process terms in targets of miRNAs with EGR1 ChIP-Seq hits.

GO term IDs	Enriched GO terms	P-value
GO:0009987	Cellular process	5.06E − 08
GO:0048513	Organ development	9.35E − 06
GO:0006357	Regulation of transcription from RNA polymerase II promoter	1.09E − 05
GO:0060056	Mammary gland involution	1.29E − 05
GO:0044237	Cellular metabolic process	3.30E − 05
GO:0060443	Mammary gland morphogenesis	3.65E − 05
GO:0022612	Gland morphogenesis	6.49E − 05
GO:0048731	System development	1.26E − 04
GO:0030097	Hemopoiesis	1.32E − 04
GO:0002520	Immune system development	1.43E − 04
GO:0045944	Positive regulation of transcription from RNA polymerase II promoter	1.45E − 04
GO:0030154	Cell differentiation	1.53E − 04
GO:0030155	Regulation of cell adhesion	1.80E − 04
GO:0048856	Anatomical structure development	1.85E − 04
GO:0042127	Regulation of cell proliferation	2.03E − 04
GO:0009888	Tissue development	3.87E − 04
GO:0048869	Cellular developmental process	3.87E − 04
GO:0048534	Hemopoietic or lymphoid organ development	3.87E − 04
GO:0010468	Regulation of gene expression	4.29E − 04
GO:0008285	Negative regulation of cell proliferation	5.44E − 04
GO:0006928	Cellular component movement	6.05E − 04
GO:0007162	Negative regulation of cell adhesion	6.10E − 04
GO:0042060	Wound healing	6.14E − 04
GO:0014068	Positive regulation of phosphoinositide 3-kinase cascade	8.61E − 04
GO:0032502	Developmental process	9.71E − 04

Table 3

Enrichment of Gene Ontology molecular function terms in targets of miRNAs with EGR1 ChIP-Seq hits.

GO term IDs	Enriched GO terms	P-value
GO:0005515	Protein binding	2.93E − 05
GO:0030528	Transcription regulator activity	1.89E − 04
GO:0051879	Hsp90 protein binding	2.64E − 04
GO:0003700	Transcription factor activity	8.36E − 04
GO:0004722	Protein serine/threonine phosphatase activity	9.73E − 04
GO:0004716	Receptor signaling protein tyrosine kinase activity	1.53E − 03
GO:0030618	Transforming growth factor beta receptor, pathway-specific cytoplasmic mediator activity	1.62E − 03
GO:0005488	Binding	2.74E − 03
GO:0003702	RNA polymerase II transcription factor activity	3.38E − 03
GO:0033613	Transcription activator binding	6.08E − 03
GO:0016563	Transcription activator activity	7.22E − 03
GO:0005072	Transforming growth factor beta receptor, cytoplasmic mediator activity	8.44E − 03
GO:0034713	Type I transforming growth factor beta receptor binding	8.44E − 03

3.4. An Evaluation of the Predicted miRNA Promoters Using ChIP-Seq

All the predicted miRNA promoters were also scanned by Egr1 canonical consensus sequence (degenerate pattern: 5′-G[T]C[A]GG[T]GGGCG[A,T]-3′ [17, 25]) and 77 predicted promoters (refer to 74 miRNA genes) were found canonical Egr1 motifs. Among the 77 promoters with motif, 33 (42.9%, refer to 31 miRNA genes) have real Egr1 binding signals based on our ChIP-Seq experiment and they are listed in Table 4. According to the data, we found that hsa-mir-142 has 45 consensus sites dispersed among its promoters and correspondingly, 10 Egr1 binding regions identified by ChIP-Seq were enriched with reads count above nine. miR-142 is one of five miRNAs (miR-142, miR-144, miR-150, miR-155 and miR-223) which are highly specific for hematopoietic cells according expression comparison [35].

Table 4

The list of miRNA promoters with both Egr1 motif and ChIP-Seq peak. mir_TSS_start and mir_TSS_end is the start and stop of the predicted promoters. score is the cumulative score for each predicted promoter. motif_num is the count of the Egr1 motifs in the promoter regions of miRNAs.

mir_name	mir_TSS_start	mir_TSS_end	score	strand	chr	reads_count	peak_summit	motif_num
hsa-mir-200b	1088265	1090140	10	+	chr1	11	1091461	2
hsa-mir-200b	1088265	1090140	10	+	chr1	10	1094077	2
hsa-mir-200b	1088265	1090140	10	+	chr1	9	1086834	2
hsa-mir-30e/30c-1(cluster-hsa-mir-30a)	40929851	40930051	19	+	chr1	25	40930041	1
hsa-mir-92b	153429179	153434335	10	+	chr1	13	153430306	3
hsa-mir-92b	153429179	153434335	10	+	chr1	11	153428517	3
hsa-mir-146b	104184797	104186757	10	+	chr10	23	104185405	2
hsa-mir-132/212 (cluster-hsa-mir-132)	1899412	1901670	10	−	chr17	23	1900085	4
hsa-mir-132/212 (cluster-hsa-mir-132)	1899412	1901670	10	−	chr17	23	1903969	4
hsa-mir-10a	44014310	44014709	15	−	chr17	51	44010351	1
hsa-mir-10a	44014310	44014709	15	−	chr17	9	44014535	1
hsa-mir-142	53379405	53785839	10	−	chr17	14	53757151	45
hsa-mir-142	53379405	53785839	10	−	chr17	14	53516188	45
hsa-mir-142	53379405	53785839	10	−	chr17	13	53765628	45
hsa-mir-142	53379405	53785839	10	−	chr17	11	53760775	45
hsa-mir-142	53379405	53785839	10	−	chr17	11	53383253	45
hsa-mir-142	53379405	53785839	10	−	chr17	10	53740717	45
hsa-mir-142	53379405	53785839	10	−	chr17	10	53671689	45
hsa-mir-142	53379405	53785839	10	−	chr17	9	53739084	45
hsa-mir-142	53379405	53785839	10	−	chr17	9	53711312	45
hsa-mir-142	53379405	53785839	10	−	chr17	9	53747880	45
hsa-mir-24-2/27a/23a (cluster-has-mir-23a)	13818427	13819944	1	−	chr19	10	13814707	2
hsa-mir-642	50863241	50863441	20	+	chr19	18	50863778	1
hsa-mir-99b/let-7e/125a (cluster-has-mir-99b)	56883717	56886813	10	+	chr19	17	56886253	5
hsa-mir-99b/let-7e/125a (cluster-has-mir-99b)	56883717	56886813	10	+	chr19	10	56884005	5
hsa-mir-149	241040800	241041000	9	+	chr2	10	241042627	1
hsa-mir-149	241043770	241044582	10	+	chr2	10	241042627	1
hsa-mir-124-3	61276703	61277978	9	+	chr20	33	61280947	2
hsa-mir-124-3	61276703	61277978	9	+	chr20	10	61274186	2
hsa-mir-124-3	61278418	61280978	10	+	chr20	33	61280947	3
hsa-mir-124-3	61278418	61280978	10	+	chr20	10	61274186	3
hsa-mir-185	18388530	18388730	20	+	chr22	38	18388406	1
hsa-mir-185	18388530	18388730	20	+	chr22	36	18384025	1
hsa-mir-185	18388530	18388730	20	+	chr22	10	18392254	1
hsa-mir-301b/130b (cluster-hsa-mir-130b)	20335495	20337979	10	+	chr22	20	20338454	1
hsa-mir-301b/130b (cluster-hsa-mir-130b)	20335495	20337979	10	+	chr22	19	20341434	1
hsa-mir-301b/130b (cluster-hsa-mir-130b)	20335495	20337979	10	+	chr22	9	20332712	1
hsa-mir-301b/130b (cluster-hsa-mir-130b)	20335495	20337979	10	+	chr22	9	20330784	1
hsa-mir-425/191 (cluster-hsa-mir-191)	49029742	49035883	10	−	chr3	11	49033588	3
hsa-mir-148a	25955148	25959801	10	−	chr7	12	25957595	2
hsa-mir-96/183 (cluster-hsa-mir-96)	129204548	129210298	10	−	chr7	23	129207821	7
hsa-mir-96/183 (cluster-hsa-mir-96)	129204548	129210298	10	−	chr7	9	129200022	7
hsa-mir-335	129919123	129919323	10	+	chr7	13	129919298	2
hsa-mir-596	1752393	1753192	10	+	chr8	13	1752487	1
hsa-mir-32	110921946	110922146	19	−	chr9	9	110922058	1
hsa-mir-455	115957951	115958151	20	+	chr9	14	115956609	1
hsa-mir-219-2	130194040	130195803	10	-	chr9	9	130196326	2

3.5. Egr1 Bound-miRNAs' Expression in PMA-Induced K562 Cells

The set of Egr1 bound-miRNA genes identified by ChIP-Seq do not necessarily mean that they are regulated by the transcription factor Egr1. However, if the expression level of Egr1 bound-miRNAs fluctuates concomitantly with the expression changes of Egr1 during MK differentiation induced by PMA, it may suggest the function relevance of their interactions. Recently, Navarro et al. compared 286 miRNAs expression at different time point (0 days, 2 days, and 4 days) after PMA treatment of K562 cells by miRNA microarray [11]. The expression level of most miRNA genes increased with PMA induction according to the microarray analysis. Since the K562 cells used in our ChIP experiment were also treated by PMA and the snapshot of interactions between Egr1 and its target DNAs were captured after 1.5 h PMA treatment, we directly applied this expression profile of miRNAs to investigate if the binding by Egr1 contributes to the regulation of miRNA expression. According to the microarray analysis of protein-coding genes (data not shown), we found that the expression of Egr1 in PMA-treated K562 cells increased to 7.5-fold of that in the untreated cells. We found that among those miRNAs with Egr1 binding sites, the expression level of 46 mature miRNAs had increased at least two fold after two days PMA treatment and 16 of them showed significant expression changes with at least five fold increase based on Navarro's miRNA microarray analysis [11]. These 16 mature miRNAs were hsa-miR-181a, hsa-miR-181b, hsa-miR-212, hsa-miR-132, hsa-miR-135b, hsa-miR-141, hsa-miR-152, hsa-miR-199b, hsa-let-7g, hsa-let-7b, hsa-miR-149, hsa-miR-153, hsa-miR-346, hsa-miR-375, hsa-miR-200a, and hsa-miR-200b. miR-181b and miR-181a are members of miR-181 cluster, which increased 133-fold and 19-fold, respectively. miR-132 and miR-212, in the miRNA miR-132/miR-212 cluster, share a promoter, whose region has two Egr1 binding sites. They were also highly up-regulated with the former increasing 10-fold and the latter increasing 6-fold. These data may suggest that due to PMA stimulating, coupling to Egr1 up-regulation during megakaryocytic differentiation, the expression of miRNA genes, which are bound and activated by the transcription factor Egr1, also increased.

3.6. Analyzing Egr1-miRNA Binding Events by ChIP-PCR

ChIP-PCR assay is a conventional and reliable measure of whether a locus is a true binding target of a protein of interest. To experimentally validate the enriched regions from ChIP-Seq, we performed ChIP-PCR assays for 12 miRNA genes. As showed in Figure 3, significant binding signals of the transcription factor Egr1 can be observed for 12 miRNA genes within the PMA-induced K562 cells, while for the input DNA, there is no signal can be observed (see Figures 3(a) and 3(b)). The results here can partially verified the ChIP-Seq experiment, however, it cannot be concluded that the sensitivity of this ChIP-Seq is as high as 100%. It can be easily envisioned if more ChIP-Seq found targets were subjected to ChIP-PCR, false positive rate may arise.

To further check that if the expression or the increase of expression levels of miRNAs can partially attribute to the binding of the transcription factor Egr1, we also performed ChIP-PCR analysis for 12 miRNA genes before PMA treatment (see Figure 3) and made a qualitative comparison for their Egr1-miRNA binding levels between these two cell conditions. Among the 12 miRNAs, 5 miRNA genes (the right column in Figure 3) showed at least 5-fold expression changes before and after PMA-treatment according to Navarro's microarray analysis for microRNA genes [11]. The ChIP-PCR results for two cell conditions showed that some Egr1-miRNA binding levels significantly changed after PMA treatment as well as those showed similar binding signals for different cell conditions. Among the 5 miRNAs whose expression levels significantly increased after PMA-treatment, miR-135b, miR-141, miR-152 showed no Egr1 binding signals before PMA treatment (see the right ChIP-PCR results in Figure 3), however, the PMA induced expression increase of transcription factor Egr1 subsequently turned on the Egr1 binding switch to these three miRNA genes. As analyzed above, miR-135b, miR-141, miR-152 were among the 16 mature miRNAs with at least five fold expression increase after PMA treatment, so before PMA induced, they still express at a relatively low level due to other factors although without the factor Egr1's transcription. We suspect that like general housekeeping genes which have low levels of genic transcription whereby gene products are constitutively produced at low levels [36], the low level expression of miRNA genes may also functionalize as housekeeping miRNAs. Once bound by the key transcription factors and transcriptionally stimulated, the explosion of their expression is likely to functionally correlate with the cell fate. For miRNA genes let-7g and miR-199b, we can observe an Egr1 ChIP binding signal enhancement before and after PMA treatment, which may mainly be caused by up-regulation of Egr1. As let-7g and miR-199b are also among the 16 mature miRNAs with significant fold changes, we can conclude that Egr1 may actively involve in these miRNA genes transcription after PMA induction.

For the left seven miRNA genes, based on Navrro' microarray analysis, miR-10a, miR-95, miR-326, miR-23b, miR-130b also showed two fold expression increase, despite not as high as that of the above five miRNA genes, and the expression level of miR-25 and miR-17 did not show any changes. The ChIP-PCR results of six of them within the nontreatment K562 cells showed positive binding signals. The similar binding signal before and after PMA treatment of K562 cell between the transcription factor Egr1 and these, or may be not limited to, six miRNA genes, suggests that the great expression increases of Egr1 may not be able to increase all the Egr1-miRNA bindings and therefore stimulate their transcriptions. It also suggests that only a part of miRNA genes, which form a regulatory pattern, contributes to this PMA-induced differentiation process. As the gene expression is such a complicated cellular process, different regulatory patterns may exist at distinct cellular processes.

4. Discussion

Cellular activities can be considered as a succession of hierarchically acting regulatory states [37]. Identification of miRNA's targets is necessarily of great significance, but it is also significant to detect what factors determine miRNAs expression. In this report, we have taken the approach of ChIP-Seq to identify all the miRNAs with known predicted promoters which are bound by the transcription factor Egr1. We found that the promoters of 124 distinct miRNAs have at least one Egr1 binding site. Besides, EGR1 binding signal was also found near extra 63 premiRNAs. Combing with the miRNA expression profile of PMA-induced K562 cells, we were able to find that the binding of the transcription factor Egr1 may be the reason for expression changes of parts of miRNAs. We believe that these will be instrumental in unraveling the mechanism of transcriptional regulation of miRNA genes.

Lineage specification is critical in developmental biology and evidences showed that miRNAs mediate control of cell fate in megakaryocy-erythrocyte progenitors [15]. As accumulative evidence showed that miRNAs expression differs between PMA- and hemin-induced K562 differentiation, the possible driven power is under revealing. EGR1 is a zinc finger family TF who is known as “immediate-early response factor” and has been implied in megakaryocytic differentiation of K562 cells. Evidences were also found that Egr1 activation is correlated with downregulation of erythroid-specific genes and up-regulation of the megakaryocyte-spcefic gene CD41a [8]. The identification of candidate EGR1-targeted miRNAs through this high-throughput approach significantly broadens the repertoire of biological effects attributable to the TF. As many of our annotated miRNAs are known to play crucial rules in hematopoietic differentiation (miR-27a, miR-10a, etc.), the role of Egr1 in hematopoietic differentiation is also implicated here.

During K562 cell erythroid/MK differentiation, there are miRNAs expression signatures specific to each cell lineage. For example, a subset of 12 highly up-regulated miRNAs (miR-34a, miR-181b, miR-299-5p, miR-134, miR-375, miR-181a, miR-139, miR-222*, miR-409-3p, miR-221, miR-212, miR-132) from Navarro's miRNA microarray in TPA-treated K562 cells was found to be specific to MK differentiations, while none of them was up-regulated during hemin-induced erythroid differentiation by Northern blot assay [11]. miR-223 was also up-regulated during PMA-induced megakaryocytic differentiation but downregulated during hemin-induced erythroid differentiation [12]. Induction of megakaryocytic differentiation in K562 cells by PMA markedly increased miR-27a expression; however, downregulation of miR-27a expression occurs upon induction of K562 cells toward erythroid differentiation. Runx1 was the transcription factor which plays a role in the erythroid/megakaryocytic lineage determination through mediating regulation of miR-27a [10]. PMA-induced activation of Egr1 expression has been proved to involve in megakaryocytic differentiation of K562 cells but the regulation pathway has not been investigated. We therefore investigate several Egr1-bound miRNAs which behave differently in erythroid/megakaryocytic cell lineage to see if Egr1 plays a role in erythroid/megakaryocytic lineage determination of K562 cells. The promoters of miRNA cluster hsa-mir-181a (hsa-mir-181a/181b) and cluster hsa-mir-132 (hsa-mir-132/212) cluster were found Egr1 binding sites according to our ChIP-Seq experiment. These two miRNA clusters are possible candidates through which Egr1 acts as a lineage determinator. Our ChIP-Seq data also showed that near the two predicted promoters of miR-27a, each had an Egr1 binding sites, being enriched by 13 reads and 10 reads, respectively. DNA sequence of the second ChIP-Seq enriched regions was found canonical Egr1 motif (degenerate pattern: 5′-G[T]C[A]GG[T]GGGCG[A,T]-3′) and we found three sites with motif, TAGGGGGCG, TCGTGGGCG, and TCGGGGGCA. So, we suspect that except that Runx1 exerts transcriptional stimulation on miR-27a expression, the binding of Egr1 may also regulate miR-27a expression level and therefore contribute to lineage determination of K562 cells.

Yang et al. characterized and validated that miRNAs of miR-103 and miR-10a exhibited downregulation after hemin induction [13], while Navarro's miRNA microarray analysis showed that miR-103 and miR-10a increased two fold after PMA induction [11]. These differences imply that like miR-27a, the two miRNAs may also participate in lineage determination of K562 cells. Meanwhile, our Egr1 ChIP-Seq results supported that Egr1 binds to the transcription initiation regions of hsa-mir-103-1 (reads count:14, canonical motif: TAGGGGGCG), hsa-mir-103-2 (reads count:25), hsa-mir-10a (two binding sites, reads count:51; canonical motifs: GCGGGGGCA, GCGTGGGCG, TCGGGGGCA, GCGTGGGCA; reads count:9; canonical motifs, GAGGGGGCG), so Egr1 would be the candidate transcription factor which regulates the expression changes of miR-103 and miR-10a. Our ChIP-PCR assay confirmed the Egr1-mir-10a binding. All these data provides key clues for establishment of regulation network of K562 cell differentiation. The relationship between Egr1 and its targeting miRNAs discussed above is summarized in Figure 4.

Figure 4

Schematic diagram of miRNA related Egr1's contribution to erythroid/megakaryocytic lineage determination of K562 cells. mir-27a, mir-10a, mir-103-1 and mir-103-2 were found Egr1 binding signals in this study. Based on literature data, when K562 cells were induced by PMA, the expression levels of these miRNAs increased. On the contrary, their expression level were found decreased when K562 cells were induced to erythroid differentiation.

5. Conclusion

By ChIP-Seq, we provide a global binding profile between the transcription factor Egr1 and its targeting miRNA genes in PMA-treated K562 cells and found that EGR1 binds to promoters of 124 miRNAs and precusors of 125 miRNAs. Among the miRNAs which are bound by Egr1, miR-10a, miR-25, miR-23b, miR-135b, miR-130b, miR-326, miR-141, miR-152, miR-199b, miR-95, and let-7g are validated according to ChIP-PCR.

Figure 1: The GOEAST graphical output of enriched GO terms in the biological process category for the targets of 38 miRNAs. Boxes represent GO terms, labeled by its GO ID, term definition, and detailed information, organized as ‘q/m|t/k (P-value)'. Significantly enriched GO terms are marked yellow. The degree of color saturation of each node is positively correlated with the enrichment significance of the corresponding GO term. Arrows represent connections between different GO terms. Red arrows represent relationships between two enriched GO terms, black solid arrows represent relationships between enriched and unenriched terms and black dashed arrows represent relationships between two unenriched GO terms.

Figure 2: The GOEAST graphical output of enriched GO terms in the cellular component category for the targets of 38 miRNAs.

Figure 3: The GOEAST graphical output of enriched GO terms in the molecular function category for the targets of 38 miRNAs.

Talbe 1: The list of ChIP-PCR primers for 12 miRNAs.

Table 2: miRNA promoter and pre-miRNAs annotion of Egr1 binding sites.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.