The transcriptional network that controls growth arrest and differentiation in a human myeloid leukemia cell line

Outline of the analysis strategy

In most cell line models, only a subset of cells undergoes growth arrest and differentiation. To maximize the sensitivity in this study, we identified a subclone of THP-1 cells in which the large majority of cells became adherent in response to PMA (Supplementary Fig. 1 online). Our strategy began with deepCAGE, which identified active TSSs at single-base-pair resolution, and simultaneously measured their time-dependent expression (using normalized tag frequency) as cells differentiated in response to PMA. The same RNA was subjected to cDNA microarray analysis on an Illumina platform. The differentiation of the cells was evident from the large increase in expression of macrophage-specific genes such as CD14 and CSF1R detected by both deepCAGE and microarray in all replicates (Supplementary Fig. 2 online).

Figure 1 summarizes our Motif Activity Response Analysis (MARA) strategy. Promoters were defined as local clusters of coexpressed TSSs and promoter regions as their immediate flanking sequences (Fig. 1a,b). To reconstruct transcription regulatory dynamics we refined earlier computational methods^6–8 by incorporating comparative genomic information and each TF’s positional preferences relative to the TSS in the prediction of regulatory sites. Binding sites for a comprehensive and unbiased collection of mammalian regulatory motifs were predicted in all proximal promoter regions (Fig. 1c) and the observed promoter expression profiles (Fig. 1d) were combined with the predicted site-counts (Fig. 1e) to infer time-dependent activity profiles of regulatory motifs (Fig. 1f). We inferred individual regulatory interactions (edges) between motifs and promoters by comparing the promoter expression and motif activity profiles (Fig. 1g). Rigorous Bayesian probabilistic methods were developed for all steps of the computational analysis. Finally, a core network was constructed by selecting the motifs that explained the greatest proportion of the expression variance, obtaining all predicted regulatory edges between TFs corresponding to these motifs and selecting those regulatory edges that had independent experimental support. Using this approach, we reconstructed the transcriptional regulatory dynamics associated with cellular differentiation in human THP-1 cells, and validated a subset of predicted regulatory interactions.

DeepCAGE quantification of dynamic TSS usage

CAGE tags generated from mRNA harvested at each time point were mapped to the human genome. Promoters were defined as clusters of nearby TSSs that showed identical expression profiles (within measurement noise) and were substantially expressed in at least one time point (Fig. 1a,b). Using these criteria we identified 29,857 promoters expressed in THP-1 cells containing 381,145 unique TSS positions (which is a subset of the nearly 2 million TSSs detected at least once in THP-1). These promoters were contained within 14,607 promoter regions (separated by at least 400 bp; Methods and Supplementary Fig. 3 online). The deepCAGE data was validated using genome tiling-array ChIP for markers of active transcription. Of the promoters identified, 79% and 78% were associated with H3K9Ac and RNA polymerase II, respectively (both markers of active transcription^13,14), compared to 18% and 27% for inactive promoters (Supplementary Note online).

Among the identified promoters 84% (24,984) were within 1 kb of the starts of known transcripts and 81% (24,327) could be associated with 9,452 Entrez genes. Approximately half of the remaining promoters were more than 1 kb away from the loci of known genes (Supplementary Fig. 4 online). These newly identified promoters are conserved across mammals, suggesting that they are true transcription starts of currently unknown transcripts (Supplementary Fig. 5 online). The association of 24,327 promoters with 9,452 Entrez genes extends previous evidence of alternative promoter usage¹¹—in this case even within a single cell type (Supplementary Table 1 online)—and demonstrates that promoter regions frequently contain multiple promoters with distinguishable expression profiles (Supplementary Table 2 online). In addition, for genes with known multiple promoters deepCAGE frequently identified only one promoter to be active in the THP-1 samples (Supplementary Fig. 6 online). Hence, deepCAGE samples a distinct aspect of transcriptional activity that can and does vary independently of mRNA abundances as measured by hybridization to representative microarray probes.

Promoter expression

Using the normalized tags per million (tpm) counts assigned to the promoters, we tested reproducibility among the three biological replicates and compared the outcome to the Illumina array from the same samples (Supplementary Fig. 7 online). DeepCAGE expression measurements were comparatively noisy (Supplementary Fig. 7a). Nevertheless, the median Pearson correlation between the replicate-averaged expression profiles of CAGE and microarray was around 0.72 (Supplementary Fig. 7b), which is comparable to that observed with other deep transcriptome sequencing datasets¹⁵. As predicted, the correlation is lower for genes with multiple promoter regions (Supplementary Fig. 7b and discussed further in Supplementary Note).

Comprehensive regulatory site prediction

Known binding sites from the JASPAR and TRANSFAC databases^16,17 were used to construct a set of 201 regulatory motifs (position-specific weight matrices, WMs), which represent the DNA binding specificities of 342 human TFs. We predicted transcription factor binding sites (TFBSs) for all motifs within the proximal promoter regions (−300 to + 100 bps) of all CAGE-defined promoters. Extending the proximal promoter regions beyond the −300 to +100 window decreased the quality of the fitted model described below (data not shown). In contrast to previous approaches that used simple WM scanning⁶, we incorporated information from orthologous sequences in six other mammals and used a Bayesian regulatory-site prediction algorithm that uses explicit models for the evolution of regulatory sites^18,19 (Fig. 1c and Methods). Notably, different motifs had distinct and highly specific positional preferences with respect to TSS (Supplementary Fig. 8 online), extending a previous genome-scale analysis²⁰. Positional preferences were incorporated in the TFBS prediction by assigning each site a probability that it is under selection and correctly positioned. This analysis generated approximately 245,000 predicted TFBSs for the 201 motifs genome-wide. For each promoter–motif combination, the TFBS prediction was summarized by a count N_pm, which represents the estimated total number of functional TFBSs for motif m in promoter p. The TFBS predictions were compared with published high-throughput protein–DNA interaction datasets (ChlP-chip) and predicted target genes were significantly (P values ranged from 0.02 for ETS1 to 6.60E–263 for GABPA) enriched among genes for which binding was observed (Supplementary Table 3 online).

Inferring key TFs and their time-dependent activities

The details of our Motif Activity Response Analysis (MARA) are described in Methods. Briefly, for each motif m and each time point t, there is an (unknown) motif activity A_mt, which represents the time-dependent nuclear activity of positive and negative regulatory factors that bind to the sites of the motif (for example, the E2F activity will depend on nuclear E2F1-8, and DP1-2 levels, as well as RB1 phosphorylation status). As in previous work^6–8,21, motif activities were inferred by assuming that the expression e_pt of promoter p at time t is a linear function of the activities A_mt of those motifs that have predicted sites in p. Additionally, the effect of motif m on the expression of promoter p is assumed to be proportional to the predicted number of functional sites N_pm. Assuming that the deviations of the predicted expression levels $e_{pt}^{\text{theo}}=\text{constant}+{\Sigma }_m{N}_{pm}{A}_{mt}$ from the observed levels e_pt are Gaussian distributed, and using a Gaussian prior on the activities, we determine fitted activities $A_{mt}^{\ast }$ that have maximal posterior probability (Methods).

The inferred motif activities were validated using a number of internal tests. First, our Bayesian procedure quantifies both the significance of each motif in explaining the observed expression variation as well as the reproducibility of its activity across replicates (Fig. 2 and Supplementary Table 4 online). The activity profiles of the top motifs are extremely reproducible across replicates and different measurement technologies (Figs. 2 and 3a and Supplementary Fig. 9 online). It should be stressed that, although motif activities are inferred by fitting the expression profiles of all promoters, the model cannot be expected to predict expression profiles of individual genes from the predicted TFBS in proximal promoters alone. The effects of chromatin structure, distal regulatory sites, nonlinear interactions between regulatory sites, and the contribution of the large numbers of human TFs for which no motif is known, are not considered. Furthermore, especially for genes that are dynamically regulated, mature mRNA abundance can be dynamically regulated independently of transcription initiation and promoter activity through selective mRNA elongation, processing and degradation. Our aim is not to predict expression profiles of individual genes but rather to predict the key regulators and their time-dependent activities, which can be inferred from integration of global expression information in a system undergoing dynamic change. We validated the significance of the inferred activity profiles by comparing the fraction of the ‘expression signal’ (expression variance minus replicate noise) that is explained by the model, compared to randomized versions, and under a tenfold cross-validation test (Supplementary Fig. 10 online). The explained expression signal is highly significant and this significance is maintained under tenfold cross-validation (Methods). In addition, the highly peaked positional profiles of TFBSs (Supplementary Fig. 8) suggest that knowing the exact TSS is important for accurate TFBS prediction. Indeed, the predicted TFBSs from CAGE promoters explain substantially more of the expression signal in microarrays than predicted TFBSs of the associated RefSeq promoters (Supplementary Fig. 10). We observe that the model better predicts the expression profiles of those promoters that are more strongly expressed, more reproducible across replicates, and have higher expression variance (Supplementary Fig. 11 online). Similarly, samples at the start and end of the differentiation time course are better predicted than those at intermediate time points (Supplementary Fig. 12 online), possibly because individual cells differentiate at different rates and leave the cell populations less homogeneous at intermediate time points.

Motif activities that were independently inferred from all 11,995 expressed microarray probes were combined with the inferred motif activities from all CAGE and microarray replicates into a final set of time-dependent motif activities (Methods). From these, we selected 30 ‘core’ motifs that contribute most to explaining the expression variation (red dots in Fig. 2) and segregated their activity profiles using a Bayesian procedure into nine clusters (Fig. 3b and Methods), including three clusters of upregulated motifs, three clusters of down-regulated motifs and three clusters containing single motifs with profiles involving different transient dynamics. The genome-wide set of target promoters for each of the motifs was determined as described in Methods. The significance of each regulatory ‘edge’ from a motif to a putative target promoter (containing a predicted TFBS) was quantified by the z value of the correlation between the motif’s activity profile and the promoter’s expression profile (Fig. 1e).

Core transcriptional regulatory network

The final aim in reconstructing transcriptional regulatory networks is to infer not only the key regulators and their target gene sets, but also the way in which the actions of these key regulators are coordinated. For this purpose, we collected all 199 predicted regulatory edges (z value ≥1.5) between the 30 core motifs. Recognizing that the prediction of individual regulatory edges is still prone to error, we constructed a core regulatory network (Fig. 4) of 55 highly trusted edges by filtering the predicted edges according to experimental validation, either within our data or in existing literature (Supplementary Table 5 online). In addition, for each core motif we extracted the set of predicted target genes (z value ≥1.5) and checked for enrichment of gene ontology terms. A selection of significantly enriched terms is shown as oval nodes in Figure 4 (full set of GO enrichments are available as Supplementary Table 6 online).

Whereas our method infers the key regulators ab initio, the majority of factors within this core network are known to be important in the monocyte-macrophage lineage, thereby validating the method. In addition the predicted targets of these motifs are enriched for biological processes known to be involved in differentiation of the monocytic lineage.

The gene ontology enrichments can broadly be divided into four groups. Downregulated motifs E2F1-5, NFYA,B,C and MYB are associated with cell cycle–related terms, consistent with the growth arrest observed during PMA-induced differentiation and the specific downregulation of numerous genes required for DNA synthesis and cell cycle progression within 24 h of PMA addition. Notably, MYB targets are also enriched specifically for microtubule-cytoskeleton–associated genes. Conversely, targets of upregulated motifs are associated with the terms immune response, cell adhesion, plasma membrane, vacuole and lysosome, all of which are consistent with differentiation into an adherent monocyte-like cell. The targeting of lysosomal genes by cholesterol-regulated SREBFs (sterol regulatory element-binding transcription factors) is of note, as lipid homeostasis is important in the macrophage in atherosclerosis and lysosomal storage diseases²². We also saw enrichment of signal transduction genes among targets of the early induced motifs EGR1-3 and TBP. Finally, there is a set of motifs whose targets are enriched in TFs. These motifs correspond to the transiently induced/repressed motifs, ATF5_CREB3, FOXO1,3,4 and SRF, and the repressed pair of OCT4 and FOXI1,J2 motifs.

Validation of edge predictions

THP-1 cells, even in an ‘undifferentiated’ state, are clearly a myeloid cell line. In seeking to validate the transcriptional network, we noted that there was a large set of TF genes expressed constitutively in the cells that were rapidly downregulated in response to PMA, of which MYB is an example, and another set that was expressed but further upregulated during differentiation. It is technically difficult to apply siRNA knockdown to genes that are only expressed later in the differentiation. To validate predicted edges empirically, we therefore chose to carry out siRNA knockdowns in undifferentiated THP-1 cells for genes encoding 28 TFs that are expressed in the undifferentiated state and for which we have associated motifs. To assess whether siRNA knockdown carried out in the undifferentiated state is appropriate to address factors that increase expression during the time course, we carried out the technically more difficult experiment of siRNA knockdown combined with PMA treatment for SPI1 (more commonly known in the literature as PU.1). All knockdowns were carried out in biological triplicate and qRT-PCR was used to confirm RNA-level knockdown, which in most cases was greater than 80% (Supplementary Table 7 online; in addition, protein-level knockdown was confirmed by protein blot for 14 siRNAs, see Supplementary Fig. 13 online). Changes in gene expression caused by TF knockdown were measured by Illumina microarrays. For each knocked-down TF gene, we obtained the list of predicted regulatory targets for the associated motif and divided the microarray probes into predicted targets and nontargets for a range of z-value thresholds. Higher-confidence targets in general show greater expression changes upon knockdown (Fig. 5a shows the example TF genes MYB, SNAI3, EGR1 and RUNX1; additional examples are shown in Supplementary Fig. 14 online). For SPI1, even in the absence of PMA treatment siRNA knockdown caused significant downregulation of predicted SPI1 targets, but the effects were much stronger when knockdown was combined with 1 h or 24 h of PMA treatment (Fig. 5b), confirming that PMA causes upregulation of SPI1 activity. A good correlation between target confidence (z-value cut-off) and average log expression ratio was observed for the large majority of experiments (Fig. 5c). For an intermediate cut-off of z = 1.5 we quantified the difference in log expression ratio of predicted targets and nontargets (Fig. 5d) and found significant changes (z-value larger than 2) for 23 of 33 cases with SPI1 knockdown combined with 24 h of PMA treatment and MYB knockdown being the most significant (Supplementary Fig. 15 online shows the entire distribution of log expression ratios of targets and nontargets for eight example TFs). Notably, for the TF genes LMO2, MXI1 and SP1, the knockdown led to a significant upregulation of their targets, suggesting that the three encoded TFs act primarily as repressors in undifferentiated THP-1 cells (Fig. 5d, also see Supplementary Fig. 14a). Together these results provide compelling experimental validation of our predicted regulatory edges.

Single TF knockdowns affect multiple motif activities

Besides validating predicted targets, the siRNA knockdowns can also be used to assess the effects of the knockdown of one TF gene on the motif activities of other TFs. In addition to the 28 TFs perturbed above, we included a further 24 TFs that lacked motifs but were naturally repressed during PMA differentiation, or had been reported to have a role in myeloid differentiation or leukemia (Supplementary Table 8 online).

The motif activity inference method was used to determine the changes in activities of all motifs upon knockdown of each TF gene. To assess the role of each TF in differentiation, we defined the differentiative overlap between a TF gene knockdown and the PMA time course as the fraction of all motifs that significantly changed their activity in the same direction upon TF gene knockdown as in the PMA differentiation (Methods). By far the largest differentiative overlap (69%) was observed for the MYB knockdown, which not only affected MYB motif activity, but also the activity of most motifs in the core network, with the most significant activity changes all in the same direction as in the PMA time course (Fig. 6a). Knockdown of 13 other TF genes generated an overlap greater than the negative control (Supplementary Table 9 online), and Figure 6 shows three further examples (E2F1, HOXA9 and CEBPG).

As for MYB, E2F1 knockdown reproduced some of the down-regulation of MYB and E2F activity observed upon PMA stimulation, but it failed to reproduce the upregulation of SREBF1,2, PU.1, NFATC1-3 and FOS,B,L1_JUNB,D activity (Fig. 6b). Similarly, the activity changes that HOXA9 knockdown induced were mostly in the same direction as in the PMA differentiation; however, the SNAI1-3 and IRF1,2 motif activities failed to be induced and the GATA4 and TBX4,5 motif activities failed to be downregulated (Fig. 6c). Notably, knockdown of CEBPG, encoding one of the PMA-downregulated factors, for which we do not have a motif, also generated activity changes that significantly overlapped those observed in response to PMA (Fig. 6d). Finally, instead of comparing the motif activity changes that different knockdowns induced, we can also directly compare the expression changes of all genes with the expression changes observed in the PMA time course. We found that MYB, HOXA9, CEBPG, GFI1, CEBPA, FLI1 and MLLT3 knockdowns all generated changes in gene expression that reiterated some of those observed with PMA treatment (Supplementary Table 8). MYB knockdown was exceptional, as it induced 35% (340/967) and repressed 19% (172/916) of the genes upregulated and downregulated with PMA, respectively. In addition the cells became adherent (Supplementary Fig. 16 online) and began to express the monocytic markers CD11B (ITGAM), CD54 (ICAM1), CD14, APOE and CSF1R (Supplementary Fig. 2), three of which we confirmed by flow cytometry (Supplementary Table 10 online). This development of adherence could be linked to the GO enrichment for cytoskeleton-associated genes among MYB targets noted above. Given these observations one might wonder whether MYB is a master regulator of the differentiation process and whether stronger and longer knockdown would have reproduced the complete differentiation observed under PMA treatment. Several observations argue strongly against this. First, the gene sets perturbed by MYB and by the other pro-differentiative TFs overlap only partially (Supplementary Table 11 online). Second, of the six other pro-differentiative TF genes only two (CEBPG and GFI1) are affected by MYB knockdown. Both these facts indicate that the other pro-differentiative TF genes are not simply downstream of MYB. Third, MYB downregulation does not occur until after the second hour of the PMA time course (Fig. 3b), which is at odds with the idea of MYB sitting at the top of the regulatory hierarchy. It is also worth noting that THP-1 cells harbor a leukemogenic fusion²³ between MLL (mixed-lineage leukemia) and MLLT3 (MLL translocation partner 3) and that the MLLT3 siRNA targets this leukemogenic fusion (note that full-length MLLT3 does not seem to be expressed in THP-1 as there is no CAGE 5′ signal for this gene). Our data indicate that this fusion interferes with differentiation and that neither PMA treatment nor MYB knockdown affects MLL-MLLT3 levels, suggesting these stimuli can bypass the differentiative block. Conversely, MLLT3 knockdown had no effect on MYB levels. These results agree with previous RNAi studies that conclude that down-regulation of MLL leukemogenic fusion proteins can promote growth arrest but is not required for terminal differentiation^24,25. Thus, individual TF gene knockdowns affect the activities of multiple motifs and elicit different, but overlapping, subsets of the regulatory changes observed in the PMA time course. Taken together, the data indicate that the independent perturbation of expression of multiple TFs in response to PMA is both necessary and sufficient to initiate partial differentiation.

Many TFs are involved in the differentiation process

The network predictions and the siRNA results above suggest that upregulation and downregulation of the activities of multiple co-operating TFs is required for differentiation. Of a curated list²⁶ of 1,322 human TFs, 610 were detected by both CAGE and microarray in at least one time point (Supplementary Table 12 online); however, only 155 of these are covered by weight matrices, suggesting that other factors may well be important in these cells. Of the 610 expressed TFs 64 were most highly expressed in the undifferentiated and 34 in the differentiated state. In addition, 101 TFs were transiently induced or repressed during differentiation. To elucidate the connection of these TFs to the inferred network, we compared the predicted regulatory inputs of co-regulated subsets of TFs with the predicted regulatory inputs of the set of all 610 expressed TFs.

Whereas no motifs are overrepresented among inputs of statically expressed TFs, inputs of dynamically expressed TFs showed enrichment for a subset of motifs. TFs downregulated from 0 to 96 h PMA were most enriched for three downregulated motifs of the core network: OCT4 (3.4×), GATA4 (3.3×) and NFYA,B,C (2.2×) (Supplementary Table 13a online). Similarly, TFs upregulated from 0 to 96 h were most enriched for core network motifs that increase activity during differentiation: SNAI1-3 (4.6×) and TBP (5.2×) (Supplementary Table 13b). Finally, transiently regulated TFs were enriched for the SRF (3.5×) and NHLH1,2 (3×) motifs (Supplementary Table 13c).

Notably, TFs that are predicted targets of SRF are mostly induced in the first hour of PMA-induced differentiation. During this first hour 55 of the 57 genes whose expression was perturbed are induced and 30% encode TFs (Supplementary Fig. 17a online). The regulatory inputs of these early-induced TFs are enriched for the motifs SRF, TBP and FOSL2 (Supplementary Table 13d), which all correspond to known PMA-responsive TFs^27–30. Among the early-induced TFs, five correspond to upregulated core network motifs themselves (FOSB, EGR1-3 and SNAI1) and two (MAFB and EGR1) are known to induce pro-differentiative changes^31,32. It is also worth noting that significant downregulation did not occur until the second hour, and this may require both early induction of transcriptional repressors and the RNA degradation proteins BTG2 and ZFP36 (tristetraprolin)^33,34 (Supplementary Fig. 17b). Together, these results suggest that induction of SRF target genes in the first hour is critical to establishing the differentiative program and is required before factors maintaining the undifferentiated state are downregulated (Supplementary Fig. 17b,c).

Web interface to data and analysis results

To facilitate the use of the data and analysis of results amassed here, we provide an online tool, EdgeExpressDB, as part of the FANTOM4 web resource, which allows users to explore our annotations of the structure, expression and regulation of promoters genome-wide. It also integrates published TF–promoter interactions, the siRNA perturbations and genome-wide chromatin immunoprecipitation experiments. Our complete set of regulatory-interaction predictions provides a large collection of hypotheses that can be targeted for validation, for example, through chromatin immunoprecipitation, gel shift assays or reporter assays. The value of this resource is illustrated by detailed examination of individual loci. For example, the osteo-pontin gene (SPP1) is massively induced from 12 h of differentiation (Supplementary Note). Our predictions confirm RUNX and PU.1 as regulators and support a previous analysis in mouse implicating the TGIF1 factor. In addition our analysis identifies NFAT, STAT, NKX6.2 and LIM domain and homeobox proteins as candidates for further testing.

Finally, our set of human promoters, TF motifs, genome-wide annotation of TF-binding sites and their predicted effects on the expression of the target promoters are available through the Swiss-Regulon website. A web interface, allowing researchers to automatically perform Motif Activity Response Analysis (MARA) of their own expression data in terms of our genome-wide predictions of TFBSs, is also available at SwissRegulon.

URLs.

FANTOM4 web resource, http://fantom.gsc.riken.jp/4/; SwissRegulon, http://www.swissregulon.unibas.ch. All methods are described in Supplementary Methods.

Accession codes.

Accession numbers of the data sets, deposited in public databases, are shown in Supplementary Table 15 online.