Glucocorticoid receptor coordinates transcription factor-dominated regulatory network in macrophages

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 $Z=\left(X_{\mathit{\text{gene}}}-\stackrel{-}{X}\right)/\sigma$ where X_gene is an expression value for a given gene at a given condition, $\stackrel{-}{X}$ is an average expression value across conditions and σ is a standard deviation of $\stackrel{-}{X}$ 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).

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, χ² = 18.33, p = 0.0001, Figure 2B) or LPS-repressed genes (Module 3, χ² = 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, χ² = 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, R ²  = 0.9225), however, at the later time, when the contribution of degradation rates becomes significant, the quality of fit decreases (R ²  = 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 (R ²  = 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].

Glucocorticoids- and LPS-regulated gene expression programs

GR is a ubiquitous ligand-dependent TF capable of eliciting highly divergent transcription programs with up to a third of protein-coding genes differentially expressed following a 24-h glucocorticoid treatment [7]. Establishing the hierarchy of regulatory events upon prolonged hormonal exposure in individual cell types is challenging, which complicates both accurate mechanistic predictions and clinical decisions. Multiple GR ligands have been designed in an attempt to create a highly specific compound that selectively regulates desired subsets of genes. Mechanistic analyses of these ligands usually focus on a specific group of disease-relevant genes and often involve long-term treatments, which obscure primary and transient responses to GR activation by a plethora of secondary pathways. In the context of inflammation, both immediate and delayed regulatory events are clinically relevant as they reflect typical glucocorticoid treatment modalities. We reasoned that by analyzing early transcriptomes elicited by the inflammatory and glucocorticoid exposure in MФ, a clinically relevant cell type, we will be able to isolate a key set of immediate GR targets responsible for the delayed gene expression patterns. Our results indicate that early glucocorticoid- and LPS-dependent changes establish a highly organized program of gene expression with distinct groups of genes following cooperative and antagonistic regulation. As expected from previous work [22, 34, 45] a large group of LPS-induced genes that included among others inflammatory cytokines, was rapidly downregulated by glucocorticoids. Another group encompassing glucocorticoid-induced genes, some of which encode TFs (Cited2, Nfil3, Jdp2) or signaling proteins (Dusp1, Tsc22d3), are involved in curbing inflammatory signaling [46]. We identified several previously unreported glucocorticoid-induced genes whose products are involved in signaling (Ccl17, Il15ra and Fzd4), regulation of transcription (Klf2, Jdp2, Ncoa5 and Tiparp) and mRNA stability (Zfand5). Several of these genes add to the arsenal of anti-inflammatory mediators regulated by GR. For example, KLF2 interferes with AP1- and NFκB-mediated transcription of Tnf and several chemokines including Ccl2, 3 and 4 [47]. Furthermore, Klf2 haploinsufficiency in mice results in an exaggerated inflammatory response and more severe disease in arthritis models [48]. CCL17, another previously unreported glucocorticoid-induced chemokine, is a marker and promoter of the polarization of ‘alternatively activated’ M2 MФ, which are considered anti-inflammatory and mediate tissue repair and wound healing [49, 50]. In addition to repressing cytokine gene transcription, glucocorticoids downregulate expression of several TFs including Atf3, Junb, Irf1, Bcl3, Tgif1, some (e.g., Rela, Nfkb2, Myc, Ets1) in the context of LPS induction, and, unexpectedly, Pparg An enrichment in positive regulators of inflammation and cell proliferation among Dex-downregulated TFs is consistent with the anti-inflammatory and anti-proliferative effects of glucocorticoids. The role of GR in repression of the Pparg gene in MФ has not been previously reported, the effect might be indirect and mediated by a well-established GR target GILZ [51], which may also account for the delay (Figure 5A). Finally, we described a previously overlooked group of LPS-downregulated genes encoding proteins with the C2H2 Zinc-fingers adjacent to the KRAB domain. Despite being one of the largest TF family, KRAB proteins remain poorly characterized. Among those whose functions were described, several are involved in transcriptional regulation, RNA and DNA binding and splicing [42–44, 52]. KRAB domains interact with a scaffolding co-repressor TRIM28 (KAP1, TIF1β) which in turn binds the heterochromatin protein 1, chromatin remodeler MI2A and H3K9-specific methyltransferase [53]. Indeed, some KRAB proteins reportedly repress transcription by heterochromatin spreading [52]. Interestingly, several KRAB proteins have been linked to NR actions [54, 55]. The role of KRAB proteins in inflammation is essentially unknown; however, genomic studies indicate that inflammatory signaling increases accessibility of large sections of the genome [56]. It is tempting to speculate that a broad downregulation of proteins involved in heterochromatin maintenance and spreading serves to increase DNA accessibility and inflammatory gene transcription.

The dynamic response to GR activation is consistent with feed forward logic

Functional relationships between GR and its targets are often classified as “direct”, that involve GR recruitment to genomic binding sites associated with regulated genes, and “indirect”, whereby primary GR-regulated factors, rather than GR itself, are responsible for activation of the downstream targets. Thus, the activation of these secondary targets is often described as sequential or delayed. Such a model, however, cannot explain many instances of non-monotonous expression dynamics (see Figure 5) and non-linear response to varying hormone concentration of many GRE-driven genes [57]. The large number of shared neighbors, overrepresentation of TFs and their high interconnectivity in GR regulatory networks (Figure 2) are consistent with more intricate regulatory modalities such as FFL. Variations in kinetic parameters for participating TFs, target gene structure and activation/repression thresholds often lead to paradoxical responses to stimulation of the master TF with profound functional implications. C-FFLs serve as delayed response organizers that detect the duration/strength of a signal that activates the initiating TF [11, 12]. Interestingly, the dynamics of Fkbp5 induction by Dex, characterized by a substantial post-exposure delay followed by a robust expression (Figure 5A), is reminiscent of the C-FFL in which the jointly regulated gene is activated by both the master and intermediate TFs [12]. Although additional experiments are required to establish the precise mode of Fkbp5 regulation, this gene is a known direct GR target that recruits GR to several GREs (Additional file 3: Figure S5).

Incoherent loops are responsible for negative and positive pulse generation, accelerated response and fold change sensing [11, 13, 14]. Here, we observed that several GR target genes exhibit both positive (Klf2, Tiparp, Tgfb3 and Mt2) and negative (Tgif1, Junb and Bcl3) pulse-like dynamics consistent with the I-FFL. In keeping with the role of a potential master regulator, GR binds to the GREs in regulatory regions of many of these genes (Figure 6). Furthermore, using a system of ordinary differential equations which describe FFLs in the “fold sensors” model [13], we showed that Klf2 expression is consistent with that of a gene under joint control of GR and a strong GR-activated repressor (Figure 5C). Several GR-activated genes are either known transcription repressors (e.g., Klf4, 9, Nfil3, Per1 and Jdp2) or may downregulate gene expression by destabilizing RNA transcripts (Zfand5, [58]). Curiously, the expression dynamics of Klf9 fits closely with the computational prediction of an intermediate repressor in the GR-R-Klf2 I-FFL (Additional file 3: Figure S3). GR is recruited to the Klf9 and Klf2 GREs as early as 40 min of Dex treatment. Both Klf9 and Klf2 regulatory regions also contain functional GAGGCGTGG KLF sites ([36], Figure 7B) which can be occupied by various TFs of the KLF family [59] including KLF9. Finally, in KLF9-KO macrophages, the induction profile of Klf2 loses the early peak followed by a decrease and acquires monotonous kinetics (Figure 5D) strongly suggesting a collapse of the I1-FFL to simple GR-dependent activation. Interestingly, KLF binding sites are overrepresented in glucocorticoid-regulated genes and are located near GREs in several bona fide GR target genes suggesting that these factors may co-regulate a number of GR targets.

KLF proteins in inflammation

Both KLF2 and KLF4 have been implicated in myeloid cell biology. KLF2 inhibits monocyte activation by inhibiting NFκB activity, which correlates with decreased expression of multiple cytokines and HIF1α, a TF that regulates myeloid cell response to bacterial infection and reactive oxygen species [60]. Consistent with the anti-inflammatory role of KLF2, mice hemizygous for Klf2 have elevated levels of inflammatory mediators, such as CCL2 and PTGS2 (COX-2). By extension, in KLF2-deficient mice, the manifestations of both Me-BSA- and IL1β-experimentally induced inflammatory arthritis are more severe [48].

KLF4 is involved in inflammatory monocyte differentiation [61, 62] and in MФ polarization toward the M2 anti-inflammatory phenotype [63]. KLF9 can act as either a transcriptional activator or a repressor [64, 65], however its role, if any, in inflammation has not been described. We showed here that GR regulates Klf genes with distinct temporal dynamics and proposed that KLF9 may act as a GR-induced Klf2 repressor. Thus, it is tempting to speculate that GR anti-inflammatory activities rely in part on the activation of Klf genes whose products regulate transcription of additional targets in concert with GR. Indeed, glucocorticoids and KLF4 regulate partially overlapping set of genes during epidermal barrier establishment in embryogenesis [66]. The proximity of the GREs and KLF binding sites in the genome suggests an intriguing possibility that GR and KLFs interact functionally or physically. Curiously, a functional interaction with the I-FFL logic has been reported for GR and another member of KLF family, KLF15 [67]. Although KLF15 is not expressed in MФ, our studies strongly suggest extensive crosstalk between GR and other KLF family members in the innate immune cells.

Mice and macrophage cultures

C57BL/6 mice (NCI, Charles River Laboratories) were maintained in the Hospital for Special Surgery Animal Facility in full compliance with institutional guidelines approved by the HSS Animal Care and Use Committee. Klf9-KO mice [30] were generously provided by Dr. R. Simmen (U. of Arkansas). BMMФ were prepared from 8–12 weeks old mice as in [23]. For RNA-seq, BMMФ from two independent mice were treated with vehicle, Dex (100 nM), LPS (10 ng/ml) or LPS + Dex for 1 h. For qPCR and ChIP analyses, BMMФ were treated as above for time indicated in Figure Legends.

RNA isolation, RT-qPCR and RNA-seq

Total RNA was isolated using the RNeasy kit (Qiagen). 0.25 μg of total RNA was used for random primed cDNA synthesis which was performed with M-MuLV reverse transcriptase (NEB) according to the manufacturer’s recommendations. Quantitative PCR (qPCR) was performed using Maxima Sybr Green/ROX/ 2x master mix (Fermentas) on StepOne Plus real time PCR system (ABI) and analyzed using δδCt method as described previously [22] with Hprt or Act1 as a normalization control. Primer pairs are listed in Additional file 2: Table S3. RNA-seq is described in (Additional file 1).

ChIP-qPCR and ChIP-seq

BMMФ were incubated -/+100 nM Dex for 40 min and ChIPs were performed as in [23] using N499 [22] and sc1004 (Santa Cruz Biotechnology) anti-GR antibodies, or normal rabbit IgG as a background control. The data for each recruitment site was normalized to non-specific signals at the unrelated 28S ribosomal gene. Primer pairs are listed in (Additional file 2: Table S3). ChIP-seq is detailed in (Additional file 1).

Gene association network construction, analysis and expression data modeling

GeneMANIA algorithm was used to build gene association networks. Starting from a gene list of interest (e.g., combined list of Dex up- and downregulated genes), geneMANIA algorithm creates a consensus network and predicts gene functions based on integration of multiple prebuilt gene association networks. The detailed description of Gene association network construction, analysis and expression data modeling is in the Supplemental Information section.