Macrophages are immune cells involved in various biological processes including host defence, homeostasis, differentiation, and organogenesis. Disruption of macrophage biology has been linked to increased pathogen infection, inflammation and malignant diseases. Differential gene expression observed in monocytic differentiation is primarily regulated by interacting transcription factors (TFs). Current research suggests that microRNAs (miRNAs) degrade and repress translation of mRNA, but also may target genes involved in differentiation. We focus on getting insights into the transcriptional circuitry regulating miRNA genes expressed during monocytic differentiation.
We computationally analysed the transcriptional circuitry of miRNA genes during monocytic differentiation using in vitro time-course expression data for TFs and miRNAs. A set of TF→miRNA associations was derived from predicted TF binding sites in promoter regions of miRNA genes. Time-lagged expression correlation analysis was utilised to evaluate the TF→miRNA associations. Our analysis identified 12 TFs that potentially play a central role in regulating miRNAs throughout the differentiation process. Six of these 12 TFs (ATF2, E2F3, HOXA4, NFE2L1, SP3, and YY1) have not previously been described to be important for monocytic differentiation. The remaining six TFs are CEBPB, CREB1, ELK1, NFE2L2, RUNX1, and USF2. For several miRNAs (miR-21, miR-155, miR-424, and miR-17-92), we show how their inferred transcriptional regulation impacts monocytic differentiation.
The study demonstrates that miRNAs and their transcriptional regulatory control are integral molecular mechanisms during differentiation. Furthermore, it is the first study to decipher on a large-scale, how miRNAs are controlled by TFs during human monocytic differentiation. Subsequently, we have identified 12 candidate key controllers of miRNAs during this differentiation process.
The mononuclear phagocyte system is defined as a family of cells comprising of bone marrow progenitors and is derived from hematopoietic stem cells. Hematopoietic stem cells sequentially differentiate into monoblasts, promonocytes, monocytes and terminal macrophage cells . The human monocytic leukemic cell line, THP-1 , is an accepted model system utilised to explore molecular events surrounding monocytic differentiation. Phorbol 12-myristate 13-acetate (PMA) induces the differentiation of monocytic THP-1 cells into macrophages/mature THP-1 cells . Before inducing differentiation, PMA first inhibits cell growth and blocks THP-1 cells in G1-phase of the cell cycle by up-regulating the expression of p21WAF1/CIP1, enhancing binding of the SP1 factor to the p21WAF1/CIP1 promoter. PMA inhibition of cell growth is mediated by several signalling pathways such as MAPK and ROS-dependent Raf/MEK/ERK pathway . Human monocytic maturation incorporates lipid and protein metabolic processes together with several G-protein coupled receptors (GPCRs) .
Differential gene expression that results in human monocytic differentiation is regulated by numerous interacting transcription factors (TFs) [4–6]. Current research suggests that microRNAs (miRNAs) target several genes that are differentially expressed in the differentiation process . miRNAs are ~22 nucleotides (nt) long non-coding RNAs, which play a key role in the repression of translation and degradation of coding mRNA [8–12]. Several computational tools are available for prediction of miRNA targets [9, 13–16].
Canonical miRNA biogenesis begins with the transcription of the pri-miRNA by RNA polymerase II [17–19]. These pri-miRNAs are cleaved into 60~70 nt pre-miRNAs by the microprocessor complex Drosha (RNase II endonuclease) and DGCR8 (a double-stranded RNA binding protein) [20, 21]. Pre-miRNAs are then exported to the cytoplasm with the help of Exportin-5 and its co-factor RanGTP . Dicer, a RNase III endonuclease, cleaves 22-nucleotide from the Drosha cleavage site to yield the mature miRNA [8, 23]. The generation of pri-miRNA by RNA polymerase II suggests that miRNA genes are controlled through the same regulatory machinery as the protein coding genes.
A straightforward analysis of the transcriptional regulation of miRNA genes is difficult. Even though most miRNAs have their own transcriptional units , it is known that several miRNAs are transcribed together as a single pri-miRNA [24–26]. These clustered miRNAs are thus co-regulated. On the other hand, miRNAs can also be transcribed together with a protein-coding host gene . In addition, a mature miRNA can be produced from several locations in the genome [8, 27]. Furthermore, it is not clear how to define the regulatory regions for miRNA genes. Current research suggests that at transcription start sites (TSSs) of genes, histones are generally trimethylated at lysine 4 residues [28, 29]. This has led to a potential definition of promoter regions for miRNAs  in human embryonic stem cells using such determined TSSs as the reference points.
As the transcriptional regulation of miRNAs is not well understood, we focus our study on the analysis of transcriptional regulation of miRNAs during monocytic differentiation. Gene expression of miRNAs and TFs was measured prior to PMA stimulation and over a 96 hour time-course post-PMA stimulation. We first utilised a general method to identify miRNAs whose expression levels differed due to PMA stimulation in THP-1 cells. We extracted promoter regions for these miRNAs and computationally mapped TF binding sites (TFBSs) to the promoter sequences. We made use of a time-lagged expression correlation analysis [31, 32] to evaluate the predicted TF→miRNA associations by combining our in silico TFBS analysis with the measured in vitro expression data. This kind of a time-lagged expression correlation analysis has been utilised before to either predict or score TF→gene or gene→gene associations [33–35]. From these TF→miRNA associations we identified 12 TFs likely to play a central role in regulating miRNAs throughout the considered differentiation process. Six of these 12 TFs (ATF2, E2F3, HOXA4, NFE2L1, SP3, and YY1) have not been previously described as important for monocytic differentiation. The remaining TFs, CEBPB, CREB1, ELK1, NFE2L2, RUNX1, and USF2, although known to be involved in monocytic differentiation, were not known to play role in transcription regulation of miRNAs in this process. We concluded the analysis by highlighting several inferred regulatory networks that suggest interplay of TFs, miRNAs, and miRNA targets and that are likely to have an impact on the differentiation process.
To the best of our knowledge this research is the first large-scale study that attempts to decipher the transcriptional circuitry that regulates the expression of miRNAs during human monocytic differentiation and identifies potential new avenues for further research.
Results and Discussion
In what follows we present and discuss the main results of the study. Figure 1 gives an overview of the analysis steps. First, we analysed the miRNA expression data to identify miRNAs that are mostly affected by the PMA stimulation. We extracted promoter regions for the identified miRNAs and predicted TFBSs in these regions. Subsequently, we scored each predicted TF→miRNA association using a time-lagged expression correlation analysis to get a measure of reliability for the predicted associations. Afterwards, we statistically identified TFs that are likely to play a central role in regulating miRNAs during the monocytic differentiation process. Finally, for several miRNAs we investigated the predicted transcription regulations and their potential influence on the differentiation process.
Identification of miRNAs most influenced by PMA stimulation
We are interested in the transcriptional regulation of those miRNAs whose expression is most influenced by the PMA stimulation. Three biological replicates of miRNA expression data provided measured expression levels at nine time-points after PMA stimuli and a zero hour control prior to PMA stimulation (see Methods). We required that two criteria were met for the inclusion of a miRNA expression time-series ('expression series' in further text) in the analysis:
i/ Expression of the miRNA should be denoted as "present" in at least one time point, otherwise we assume that the expression series for the miRNA is invalid. In this manner, we identified 155, 238, and 191 miRNAs and associated expression series for the first, second, and third replicate, respectively.
ii/ For a miRNA, i/ must hold true in at least two of the three biological replicates.
The expression values of different biological replicates for a miRNA that satisfy the criteria have been averaged at each time point to generate one expression series per miRNA. This resulted in expression series for 187 miRNAs (see Methods).
In order to find the set of 'most relevant' miRNAs, we calculated for each of the 187 identified miRNAs the log2fc (fc standing for fold-change relative to time zero) at each of the measured time points (see Methods). A miRNA we considered to be influenced by PMA stimulation if its log2fc > 1 or log2fc < -1 at any measured time point post-PMA stimulation (see Figure 2). Figure 2 shows that the majority of the miRNA expression does not change significantly over time and is confined within the selected threshold values. We found a total of 81 miRNAs that satisfied this criterion. To determine those miRNAs that deviate from the baseline expression we proceeded as follows. For each time point t where log2fc > 1 or log2fc < -1 was satisfied for a miRNA, we calculated the difference dt of the expression et at time point t and its expression e0 at the zero time point. We sub-selected those miRNAs for which abs(dt) > 0.1 for at least one time point. This resulted in a set of 53 miRNAs for which we are more confident that their expression is affected by the PMA stimulation.
The fc does not take the level of expression into account. It is important to note that miRNAs that have very high expression level and change only little over time might have a strong biological impact, even though this is not reflected by variation in the expression level. Our approach, based on fc excludes such cases. On the other hand, miRNAs with very low expression levels might have high fc values that may suggest a strong biological impact, even though this may be arguable since the changes in expression levels could be very small. Hence, we introduced a second threshold for the difference in expression values of 0.1, even though no guideline exits for choosing this threshold.
TFBS analysis of miRNA promoter regions
Promoter regions of miRNAs are regions of DNA where TFs bind to regulate the transcription of miRNA genes into pri-miRNAs. A pri-miRNA can be associated to several promoter regions derived from different TSSs. The transcriptional control of TFs is towards the pri-miRNA that can be cleaved into several pre-miRNAs . Thus, we consider the miRNAs that form such clusters to be generally regulated in the same manner.
Marson et al.  defined promoter regions of miRNAs using TSSs determined based on trimethylated histones. We chose to analyze these promoter regions. For 34 of the 53 earlier identified mature miRNAs we were able to extract 38 promoter regions for 37 associated miRNAs (see Methods and Additional Files 1).
To map TFBSs to the 38 promoters we utilised TRANSFAC Professional database (version 11.4) [37, 38]. TRANSFAC's 522 mammalian minimum false positive matrix profiles of binding sites were mapped to the promoter regions (see Methods). These matrices, which correspond to the predicted TFBSs, are associated with TFs that possibly bind these TFBSs (see Methods). By mapping the matrices to their corresponding TFs, we obtained 5,788 unique TF→miRNA associations for 673 TFs and 37 miRNAs.
Evaluation of predicted TF→miRNA associations
Each predicted TF→miRNA association has been evaluated to get the most accurate picture of miRNA gene regulation during human monocytic differentiation. The result of this evaluation relates to our confidence that we are dealing with a genuine TF→miRNA association. The evaluation was based on time-lagged expression correlation between the gene expression series of the TF and that of the mature miRNA (see Methods). Expressions for miRNAs and TFs have been measured in human THP-1 cells prior PMA stimulus at one time point and post-PMA stimulus at non-equidistant time points up to 96 hours.
We interpolated the expression series for each of the 34 mature miRNAs using half an hour steps (see Additional Files 2). In concordance with the miRNA expression data, we averaged the TF qRT-PCR expression series over the two biological replicates at the same time points and interpolated each expression series using half an hour steps (see Methods). In this manner, we derived expression series for 2,197 TFs (see Methods).
The TF→miRNA associations were inferred from TFBS analysis of promoter regions of miRNA genes. From the predicted 5,788 TF→miRNA associations, we discarded all associations for which we do not have expression data for the TF in the above mentioned averaged expression set. After calculating Pearson's correlation coefficient (PCC) for each TF→miRNA associations using a time-lagged expression correlation analysis and the interpolated expression data for TFs and mature miRNAs, we finally derived a set of 1,989 TF→miRNA associations (see Additional Files 3) for 37 miRNAs and 258 TFs (see Additional Files 4), each associated with a PCC value (see Methods). In Figure 3A we show the number of TF→miRNA associations that have PCCs equal to or greater than selected thresholds. As expected, the number of associations steadily decreases with increasingly stringent PCC thresholds.
Previous research demonstrated that the regulatory effects of a TF on its target genes is not instantaneous but with a time-lag [39–41]. Unfortunately, the correct time-shifts are undetermined. In our analyses, we incorporated time-shifts in a range from 0.5 hours to six hours to allow for a sufficient time-delay for the regulation by the TF to exert an effect on the transcription of its target miRNA genes. We calculated for each of the 1,989 TF→miRNA associations the most favourable time-shift and with this, the time-lagged PCC of expression as the score for the association (see Methods). The higher the absolute value of the PCC for an association, the more confidence we have that the association is genuine and could play an important role in the differentiation process. For each miRNA/miRNA-cluster and its regulating TFs, the maximum PCC s were calculated individually (see Methods). Other approaches considered all TFs that regulate a gene to extract a common time-shift for all TFs and the gene  or compute the best time-shift depending on known examples of regulation . Up to now, too few experimentally verified examples of TFs that regulate miRNAs are known, thus a model to introduce the "correct" time-shift could not be inferred. Furthermore, certain miRNAs were predicted to be clustered and share common promoter regions. Hence, a time-shift common to all miRNAs in a cluster was calculated for each of the associated TFs. As a criterion, common time-shifts were only taken into account if all PCC s between the TF and all miRNAs that form the cluster had the same sign (e.g. all positive or all negative) to avoid contradicting effects of the same TF on different miRNAs of the cluster. TF→miRNA associations where all considered time-shifts were discarded (because of sign disagreement) were excluded from further analysis.
Identification of TFs central to regulation of miRNA genes
In order to find the TFs that have the most influence on miRNAs during the differentiation process, we analysed TFs corresponding to TF→miRNA associations having the highest absolute PCC. We ranked 1,989 TF→miRNA associations according to the absolute value of their corresponding PCC s. From the ranked associations we selected the upper quartile (with the highest absolute PCCs). In this manner, we obtained 498 associations with an absolute PCC greater than 0.775 (see Figure 3B). The 498 associations are formed by 111 unique TFs and 35 unique miRNAs. TFs that appear significantly more often in the upper quartile of associations are assumed to more likely play a central role in regulating miRNAs during the differentiation process. We utilised the one-sided Fisher's exact test to calculate the Bonferroni-corrected p-value for enrichment of each TF in the subset of 498 associations, in contrast to the remaining set of 1,491 associations. The correction factor is the number of unique TFs (258) in the complete set of all associations (1,989). In this manner, we found that 12 TFs are statistically significantly enriched in the set of 498 associations with a corrected p-value smaller than 0.01 (see Table 1). Six of these 12 TFs (ATF2, E2F3, HOXA4, NFE2L1, SP3, and YY1) have not been previously described as important for monocytic differentiation. The remaining TFs (namely, CEBPB , CREB1 , ELK1 , NFE2L2 , RUNX1 , and USF2 ) are known to play role in monocytic differentiation, but not explicitly as regulators of miRNAs in this process.
Our approach attempts to identify the most dominant TFs that putatively regulate miRNAs from the selected subset of TF→miRNA associations with highest PCCs. The complete set of 1,989 TF→miRNA associations consists of many associations with a low PCC (see Figure 3). In order to be able to focus on associations that are most likely genuine, we sub-selected the associations with the highest PCCs. At the same time we did not want to restrict the analysis to too few associations, so as to be able to deduce the general participants in the transcriptional regulation process of miRNAs. Consequently, we selected the upper quartile of TF→miRNA associations ranked based on decreasing absolute values of PCC as a reasonable compromise between sensitivity and specificity.
Transcriptional circuitry of miRNAs during monocytic differentiation
To shed light on a portion of the molecular underpinnings of monocytic differentiation we will discuss the TF→miRNA associations for miRNAs that have been described earlier to be affected by PMA stimulation. In this manner, we can confer whether or not our findings correspond to the published scientific findings and further introduce novel TF→miRNA associations. An overview of the regulatory effects of the TF subset (defined above) on the miRNAs is presented in Figure 4. The figure shows each association, from within the subset of the upper quartile of associations, in form of a coloured dot in a heat-map style of format using the TIGR Multiexperiment Viewer (version 4.3) (TMEV, ). We can observe certain clusters of miRNAs that are regulated by the same set of TFs. In the following discussion, we mainly focused on the upper quartile of TF→miRNA associations and on the TFs illustrated in Figure 4 that we have identified to be central to monocytic differentiation. For the sake of completeness, we also discuss several TFs that are known to be regulators of certain miRNAs, even though they might not appear in our set of "best" TF→miRNA associations. Subsets of miRNAs that have support through literature establishing their expression during PMA-induced differentiation are discussed. All network graphics in the following figures have been produced with the help of Cytoscape  and all pathway analyses were based on KEGG  using DAVID .
Fugita et al. demonstrated that mir-21 is expressed during PMA-induced differentiation in the human promyelocytic leukaemia cell line, HL-60 . Our expression data demonstrate that miR-21 is up-regulated during the differentiation process (see Figure 5C). Our correlation data suggest that several of the 12 TFs (see above), which we identified as being central to the considered differentiation process bind in the promoter region of miR-21 (YY1, NFE2L2, ATF2 and NFE2L1, see Figure 4). Additionally, the binding of TFs, AP-1/c-jun, and c-fos to the promoter region of mir-21 has been demonstrated via chromatin immunoprecipitation (ChIP) in the human promyelocytic leukaemia cell line, HL-60 after 4 h PMA induction . Our TFBS analysis results suggests the binding of several members of the JUN-FOS family (JUN, JUNB, JUND, FOS, FOSB, FOSL1, and FOSL2) to the promoter region of mir-21, even though they do not appear in the upper quartile of TF→miRNA associations. The expression data for the JUN family members displayed continued up-regulation for 96 hours, whereas FOS family members, with exception of FOSL1, were down-regulated after 4 hours (see Figure 5B). AP-1/c-jun form a complex with the JUN-FOS family members during transcription, and AP-1/c-jun is known to be activated by PMA induction which is supported by our findings (data not shown) . Fugita et al. also demonstrated that AP-1 and SPI1 synergistically mediate the transcriptional process . Our method predicted a SPI1 binding site in the promoter region of the mir-21 gene. The time-lagged expression correlation analysis demonstrated that SPI1 is highly correlated to miR-21 (PCC = 0.798; see Figures 5B and 5C).
miR-21 has been found to display anti-apoptotic functioning and targets tumour suppressor genes, like the PTEN gene in human hepatocellular cancer cells  and the tropomyosin 1 (TPM1), PDCD4, and maspin gene in the human breast cancer cell line, MDA-MB-231 . The miR-21's predicted targets (see Methods) were found to be primarily involved in pathways such as TGF-β signalling pathway, MAPK signalling pathway and the JAK-STAT signalling pathway (see Figure 5A). The TGF-β signalling pathway and MAPK signalling pathway is primarily involved in differentiation, proliferation, apoptosis and developmental processes, while the JAK-STAT signalling pathway is involved in immune responses. We found that several TFs such as ATF2, FOS, JUN and JUND included in the predicted TF→mir-21 associations are involved in the MAPK signalling pathway (see Figure 5A).
Time-lagged expression correlation analysis demonstrated that NFE2L1 and SPI1 are highly correlated to miR-21 as opposed to YY1, NFE2L2, and ATF2, which have negative PCC s (see Figure 4). Besides JUN-FOS family members and SPI1 that are known to regulate the miR-21, our results suggest a novel NFE2L1→miR-21 association, which seems to play an important role in monocytic differentiation (see Figure 5A).
Rosa et al. reported that mir-424 is expressed during PMA-induced differentiation and that mir-424 is transcribed by SPI1 in the CD34+ human cord blood cells and CEBPA (C/EBPα) blocks SPI1 induced dendritic cell development from CD34+ human cord blood cells by displacing the co-activator c-Jun [55, 56]. The up-regulation of miR-424 (see Figure 6C) leads to the repression of NFIA which allows for the activation of differentiation specific genes such as M-CSFr (CSF1R) . Furthermore, the pre-mir-424 is transcribed together with pre-mir-503 and pre-mir-542 as one transcript. These pre-miRNAs form the mature miRNAs miR-424, miR-503, miR-542-5p, and miR-542-3p. Our data suggest that several of the 12 TFs (see above), which we identified as being central to the considered differentiation process bind in the promoter region of miR-424 (RUNX1, E2F3, SP3, YY1, NFE2L2, CREB1, ATF2, USF2, ELK1, CEBPB and HOXA4; see Figure 4). Figure 4 shows that mir-424 and mir-542 are regulated by the same TFs and are thus as well clustered in the heat-map. However, mir-503, part of the same cluster and thus subject to the same regulations, is not displayed in Figure 4. This is a consequence of the expression data obtained for miR-503 causing the PCC s for the TF→miRNA associations to decrease and thus not being part of the top quartile of associations (see above). We further predicted a SPI1 and CEBPA binding site in the promoter region of these clustered miRNAs, which corresponds to findings reported by Rosa et al. . SPI1 is positively correlated to miR-424 and CEBPA negatively. Furthermore, both associations are not within the top quartile of associations with highest PCCs. Nevertheless, these observations indicate that SPI1 enhances the expression of the mir-424 cluster and might work in conjunction with the other identified TFs to influence the miRNA's transcription.
The predicted targets of miR-424 were found to be involved in the same pathways as the targets of miR-21; the TGF-β signalling pathway, MAPK signalling pathway and JAK-STAT signalling pathway with additional pathways such as acute myeloid leukaemia and antigen processing and presentation, the p53 signalling pathway and SNARE interactions in vesicular transport. We found that several TFs included in the predicted TF→mir-424 associations, are involved in the MAPK signalling pathway (ELK1, ATF2), acute myeloid leukaemia (E2F3, RUNX1) and antigen processing and presentation (CREB1) (see Figure 6A).
The time-lagged expression correlation analysis demonstrated that of the 12 TFs (see above) only ELK1, USF2, CEBPB and HOXA4 were positively correlated to the expression of miR-424 (see Figure 4 and Figures 6B and 6C). Besides the earlier mentioned involvement of SPI1 in regulating mir-424 , our analysis suggests that ELK1, USF2, CEBPB and HOXA4 may be the TFs most likely responsible for the expression of mir-424 in monocytic differentiation (see Figure 6A).
Chen et al. reported that mir-155 is expressed during PMA-induced differentiation in the human promyelocytic leukaemia cell line, HL-60 . Our expression data demonstrate that miR-155 is up-regulated during the differentiation process (see Figure 7C). Our TFBS analysis data suggest that several of the 12 TFs (see above), which we identified as being central to the considered differentiation process, bind in the promoter region of miR-155 (SP3, NFE2L2, CREB1, NFE2L1 and ELK1; see Figure 4). Zeller et al. demonstrated binding of MYC to the promoter region of mir-155 in the human burkitt lymphoma cell line (P493-6) . Also, Yin et al. demonstrated binding of FOSB and JUNB to the promoter region of mir-155 using chromatin immunoprecipitation (ChIP) in the human B-cell line . miR-155 has been linked to Epstein-Barr virus (EBV) related diseases that are associated with latency during which only a subset of viral genes are transcribed with a set of EBV-encoded microRNAs. One such EBV gene is LMP1 which is a known oncogene that induces miR-155 in DeFew cells . Gatto et al. demonstrated the positive expression of miR-155 in DeFew cells induced with PMA and that the promoter region has two NF-κB (NFKB1) binding sites . Once again, our results predict the binding of several members of the JUN-FOS family to the promoter region of mir-155 but neither MYC nor NF-κB, this may be a consequent of the extracted regulatory region for mir-155, being incomplete. The expression data demonstrated the up-regulation of JUN-FOS (see Figure 5B) family members and NF-κB but a down-regulation of MYC (data not shown). Our observations indicate that JUN-FOS family enhances the expression of the miR-155 even though the predicted associations are not within the upper quartile of associations with highest PCCs.
MiR-155s predicted targets were found to be involved in the same pathways as the targets of miR-21 and miR-424; the TGF-β signalling pathway, MAPK signalling pathway and JAK-STAT signalling pathway with additional pathways such as acute myeloid leukaemia and Wnt signalling pathway (see Figure 7A). We found that several TFs such as ATF2 and ELK1, included in the predicted TF→mir-155 associations, are involved in the MAPK signalling pathway and CREB1 was found to be involved in antigen processing and presentation (see Figure 7A).
The time-lagged expression correlation analysis demonstrated that of the 12 TFs (see above) only NFE2L1 and ELK1 had TFBSs predicted within the promoter of miR-155 and were positively correlated to miR-155 (see Figure 4 and Figure 7B) and thus our findings propose that the NFE2L1→mir-155 and the ELK1→mir-155 associations are likely to be important to the monocytic differentiation process.
Members of the miRNA cluster mir-17-92 are known to be down regulated in the HL-60 cell line after PMA stimulation . The miRNA cluster on chromosome 13 contains several miRNAs (hsa-mir-17, hsa-mir-18a, hsa-mir-19a, hsa-mir-20a, hsa-mir-19b-1, and hsa-mir-92-1 (hsa-mir-92-1 excluded from analysis, due to ambiguous nomenclature)) that are transcribed as a single transcript. Our data shows that members of miR-17-92 are indeed down regulated after PMA stimulation and furthermore, that the lowest PCC between the expression series of the miRNA cluster members is ~0.86, which supports the cluster membership. Even though the function of miR-17-92 is largely unknown, lymphomas that express these miRNAs at a high level have reduced apoptosis [61, 62] and the miRNAs target multiple cell cycle regulators and promote G1→S phase transition . Expression of miR-17-92 is high in proliferating cells and is positively regulated, in part, by MYC (c-Myc) . E2F1, an activator of MYC, is itself a target of miR-17 and miR-20a  indicating that both MYC and E2F1 are under the control of a feedback loop. It has been experimentally shown that E2F3 activates the transcription of the miR-17-92 cluster [62, 36]. A model has been proposed that miR-17-92 promotes cell proliferation by targeting pro-apoptotic E2F1 and thereby favouring proliferation through E2F3 mediated pathways . Additionally, E2F3 is shown to be a predominant isoform that regulates miR-17-92 transcription . We show that after ranking PCC s of gene expression between miRNAs and putative TFs, E2F3 is the only TF appearing significantly associated with miR-17-92 within the upper quartile of TF→miRNA associations (see Figure 4).
Amongst a small set of eight predicted regulators (E2F1, E2F3, E2F4, TFAP2A, TFAP2B, TFAP2C, TFDP1, SP1), TFDP1 is known to form a heterodimer with another putative TF, E2F1 , implicating TFDP1/E2F1 complex as a regulator of miR-17-92 transcription.
In Figure 8A we show the putative regulation of miR-17-92 and its known effects in proliferation, differentiation and apoptotic pathways. Specifically, we predict E2F1 and E2F3 to regulate the miR-17-92 cluster. Figure 8B shows that expression of miR-17-92 members are correlated to E2F3 with a minimum PCC of ~0.9. Conversely, miR-17-92 members are correlated with E2F1 by a maximum PCC of ~-0.65. A disproportionately high PCC of E2F3 gene expression to miR-17-92 as compared to other TFs seems to support the claims made by Woods et. al. that E2F3 is indeed the predominant TF in this regulatory context . In addition, Cloonan et al. demonstrated that the pri-miRNA is cell cycle regulated, which supports the claim that the cluster is under the control of E2F family members, which are master regulators of the cell cycle . On inspection of the log2fc of TF gene expression over time (see Figure 8C) we observed that E2F3 is sharply up-regulated at 6 hours by ~2 fold, whilst its closely related and pro-apoptotic family member, E2F1, is down-regulated by a factor of ~5.7. After ~70 hours E2F3 and E2F1 gene expression levels return near to baseline, this corresponds to a progression towards a differentiated state before 96 hours post-PMA stimulation. Yet, regardless of the high PCC between E2F3 gene expression and the miR-17-92 cluster, the miRNA cluster is generally down-regulated (see Figure 8D). Acknowledging that the miRNA cluster targets and inhibits a well known RUNX1 (AML1) induced differentiation and proliferation pathway , these results strongly suggest that PMA stimulation disfavours both E2F1 induced proliferative and E2F1 induced apoptotic pathways. Whilst, equally, given that both ETS1 and ETS2, components of the above mentioned RUNX1 differentiation and proliferation pathway, are up-regulated (data not shown), these results indicate that PMA-treated monocytes up-regulate members of differentiation pathways. In light of the above findings we hypothesize, that since members of the AP-1 complex are concurrently up-regulated in the early stages after PMA stimulation, that monocytic differentiation is mediated by the M-CSF receptor-ligand RAS signalling pathway and indirectly controlled by miR-17-92 through the E2F TF family members E2F1 and E2F3. Generally, this hypothesis seems to be plausible, since RUNX1 is also an inhibitor of miR-17-92  indicating its dual role to both suppress transcription of the pro-proliferative miRNA cluster miR-17-92, and to mediate an M-CSF receptor differentiation pathway. Additionally, patterns of expression observed for miR-17-92 during monocytic differentiation resemble a previous analysis of miR-17-92 expression levels during lung development  supporting the general involvement of miR-17-92 amongst differentiation pathways.
TFAP2A (AP-2) and SP1 are two TFs predicted to regulate the miR-17-92 cluster and are notably up-regulated along with the cluster in the first 20 hours post-PMA stimulation. TFAP2A and SP1 are known to activate transcription of an enzyme involved in the sphingolipid metabolism consisting of several metabolites known to affect cellular proliferation . TFAP2A and SP1 transcribe sphingomyelin phosphodiesterase 1 (SMPD1) during monocytic differentiation in THP-1 cells after PMA stimuli . SMPD1 is required for the cleavage of sphingomyelin to phosphocholine and ceramide. As ceramide is a known inhibitor of proliferation , it seems reasonable that TFs of SMPD1 are up-regulated during differentiation. However, ceramide is also a substrate for several other enzymes whose products have not been implicated in proliferation, apoptosis or differentiation. Interestingly, miR-19a and miR-19b (part of the miR-17-92 cluster), are predicted to target sphingosine kinase 2 (SPHK2) mRNA in four independent databases (see Methods). SPHK2 is an enzyme that metabolizes downstream ceramide products. In the sphingolipid metabolism, SPHK2 has two functions. First, it catalyses the production of sphingosine 1-phosphate from sphingosine, which is produced from ceramides; and second, it catalyses the production of sphinganine 1-phosphate from sphinganine . Sphinganine and sphinganine 1-phosphate have been shown to inhibit and promote cell growth, respectively . Thus, we note that the predicted targeting and down-regulation of SPHK2 by miR-19a and miR-19b in the first 20 hours post-PMA stimulation could prevent the metabolism of two anti-proliferative metabolites simultaneously, thereby inhibiting proliferation. It is known that PMA stimulation can block proliferation of THP-1 cells up to 24 hrs . Thus, we propose an additional regulatory effect of TFAP2A and SP1 on the sphingolipid metabolism via the miRNA cluster miR-17-92. TFAP2A/SP1 mediated transcription of SMPD1 alone might not be enough to maintain an anti-proliferative ceramide signal, as ceramide is metabolized by other factors. On the other hand, TFAP2A/SP1 co-transcription of miRNAs targeting SPHK2 could provide an efficient and succinct means to retaining the ceramide signal.
We have computationally analysed the regulatory machinery that potentially affects transcription of miRNA genes during monocytic differentiation. Our methodology included the extraction of promoter regions for miRNA genes defined by trimethylated histones, computational prediction of TFBSs to establish TF→miRNA associations, and the use of time-course expression data for TFs and miRNAs measured during monocytic differentiation to assess reliability of the predicted TF→miRNA associations via time-lagged expression correlation analysis.
Several TFs (CEBPB, CREB1, ELK1, NFE2L2, RUNX1, and USF2), which are known to play a role in monocytic differentiation, have been identified. Our analysis suggests that their role in the differentiation process could be further expanded through consideration of the transcriptional regulation of miRNAs they affect. In addition, we propose several TFs (NFE2L1, E2F3, ATF2, HOXA4, SP3, and YY1) to have a central role in the regulation of miRNA transcription during the differentiation process. We have shown for several miRNAs (miR-21, miR-155, miR-424, and miR-17-92) how their predicted transcriptional regulation could impact the differentiation process.
The process of identifying a complete list of TF→miRNA associations is hampered by the correct definition of promoter/regulatory regions being an unresolved issue that has a great impact on all studies that deal with gene regulation. We utilised a recent set of promoters defined based on the observation that histones are generally trimethylated at lysine 4 residues at TSSs of genes. Due to the employed definition of promoters by Marson et al. we find that for several miRNAs we were not able to extract regulatory regions. Furthermore, we note that the here utilised promoter regions defined by Marson et al. range in length between 200 and ~4,700 bp with 60 percent of the utilised promoter regions being of length below 202 bp. Consequently, the promoter set defined by Marson et al. allows us to mostly analyze regulatory elements proximal to the TSS. Nevertheless, it has been well documented [70, 71] that proximal regulatory elements such as the TATA box play an important role in type II polymerase gene transcription. However, the utilised promoter set in this study represents one of the first sets of regulatory regions for miRNA genes.
It is important to note that the transcriptional circuitry described in our results is biased towards monocytic differentiation expression data, as several of TF→miRNA associations were discarded due to missing/incomplete expression data for either TF or miRNA. Furthermore, the expression based approach is limited in so far, as mature miRNAs are not the direct product of the TFs-mediated regulation but can undergo post-transcriptional regulation on pri- and pre-miRNA level . Thus, it is possible that miRNAs that are transcribed together as one primary transcript, show different expression profiles on the mature miRNA level. The three main reasons that constrained the set of TF→miRNA associations we determined in this study are as follows: 1/ An incomplete promoter set for miRNA genes. 2/ An incomplete/inaccurate motif set for the prediction of TFBSs. 3/ An incomplete expression set for TFs and miRNAs. Each of the reasons impacts on the accuracy of the predicted TF→miRNA associations.
Nevertheless, our analysis provides the first large-scale insights into the transcriptional circuitry of miRNA genes in monocytic differentiation. Taken together, our results suggest important regulatory functions of several TFs on the transcriptional regulation of miRNAs. The regulatory networks discussed here form only the starting point for an in-depth analysis of the regulatory mechanisms involved. The predicted TF→miRNA associations and their corresponding PCC s can provide the basis for a more detailed experimental analysis of miRNA regulation during monocytic differentiation.
We have computationally analysed the regulatory machinery that potentially controls the transcription of miRNA genes during monocytic differentiation. We made use of TFBS predictions in promoter regions of miRNA genes to associate TFs to miRNAs that they potentially regulate. With the help of time-course expression data for miRNAs and TFs during monocytic differentiation we evaluated each predicted association using a time-lagged expression correlation analysis. In this manner we derived a putative picture of the transcriptional circuitry that regulates miRNAs involved in human monocytic differentiation and determined potential key transcriptional regulators of miRNAs for this differentiation process.
miRNA time-course expression data
The miRNA expression profiles were obtained using Agilent's Human miRNA microarrays as described in . Three biological replicates have been measured before PMA stimulation and post-PMA stimulation at nine time points ranging from 1-96 hrs (1 hr, 2 hr, 4 hr, 6 hr, 12 hr, 24 hr, 48 hr, 72 hr, 96 hr). We required that two criteria were met for the inclusion of a miRNA expression time-series in the analysis:
i/ Expression of each miRNA should be denoted as "present" in at least one time point. Otherwise we assume that the expression series for the miRNA is insignificant.
ii/ For a miRNA, i/ must hold true in at least two of the three biological replicates.
The expression values of different biological replicates for a miRNA that satisfy the criteria have been averaged at each time point to generate one expression series per miRNA. Finally, each expression series was interpolated using piecewise cubic hermite interpolation [74, 75] with half an hour steps. In this manner, we obtained 193 (0-96 hrs) expression values for each individual miRNA expression series.
Identification of miRNAs showing differential gene expression
We calculate the log2fc by dividing each expression value of a miRNA by its expression value at zero hour (control) and taking the logarithm of base two (log2) of that ratio. A miRNA is considered to be influenced by the PMA stimulation in the differentiation process, if
i/ In at least one time point t its log2fc > 1 or log2fc < -1.
ii/ At any time point t where i/ holds true, the absolute difference dt in expression et at time point t and the expression e0 at zero hours must be greater than 0.1.
Transcription factor time-course expression data
The TF expression profiles were obtained using qRT-PCR as described in [6, 76]. Two biological replicates have been measured prior to PMA stimulation and in nine time points post-PMA simulation (1 hr, 2 hr, 4 hr, 6 hr, 12 hr, 24 hr, 48 hr, 72 hr, 96 hr). Primer design, RNA preparation, and cDNA synthesis have been performed analogously to . Normalization of the expression data of both replicates have been done as described in [6, 77]. All expression series for a TF that had available expression data within two biological replicates have been averaged over the respective biological replicates to produce one series of expression values per TF. Finally, each expression series was interpolated in half an hour steps using piecewise cubic hermite interpolation. Thus, we obtain 193 (0-96 hrs) expression values for each individual TF expression series.
Defining the promoter regions for miRNAs
We adopted the definition of miRNA promoters from . Each of the promoter regions had a score associated (as defined in ) that represents the confidence of dealing with a genuine regulatory region. We extracted all promoter regions with a score greater or equal to zero. The coordinates of the promoter regions were translated from the Human genome build 17 (hg17) to the Human genome build 18 (hg18)  using the UCSC liftover program .
TFBS analysis of miRNA promoter regions
TFBSs were mapped to the promoter region of the miRNAs with the MATCH™ program  using 522 mammalian matrices of TRANSFAC Professional Database (version 11.4) with their corresponding minimum false positive threshold profiles. Since TRANSFAC matrices are frequently associated with several TFs whose binding sites were used in building these matrices, we associated to each matrix all respective TFs (that have an Entrez Gene identifier associated). For example, we can associate several members of the JUN-FOS family (JUN, JUNB, JUND, FOS, FOSB, etc.) to matrix M00517. Binding sites of these TFs have been utilised to create this matrix. Thus, all of the TFs might be able to bind the TFBS predicted by the matrix.
Weighting associations using Pearson correlation
For each of the predicted TF→miRNA associations, scores (PCCs) were calculated as an indicator of how reliable the predicted association is, and as a measure of the strength of the association within the context of monocytic differentiation. The expression data for TFs and mature miRNAs during monocytic differentiation were utilised to calculate the best time-lagged expression correlation for a TF→miRNA association. The time-lagged expression correlation analysis calculates PCC between the TF expression and the time-shifted mature miRNA expression at different time-delays in order to take the influence of the TF on the miRNA transcription over time into account. We find the time-delay that maximizes the absolute value of PCC between the expression of the TF and that of the mature miRNA. The associations between pre-miRNA and the mature miRNA have been extracted using miRBase sequence database (version 10.1) [81, 13, 14].
For each predicted TF→miRNA association, where the miRNA does not share the same promoter with other miRNAs (i.e. not in a cluster), we calculate the PCC as follows:
i/ Identify the time-shift st. This is the time-shift where the absolute value of the PCC between the expression of the TF and the respective mature miRNA is maximal. We calculate the PCC for time-shifts ranging from 0.5 hour to six hours in intervals of half an hour.
ii/ The PCC for the association is calculated as PCC of the expression of TF and mature miRNA at the time-shift stfound in i/.
If a miRNA appears in a cluster with other miRNAs on the genome, then the predicted TF in the promoter of that cluster is associated to each of the respective miRNAs. Since the cluster is transcribed as one primary transcript we assume that a TF regulates each miRNA within the cluster with the same time-shift. Thus, we calculate one common time-shift st for the considered TF and all miRNAs within the cluster. The time-shift st is calculated as follows:
i/ The PCC of expression between the TF and each miRNA in the cluster is calculated for each considered time-shift (0.5 hour to six hours).
ii/ The average of all PCC s derived in i/ was calculated for each time-shift (0.5 hour to six hours). As a criterion for inclusion, the calculated PCC s for all associations should to have the same sign.
iii/ If ii/ could not be calculated at any time-shift (due to the sign rule), we did not assume that the TF X regulates any miRNA in that cluster and all X→miRNA associations of that cluster were discarded.
iv/ If not iii/, then the time-shift st is determined as the time-shift that maximizes the average calculated in ii/.
PCC of one TF→miRNA association where the miRNA is part of a cluster forms the PCC of expression of the TF and the respective mature miRNA at the determined time-shift st for the TF and the cluster. If a pre-miRNA is associated to more than one mature miRNA from its 5' and 3' arm, then the PCC is calculated independently for each mature miRNA and the maximum PCC is chosen.
Target predictions of miRNAs
The target gene predictions of human miRNAs have been gathered from four public available databases for miRNA target predictions (microRNA.org version 4 , TargetScan version 4.2 , miRBase version 5 [13, 14], and EIMMO2  with a cut-off value greater than 0.5). All target gene identifiers utilised in the respective databases have been converted to Entrez Gene identifiers using BioMart . If this was not possible the prediction has been discarded. We considered only predictions that are present in at least three out of the four databases.
transcription factor binding site
transcription start site
Pearson correlation coefficient
van Furth R, Cohn ZA, Hirsch JG, et al: The mononuclear phagocyte system: a new classification of macrophages, monocytes, and their precursor cells. Bull World Health Organ. 1972, 46: 845-52.
Traore K, Trush MA, George M, et al: Signal transduction of phorbol 12-myristate 13-acetate (PMA)-induced growth inhibition of human monocytic leukemia THP-1 cells is reactive oxygen dependent. Leuk Res. 2005, 29: 863-79. 10.1016/j.leukres.2004.12.011.
Martinez FO, Gordon S, Locati M, Mantovani A: Transcriptional profiling of the human monocyte-to-macrophage differentiation and polarization: new molecules and patterns of gene expression. J Immunol. 2006, 177: 7303-11.
Suzuki H, Forrest ARR, van Nimwegen E, et al: The transcriptional network that controls growth arrest and differentiation in a human myeloid leukemia cell line. Nat Genet. 2009, 41 (5): 553-62. 10.1038/ng.375.
Gavin IM, Glesne D, Zhao Y, Kubera C, Huberman E: Spermine acts as a negative regulator of macrophage differentiation in human myeloid leukemia cells. Cancer Res. 2004, 64: 7432-8. 10.1158/0008-5472.CAN-04-0051.
Mollinedo F, Gajate C, Tugores A, Flores I, Naranjo JR: Differences in expression of transcription factor AP-1 in human promyelocytic HL-60 cells during differentiation towards macrophages versus granulocytes. Biochem J. 1993, 294 (Pt 1): 137-44.
Rosa A, Ballarino M, Sorrentino A, et al: The interplay between the master transcription factor PU.1 and miR-424 regulates human monocyte/macrophage differentiation. Proc Natl Acad Sci USA. 2007, 104: 19849-54. 10.1073/pnas.0706963104.
Reddy VA, Iwama A, Iotzova G, et al: Granulocyte inducer C/EBPalpha inactivates the myeloid master regulator PU.1: possible role in lineage commitment decisions. Blood. 2002, 100: 483-90. 10.1182/blood.V100.2.483.
Chen A, Luo M, Yuan G, et al: Complementary analysis of microRNA and mRNA expression during phorbol 12-myristate 13-acetate (TPA)-induced differentiation of HL-60 cells. Biotechnol Lett. 2008, 30: 2045-52. 10.1007/s10529-008-9800-8.
Lu Y, Thomson JM, Wong HYF, Hammond SM, Hogan BL: Transgenic over-expression of the microRNA miR-17-92 cluster promotes proliferation and inhibits differentiation of lung epithelial progenitor cells. Developmental Biology. 2007, 310: 442-453. 10.1016/j.ydbio.2007.08.007.
Suzuki H, Okunishi R, Hashizume W, et al: Identification of region-specific transcription factor genes in the adult mouse brain by medium-scale real-time RT-PCR. FEBS Lett. 2004, 573: 214-8. 10.1016/j.febslet.2004.07.068.
Durinck S, Moreau Y, Kasprzyk A, et al: BioMart and Bioconductor: a powerful link between biological databases and microarray data analysis. Bioinformatics. 2005, 21: 3439-40. 10.1093/bioinformatics/bti525.
We thank Alistair Forrest and Mutsumi Kanamori for providing the normalized miRNA array data. SS, US and VBB were partly supported by the National Bioinformatics Network grants; ME has been supported by a Scarce Skills Scholarship from the National Research Foundation (SFH2007080600011); MK has been supported by the postdoctoral fellowship from the Claude Leon Foundation, South Africa; VBB has been partly supported by the National Research Foundation grants (61070 and 62302); CRM and VBB were partly supported by the DSTNRF Research Chair grant (64751). This study was also partly supported by a grant of the Genome Network Project from the Ministry of Education, Culture, Sports, Science and Technology, Japan (MEXT) to YH and a Research Grant for RIKEN Omics Science Center from MEXT to YH.
Authors and Affiliations
South African National Bioinformatics Institute, University of the Western Cape, Modderdam Road, Bellville, South Africa
Sebastian Schmeier, Cameron R MacPherson, Magbubah Essack, Mandeep Kaur, Ulf Schaefer & Vladimir B Bajic
RIKEN Omics Science Center, RIKEN Yokohama Institute, 1-7-22 Suehiro-cho, Tsurumi-ku Yokohama, Kanagawa, 230-0045, Japan
Harukazu Suzuki & Yoshihide Hayashizaki
Computational Bioscience Research Center (CBRC), King Abdullah University of Science and Technology (KAUST), Thuwal, 23955-6900, Kingdom of Saudi Arabia
Sebastian Schmeier, Mandeep Kaur, Ulf Schaefer & Vladimir B Bajic
SS and VBB conceptualized the study and wrote the paper; SS performed the analysis; CRM, ME, MK, and US performed parts of the analysis; HS and YH contributed experimental data. All authors have read and approved the final version of the manuscript.
Additional file 1: Promoter regions miRNAs. The first column contains the pre-miRNA identifier. The second column contains the chromosomal positions of the associated promoter regions separated by comma. Note that several pre-miRNAs can have the same promoter region associated, because they are transcribed together in form of a cluster. (XLS 23 KB)
Additional file 2: Interpolated expression data for 34 mature miRNA. The file consists of 195 columns. The first column contains the mature miRNA identifier. The second column contains the associated pre-miRNA identifier(s). Column 3-194 contain the interpolated expression values ranging in half an hour steps from 0 to 96 hours. (XLS 262 KB)
Additional file 3: Predicted TF→miRNA associations and their inferredPCCvalues. The file consists of three columns. The first column contains the TF. An identifier consists of the Gene Symbol separated by an underscore with the Entrez Gene id. The second column contains the pre-miRNA identifier of the miRNA that forms an association with the TF of the first column. The third column contains the inferred PCC for the association, which is based on the expression data of the TF and the mature miRNA associated to the pre-miRNA(s). In total the file contains 1,989 TF→miRNA associations. (XLS 319 KB)
Additional file 4: Interpolated expression data for 258 TFs. The file consists of interpolated expression data for 258 TFs that are present in the predicted TF→miRNA associations. Furthermore, the file consists of 194 columns. The first column is the TF identifier (Entrez Gene Id). Column 2-194 contain the interpolated expression values ranging in half an hour steps from 0 to 96 hours. (XLS 2 MB)
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
This article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License (
), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Schmeier, S., MacPherson, C.R., Essack, M. et al. Deciphering the transcriptional circuitry of microRNA genes expressed during human monocytic differentiation.
BMC Genomics10, 595 (2009). https://doi.org/10.1186/1471-2164-10-595