Transcriptome analysis of mouse macrophages exposed to acute glucocorticoid and LPS stimulation
To analyze early regulatory events initiated by glucocorticoids and inflammatory stimuli we treated BMMФ with either ethanol vehicle (U), LPS (L), Dex (D), or a combination of the two (L + D) for 1 h, isolated and sequenced PolyA-enriched RNA as described in (Additional file 1). The sequencing results are summarized in Additional file 2: Table S1.
To uncover the regulatory patterns in gene expression data, we performed k-mean cluster analysis of ANOVA-filtered differentially expressed genes (see Additional file 1). To minimize magnitude-based clustering, the log2-transformed expression values were first converted to Z-scores as in where Xgene is an expression value for a given gene at a given condition, is an average expression value across conditions and σ is a standard deviation of across conditions. The optimal number of clusters was determined using the “elbow” method by plotting within-cluster variance vs. the number of clusters. The analysis was performed using Euclidian distance with progressively increasing number of clusters from 8 to 14 to determine a stable configuration using the cluster analysis module of STATISTICA 8.0. We grouped ANOVA-filtered data into 12 clusters of co-regulated genes (Additional file 2: Table S2). We further evaluated the significance of differences in gene expression within clusters by the Mann–Whitney test and validated the results for a limited number of genes by RT-qPCR. As many glucocorticoid- and LPS-responsive genes have been previously characterized by us and others [20–22], we selected for validation either poorly characterized or novel target genes. Based on patterns of co-regulation, we grouped these clusters into four larger categories (Figure 1).
I) Genes activated by either LPS (Cluster 2) or glucocorticoids (Cluster 7)
The majority of genes in these two clusters are upregulated by either LPS or Dex independently, with little evidence for inter-treatment interactions at 1 h. The LPS-induced Cluster 2 (Figure 1, Mann–Whitney, PU-L = 4.11*10-13) contains genes encoding pro- and anti-inflammatory cytokines and chemokines (Il10, Cxcl1, 3, 5 and 7, Ccl7 and Tnfsf9), TFs involved in stress response (Maff, Ets2, Fosl2 and Kdm6b) and proteins involved in TLR signaling (Tlr2, Cd14 and Cd40) and signal transduction (Itpkc, Rabgef1, Gbp5a).
Cluster 7 contains Dex-induced genes (PU-D = 0.0013), including several well-characterized GR targets such as TFs Per1 and Klf9, immunophilin Fkbp5, potassium channel Kcnk6. In addition, this cluster includes several genes whose regulation by Dex has not been previously reported: Interleukin 15 receptor alpha (Il15ra), the Wnt pathway receptor Fzd4, the TF Klf2 and chemokine Ccl17.
II) Genes co-activated by LPS and Dex
These genes display either predominantly additive (cluster 1) or synergistic (cluster 6) activation by LPS and Dex. Several genes in these clusters are previously characterized GR targets including Dusp1, Nfil3 and Cited2 (cluster 1) and Mt2 and Pfkfb3 (cluster 6), whereas the glucocorticoid and LPS responsiveness of several others, such as histone chaperon Jdp2 (cluster 1), has not been reported previously.
III) Genes induced by LPS and repressed by Dex
This large group of genes is represented by clusters 3, 4 and 5. Cluster 3 contains LPS-induced genes (PU-L = 1.95*10-11) expressed at relatively high level in resting BMMФ. The basal expression of these genes is significantly more sensitive to hormonal treatment (PU-D = 0.0107) than their LPS-induced expression. This cluster encompasses a number of inflammatory cytokines (Ccl2, 3 and 4, Tnf, Tnfaip2), TFs (Ier5, Junb, Bcl6, Prdm1 and Irf1) and proteins involved in signal transduction (Gadd45b, Dusp5, Rasgef1b). Interestingly, several genes in this cluster (Ccl2, 3 and 4, Tnf) are characterized by the presence of the stalled RNA Pol II near the transcription start site in uninduced conditions and are activated primarily at the level of the Pol II pause release during early elongation [23–25]. Cluster 4 combines a heterogeneous group of genes with low basal expression (PUcluster 4<Ucluster3 = 0.0003) that are strongly induced by LPS and includes inflammatory cytokines (Il1a, Il1b, Cxcl10, Ifnb1, Tnfsf4 and Il1f9) and other direct mediators of inflammation (Ptgs2, mIR-155 host gene, Hsp1a), TFs (Egr3) and signaling molecules (Gpr84, Areg). Cluster 5 contains many genes whose LPS induction is strongly inhibited by Dex treatment (PL-(L+D) = 0.02) including several cytokines and chemokines (Il12b, Lif, Il1rn and Il17ra ) and TFs (Mxd1, Etv3, Klf7 and Ets1).
Cluster 8 contains statistically heterogeneous previously reported (Atf3, Egr2 and Ier3) as well as novel (Enc1 and Bhlhe40) targets for GR-mediated repression that are largely unaffected by LPS treatment. Thus, this cluster is formally outside of group III.
IV) Genes downregulated by LPS
LPS-repressed genes were separated into 3 clusters based on the combined effect of Dex and LPS on gene expression. The LPS-mediated downregulation is either weakly potentiated (cluster 12) or antagonized (cluster 11) by Dex. The functions of the majority of these genes are poorly understood or unknown, however, 32% of genes in cluster 11 and 20% in cluster 12 encode uncharacterized C2H2 Zn-finger proteins implicated in transcription and chromosome maintenance.
Expression of genes in cluster 10 is downregulated by Dex and LPS in an additive manner. At least one gene in this cluster (Angptl4) is a known GR target. Several genes encode regulators of immune cells activities including Rit1 and Cd300lb. Cluster 9 contains Dex-induced genes that are weakly repressed by LPS including previously reported Ddit4, Arl4d and Sik1.
Cluster validation by RT-qPCR
Two representative genes in each cluster were chosen for validation. Total RNA was isolated from treated (D, L, or L + D for 1 h) and control BMMФ and transcript levels for indicated genes (Figure 1) were determined by RT-qPCR using Act1 or Hprt as housekeeping control genes. As several LPS-induced, Dex-repressed genes coding for various cytokines that we found in clusters 3 (Tnf, Ccl2, 3 and 4), 4 (Il1a and b) and 5 (Lif, Niacr1, Mmp13) have been previously characterized by us and others [20–22], we focused on uncharacterized LPS/Dex targets. The expression patterns determined by RT-qPCR closely resembled those determined by RNA-seq, with the exception of several weakly expressed genes that were not detectably repressed by Dex following a 1-h treatment (e.g., Bcl2, Med21; data not shown).
Glucocorticoid-regulated genes form a highly interconnected association network with distinct response-specific modules
To determine the prevalent functional patterns in groups of co-regulated genes, we performed gene enrichment analysis and visualization using GeneMANIA plugin [26] for Cytoscape 2.8 and Exploratory Gene Association Networks (EGAN) software [27].
Using a list of genes that were up- (clusters 1, 6, 7 and 9) and down- (clusters 3, 4, 5, 8 and 10) regulated by Dex, we generated a consensus association network consisting of 333 nodes and 8296 edges. To discern underlying data structure in the Dex-regulated network, we decomposed this network using the Newman-Girvan community clustering algorithm [28], a divisive procedure that iteratively removes network edges with largest “edge betweenness”, recalculates this metric for a novel network and repeats the procedure until the network is split into several groups. The Newman-Girvan algorithm disregards edge weights and uses only connectivity to define communities. Community clustering partitioned Dex-regulated network into three unequal modules that were significantly enriched with Dex-repressed genes (Module 1, Figure 2A), Dex-induced genes (Module 2, χ2 = 18.33, p = 0.0001, Figure 2B) or LPS-repressed genes (Module 3, χ2 = 15.347, p = 0.00046, Figure 2C).Network topology analysis of Dex-responsive modules revealed significantly greater network densities and clustering coefficients in Modules 1 and 2 compared to networks generated from the same number of non-expressed genes extracted from the same BMMФ experiments (Figure 2D). Similarly, the average neighborhood connectivities and average clustering coefficients for Modules 1 and 2 were considerably greater than the values for non-expressors (Figure 2E), indicating that nodes in these modules form tight interconnected local groups. A broader shared neighbors distribution in Modules 1 and 2 indicates high prevalence of shared nodes, suggesting an enrichment for multiple input motifs. Interestingly, Module 2 (enriched with Dex-induced genes) was more structured than Module 1 (predominantly Dex-repressed genes) as evidenced by consistently higher average neighborhood connectivities and average clustering coefficients (Figure 2E). Although Module 3 has a non-random composition, with the exception of network heterogeneity, all other analyzed topological metrics were typical of networks composed of randomly selected non-expressed genes (Figure 2D). Therefore, we focused the rest of the analysis of Modules 1 and 2. Relatively high heterogeneity for both modules (0.506 and 0.581; Figure 2D) indicates the presence of network hubs - nodes with high degree of connectivity. Indeed, 10 most connected nodes in Modules 1 and 2 account for 36 and 32% of all edges, respectively (Figure 2F). Interestingly, eight out of 10 most connected nodes in Module 2 are sequence-specific DNA-binding TFs (bold in Figure 2F).
Glucocorticoid response-specific modules are functionally distinct
Using gene ontology (GO) analysis to determine enriched gene categories in subsets of functionally related genes, we identified 425 enriched GO categories for Module 1, 285 for Module 2 (FDR corrected p < 10-3) and 77 for Module 3 (FDR corrected p < 10-2). Only 115 GO categories overlapped between Modules 1 and 2 (Figure 3A, 3C). To facilitate visualization and interpretation of these results and compare enriched functional categories among groups of Dex-regulated genes, we generated GO terms similarity networks using Gene Set Enrichment Mapping Cytoscape plug-in. Multiple GO categories related to regulation of metabolic processes, embryonic and post-embryonic development and regulation of apoptosis and signaling are enriched in Module 2 (Figure 3A) that contains a large number of Dex-upregulated genes. Notably, 32/285 GO categories enriched in Module 2 were related to regulation of gene expression, regulation of transcription, sequence-specific DNA binding transcription factors. For example, negative regulation of gene expression (GO:10629), negative regulation of transcription (GO:16481), negative regulation of transcription - DNA-dependent (GO:45892), sequence-specific DNA binding (GO:43564), negative regulation of transcription from RNA polymerase II promoter (GO:122) and transcription regulator activity (GO: 30528) were all enriched in Module 2, but not Module 1 (Figure 3C). Overrepresentation of genes coding for regulators of gene expression in the early Dex-responsive transcriptome suggests that GR initiates a transcriptional program that relies on the step-wise activation of multiple TFs. Only a few categories related to immune/inflammatory responses have been found in Module 2 (Figure 3A).
Conversely, the majority of enriched GO categories in Module 1, which contains predominantly Dex-repressed genes, are related to immune and inflammatory responses, signaling and regulation of signal transduction and metabolic regulation including immune response (GO:6955), immune system process (GO:2376), inflammatory response (GO:6954) and regulation of cytokine production (GO:1817) (Figure 3B, 3C). The fraction of gene expression-related GO categories in this module is significantly smaller (22/425, χ2 = 8.05, df = 1, p = 0.00455) than in Module 2.
Dex-responsive transcription regulators
We identified 37 Dex-responsive genes whose products are involved in the regulation of gene expression (Figure 4). 12 of these genes including TFs Klf2, 4 and 9, Per1, Jdp2, Cited2, Nfil3 and Tiparp are upregulated by Dex; 12 others including TFs Junb, Atf3, Tgif1, Irf1 and Bcl3 are downregulated; for 13 regulatory proteins (e.g., Zfp131, Zfp36, Nr4a3, Rela, Nfkb2, Klf7 and Ets1), the inhibitory effect of Dex is apparent only in LPS-induced MФ. For all Dex-induced and a subset of Dex-repressed TFs, we have independently confirmed RNA-seq data by RT-qPCR (Additional file 3: Figure S1; also see Figure 1 for Cited2, Irf1, Etv3 and Atf3).To uncover potential functional interactions between Dex-regulated TFs, we treated BMMФ with Dex for up to 9 h and determined expression levels of a subset of Dex-regulated genes by RT-qPCR. We observed a striking difference in expression patterns over time. Nfil3, Cited2, Jdp2 and Per1 (Figure 5A) are characterized by an accelerated burst phase, with the mRNA level reaching maximum within 30–60 min and then remaining constant (Nfil3 and Cited2) or slowly declining over time (Per1 and Ncoa5). Klf4, Klf9, Tsc22d3 and Ddit4 are strongly induced within the first 2 h, and their RNA levels continue to increase for the next 6 h (Figure 5A). The expression profile of Fkbp5 exhibits an initial delay, followed by a robust and sustained induction (Figure 5A). Conversely, Klf2 and Tiparp displayed pulse-like rapid upregulation within 1–3 h followed by a decline in transcript level, which in the case of Klf2 reaches baseline (Figure 5A); a similar biphasic pattern of expression was observed for Tgfb3, Il15ra and Mt2 (Figure 5A). Interestingly, Bcl3, Junb and Tgif1 responded with rapid pulse-like downregulation followed by a slow return to basal expression level, whereas Atf3 was rapidly downregulated within the first hour and remained repressed throughout the time course (Figure 5A). Unexpectedly, the Pparg expression was only modestly induced by Dex at the early time points, then decreased dramatically by 3 h and remained low for up to 9 h (Figure 5A).
The dynamics of expression for several Dex-regulated TFs suggests that they are under combinatorial controls that involve GR and additional GR targets which either cooperate with or antagonize GR actions. As such a model implies transcription/protein production of these putative GR targets, we first examined the expression of Dex-regulated genes in the presence of a protein synthesis inhibitor cycloheximide (Chx). Treatment with Chx up for to 3 h had no dramatic effects on Dex-mediated regulation of Klf9, Tcs22d3, Tgfb3 and Bcl3 (Figure 5B), suggesting a direct regulation by GR that does not rely on synthesis of additional proteins. Conversely, the expression of Atf3 became refractory to Dex in the presence of Chx (Figure 5B) suggesting that additional proteins induced by Dex rather than GR itself are likely to directly regulate this gene. For several Dex-responsive genes (Klf2, Klf4, Nfil3, Tiparp, Tgif1), however, treatment with Chx dramatically upregulated their basal expression, complicating the assessment of the relative contribution of direct vs. indirect effects of GR to their regulation (Additional file 3: Figure S2) and necessitating an alternative approach.
Temporal dynamics of hormone-regulated gene expression is consistent with feed forward logic
The dynamics of the transcriptional response of several genes to Dex imply the existence of some feedback mechanism that limits activation by GR yet, at the same time, is GR-dependent. Because Chx elicits many off-target effects and does not enable us to discriminate between the secondary targets of GR and those jointly controlled by GR and a GR-regulated TF, we performed dynamic modeling of expression data in an attempt to discern specific regulatory patterns. Several mechanisms including positive and negative autoregulation, positive and negative feedback and feed forward loops (FFL) could account for deviations from a simple model with a single TF regulating gene expression via a single DNA binding site [12]. The kinetics of GR expression following Dex stimulation is not consistent with auto-regulatory models (Figure 5A). Depending on the overall regulatory outputs and activities of individual FFL components, two types of FFL have been recognized – coherent (C-FFL) and incoherent (I-FFL). In type 1 I-FFL (I1-FFL), the activating master TF and a repressing intermediate regulator have opposite effect on a jointly regulated gene. Dynamic modeling and experimental studies of I1-FFL dynamics have demonstrated several properties of this network motif, including its ability to produce sharp pulse-like activation of a jointly regulated gene with a fast relaxation time [11, 13]. Several GR-regulated genes, including Klf2, Tiparp, Tgfb3, Mt2 exhibit pulse-like kinetics at constant Dex exposure (Figure 5A). Near-baseline relaxation of Klf2 expression suggests that this gene is under control of both GR and a strong Dex-induced repressor. Because GR is largely inactive in the absence of ligand, glucocorticoids act as a low-latency on-off switch eliminating the need to correct for a baseline activity of GR. The dynamics of Klf2 and repressor (R) expression is described by the ordinary differential equations (2) and (3) (see Additional file 1) [11, 13].
Assuming equal degradation rates of Klf2 and R, these equations can be solved analytically (equation (1) and [13]) and used to fit the expression data for Klf2. When limited to the early data points (up to 4 h), the expression data fit very well to the predicted expression pattern (Figure 5C, R2 = 0.9225), however, at the later time, when the contribution of degradation rates becomes significant, the quality of fit decreases (R2 = 0.5491). Using parameter estimates derived from equation (1) fitting, we solved equations (2) and (3) numerically. Figure 5C shows a good concordance between the experimental data for Klf2 expression as determined by RT-qPCR and the numerical solution (R2 = 0.8317) that describes the dynamics of the I1-FFL. This result strongly suggests that Klf2 and, potentially, several other GR targets that exhibit similar expression dynamics (Tiparp, Tgfb3 and Mt2) are jointly regulated by GR and a GR-induced repressor via the I1-FFL network motifs.
Numerical solutions of equations 2 and 3 also provide a theoretical prediction of an intermediate repressor dynamics in the GR-R-Klf2 I-FFL. Interestingly, among several known transcription repressors activated by GR, the expression kinetics of Klf9 fits the best to the predicted model (Additional file 3: Figure S3). The dynamics of the GR-R-KLF2 FFL can be tested by perturbing the concentration of a hypothetical intermediate repressor which should uncouple the FFL thus shifting peak-like FFL-mediated kinetics to simple monotonous kinetics eventually converging to a steady-state level. To test the role of Klf9 as a potential GR-activated repressor of Klf2 transcription, we derived M from Klf9-KO mice [30] kindly provided by Dr. Simmen and treated them with 100 nM Dex as above. Interestingly, in Klf9 null BMMФ the peak-like Klf2 induction profile “degenerated” to monotonous activation kinetics that plateaued at a steady state level by 3 h (Figure 5D), replicating the profiles of previously reported uncoupled experimental I-FFLs [11, 31], consistent with the proposed role of KLF9 as an intermediate repressor in the GR-KLF9-KLF2 I-FFL. At the same time, deletion of Klf9 did not affect the expression dynamics of either GR itself or the GR target gene Tsc22d3 with a simple monotonous activation profile (Figure 5D).
GR is recruited to binding sites associated with Dex-regulated genes
The FFL gene regulatory circuitry predicts that the master TF binds DNA to regulate transcription of both FFL nodes. Using several published mouse ChIP-seq datasets of acute GR recruitment [32, 33] we interrogated Dex-induced genes for the presence of GR binding sites within the gene and 15 Kb upstream and downstream of the gene. With the exception of Sox4 and Klf4, all genes encoding Dex-induced TFs contained at least one peak within the analyzed intervals (Figure 6A, Figure 4). To compare these peaks to those in MФ, we analyzed GR recruitment by ChIP-seq using untreated MФ as a control and identified 16,657 peaks induced by a 40-min Dex exposure at 2% FDR. Selective comparison of binding site distributions revealed a high level of concordance between Dex-induced peaks in MФ and those previously described in adipocytes [32] (Figures 6A, 7B and Additional file 3: Figure S4) and a partial overlap with a GR cistrome in MФ polarized with high dose long-term glucocorticoid exposure [34]. By ChIP-qPCR, we detected GR recruitment as early as 40 min post Dex treatment at multiple putative GR binding sites, including those at Per1, Cited2, Klf2, Klf9, Nfil3, Jdp2, Tiparp and Ncoa5 (Figure 6B). These observations correlate well with the expression data (Figures 1 and 5A). Although Klf4 was strongly induced by Dex, no glucocorticoid response elements (GREs) near the gene has been previously reported. We performed a scanning ChIP with evenly spaced primers within the Klf4 gene and several primers flanking potential GR binding sites (Figure 6C). Two of the putative GREs located at -3830 bp (gGcACAgcaTGTaTC) and +5896 bp (aGaACAgaaTGTagttc) relative to the Klf4 transcription start site recruited GR following a 40-min treatment with Dex (Figure 6C), consistent with the notion that, similar to genes shown in Figure 6B, GR is likely to regulate Klf4 directly.
We then used a set of Dex/LPS-regulated genes to analyze the distribution and enrichment of binding sites for TFs that were present in Chip-Enrichment Analysis (ChEA) database [35]. The ChEA database contains curated published genome-wide datasets of TF binding sites in human, mouse and rat. After filtering out TFs that were not expressed in MФ (RPKM < 1), we noted that binding sites for several Dex-responsive TFs, such as KLF2, KLF4, ATF3, EGR1, CEBPβ and IRF1 are enriched among Dex/LPS-regulated genes. Interestingly, binding sites for PPARγ, whose expression was inhibited upon prolonged Dex treatment, were found near the majority (22/37) of Dex-responsive gene expression regulators (Figure 4) and highly enriched among Dex/LPS-regulated genes in general.
We then determined the frequency of genes associated with binding sites for several TFs identified by ChEA and induced by Dex in individual clusters of Dex/LPS-regulated genes (Figure 1). We defined a binding site as ‘gene-associated’ if its genomic intervals overlapped with genomic intervals encompassing all mouse genes annotated in mm9 +/- 15 Kb by at least one nucleotide. In good correlation with the RNA-seq data, acute GR recruitment peaks previously identified by ChIP-seq in Dex-treated adipocytes [32] were enriched in Dex-induced clusters 1, 6, 7 and 9 (Figure 7A, arrow up). Among Dex-regulated TFs, mouse genome-wide binding datasets are currently available for KLF4 (ChIP-seq), KLF2 (chip-on-chip), PPARγ (ChIP-seq) and NFIL3 (chip-on-chip) [36–40]. Although KLF4 sites are enriched in the entire Dex/LPS dataset compared to non-expressing genes, only in cluster 8 (Dex-repressed genes) the enrichment level attains significance (Figure 7A). Conversely, KLF4 binding sites are underrepresented in the LPS-induced cluster 2. Although the KLF2 chip-on-chip dataset available was relatively small, several KLF2 targets were regulated by Dex and contained GR binding sites, including Klf2 itself, Klf4, Klf9 and Tgfb3 [37].
Two mouse PPARγ datasets, one from differentiated adipocytes and one from resting MФ are currently available [38, 39]. Consistent with the previously reported role of PPARγ in repression of inflammatory genes, PPARγ binding sites are overrepresented in several LPS-induced clusters (clusters 2, 3 and 5; Figure 7A) in both datasets. In addition, in the MФ dataset, PPARγ binding sites were enriched in Dex-induced and -repressed clusters, 7 and 8, respectively (Figure 7A).
The only available genome-wide dataset of Nfil3 binding was acquired in a neuronal cell line [40]. Among Nfil3-occupied genes identified in that study, only four genes overlap with Dex/LPS-regulated dataset reported here; however, one of them, Tsc22d3 (GILZ), is a well-characterized GR target.To identify genes that might be under combinatorial control by GR and another Dex-responsive TF, we searched for loci that contained GR, KLF and/or PPARγ binding sites located close to each other within a gene +/- 15 Kb. Intriguingly, several GR targets including Klf2, Klf9, Cited2 and Mt2 contained tight clusters of binding sites for GR, KLF(4) and PPARγ (Figure 7B) suggesting that these TFs may interact functionally or physically at a subset of GR-regulated genes.
LPS downregulates a unique class of genes encoding the C2H2-KRAB gene expression regulators
Even a brief LPS treatment results in a marked downregulation of a large number of genes confined to clusters 10–12. GO overrepresentation analysis of LPS-repressed genes revealed that many of them participate in the regulation of nucleic acid metabolism, gene expression and transcription however, detailed information on specific functions of many of these proteins is lacking. Interestingly, 33 proteins in these clusters contained tandem zinc-finger motifs (COG:5048, p = 3.49*10-29, Figure 8A). We confirmed that the expression of 10 of these Zn-finger proteins is rapidly downregulated by LPS in MФ, and remains suppressed for up to 6 h of treatment (Figure 8B). Further analysis of domain architecture revealed that in the majority of these proteins, tandem Zn-fingers co-occur with domains such as Kruppel-Associated Box (KRAB, Pfam01352, p = 1.029*10-19), BTB (Pfam00651, p = 5.68*10-5) and SCAN (Pfam 02023, p = 2.1*10-4, Figure 8A). KRAB is a Tetrapoda-specific domain that defines one of the largest sub-families of Zn-finger proteins [41] which are involved in nucleic acid binding and regulation of gene expression. Although the specific functions of the majority of KRAB proteins with respect to innate immunity are not well studied, in a few characterized cases KRAB proteins have been associated with transcriptional repression, establishing reversible patterns of histone and DNA methylation and reversible heterochromatization [42–44].