Deciphering the transcriptional circuitry of microRNA genes expressed during human monocytic differentiation

Background 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. Results 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. Conclusions 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.


Background
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 [1]. The human monocytic leukemic cell line, THP-1 [2], 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 [3]. 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 p21 WAF1/ CIP1 , enhancing binding of the SP1 factor to the p21 WAF1/ CIP1 promoter. PMA inhibition of cell growth is mediated by several signalling pathways such as MAPK and ROSdependent Raf/MEK/ERK pathway [4]. Human monocytic maturation incorporates lipid and protein metabolic processes together with several G-protein coupled receptors (GPCRs) [5].
Canonical miRNA biogenesis begins with the transcription of the pri-miRNA by RNA polymerase II [17][18][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 [22]. Dicer, a RNase III endonuclease, cleaves 22nucleotide 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 miR-NAs have their own transcriptional units [8], it is known that several miRNAs are transcribed together as a single pri-miRNA [24][25][26]. These clustered miRNAs are thus coregulated. On the other hand, miRNAs can also be transcribed together with a protein-coding host gene [8]. 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 [30] 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][34][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 largescale 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. After-wards, 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 log 2 fc (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 log 2 fc > 1 or log 2 fc < -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 log 2 fc > 1 or log 2 fc < -1 was satisfied for a miRNA, we calculated the difference d t of the expression e t at time point t and its expression e 0 at the zero time point. We sub-selected those miRNAs for which abs(d t ) > 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 [36]. Thus, we consider the miRNAs that form such clusters to be generally regulated in the same manner.
Marson et al. [30] 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 Overview of the analysis Figure 1 Overview of the analysis. The figure illustrates the analysis steps (blue/green boxes). In addition, the figure shows the data (red boxes) that has been utilised within individual analysis steps. To map TFBSs to the 38 promoters we utilised TRANSFAC Professional database (version 11.4) [37,38]. TRANS-FAC'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 miR-NAs.

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 miR-NAs 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][40][41]. Unfortunately, the correct time-shifts are undetermined. In our analyses, we incorporated timeshifts 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 PCCs 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 [33] or compute the best time-shift depending on known examples of regulation [31]. Up to now, too few experimentally verified examples of TFs that regulate miRNAs are known, thus a model to introduce the "correct" timeshift 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 PCCs 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 PCCs. 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 [42], CREB1 [43], ELK1 [44], NFE2L2 [45], RUNX1 [42], and USF2 [46]) 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, [47]). 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 [48] and all pathway analyses were based on KEGG [49] using DAVID [50].

miR-21
Fugita et al. demonstrated that mir-21 is expressed during PMA-induced differentiation in the human promyelocytic leukaemia cell line, HL-60 [51]. 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 [51]. Our TF→miRNA associations and their inferred Pearson correlation coefficients. A/ Depicted is the number of TF→miRNA associations that have a score equal or greater then specific PCCs. The blue blocks indicate the number of associations that have a positive PCC greater or equal to the positive value indicated on the x-axis. The red blocks indicate the number of associations with a negative PCC smaller or equal to the negative value indicated on the x-axis. As expected, the number of associations steadily decreases with increasing absolute PCC. B/ Depicted is the distribution of the absolute value of the calculated PCCs for all 1,989 TF→miRNA associations. The red line indicates the cut-off value that was utilised to select the top quartile of the associations. The distribution is not normal distributed, but skewed towards higher PCCs resulting from the chosen method of time-shifts, which favours higher PCCs over lower ones.
known to be activated by PMA induction which is supported by our findings (data not shown) [52]. Fugita et al. also demonstrated that AP-1 and SPI1 synergistically mediate the transcriptional process [51]. 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 [53] and the tropomyosin 1 (TPM1), PDCD4, and maspin gene in the human breast cancer cell line, MDA-MB-231 [54]. 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 PCCs (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).

miR-424
Rosa et al. reported that mir-424 is expressed during PMAinduced 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 coactivator 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) [55]. 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 miR-NAs 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 PCCs 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. [55]. 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 [55], 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).

miR-155
Chen et al. reported that mir-155 is expressed during PMA-induced differentiation in the human promyelocytic leukaemia cell line, HL-60 [57]. 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 Overview of 12 TFs and their regulatory effect on miRNA  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).

Involvement of miR-21 in monocytic differentiation
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.

miR-17-92
Members of the miRNA cluster mir-17-92 are known to be down regulated in the HL-60 cell line after PMA stimula-tion [57]. 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 [63]. Expression of miR-17-92 is high in proliferating cells and is positively regulated, in part, by MYC (c-Myc) [64]. E2F1, an activator of MYC, is itself a target of miR-17 and miR-20a [61] 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 [36]. Additionally, E2F3 is shown to be a predominant isoform that regulates miR-17-92 transcription [36]. We show that after ranking PCCs 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).
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  [36]. 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 [63]. On inspection of the log 2 fc 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 [66], 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 PMAtreated 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 [66] 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 [67] 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 [68]. TFAP2A and SP1 transcribe sphingomyelin phosphodiesterase 1 (SMPD1) during monocytic differentiation in THP-1 cells after PMA stimuli [68]. SMPD1 is required for the cleavage of sphingomyelin to phosphocholine and ceramide. As ceramide is a known inhibitor of proliferation [69], 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 1phosphate from sphingosine, which is produced from ceramides; and second, it catalyses the production of sphinganine 1-phosphate from sphinganine [69]. Sphinganine and sphinganine 1-phosphate have been shown to inhibit and promote cell growth, respectively [69]. 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 [4]. 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.

Summary
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 dif-ferentiation, 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 trimethyl-   [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.

Involvement of miR-155 in monocytic differentiation
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 [72]. 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/inaccu-

Conclusions
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 [73]. 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 timeseries 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 log 2 fc by dividing each expression value of a miRNA by its expression value at zero hour (control) and taking the logarithm of base two (log 2 ) 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 log 2 fc > 1 or log 2 fc < -1.
ii/ At any time point t where i/ holds true, the absolute difference d t in expression e t at time point t and the expression e 0 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 [76]. 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 [30]. Each of the promoter regions had a score associated (as defined in [30]) 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) [78] using the UCSC liftover program [79].

TFBS analysis of miRNA promoter regions
TFBSs were mapped to the promoter region of the miR-NAs with the MATCH™ program [80] 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 s t . 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 s t found 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 s t for the considered TF and all miRNAs within the cluster. The time-shift s t 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 PCCs derived in i/ was calculated for each time-shift (0.5 hour to six hours). As a criterion for inclusion, the calculated PCCs 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 s t 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 s t 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 [15], TargetScan version 4.2 [9], miRBase version 5 [13,14], and EIMMO2 [16] 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 [82]. 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.