Transcriptional landscape of Mycobacterium tuberculosis infection in macrophages

Mycobacterium tuberculosis (Mtb) infection reveals complex and dynamic host-pathogen interactions, leading to host protection or pathogenesis. Using a unique transcriptome technology (CAGE), we investigated the promoter-based transcriptional landscape of IFNγ (M1) or IL-4/IL-13 (M2) stimulated macrophages during Mtb infection in a time-kinetic manner. Mtb infection widely and drastically altered macrophage-specific gene expression, which is far larger than that of M1 or M2 activations. Gene Ontology enrichment analysis for Mtb-induced differentially expressed genes revealed various terms, related to host-protection and inflammation, enriched in up-regulated genes. On the other hand, terms related to dis-regulation of cellular functions were enriched in down-regulated genes. Differential expression analysis revealed known as well as novel transcription factor genes in Mtb infection, many of them significantly down-regulated. IFNγ or IL-4/IL-13 pre-stimulation induce additional differentially expressed genes in Mtb-infected macrophages. Cluster analysis uncovered significant numbers, prolonging their expressional changes. Furthermore, Mtb infection augmented cytokine-mediated M1 and M2 pre-activations. In addition, we identified unique transcriptional features of Mtb-mediated differentially expressed lncRNAs. In summary we provide a comprehensive in depth gene expression/regulation profile in Mtb-infected macrophages, an important step forward for a better understanding of host-pathogen interaction dynamics in Mtb infection.


Results
Construction of time-course promoter activity profiles for Mycobacterium tuberculosis-infected macrophages. Macrophages are innate immune cells that can be temporally activated to classically (M1) or alternatively (M2) by IFNγ or IL-4 and IL-13, respectively. Macrophages are also the primary (obligatory) target cells for Mycobacterium tuberculosis (Mtb) infection. Using in vitro cultured primary mouse bone marrow-derived macrophages (BMDMs), we designed and characterized the transcriptional and kinetic (0, 4, 12, 24, 48 hours) landscape of macrophages during Mtb infection using hyper-virulent strain Bejing Mtb HN878. Differentially expressed genes in each condition/time point were compared among M1-and M2-stimulated BMDMs and non-activated Mtb-infected BMDM's ( Fig. 1a, see comparison I). Further, we compared differentially expressed genes of M1-and M2-pre-activated Mtb-infected BMDM's (IFNγ_Mtb and IL-4/IL-13_Mtb) with those of non-pre-activated Mtb-infected BMDM's ( Fig. 1a, see comparison II), which was designed to explore the differential effects of macrophage M1 and M2 pre-activation in Mtb infection. Total RNA was extracted from harvested cells and subjected to non-amplified deepCAGE transcriptomics using single molecule Helicos Scientific RepORTS | (2018) 8:6758 | DOI: 10.1038/s41598-018-24509-6 sequencers. A flow chart describing the CAGE analysis pipeline was shown for the better understanding of deep-CAGE transcriptome analysis ( Supplementary Fig. S1). The CAGE tags, computationally mapped on mouse genome (version mm9), ranged from 460,316 to 18,969,698 with a median of 2,655,669 mapped tags among all these libraries (Supplementary Table S1, also http://fantom.gsc.riken.jp/5/). Tags were clustered to identify transcription start site (TSS) regions and create promoter activity profiles (see Materials and Methods). The profiles of three biological replicates showed high correlation in each data point ( Supplementary Fig. S2). To verify the quality, we explored the induction of typical and well-known marker genes for M1/M2 activations and Mtb infection; as expected, expression of these marker genes were drastically up-regulated during M1/M2 activation ( Supplementary Fig. S3) as well as during Mtb infection ( Supplementary Fig. S4), suggesting that both activation and infection were successful. Further, we confirmed the CAGE expression profile of typical marker genes (such as Tnf, Il10, Il1b, Arg1) in Mtb, IFNγ_Mtb, IL-4/IL-13_Mtb by qRT-PCR ( Supplementary Fig. S5).
For further examining robustness of the genome-wide transcriptomics experiments, promoter level expression data was subjected to the principal component analysis (PCA) in each comparison group (Fig. 1b,c,  Supplementary Fig. S6a and b). In comparison I, PCA results suggested that data for Mtb samples were largely separated from IFNγ and IL-4/IL-13 samples during the time kinetic ( Fig. 1b and Supplementary Fig. S6a). This suggests that Mtb infection is clearly distinct from macrophage activation, as biologically expected. PCA analysis was also performed for the comparison II among Mtb, IFNγ_Mtb and IL-4/IL-13_Mtb. The results indicate that these three group data were separated in a time dependent manner, with individual time-points of different treatment clustering together ( Fig. 1c and Supplementary Fig. S6b). This suggests that differently activated macrophages may react similar during Mtb infection, or that the effects of the much stronger Mtb stimulus are overshadowing the effects of activation.
Redefining the transcriptional regulatory dynamics of Mtb-infected macrophages. In order to characterize the transcriptome of Mtb-infected macrophages, we extracted differentially expressed genes (>2 fold change, FDR < 0.05) at each time point (Table 1 and Supplementary Table S2a and b). Mtb infection drastically increased the number of differentially expressed gene transcripts (>10 fold) in comparison to the transcriptional activity of M1 or M2 activation (Supplementary Table S2c and S2d Table), observed until 48 hours post Mtb infection. In total, 1394 and 915 genes were significantly up-and down-regulated during Mtb infection, respectively. This amounted to 25.4% of the total 9,052 expressed genes of non-stimulated macrophage at 0 h. Of note, Mtb infection also influenced the transcriptional activity of down-regulated genes, particular at 4 and 12 h with 38.3% and 41.7% of differentially expressed genes, respectively, when compared with IFN-or IL-4/IL-13-stimulation at the same time points (11.1~26.3%) ( Table 1). For more insights into the differential expression between Mtb infection and macrophage M1 or M2 activation, we explored commonly altered genes. Of interest, IFNγ-induced M1 genes were particularly highly overlapping with Mtb-induced genes with 71 out of 92 (75%), see parenthesis in Table 1. The overlapping genes included Nos2, Tnf, Cxcl9, Cxcl10, Irg1, among others, consistent with our previous reports (and others) that Mtb infection does induce interferon-related inflammatory responses in macrophages 26,38,39 . Taken together, Mtb infection largely affects macrophage gene expression including inflammatory response genes.
To further explore the global effect of Mtb infection in macrophages, differentially expressed genes were subjected to gene ontology analysis at each time point (Fig. 2a,b and Supplementary Table S3a and b). As might be expected from above, Mtb infection induced up-regulated genes revealing significant enrichment of the ontology term "immune response" (GO: 0006955), particularly at 12 and 24 h (Fig. 2a), which was also enriched in IFNγ up-regulated genes at 4 h (Supplementary Table S3c Table S2a). The ontology term "response to wounding" (GO: 0009611) was also highly enriched for Mtb-up-regulated genes during 4 to 24 h. At earlier time point from 4 and 12 h, apoptosis-related ontology terms, "regulation of apoptosis" (GO: 0042981), "regulation of programmed cell death" (GO: 0043067) and "regulation of cell death" (GO: 0010941), were enriched for the up-regulated genes and disappeared at 24 h. In addition to the above ontology terms, a time-dependent shift of other ontology terms were observed; those were nucleotide binding related terms at 4 h, "nuclectide binding" (GO: 000166), "purinenucleotide binding" (GO: 0017076), "ribonucleotide binding" (GO: 0032553), "purine ribonucleotide binding" (GO: 0032555), and terms related to inflammation regulation at 12 and/or 24 h, "regulation of cytokine production" (GO: 0001817), "negative regulation of molecular function" (GO: 0044092), "cytokine activity" (GO: 0005125), "inflammatory response" (GO: 0006954), "regulation of leukocyte activation" (GO: 0002694), and "regulation of cell activation" (GO: 0050865). Interestingly, term related to inflammation regulation was further enriched at 48 h such as "positive regulation of inflammatory response" (GO: 0050729), "inflammatory response"(GO: 0006954), "negative regulation of cell proliferation" (GO: 0008285), "chemokine-mediated signaling pathway" (GO: 0070098), "positive regulation of ERK1 and ERK2 cascade" (GO: 0070374). Time-dependent shift of gene ontology terms were also observed for Mtb down-regulated genes ( Fig. 2b and Supplementary Table S3b). In particular at the late time point at 24 h, Mtb down-regulated genes were highly significant enrichment for ontology terms such as "nucleosome assembly" (GO: 0006334), "chromatin assembly" (GO: 0031497), "protein-DNA complex assembly" (GO: 0065004), "DNA packaging" (GO: 0006323), which may imply that possible dis-regulation of basic cellular functions may occur in macrophages around 24 h post Mtb-infection. Together, the marked transcriptional regulatory dynamics of Mtb-infected macrophages seen might be the result of the cellular war between host immune defense and pathogen evasion responses.

Differentially expressed TF genes in Mtb infection. The gene expression dynamics during Mtb infec-
tion underlays global changes of transcription factor (TF) gene expression. Hence, we explored differentially expressed TF genes in Mtb-infected macrophages. Ninety-nine and sixty TF genes were significantly up-and down-regulated in response to Mtb infection (Supplementary Table S2a and b). Considering 953 TF genes were expressed in BMDM at time 0 h, 16.6% (159/953) of TF gene expression was significantly altered in Mtb infection. The average expression features of up-or down-regulated TF genes revealed a rapid alteration at 4 h and quickly return to the original expression level in the following time point (Fig. 2c,d), consisting with the up-and down-regulated gene expression as this time point.
Next we examined how M1/M2 up-regulated TF genes in this study are affected by Mtb infection. We found that 5 out of 7 M1 up-regulated TF genes (Batf2, Irf1, Stat1, Zfp281 and Zfp800) and 3 out of 4 M2 up-regulated TF genes (Bhlhe40, Egr2 and Tfec) were also up-regulated during Mtb infection (Supplementary Table S2a

Pre-stimulation with IFNγ and IL-4/IL-13 showed global effects in transcriptomics of Mtb-infected macrophages.
Next, we explored how Mtb infection with IFNγ or IL-4/IL-13 pre-stimulation affects gene expression changes in macrophage (comparison II in Fig. 1a). As shown in the comparison I (Table 1), IFNγ-and IL-4/ IL-13-mediated gene expression changes were temporal, as a consequence few differentially expressed genes were observed at 24 h post stimulation, the starting time for Mtb infection. Differentially expressed genes were extracted at each time point after Mtb infection with IFNγ or IL-4/IL-13 pre-stimulation (IFNγ_Mtb or IL-4/IL-13_Mtb, respectively). We obtained 3162 and 2935 differentially expressed genes for IFNγ_Mtb and IL-4/IL-13_Mtb, respectively (Supplementary Table S4a to d). Including the differentially expressed genes in non-stimulated and pre-stimulated Mtb, a total of 3806 differentially expressed genes were found overall. More than half of the genes (1884 genes) were commonly altered in three conditions (see Venn diagram in Fig. 3a). In non-stimulated Mtb-infected macrophages, the number of differentially expressed genes was small (75 genes), while in IFNγ_Mtb or IL-4/IL-13_Mtb we found that the number was several fold higher with 565 and 464 differentially expressed genes, respectively. This suggests that M1 and M2 pre-stimulation promotes distinct and additional differentially expressed genes in Mtb infection. We next explored time kinetics of differentially expressed genes of the three conditions (Fig. 3b,c). The early time point at 4 h, commonly up-and down-regulated genes were dominating, but rapidly decreasing with time, with small number of common genes expressed at 48 h. The differential genes in non-stimulated Mtb-infected macrophages were mainly found in 4 and 12 h post infection, reclining drastically thereafter.
This scenario was very different in up-and down-regulated differentially expressed genes for IFNγ_Mtb and IL-4/IL-13_Mtb. As shown in Fig. 3b in a kinetic expression, up-and down-regulated differential genes in IFNγ_Mtb raised within 12 h to over 500 genes and stayed during the 48 h kinetic. A slightly different dynamic was observed in IL-4/IL-13_Mtb. Indeed, differentially up-regulated genes showed a similar kinetic dynamic as IFNγ_Mtb, raising within 12 h close to over 500 genes and only reclining slowly during the 48 h. Together, Scientific RepORTS | (2018) 8:6758 | DOI:10.1038/s41598-018-24509-6 this demonstrates that Mtb infection within M1 and M2 has drastic transcriptional changes, compared to Mtb-infected macrophages.
To identify groups of genes that share common patterns of expression change within the duration of the time courses, we clustered the 3806 differentially expressed genes from Fig. 3a by k-means clustering. We clustered the differentially expressed genes into 7 clusters based on their expression profile (Fig. 3d). Cluster 1 showed progressive up-regulation, clusters 2, 3 and 7 showed sustained up-and down-regulation during time, and Cluster 4 to 6 transient up-or down-regulation. The number of genes in each cluster shows that Mtb has small populations of progressively upand sustained down-regulated genes (clusters 1 and 7, respectively) in comparison with IFNγ_Mtb or IL-4/IL-13_Mtb (Supplementary Table S5). The contents analysis revealed that genes involved in clusters 1 and 7 of IFNγ_Mtb and IL-4/IL-13_Mtb consist of higher population of un-differentially expressed genes in Mtb (Table 2). Further, significant population of sustained up-regulated genes in cluster 2 of IFNγ_Mtb and IL-4/IL-13_Mtb were transiently up-regulated genes in cluster 4 of Mtb. Similarly, significant population of sustained down-regulated genes in clusters 7 of IFNγ_Mtb and IL-4/IL-13_Mtb are transiently down-regulated genes in clusters 6 of Mtb. The shifts indicate that pre-stimulation prolonged differential expression in a set of genes. On the other hand, there are few genes that shift opposite direction. Taken together, the results demonstrate global effect that IFNγ_Mtb and IL-4/IL-13_Mtb reveal differentially expressed genes for long time which is in part due to prolonged gene expression change.
Finally, we carried out gene ontology enrichment analysis using exclusively differentially expressed genes for IFNγ_Mtb or IL-4/IL-13_Mtb (565 and 464 in Fig. 3a, respectively). We found that another set of genes involved in the ontology term "immune system process" (GO: 0002376), "innate immune response" (GO: 0045087) are  Table S3d and e). These results indicate that IFNγ or IL4/IL13 pre-stimulation in Mtb infection augments and prolongs M1 or M2 activations, respectively.

Differentially expressed lncRNAs in Mtb-infected macrophages. Long non-coding RNAs (lncR-
NAs) play important roles of regulation of gene expression in various ways, although the function of majority of lncRNAs is unknown. We explored differentially expressed lncRNAs at each time point in Mtb infection. We obtained a total 151 differentially expressed lncRNAs species, as detailed for Mtb (51), IFNγ_Mtb (129) or IL-4/IL-13_Mtb (91) (Fig. 6a and Supplementary Table S7a to c). It is evident that pre-stimulations distinctly enhanced the number of differentially expressed lncRNA, as well as protein-coding genes, although higher number of lncRNAs were differentially expressed in IFNγ_Mtb (129/151 = 85.4%) in comparison with IL-4/IL-13_ Mtb (91/151 = 60.3%). Interestingly, Venn diagram analysis revealed that only 33 out of 151 totally differentially expressed lncRNAs (21.9%) were commonly altered (Fig. 6a), which is in contrast with 1884 out of 3806 (49.5%) in protein-coding genes (Fig. 3a). Another interesting feature was that the majority of differentially expressed lncRNAs were down-regulated (42/51 = 82.3%, 96/129 = 74.4% and 69/91 = 75.8%) in Mtb, IFNγ_Mtb and IL-4/ IL-13_Mtb, respectively) (Supplementary Table S7). This is in contrast with that major differential expression of  protein coding genes were up-regulated in Mtb infections and M1 or M2 activations (Table 1). Next we explored the number of differentially expressed lncRNAs at each time point (Fig. 6b,c). The number slightly decreased within the three conditions, as there was no drastic transient up-regulation nor down-regulation in differentially expressed lncRNAs, which was in contrast with those in protein-coding genes (Fig. 3b,c). We also investigated expression correlation between differentially expressed lncRNAs and their nearest protein coding genes ( Fig. 6d and Supplementary Table S7a to c). It became evident that there is a shift to positive correlation on average in Mtb-infected samples, either with or without M1 or M2 pre-activations, although each correlation value is variable (Fig. 6d). Interestingly, number of lncRNAs, with their expression extreme positive correlation (r > 0.90) with that of their nearest genes, was only 1 for Mtb but drastically increased to 5 for IFNγ_Mtb and 6 for IL-4/IL-13_Mtb (Supplementary Table S7a to c).
Conclusively, Mtb infection in M1 as well as in M2 drastically induced lncRNAs compared to unstimulated macrophages.

Discussion
In this study, we comprehensively analyzed the transcriptome of Mtb-infected macrophages and the effect of IFNγ (M1) or IL-4/IL-13 (M2) pre-stimulation in a time-dependent manner using our unique CAGE technology. Our results showed that Mtb infection widely and drastically altered macrophage gene expression including induction of inflammatory response genes, which is far larger than that of M1 or M2 macrophage activation. This drastic gene expression change was mediated by expressional alteration of various known and novel Mtb-infection related TF genes including both M1 and M2 activation. We also showed that M1 or M2 pre-activation induces global effect on transcriptional landscape of Mtb-infected macrophages by augmenting M1 or M2 genes and promoting additional genes. Furthermore, we described distinct features of lncRNA expression in Mtb-infection.
Mtb-induced alteration of large number of macrophage gene expression has already been reported in several previous studies 4 . We extended this field in detail by utilizing a kinetic Mtb infection approach. One of the major role of macrophages is to produce inflammatory response, processing of antigen through MHC class I molecule and to engage in pathogen killing, which is regulated by M1 activation 13,40 . However, GO enrichment analysis of Mtb-induced up-regulated genes revealed that not only M1-related immune/inflammation response genes but also M2 genes related to cell wounding and apoptosis are involved. Actually, many of both M1 and M2 TF genes were up-regulated in Mtb infection. Thus, this Mtb-mediated wide range response may be partly responsible for alteration of large number of gene expression in macrophages. We also observed high population of Mtb-induced down-regulated genes, which is consistent with relatively high number of down-regulated TF genes. Because down-regulated genes are rare in M1 or M2 activation, this observation is one of characteristic in Mtb infection. Actually, GO terms related to dis-regulation of basic cellular functions are enriched in the down-regulated genes. However, it is elusive whether those down-regulated genes are pathogenic effect of Mtb infection, since it may be a side reaction that macrophages converge host-protective response.
We showed that M1 or M2 pre-activation in Mtb-infected macrophages augmented M1 or M2 effector genes, respectively. Although molecular mechanism how the augmentation occurs remains to be elucidated, this reaction is interesting from the view of host protection and subversiveness. Mtb infection in M1 pre-activated macrophages enhances expression of inflammatory genes, which seems a reasonable host-protective reaction 4,26 . On the other hand, because Mtb favors to reside within M2 polarized macrophages 20 , augmentation of M2 effector genes in Mtb-infected M2 pre-activated macrophages is considered to be host subversive. The later reaction may be host evasion mechanisms employed by Mtb to establish its persistence and survival within M2 polarized macrophages.
We also successfully identified differentially expressed lncRNAs in Mtb-infected macrophages. Although M1 or M2 pre-activation enhanced altered lncRNAs as well as coding transcripts, the number of commonly altered lncRNAs was not high and the majority of the alteration was indeed down-regulation 41,42 . We showed in a previous report that majority of M1-or M2-mediated differentially expressed lncRNAs was up-regulation 37 . Thus, unique feature of altered lncRNAs in Mtb infection may have significance. Recently, many lncRNAs have been demonstrated to play a role of transcriptional regulation of their neighbor genes 43 . In fact, we detected positive correlation of expressional change between altered lncRNAs and their neighbor genes (Fig. 6d). These functional perturbation analysis of newly identified lncRNAs should enable us for a deeper and better understanding of the role of these transcripts in Mtb-infected macrophages.
This comprehensive promoter-based transcriptome data of Mtb-infected macrophages, together with previously reported data, will be a valuable resource for the research community, particular in immunology, for gaining new insights into Mtb host evasion. Newly identified protein-coding and noncoding transcripts, which altered in Mtb infection in the presence and absence of M1 or M2 pre-activation, may serve as potential transcriptional biomarkers of Mtb infection. Further, those altered transcripts could be a potential source, leading to host-directed Mtb drugs, as part of them may be regulated by Mtb to be host subversive, and may be resistant to drug tolerance and prolonging current Mtb drugs.

Materials and Methods
Generation of bone marrow-derived macrophages. BALB/c mice were purchased from Jackson Laboratories and bred in South Africa. Bone marrow-derived macrophages (BMDMs) were generated from [8][9][10][11][12] week old male BALB/c mice as described previously 37 Fig. S7). Cells were treated as mentioned above with Qiazol and stored at minus 80 o C for RNA extraction. All Mtb infection experiments were performed at the Biosafety Level 3 (BSL-3) laboratory, Faculty of Health Sciences, University of Cape Town. Town Total RNA was prepared using miRNAeasy kit (Qiagen, Valencia, CA, USA) and concentration and quality of each RNA samples was checked as described previously 37 . Preparation of Helicos CAGE library, sequencing. Single molecule Helicos CAGE library was prepared, sequenced, mapped and clustered into TSS regions as described previously 36,37 . Briefly, each library was prepared by using 5 μg of total RNA, with RIN value of more than 7.5 (Supplementary Table S1). Three to four biological replicates were prepared per time point.
Construction of promoter data. To identify peaks (TSS clusters) in the CAGE profiles, decomposition peak identification (DPI) was used as described previously 36 . To calculate normalization factors for the expression of promoters we used "relative log expression (RLE)" 44 method as described previously 33,37 . Reproducibility between replicates was assessed by computing pair wise Pearson correlations between samples across TSS regions after addition of a pseudo count of one and logarithmic transformation (Supplementary Fig. S2).
Principal component analysis. Principal component analysis (PCA) was performed using the R-function 'prcomp' on a TSS level. TSS regions of triplicates were averaged prior to PCA. We only considered TSS regions with a minimum of 5 tags in at least one sample 37 .
Differential expression analysis and qualitative analysis of protein-coding genes. Reads mapped to DPI/promoter regions were counted and used as input to perform differential expression (DE) analysis as described previously using edgeR 44,45 . IFNγ-and IL-4/IL-13-stimulated BMDM samples at time-points 4, 12, 24, and 48 h post-stimulation were compared to unstimulated BMDM samples at the same time-points, except for 4 and 12 h time-points, which were compared to unstimulated 0 h samples as described previously. Furthermore, Mtb infected (Mtb), IFNγ or IL4/IL13 pre-stimulated and Mtb infected (IFNγ_Mtb or IL-4/IL-13_Mtb) samples at time-points 4, 12, 24, 48 h were compared with unstimulated samples of each time points. Differentially expressed genes were selected using an FDR < 0.05 and a log2 fold-change >1 (log2 fold <−1 in case of down-regulation) in each comparison. Finally, we filtered for robustly expressed DE genes using a 20 tags-per-million (TPM) threshold in the respective comparison. We also created differentially expressed genes for IFNγ_Mtb or IL-4/IL-13_Mtb by comparing with IFNγ or IL-4/IL-13 pre-activated BMDM, resulting in similar differentially expressed genes with the above comparison. Therefore we did not use this comparison for our analysis.
The qualitative analysis was performed using DE gene data of Mtb, IFNγ_Mtb, IL-4/IL-13_Mtb samples. Gene expression up-or down-regulation was calculated based on TPM expression fold-changes between Mtb, IFNγ_ Mtb, and IL-4/IL-13_Mtb samples. Several-fold elevated or suppressed genes were identified by comparing the expression values of IFNγ_Mtb with Mtb, IL-4/IL-13_Mtb samples. Similarly, several-fold elevated or suppressed genes in IL-4/IL-13_Mtb samples were identified by comparing them with Mtb, IFNγ_Mtb samples. We selected non-TF genes (more or less than 10-fold expressed) and TF genes (more or less than 3-fold expressed) among all significantly differentially expressed genes.
Differential expression analysis of lncRNA promoters. Differential expression analysis of lncRNAs was performed by counting all mapped reads within the DPI/promoter regions as described previously using edgeR 44 . We subselected transcripts that had at least 1 TPM expression in at least two replicates in each sample groups. Differentially expressed lncRNA promoters were selected using FDR < 0.05 (Benjamini-Hochberg FDR) and a log2 fold-change >1 (log2 fold <−1 in case of down-regulation) in each comparison.
Gene ontology analysis. Gene ontology analysis was performed using the database for Annotation, Visualization and integrated discovery (DAVID). The top gene ontology terms were selected using p value of 0.001 and FDR 0.05. K-means clustering analysis. We performed k-means clustering to identify groups of genes that share common patterns of expression change within the duration of the time course. First, for each time-course experiment we extracted the normalized expression values of differentially expressed genes. The triplicates were averaged into one expression profile. We then merged the expression tables from 3 time-courses, thus a given gene is represented in each time-course in which it was found differentially expressed. In total, we obtained 8,426 expression profiles (since most of the genes were differentially expressed in more than 1 condition). The data were log10 transformed, after adding an offset of 0.1. Subsequently, for each gene we scaled the data to the mean expression of 0 within the control time points and to standard deviation of 1. This way we calculated the expression values that were relative to the average control time points, which allowed us to cluster the time-dependent response pattern relative to the baseline. We then clustered the 8,426 expression profiles into 7 clusters using k-means algorithm as implemented in R (seed fixed at: 123, default settings).