THP-1 cells (high differentiation clone 5, from the FANTOM4 consortium, passage number 8 (P8), provided by Dr. Mark Barnett, The Roslin Institute, United Kingdom) were cultured in sterile-filtered (0.2 μm) RPMI1640 medium with HEPES modification (25 mM HEPES + NaHCO₃) (R5886 Sigma-Aldrich, Gillingham, United Kingdom), 10% heat inactivated fetal bovine serum (FBS, GE Healthcare, PAA laboratories, Pasching, Austria), 1x Glutamax (35050-038, Gibco), 1 mM sodium pyruvate (11360-039, Gibco), 1x MEM non-essential amino acids solution (11140-035, Gibco), 50 μM 2-mercaptoethanol (31350-038, Gibco) and penicillin-streptomycin (16 U/ml and 16 μg/ml, Gibco) (THP-1 medium). Cells were incubated at 37°C, 5% CO₂ and their concentration was maintained between 2 × 10⁵ – 8 × 10⁵ cells per ml by splitting in half three times a week. For the THP-1 differentiation time course, THP-1 cells were grown up at a low passage number (typically 2 or 3 passages from the P8 stock). The day before the start of the differentiation, cells were counted by hemocytometer and for each time point between 3 × 10⁶ and 5 × 10⁶ cells were pelleted and resuspended in 10 ml fresh medium. THP-1 cells were then differentiated by adding 30 ng/ml (48.6 nM) phorbol 12-myristate 13-acetate (PMA; P1585, Sigma-Aldrich) in DMSO. The cells were plated on a Sterilin bacteriological plate and at appropriate time points, the cells were lifted off by flushing with a blunt-end needle syringe, resuspended in 1 ml of RNA Bee (AMS Biosciences, Frienswood, TX, United States) and frozen at −80°C for later extraction. Medium was not changed for the duration of the experiment to avoid initiation of a serum-specific response during the time course. For each experiment, one culture was used to set up two batches of cells. The first batch was treated with PMA on the morning of the first day. Time points between 0 and 12 h were taken from this batch. The second batch from the same culture was treated with PMA in the evening of the first day (which is the 0 h time point from this batch, taken at the same time as the 12 h time point of the first batch). All cells were then incubated overnight. On the morning of the second day, the 24 h time point of the first batch and the 12 h time point of the second batch were taken, followed by time points between 12 and 24 h taken from the second batch. Batches had several overlapping time points (0, 12, 24, and 36 h) to ensure that cells were at similar state of differentiation. RNA was extracted from the RNA Bee lysate according to manufacturer’s instructions. Extracted RNA was treated with DNase I (Ambion DNase kit AM1906, Thermo Fischer Scientific) and the quantity and quality assessed using a NanoDrop spectrophotometer ND-1000 (Nanodrop technologies, Wilmington, DE, United States) and Agilent RNA ScreenTape System (Agilent Technologies, Santa Clara, CA, United States) according to manufacturers’ instructions. Only RNA samples with RIN^e values (Imbeaud et al., 2005) around or above 8 (out of 10) were used for CAGE sequencing.

Time points chosen for CAGE sequencing were hourly from 0 to 12 h, two hourly from 12 to 18 h, and then at 24, 36, 48, and 96 h. At least three replicates of each time point were included. Eight samples of 5 μg RNA each were made into one CAGE library following the protocol adapted from Takahashi et al. (2012). All the primer sequences were taken from Takahashi et al. (2012) and primers were synthesized by IDT (Coralville, IA, United States). A detailed description of the CAGE protocol and associated bioinformatic analysis is provided in Supplementary File S1. Results for 66 samples were processed using the R/Bioconductor package CAGEr (Haberle et al., 2015). All 66 samples were normalized to a power law distribution (Supplementary File S1), with expression levels presented as tags per million (TPM). The full quality assurance process for CAGE sequencing is outlined in Supplementary File S1. TSS were mapped to the human genome, version GRCh38.p9. The normalized and aggregated groups of TSS (clusters of transcription start sites; CTSS) were then annotated with gene names as described in Supplementary File S1.

Patterns of gene expression during differentiation were visualized using the network analysis tool Graphia Professional version 2.0¹. Mean expression values across replicates for CTSS with at least 10 TPM in at least one sample were uploaded into the Graphia Professional software. Sample-to-sample analysis was performed using a correlation coefficient threshold of r ≥ 0.95, as indicated in the Results. Two samples, both at the 4 h time point, were found to be outliers, probably due to low level of mapped reads (2.7 and 6.2 × 10⁶ reads). These were removed from subsequent analysis. Gene-to-gene analysis was performed using a correlation coefficient of r ≥ 0.75. Clusters of genes with similar expression patterns were determined with the Markov Clustering Algorithm (MCL) at an inflation value (determining granularity of the clusters) of 1.7. A list of genes associated with cell proliferation was obtained from Giotti et al. (2018) and a list of human transcription factors from Lambert et al. (2018). Further network analysis was performed on these subsets of genes, with correlation coefficients and inflation values as described in section “Results”. DAVID (Database for Annotation, Visualization and Integrated Discovery²) was used to identify functionally enriched biological themes based on gene ontology (GO) annotations, using the Functional Annotation Clustering tool (Huang et al., 2007). High annotation scores indicate significant enrichment for groups of related GO terms.

Activities of transcription factor binding motifs were estimated using MARA (Suzuki et al., 2009) using all CTSS with at least 10 TPM in at least one sample, excluding the mitochondrial chromosome and considering conserved motifs in the region −300 to +100 bp of the midpoint of the CTSS. Motif names are based on the set of consolidated position weight matrices for human transcription factors curated for the earlier project (Suzuki et al., 2009), available at http://swissregulon.unibas.ch/sr/fantom4.

Candidate enhancers were identified as bidirectional TSS as described (Andersson et al., 2014) at two levels of stringency. The more stringent set (∼1000 loci) required at least two libraries with more than one CAGE tag for each CTSS, and for the enhancers an absolute directionality score of less than 0.8, and bidirectional tags in more than two samples. For the less stringent set (∼2500 loci), the criteria were >1 tag per CTSS in at least one sample, and bidirectional tags for each enhancer in at least one sample. All bidirectional TSS regions closer than 500 bp to an annotated TSS and closer than 100 bp to an exon were removed to avoid detection of promoters of protein-coding genes. The analysis was done with use of the CAGE analysis package CAGEfightR (Thodberg et al., 2019), and the ported enhancer detection scripts from Andersson et al. (2014)³.

Gene expression levels detected by CAGE were validated for representative genes using quantitative reverse transcriptase polymerase chain reaction (qRT-PCR) on the same RNA samples. cDNA was synthesized using Superscript III Reverse Transcriptase (Invitrogen) according to manufacturer’s instructions, using 500 ng of total RNA and 100 ng of random primers. The ribonuclease inhibitor used was RNasin Plus (Promega, Madison, WI, United States; 20 units per reaction) and extension at 50°C was for 60 min. The resulting cDNA was diluted 1:1 with water. Standard curves were established using four concentrations of starting RNA: 12.5, 6.25, 3.13, and 1.56 ng/μl. qRT-PCR was carried out using SYBRGreen (Roche, Mannheim, Germany) in a Roche Light Cycler 480 according to manufacturer’s instructions. Samples were pre-incubated at 95°C for 5 min (ramp rate 4.40°C/s), then amplification steps were repeated for 45 cycles. Amplification steps were 95°C for 10s (ramp rate 4.40°C/s), 60°C for 15 s (ramp rate 2.20°C/s) and 70°C for 30 s (ramp rate 4.40°C/s). Afterward, the melting curve for primers was measured by incubating at 95°C for 5s (ramp rate 4.40°C/s), 65°C for 1 min (ramp rate 2.20°C/s) and then the temperature was increased to 97°C by 0.11 C°/s. At the end, the plate was cooled for 30 s at 40°C. Calculations were performed using the Roche Light Cycler 480 software, with the Advanced Quantification setting. ΔCt was calculated with previously established values of primer efficiencies from the standard curves (calculated using the same software, by using Abs Quant/2nd Derivative Max setting). Primer sequences and efficiencies are available in Supplementary File S1.

Propidium iodide staining was used to examine the cell cycle. 1 × 10⁶ cells were pelleted at 400 g for 5 min and resuspended in 1 ml of ice-cold 70% ethanol and then left at 4°C for at least 24 h (up to 2 months). On the day of the flow cytometry, cells were spun down at 16,000 g for 5 min and then washed twice with PBS. 300 μl of the Cell Cycle Staining Solution [38 mM sodium citrate, 10 μg/ml pancreatic RNase A (Signa-Aldrich), 68 μM propidium iodide (Sigma Aldrich)] was added to the pellet. Samples were incubated in the dark for 1–2 h and 300 μl of PBS was added right before flow cytometry. Samples were analyzed on a Fortessa cytometer (BD Biosciences, Franklin Lakes, NJ, United States) using a 610 nm laser. Results were analyzed using FlowJo software version 10.0.8r1 (Flow Jo LLC, Ashland, OR, United States) by gating single viable cells by excluding doublets. The final window showed the histogram of the area of propidium iodide channel, which was examined using the cell cycle option of FlowJo software by choosing the Watson (pragmatic) type of analysis (Watson et al., 1987).

PMA treatment of THP-1 high differentiation clone 5 cells was carried out multiple times to harvest RNA for CAGE sequencing. To ensure that differentiation from monocyte to macrophage phenotype had occurred, changes in morphology and attachment to the culture vessel surface were examined for each replicate. Undifferentiated THP-1 monocyte-like cells grew in suspension, with most cells not attached to the bottom surface of the vessel. The THP-1 cells started to adhere around 6–8 h after addition of 30 ng/ml of PMA. They then began to develop more macrophage-like morphology, with the most dramatic difference observed at 48–96 h. Their shape changed from being round to more variable and the cytoplasm contained more phagocytic vacuoles. At each time point there were a few cells where the morphology indicated they were more macrophage-like than the majority of cells. More than 99% of the cells were adherent and had similar morphology by the end of the time course at 96 h, indicating that this culture of THP-1 clone 5 had retained its high differentiation phenotype.

Prior to addition of PMA, the majority of the THP-1 cells were in G1 (growth) phase with approximately one third in S (DNA synthesis) phase (Figure 1A). After 10 h, the profile of DNA content detected by FACS was clearly changed. The number of cells in S phase was dramatically reduced and by 24 h there were <10% of cells with intermediate DNA content between 2N (G1 phase) and 4N (G2 phase) (Figures 1B,C). This low proportion was maintained for the rest of the time course (Figures 1D,E). The proportion of cells with 4N DNA content (G2) actually increased by 24 h as cells that had already entered S phase completed DNA synthesis. The proportion of cells in G2 phase gradually declined thereafter. These results indicate that differentiation of THP-1 cells to a macrophage morphology was preceded by a rapid decline in proliferation due to blocks at both entry to S phase and the well-documented G2-M phase (cytokinesis/cell division) check point, most commonly associated with DNA damage (Stark and Taylor, 2004).

To ensure that the differentiation of THP-1 cells from pro-monocytic to macrophage-like phenotype was supported by changes in gene transcription, gene expression levels derived from qPCR were assessed for two relevant markers, MYB which encodes a transcription factor associated with proliferating hematopoietic cells (reviewed by Zhou and Ness, 2011) and monocyte differentiation antigen gene CD14 (Bazil et al., 1982; Wright et al., 1990). MYB was down-regulated during THP-1 differentiation in the earlier FANTOM4 study (Suzuki et al., 2009). In the present study, MYB expression was also found to decrease after PMA stimulation in all biological replicates, with the lowest level at around 2:30–3 h (Supplementary Figure S1). CD14, which was not detected in undifferentiated THP-1 clone 5 cells, consistent with FACS analysis of the parent THP-1 line (Bosshart and Heinzelmann, 2016), was not up-regulated until the last few points of the time course, from 24 h (Supplementary Figure S1).

The time course of THP-1 differentiation for all of the replicates is available on the ZENBU genome browser⁴. The browser enables gene-specific search and includes the FANTOM5 promoter and enhancer tracks and the UCSC CpG island track. To compare the detailed time course with the earlier study (Suzuki et al., 2009) that had fewer time points, CAGE sequencing results were first analyzed for index genes highlighted as markers in that study. Figure 2A shows the time courses of representative macrophage differentiation marker genes. As with the qPCR results, CD14 increased late in the time course. Consistent with the previous data, CD14 was most strongly induced between 48 and 96 h. ICAM1 was induced maximally by around 24 h, ITGAM (encoding CD11B) increased after 12 h and APOE increased from around 18 h. By contrast to the earlier study, but in keeping with other results (Datta et al., 1992), the transcript encoding the growth factor receptor CSF1R was highly expressed by THP-1 cells without differentiation (200 TPM) but slowly induced a further 7-fold by PMA, reaching a maximum around 96 h. CSF1R is set as the default view on the ZENBU browser. The time course of CSF1R induction resembled the induction of CD36, proposed as a marker of differentiation in this system (Maess et al., 2010). CD36 mRNA was induced by 1 h and continued to increase throughout the time course. THP-1 cells respond to CSF1R ligand, CSF1 (Datta et al., 1992) and were previously shown to express CSF1 mRNA (Rambaldi et al., 1988). In the present study, CSF1 mRNA was induced with a very similar time course to CSF1R suggesting the likely involvement of an autocrine stimulatory loop in differentiation. CSF1 mRNA is also induced during the differentiation of primary human monocytes (Horiguchi et al., 1986) and indeed the highest level of expression in the FANTOM5 dataset is in monocyte-derived macrophages.

The transcriptional regulation of SPP1 (encoding osteopontin), representative of a cluster of genes induced between 12 and 24 h and associated with adherence, was dissected in detail previously (Suzuki et al., 2009). SPP1 was also massively induced by PMA in each time course replicate generated for this study. The timing of induction varied between replicates, so that the expression was highly variable between 6 and 24 h. Figure 2A shows mean data, averaged across replicates. The variation suggested that osteopontin is induced during the process of adhesion and that cluster analysis of the replicates independently might identify other genes that are involved in this key event.

A key feature of the PMA time course also dissected previously was the rapid down-regulation of genes for transcription factors involved in cell cycle regulation, notably MYB and E2F1, associated with growth arrest and induction of differentiation. Knockdown of some of these factors was sufficient to induce growth arrest and partial differentiation (Suzuki et al., 2009). Figure 2B shows a subset of these transcription factors. CAGE sequencing showed that MYB was expressed from two promoters, both of which were down-regulated rapidly in response to PMA. The related transcription factor gene, MYBL2 was highly expressed at the start of the time course (16 fold higher than the combined expression of the two MYB promoters) and declined later than MYB. MYBL2 was not assessed in the previous analysis. E2F1 was maximally repressed by around 12 h. GFI1 has been identified as a transcriptional repressor of target genes of the macrophage transcription factor, PU.1 (encoded by SPI1) (Wilson et al., 2010; Barth et al., 2019). In THP-1 cells, GFI1 was expressed from 4 separate promoters and rapidly down-regulated by PMA (Figure 2B). Consistent with repression of entry into the cell division cycle, the G1/S cyclin gene CCNE2 was rapidly repressed to undetectable levels (Figure 2C). Cell cycle-associated transcripts for genes associated with S phase and mitosis (Giotti et al., 2018) such as PCNA, BUB1, AURKA, CDK4 and the M phase cyclin, CCNB1, declined more slowly and progressively (Figure 2C), presumably as individual cells competed the cycle (Figure 1). Many of these transcripts showed some evidence of re-expression at 48–96 h (for example E2F1, CCNE2 and PCNA in Figure 2), which may reflect expansion of the small subset of cells that is resistant to PMA-induced differentiation even in this cloned line.

The time course of the loss of positive cell cycle regulators in response to PMA suggests that this may be a consequence rather than a cause of the block to S phase entry (Figure 1). The FANTOM 4 project did not analyze the expression of inducible cell cycle repressors. Others reported that the cell cycle inhibitor CDKN1A (p21WAF) was induced by PMA at the mRNA and protein level in THP-1 cells (Akashi et al., 1999; Traore et al., 2005). The induction of CDKN1A was associated with the loss of phosphorylated CDK2 protein, which was linked in turn to selective inhibition of phosphorylation of RB (retinoblastoma protein). Figure 2C shows that CDKN1A was induced from almost undetectable levels within 1 h and thereafter sustained at high levels throughout the time course.

These initial studies indicate that the current dataset broadly reproduces the temporal cascade of cell cycle arrest and monocyte-macrophage differentiation seen previously (Suzuki et al., 2009) and show that the increased temporal resolution of the present study was informative. We therefore undertook a refined network analysis of the data.

CAGE sequencing data values (means for each replicate; N = 2–6) were loaded into the network analysis tool Graphia Professional to find patterns of gene activation and repression across the time course. To look for general trends across all replicates, a sample-to-sample network at 18 different time points was generated using a correlation coefficient of ≥0.95. The sample-to-sample network (Figure 3A) highlights the global transcriptomic transition as the time after treatment increased. The transcriptome of samples at 0 h (deep red) were most different from those at 96 h (deep purple). Samples from other time points lie between these two, moving from red (1 h) through orange and yellow (2–6 h), green (7–12 h), and blue (14–48 h).

Next, CTSSs were correlated with each other based on the mean expression at 18 time points (to create a gene-to-gene network). The analysis was done with a correlation coefficient threshold of 0.85, which included 10,168 nodes out of 11,360 possible CTSS (see section “Materials and Methods”). Expression patterns of excluded CTSSs were not correlated with any other CTSS at a correlation coefficient of 0.85. A network graph comprising one large element of 10,046 nodes and 53 smaller unconnected elements of between 2 and 11 nodes was generated using Graphia Professional. Markov clustering with an inflation value of 1.7 partitioned the graph into clusters of CTSS with related expression profiles. Figure 3B shows the result for the main element. There were 92 clusters with 10 or more nodes. A number of clusters showed distinct patterns of gene expression during differentiation. The largest of these are shown in Figure 3C. The clusters in which the key genes discussed above were found are shown in Table 1. A full list of all the clusters and the mean expression patterns of those with 10 or more nodes is presented in Supplementary Table S1.

The largest cluster, Cluster001, contained genes which code for proteins associated with mature macrophage functions, including surface markers (CD14, CD163, ITGAM), lysosome activity, transmembrane location and integrin complex. The most significantly enriched GO term was lysosome (enrichment score 14.8). This cluster contained genes with expression that steadily increased throughout differentiation, peaking at the end of the time course. Cluster002 declined progressively across the time course and showed high enrichment scores for terms associated with cell cycle (enrichment score 31.7) and mitochondria (enrichment score 21.1). Three hundred and sixty seven (17%) of the genes in this cluster were contained in an annotated list of 701 cell cycle genes (Giotti et al., 2018). These two clusters are consistent with the major phenotypic impacts of PMA, growth arrest and macrophage differentiation.

Aside from these two clusters, the network analysis revealed that many transcripts were induced or repressed transiently. By 1 h after PMA treatment, early response genes with functions such as transcriptional activation were highly expressed (Cluster020 peaking at 2 h then declining rapidly and Cluster007 peaking at 4 h). Both these clusters were associated with GO terms relating to regulation of transcription (enrichment scores 2.7 and 3.9, respectively). Expression of the genes in Cluster003 began to rise after 1 h and peaked between 5 and 8 h. They were associated with both the cell cycle (enrichment score 4.9) and cell-cell adhesion (enrichment score 4.1). Expression of genes in Cluster004 slowly increased from the 2 h time point, peaking between 10 and 18 h, and this cluster was enriched for GO terms relating to lysosome (enrichment score 5.3). Images of the average expression profiles are shown in Figure 3C.

Cancer is characterized by abnormalities in cell proliferation, particularly the escape of the cancer cells from the normal controls on cell division (Sachs, 1982). In the analysis of the whole dataset, the second largest cluster (Cluster002) contained many genes associated with cell proliferation and the averaged expression profile of this cluster showed a decline during the time course with a small increase at 96 h. Note that by contrast to earlier studies (Nguyen et al., 2011; Garrison et al., 2017; Germano et al., 2017; Prange et al., 2017) (see section “Introduction”) the set of down-regulated transcripts did not include KMT2A (MLL), DOT1L, RUNX1 or ZNF521, implicated in primary leukemogenic transformation in AML. Each was expressed constitutively by THP-1 cells and unchanged upon differentiation. The cell cycle-associated cluster was explored in more detail by further clustering the expression profiles of the 2,158 genes separately, at higher stringency. A sample-to-sample analysis (r ≥ 0.9) showed that for this group of genes, the later time points were quite distinct, but the earlier time points (from 0 to 14 h) were relatively similar with a weak correlation between time point and position in the network (Figure 4A). The gene-to-gene analysis was run at a correlation coefficient threshold of 0.95, which included 1540 nodes connected by 17,446 edges. The network was then clustered with an MCL inflation value of 1.7. 618 nodes of the original Cluster002 were excluded at this threshold. We hypothesized that genes involved in different phases of the cell division cycle, or with different mRNA stability, might decay with distinct temporal profiles Consistent with that view, the annotated cell cycle genes appeared in specific subclusters with distinct profiles of decay with time (Figure 4B). For example, Subcluster002 with 99 nodes had none of the reported cell cycle genes (Giotti et al., 2018), a highly significant departure from expected (P = 7 × 10^–6). Overall expression in this cluster dropped very early in the time course. This cluster contains monocyte-associated transcription factors (IRF8, GFI1). These factors, as well as MYB and MYBL2 (which are also in the larger Cluster002 but did not form a part of a subcluster) have myeloid-specific roles in proliferation and down-regulation using siRNA can drive growth inhibition (Suzuki et al., 2009). Subcluster003 and Subcluster004 had a significant excess of cell cycle genes (P < 1 × 10^–5 and P = 3 × 10^–5 respectively). For these clusters, average expression declined slowly after 10 h and both had a small rise at 96 h, not seen in Subcluster002. Both these clusters had equal numbers of genes associated with S and G2/M phase of the cell cycle. In contrast Subcluster010 with 17 nodes had 10 G2M genes and no S phase genes (P = 5 × 10^–6). Average expression was relatively flat until 18 h when it dropped rapidly. Using the Bonferroni method to correct for multiple testing, a P-value of 0.002 would be equivalent to the accepted threshold for significance of 0.05, so these differences are highly significant. The full list of subclusters and averaged expression profiles are provided in Supplementary Table S2.

A previous study (Giotti et al., 2018) identified 701 genes that were expressed/induced in multiple cell types specifically in S phase and G2/M phases of the cell cycle. Of these, 561 (80%) were expressed in THP-1 cells, represented by 617 promoters. These 617 promoters were selected from the filtered dataset of mean values and analyzed separately. When clustered at R ≥ 0.75, 599 promoters were included in the analysis and 512 of these fell into the largest cluster (named CC-cluster01) with an average expression pattern which declined slowly until 10 h and then more steeply to 48 h, after which there was a small increase in average expression (Figure 4C). This group included the classic proliferation markers mentioned earlier. This finding provides further support for the annotation of these genes as cell cycle-related. A second cluster of 21 nodes (CC-cluster02) showed a small decline in expression from 0 h, with a large increase at 96 h. Other smaller clusters containing 2 – 11 unique genes showed similar expression patterns to the whole dataset (Figure 4C). The full CC-cluster lists and averaged expression profiles are provided in Supplementary Table S3. The set of 150 cell-cycle genes for which promoter activity was not detected in THP-1 cells is provided in Supplementary Table S4. In most cases, the apparent lack of expression is an issue of promoter annotation, and there is, in fact, expression of promoter activity in the region of the corresponding gene evident in the FANTOM5 dataset. For example, our data, and the FANTOM5 data, reveal that the transcript encoding the key regulator CDC25B, which is essential for G2/M cell cycle progression in AML (Didier et al., 2008) is expressed in THP-1 cells but is primarily driven by a cluster of three promoters that are internal to the annotated TSS in RefSeqs NM_001287516, NM_001287517, and NM_001287518. This can be seen in the data available on the ZENBU browser.

In summary the analysis of known cell cycle genes is consistent with the shut-down of entry into S phase prior to the 10 h time point as indicated by the cell cycle analysis (Figure 1) and subsequent slow decay of S phase-associated transcripts as cells complete the DNA synthesis phase of the cell cycle and are blocked in G2/M.

In cells undergoing transitions in state, transcription factor gene expression precedes the expression of the target genes being regulated (Arner et al., 2015). Of 1639 transcription factors curated recently (Lambert et al., 2018), 647 (39%) were expressed in the filtered THP-1 dataset. The expression profiles of these genes were used for further analysis. In a sample-to-sample analysis including all time points (R ≥ 0.61), the 0 h time point was an outlier, linking only with the 3 h time point (Figure 5A). A gene-to-gene analysis of the 647 transcription factor genes (796 promoters) was performed at a relatively stringent correlation coefficient of R ≥ 0.85, which included 540 nodes (transcription factor genes) and 3026 edges (connections between them at R ≥ 0.85). Clustering was performed at an inflation value of 1.7 and clusters were named as TF-cluster01 etc. The clusters are shown in Supplementary Table S5 which also contains histograms of mean activity of each cluster. The much denser time course analyzed herein, compared to the FANTOM4 project (Suzuki et al., 2009) highlights a temporal cascade in which the large majority of the transcription factor promoters exhibit regulated expression. Broadly speaking, there is a set of clusters where the transcripts are progressively down-regulated with different dynamics, a reciprocal set progressively up-regulated and a third set in which expression is regulated transiently (Figure 5B).

The transiently regulated transcription factors include the so-called immediate early genes (IEG) that have been analyzed in detail in other differentiation/activation time courses by the FANTOM 5 consortium (Vacca et al., 2018). That analysis identified a set of 42 protein-coding genes with promoters that consistently drive a rapid, but transient, wave of transcription in response to stimulus. A subset of these IEG was expressed by differentiating THP-1 cells. TF-Cluster15 contains transcription factors EGR1, EGR3, EGR4, JUNB, KLF2, KLF10 and MEF2D which were fully induced within 1 h and thereafter declined slowly. Other IEG had distinct patterns. EGR2, which has specific functions in myeloid differentiation including regulation of CSF1R transcription (Laslo et al., 2006; Krysinska et al., 2007), is part of TF-cluster10 (also including ELK3, FOXP2, JUNB, TBX3, and ZEB2) that peaks later around 8–10 h, around the time of the first appearance of adherent cells. ZEB2 was shown recently to be required for the establishment and maintenance of tissue macrophage identity (Scott et al., 2018). By contrast to other systems in which it is an IEG, and consistent with its known functions in macrophage differentiation (Shen et al., 2016), JUN was not induced rapidly but instead formed part of a progressive differentiation cluster (see below).

Mean expression of transcription factor genes in the largest cluster (TF-Cluster01) was maximal in untreated THP-1 cells and declined rapidly in response to PMA. This down-regulated cluster contains multiple transcripts associated with immature myeloid cells and maintenance of the proliferative state, including MYB, GFI1 (both shown in Figure 2B), IRF8 and RXRA. Average expression of transcripts within TF-Cluster03 declined more slowly; they include the core regulators of cell cycle gene expression, E2F1 (Figure 2B) and FOXM1 (Giotti et al., 2018). Several TF clusters show progressive up-regulation with distinct profiles. TF-Cluster04 (including JUN, MAF, NR1H3) increased across the whole time course, TF-Cluster11 (including ATF3, TCF7L2) was induced earlier, and TF-Cluster02 (including IRF1, FOXP1) was induced late, associated with terminal differentiation. One highly induced transcription factor with a unique expression profile that did not form part of a cluster was MAFB, discussed further below.

Promoters contain motifs that are recognized by DNA binding proteins. Where promoters share such motifs it is likely that they are regulated by the same proteins, and expression correlates with the number of copies of the relevant motif in the promoter region. The analysis tool MARA (Suzuki et al., 2009), which calculates motif activity based on this correlation, was used to analyze the response of THP-1 cells to PMA. The results are shown in Supplementary Table S6 and histograms of activity of selected motifs across the time course are shown in Supplementary Figure S2. A preliminary clustering analysis (r ≥ 0.6, MCL inflation 1.7) showed four patterns of activity: declining or increasing throughout the time course, transiently up-regulated and transiently down-regulated (Figure 5C). These patterns confirm the conclusions of the earlier work (Suzuki et al., 2009). The imputed transcriptional activity of MYB and E2F1.5 motifs, and of TFDP1, which binds cooperatively with E2F (Rubin et al., 2005), was down-regulated progressively over the time course. The activity of related motifs recognized by AP1 family members (FOS; FOS{B,L1}_JUN {B,D}; FOSL2; JUN) increased rapidly, as did the activity of the NFE2L2 and TGIF1 motifs (both factors induced by PMA).

There was some evidence of increased activity of the IRF1,2,7 motif, consistent with the induction of interferon target genes in Cluster032. Kubosaki et al. (2010) reported that the same genes were expressed constitutively in THP-1 cells and were targets of IRF8. Consistent with this and other studies (Phanstiel et al., 2017), IRF8 was highly expressed by THP-1 cells (the most abundant transcription factor mRNA) and repressed by PMA, but the putative target genes were low in undifferentiated cells and strongly induced. PMA also initiated clear and transient increased activity of the EGR1.3 motif. Also consistent with FANTOM4 data there was a decrease in activity of the NFY{A,B,C} and octamer (POU2F1.3) motifs, although NFYB and NFYC mRNAs (both highly expressed in THP-1 cells) were not regulated by PMA.

The greater resolution of the current dataset permitted the detection of regulated activity of additional motifs, notably the NR4A2 (also recognized by NR4A1) and TFEB motifs. Both NR4A1 and NR4A2 were amongst immediate early genes induced by PMA. The TFEB motif, recognized by members of the microphthalmia transcription factor family, is enriched in the promoters of lysosome-associated genes (Sardiello et al., 2009; Hume et al., 2010). Two members of this family, TFEC and TFE3, were constitutively expressed in THP-1 cells but two others, TFEB and MITF, were progressively induced in response to PMA in parallel with lysosomal hydrolases (e.g., ACP2, CTSB, LIPA) and components of the vacuolar ATPase (e.g., ATP6V0A1). Conversely, we did not reproduce evidence of inducible activity of SPI1 (PU.1), RUNX1.3 or SNAI1.3 motifs although transcripts encoding the related transcription factors TWIST1 and SNAI1 were induced by PMA.

The increased depth of sequencing available in this dataset also enables for the first time the detection of enhancer RNA (eRNA), which derives from the transcriptional activity of regulated enhancers, and inferred relationships with inducible promoters. The eRNAs are rapidly degraded by the exosome complex (Andersson et al., 2014) and accordingly, they are generally detected at <1 TPM. Detection is therefore subject to stochastic noise with median mapped reads of around 10 million per library. Nevertheless, integration of the entire dataset enabled detection of regulated bidirectional transcription associated with enhancers annotated in the FANTOM5 dataset and identified additional candidate enhancers. The identified enhancers were compared to those identified in FANTOM 5. 1353 of the less stringent set and 656 of the more stringent set were found in the full set of ∼63,000 enhancers (Andersson et al., 2014), indicating that this analysis has discovered some potential novel enhancers. Both sets were compared to the set of ATAC-seq peaks identified in analysis of chromatin architecture of differentiation (Phanstiel et al., 2017). Overlap was detected for 88 and 79% enhancers in the more and less stringent sets, respectively. One novel example highlighted in the chromatin analysis is the complex BCL6/TPRG1 locus. BCL6 is implicated in cell cycle regulation and interacts with TP53 and CDKN1A (Phan et al., 2005), whereas there is little known of the biology of TPRG1. The promoters of these two transcripts are separated by >200 kb (see ZENBU Browser). The TSS of TPRG1 is located upstream of the annotated RefSeq transcript NM_198485 at Chr3:188947237, downstream of a TATA box. The differentiation time course data confirm the rapid (within 1 h) and synchronous continued induction of both transcripts, and the detectable induction of enhancers with the intervening LPP locus. Figure 6 shows the profile of the activity of each of these promoters/enhancers across the time course. A similar co-regulated pattern was evident at the MITF locus. The induction of MITF with time in response to PMA was correlated with induction of the adjacent FRMD4B gene and evidence of bidirectional transcription initiation in the intervening genomic region.

To further confirm the validity of the enhancers found, they were compared to every individual FANTOM5 enhancer set. The THP-1 FANTOM5 library was the top hit for both sets, and other myeloid libraries were overrepresented in the top of the list. 219 out of 1172 FANTOM5 THP-1 enhancers were confirmed in the more stringent time course set. The results are available in Supplementary Tables S7A,B.

The macrophage-specific transcription factor, PU.1 (encoded by SPI1) is required for differentiation and controls the expression of many macrophage-enriched genes in THP-1 cells (Suzuki et al., 2009). SPI1 was expressed constitutively by THP-1 cells and was not highly regulated by PMA. As noted above, MARA did not reveal regulated activity of the PU.1/SPI1 motif across the time course. The expression of SPI1 in AML cells depends upon an enhancer around 16 kb upstream of the TSS that binds the transcription factor SATB1 (Steidl et al., 2007). The CAGE data (see ZENBU Browser) reveal that there is abundant unidirectional transcription initiation from this enhancer in THP-1 cells, likely driving an annotated lncRNA, that is marginally up-regulated by PMA. SATB1 mRNA, on the other hand, was progressively down-regulated in response to PMA.

Amongst other genes of interest, the MYB gene contains two upstream and two intronic enhancers. A candidate enhancer in the first intron was progressively down-regulated in response to PMA. This same region, conserved across species, has previously been implicated in regulated expression of Myb in mouse erythroleukemia cells (Suhasini et al., 1997). The data reveal a second major locus of regulated bidirectional transcription in intron 9 (Figure 7). Previous studies of regulated expression of CDKN1A have emphasized the impact of PMA on promoter activity and the binding of transcription factor SP1. The FANTOM5 study identified four enhancers in the region of CDKN1A. One of these, at Chr6:36667005-36667447 (around −12 kb relative to the TSS and conserved in mouse) was strongly and rapidly induced by PMA (Supplementary Figure S3). The progressive induction of the CSF1R gene in response to PMA was preceded by increased antisense transcription initiation in the first intron (Supplementary Figure S3), which contains two enhancers including the highly conserved FMS intronic regulatory element (FIRE) (Rojo et al., 2017). Finally, the FANTOM5 data identified candidate enhancers between the ITGAM (encoding CD11B) and ITGAX (CD11C) loci, both genes induced late in differentiation by PMA, and the CAGE data detect bidirectional transcription in the intergenic region (Supplementary Figure S3).