Bmal1 regulates inflammatory responses in macrophages by modulating enhancer RNA transcription

To elucidate the genome-wide function of the clock protein Bmal1 in the machinery regulating the inflammatory responses of macrophages, we first examined the effects of Arntl deletion on the transcriptome. Bone marrow-derived macrophages (BMDMs) were prepared from Arntl^−/− mice or their wild-type (WT) littermate controls. The cells were stimulated with Kdo2 lipid A (KLA), a chemically defined substructure of bacterial lipopolysaccharide (LPS) that is specifically recognized by TLR4¹⁹, after which the transcriptome was analyzed over the course of the inflammatory response through RNA-sequencing (RNA-seq). RefSeq genes with normalized counts ≥100 at one or more time points in WT macrophages (a total of 9,720 genes) were clustered into 5 groups using the Mfuzz clustering method²⁰, and the temporal changes in gene expression in each cluster were compared between WT and Arntl^−/− macrophages (Fig. 1). In all clusters, expression of many genes was affected by Arntl deletion, and expression patterns appeared to be desynchronized, particularly in clusters 1–3. Notably, the expression patterns of the KLA-induced genes in clusters 1 and 2 were strikingly disturbed in Arntl^−/− macrophages. In WT cells, genes in clusters 1 and 2 were upregulated after KLA treatment, and the expression of core genes involved in those clusters reached a peak at 6 h. In Arntl^−/− cells, by contrast, the clear induction pattern of those genes seen 6 h after KLA treatment was lost. Instead, expression of many of those genes reached peaks 24 h after KLA treatment, indicating a marked delay in the KLA response. Moreover, whereas expression of genes in cluster 3 sharply dipped after 6 h in WT cells, this nadir was absent in Arntl^−/− macrophages. Gene ontology analysis indicated that genes related to inflammatory responses and cellular metabolism were significantly enriched in clusters 1–3. In particular, the GO terms “innate immune response” and “IκB kinase/NF-κB signaling” were enriched in cluster 1 (Fig. 1). Similarly, the GO terms “NF-κB signaling pathway” and “innate immune system” were enriched in cluster 2. Cluster 3 was also associated with innate immune responses. These findings suggest that the innate immune responses mediated via NF-κB signaling were disturbed in Arntl^−/− macrophages, due to the delayed and/or desynchronized transcriptome response to TLR4 activation.

Expression of genes belonging to clusters 4 and 5 were also disturbed by Arntl deletion; however, the trends in the gene expression of the core members appeared to be less affected by Arntl deletion than in clusters 1–3. In addition, clusters 1–3 contained genes regulated by NF-κB signaling. Given these results, we will focus on the inflammatory responsive genes under the control of TLR4 signaling in the following analyses.

The results of our RNA-seq clearly indicate that TLR4-responsive transcriptome activation is delayed or desynchronized in Arntl^−/− macrophages. Based on these results, and given the role of Bmal1 as a central transcription factor within the cellular clock, we hypothesized that Bmal1 contributes to the time-dependent regulation of inflammation-related gene expression in macrophages. To test that idea, we first assessed the effect of TLR4 activation by KLA on the expression of clock genes. Expression of Arnt1 and Per1/2 was transiently downregulated 4 h after KLA treatment (Supplementary Fig. S1). Cry1/2 was upregulated after KLA treatment; its expression showed a peak at 8 h. The effects of TLR4 activation on Arntl expression are consistent with earlier reports^17,18 and demonstrate that TLR4 activation by KLA induces time-dependent changes in clock gene expression.

To characterize the effect of Arntl deletion on induction of inflammatory response-related genes, we further investigated the genes classified as cluster 1 (1,758 genes). The temporal expression profiles of these genes were evaluated by comparing the expression levels of individual genes at 0, 6 and 24 h after KLA treatment in WT and Arntl^−/− macrophages (Fig. 2a and b, Supplementary Fig. S2). While expression of the majority of cluster 1 genes was higher at 6 h than at baseline (0 h) in WT cells, expression of these genes remained unchanged in Arntl^−/− cells (Fig. 2a). Many cluster 1 genes in WT cells were slightly inhibited or unchanged at 24 h, as compared to their expression levels at 6 h. In Arntl^−/− cells, by contrast, expression levels of the majority of these genes were increased at 24 h as compared with those at 6 h (Fig. 2b). These changes were also observed by comparing gene expression levels in WT and Arntl^−/− at each time point. Expression of the majority of cluster 1 genes was higher in Arntl^−/− cells at 0 h, lower at 6 h, and higher at 24 h than in WT cells (Supplementary Fig. S2). Collectively, expression of cluster 1 genes in Arntl^−/− cells showed a trend toward basal activation and delayed and/or enhanced expression later, at 24 h.

Quantitative PCR (qPCR) analysis independently confirmed that expression of proinflammatory genes belonging to cluster 1, including Nos2 and Il1b, remained significantly elevated 24, 48 and 72 h after TLR4 activation in Arntl^−/− cells as compared to WT cells (Fig. 2c). Thus activation of these genes in Arntl^−/− macrophages was both enhanced and prolonged. Consistent with the enhanced Il1b mRNA expression, levels of secreted IL-1β protein were significantly increased in the Arntl^−/− BMDMs as compared to WT cells (Fig. 2d).

As in cluster 1, the upregulation of core members of cluster 2 genes seen at 6 h in WT cells was absent in Arntl^−/− cells (Fig. 1). For instance, our qPCR showed that in Arntl^−/− cells, expression of Csf1r remained low 6 h after LPS treatment but was higher after 24 h, suggesting induction of Csf1r was delayed and sustained in Arntl^−/− cells (Fig. 2d). qPCR also showed that expression of Hif1a, belonging to the cluster 3, was enhanced at most time points in Arntl^−/− cells. Collectively, these findings indicate that TLR4-mediated expression of a number of proinflammatory genes is enhanced and/or sustained in Arntl^−/− macrophages.

We next analyzed the mechanism by which Bmal1 affects the temporal proinflammatory gene response to TLR4 activation. To that end, ChIP-seq of Bmal1 was performed in mouse primary macrophages (thioglycolate-elicited peritoneal macrophages). As expected, Bmal1 binding to clock gene loci was observed (data not shown). Interestingly, Bmal1 binding was also observed in the regulatory regions of TLR4-inducible genes exemplified by Hif1a (cluster 1) and Csf1r (cluster 2) (Fig. 3a). In addition, de novo motif analysis identified a motif for basic helix-loop-helix (bHLH) transcription factors, including NPAS2, as the most highly enriched sequence followed by motifs for PU.1, a macrophage lineage-determining transcription factor, and AP-1 and C/EBP, two transcription factors known to be involved in macrophage activation (Fig. 3b). This suggests Bmal1 associates with enhancers bound by transcription factors important for the regulation of macrophage function, including PU.1, AP-1 and C/EBP.

Previous studies demonstrated that PU.1 is necessary for establishing macrophage-specific cistromes for signal-responsive transcription factors, including NF-κB^21,22. We found that deletion of Arntl affected expression of NF-κB-regulated genes (Fig. 1). Based on those observations, we hypothesized that Bmal1 may co-bind to its cognitive sites along with NF-κB and/or PU.1. Indeed, approximately half of the 2,026 Bmal1 peaks were co-occupied by PU.1 (1,040 sites; 51%). Moreover, 729 sites (36%) were co-occupied by p65, and 611 sites (30%) were co-occupied by both p65 and PU.1 (Fig. 3c). These results suggest that Bmal1 may directly affect PU.1-containing enhancers after TLR4 activation. Given that PU.1 acts as a macrophage lineage-determining transcription factor that is essential for establishment of macrophage-specific enhancers^13,23, the extensive co-binding of Bmal1 with PU.1 suggests Bmal1 is involved in macrophage-specific regulatory programs, such as inflammatory responses. Consistent with that idea, only a small fraction (5.5%) of macrophage Bmal1 binding peaks were shared between macrophages and the liver²⁴ (Fig. 3d), suggesting the targets regulated by Bmal1 differ between macrophages and the liver. In macrophages, moreover, the target genes containing Bmal1 binding sites predicted by GREAT (genomic regions enrichment of annotations tool)²⁵ were significantly associated with different GO terms than liver, with the exception of circadian rhythm-related terms (Supplementary Fig. S3). It thus appears that, aside circadian regulation, Bmal1 has different regulatory functions in macrophages than in the liver.

Given the role of NF-κB as a primary driver of TLR4-mediated inflammatory transcriptional activation, we hypothesized that perturbed NF-κB binding, namely delayed and prolonged association of NF-κB in response to TLR4 activation, is the mechanism responsible for the altered temporal proinflammatory gene expression in Arntl^−/− macrophages. To test that idea, we performed ChIP-seq of the p65 component of NF-κB in WT and Arntl^−/− macrophages, after which the temporal pattern of p65 binding was analyzed in WT and Arntl^−/− macrophages. Consistent with previous studies^21,26, the average p65 binding was increased at PU.1-binding sites 4 h and 8 h after KLA treatment in WT macrophages, and was then lower at 12 h and later (Fig. 4, left column). Unexpectedly, however, Arntl deletion did not affect average p65 association with PU.1 loci (Fig. 4, left column).

The unexpected finding that p65 recruitment was unaffected by loss of Arntl prompted us to analyze the enhancer landscape. To assess global enhancer activity, we performed ChIP-seq of lysine 27-acetylated histone H3 (H3K27ac), which is positively associated with transcriptional activity^27,28. Interestingly, average H3K27 acetylation levels at all PU.1 binding sites were higher in Arntl^−/− than WT macrophages at baseline and throughout the course of the response after KLA stimulation (Fig. 4, right column).

To further assess epigenetic changes, we first analyzed the density of H3K27ac in PU.1/Bmal1 co-binding regions to analyze possible direct effects of the loss of Bmal1. ChIP-seq tags for Bmal1, PU.1 and H3K27ac around the top 500 Bmal1 binding sites at baseline, were grouped into 5 clusters using the K-means clustering algorism. A heat map of Bmal1, PU.1 and H3K27ac reads around PU.1 peaks shows that while H3K27ac levels in cluster A were unchanged or decreased in Arntl^−/− cells, H3K27ac levels in clusters B, C and D were clearly higher in Arntl^−/− cells than WT cells under basal conditions (Fig. 5a). In agreement with the results of the top 500 Bmal1 sites, average H3K27ac levels around all Bmal1-PU.1 co-binding sites were increased in Arntl^−/− cells as compared to WT cells (Supplementary Fig. S4). The changes are exemplified by the KLA-inducible Arntl target genes Hif1a and Csf1r (Fig. 5b). These results indicate that loss of Bmal1 leads to increased H3K27 acetylation at most enhancers to which Bmal1 and PU.1 bind.

Although Arntl deletion globally affected the H3K27ac states at PU.1 sites (Fig. 4 right column), only a minor fraction of PU.1 sites were cobound by Bmal1 (Fig. 3c). Therefore, the differential deposition of H3K27 acetylation at PU.1-contining enhancers cannot be fully explained as a primary effect of the loss of Bmal1 binding. Accordingly, we further characterized the mechanism underlying the changes in H3K27 acetylation. TLR4 activation can both activate and repress PU.1-bound enhancer activity concomitantly with changes in H3K27 acetylation²⁶. To dissect the epigenetic effect of the loss of Bmal1 in line with our focus on KLA-activated genes in the present study, we focused on the PU.1 enhancers at which H3K27ac was increased by KLA. Furthermore, to focus on early KLA-activated enhancers, we analyzed the PU.1 sites at which H3K27ac tags were increased ≥2-fold as compared to the steady state 4 h after KLA stimulation (5,572 peaks). In WT cells, H3K27ac levels at these KLA-activated enhancers were increased at 4 h and then gradually declined (Fig. 6a). By contrast, the average deposition of H3K27ac was already high at baseline in Arntl^−/− macrophages. And while H3K27ac levels were comparable between the two genotypes 4 h after KLA, the elevation in H3K27ac persisted beyond 4 h in Arntl^−/− cells. These altered H3K27ac states are exemplified by the Irf1 locus (Fig. 6b). Moreover, the persistent activation of KLA-responsive enhancers appears to be in agreement with the delayed and prolonged mRNA response to KLA treatment for many genes in clusters 1–3 (Fig. 1). Collectively then, these findings suggest the loss of Arntl induces an increase in baseline enhancer activity at inflammatory genes, reduces the responsiveness to TLR4 activation, and prolongs inflammatory responses.

To elucidate the mechanism by which Bmal1 modulates H3K27 acetylation globally on PU.1 sites, we hypothesized that clock-controlled transcription factors are involved. Of particular interest were RevErb-α and RevErb-β, encoded by Nr1d1 and Nr1d2, respectively, because these nuclear receptors have been shown to regulate expression of genes involved in the control of metabolism and inflammatory responses in macrophages^10,29. In addition, Nr1d1 and Nr1d2 are known targets of Bmal1³⁰. Consistent with their involvement, our ChIP-seq revealed several Bmal1 binding peaks in the regulatory regions of Nr1d1 and Nr1d2 loci (Fig. 7a), and levels of Nr1d1 and Nr1d2 expression were markedly reduced in Arntl^−/− macrophages throughout the course after KLA stimulation (Fig. 7b).

RevErbs reportedly repress enhancer activity by inhibiting eRNA transcription²⁹. To test whether dysregulation of eRNA transcription modulates enhancer activity in Arntl^−/− cells, we took advantage of previously reported global run-on sequencing data (GRO-seq) from Nr1d1/Nr1d2 double knockout (DKO) macrophages²⁹. To identify eRNA-expressing enhancers, we first identified regions with significant nascent RNA signals determined with GRO-seq in WT or DKO cells. After signals detected at gene bodies were excluded, the nascent RNA signals with read counts increased ≥1.5-fold in Nr1d1/Nr1d2 DKO cells as compared to WT cells were identified. The resultant 575 regions were likely to include the regions in which eRNA expression was upregulated by DKO. Further, we identified 330 regions with PU.1-binding among the eRNA-expressing enhancers. Interestingly, the acetylation status of H3K27 around those eRNA-expressing, PU.1-binding enhancers was markedly higher in Arntl^−/− than WT macrophages (Fig. 8a), suggesting the reduced RevErb-mediated suppression of eRNA de-repressed the enhancers. Our qPCR results confirmed that the eRNA expression located at −5 kb in Mmp9 and +28 kb in Cx3cr1, both of which are targets of RevErb-mediated repression²⁹, was significantly increased in Arntl^−/− macrophages (Fig. 8b). And H3K27ac was also increased at the PU.1-containing enhancers (Fig. 8c). Collectively, these results suggest that epigenetic dysregulation of the response to inflammatory TLR4 activation in Arntl^−/− cells is at least partly mediated by dysregulation of eRNA transcription due to impaired RevErb function.

Arntl^−/− mice were generated on a C57BL/6 background as described previously³⁹. Mouse thioglycollate-elicited macrophages and bone marrow-derived macrophages were isolated from male 6- to 9-week-old C57BL/6J (Charles River laboratories) and Arntl^−/− mice. All animal experiments were approved by the Experimental Animal Care and Use Committee of Tokyo Medical and Dental University and strictly adhered to the guidelines for animal experiments of the institute.

We used BMDMs for all experiments except ChIP-seq for Bmal1. For culture of BMDMs, bone marrow was collected from mice by perfusing the medullary cavities of femurs, tibias and iliac bones, after which the bone marrow cells were cultured for 6 days in RPMI-1640 medium containing 10% FCS and 20 μg/ml M-CSF (R&D). Peritoneal macrophages were harvested by lavage 3 days after intraperitoneal injection of 3 ml of 3% thioglycollate medium (http://www.lipidmaps.org/protocols/↗). The cells were then cultured overnight and adherence selected.

Total RNA was isolated from cells, purified using RNeasy columns and treated with RNase-free DNase according to the manufacturer’s instructions (Qiagen). For quantitative PCR (qPCR) analysis, 1 nM of cDNA or ChIP DNA was used for real-time PCR with gene- or locus-specific primers. qPCR was performed on an Applied Biosystems StepOne Plus system using SYBR GreenER master mix (Invitrogen) and the following protocol: incubation for 10 min at 50 °C and 10 min at 95 °C followed by 40 cycles of 15 s at 95 °C and 30 s at 60 °C. Primer sequences for RT-qPCR analysis are listed in Supplementary Table. S1

BMDMs were stimulated with KLA (100 ng/ml) for 4 h, followed by ATP (5 mM/ml) treatment for 16 h. Levels of IL-1β in culture media were quantified using a mouse IL-1β ELISA kit (R&D, #MLB00C) according to the manufacturer’s instructions.

For ChIP assays, anti-NF-κB p65 (sc-372) antibody was purchased from Santa Cruz Biotechnology; anti-Bmal1 (ab93806) and anti-H3K27Ac (ab4729) antibodies were from Abcam. For Bmal1 ChIP, 2 × 10⁷ thioglycollate-elicited primary peritoneal macrophages were used. For NF-κB p65 and H3K27ac ChIP, 2 × 10⁷ bone marrow-derived macrophage were used. After crosslinking in 1% formaldehyde, cells were resuspended in buffer containing 10 mM HEPES/KOH (pH 7.9), 85 mM KCl, 1 mM EDTA, 0.5% IGEPAL CA-630 and protease inhibitors and incubated for 5 min. The cells were then spun down and resuspended in 500 μl of lysis buffer (50 mM Tris/HCl [pH 7.4], 1% SDS, 0.5% Empigen BB, 10 mM EDTA) with protease inhibitors, and the chromatin samples were sheared by sonication. The lysates were diluted with 750 μl of dilution buffer (20 mM Tris/HCl, 100 mM NaCl, 0.5% Triton X-100, 2 mM EDTA) and 1% was used as input DNA. With the remaining samples, immunoprecipitation was carried out overnight using Dynabeads protein G coated with specific antibody. The beads were then washed twice with wash buffer I (20 mM Tris/HCl, 150 mM NaCl, 0.1% SDS, 1% Triton X-100, 2 mM EDTA), twice with wash buffer II (10 mM Tris/HCl, 250 mM LiCl, 1% IGEPAL CA-630, 0.7% Na-deoxycholate, 1 mM EDTA), once with TE plus 0.2% triton X-100, and once with TE plus 50 mM NaCl. The beads were then eluted with elution buffer (TE, 2% SDS), after which the DNA crosslinking was reversed, and the DNA was purified using a PCR purification kit (Qiagen) according to the manufacturer’s instructions.

ChIP-seq libraries were prepared using a NEBNext Ultra DNA library prep kit for Illumina according to the manufacturer’s instructions. RNA-seq libraries were prepared using a NEBNext Ultra RNA library prep kit for Illumina according to the manufacturer’s instructions. Libraries were PCR-amplified for 12–15 cycles, and sequenced on a GAIIx, HiSeq. 1500 or HiSeq. 2500 (Illumina).

For macrophage Bmal1, p65 and PU.1 ChIP-seq, reads were mapped to the mm9 mouse genome using STAR (Dobin et al., 2013). For liver, Bmal1 reads were mapped to the mm9 using bowtie⁴⁴. Peak calling and annotation were performed using HOMER²¹. Peaks that overlapped with blacklisted regions⁴⁵ and simple repeat regions were removed. For ChIP-seq analyses, Homer-identified peaks whose scores were ≥20 were selected as high-confidence peaks, except liver Bmal1 peaks, which were used irrespective of scores in downstream analyses. For macrophage Bmal1 ChIP-seq, Bmal1 peaks were pooled from 3 experimental conditions (cells treated with KLA for 0, 4 and 24 h). The data used for liver Bmal1 ChIP-seq analysis was from ZT8 cells, as that time point had the largest the number of identified peaks in a previous report³². For GRO-seq analyses, de novo transcripts were identified using Homer’s findPeaks program with the -groseq option. De novo transcripts that overlapped with blacklisted regions and simple repeat regions were removed. Gene ontology analyses were performed using Metascape software (http://metascape.org/↗)⁴⁶. Motif analyses were performed using HOMER. PU.1 Publically available ChIP-seq, RevErb ChIP-seq and GRO-seq results were downloaded from GEO (GSE62826↗, GSE45914↗, and GSE39860↗). All RNA-seq and ChIP-seq data are available in the GEO under the accession number GSE95712↗.

Differences among more than two groups were analyzed using one-way ANOVA followed by Tukey-Kramer post-hoc tests. Values of p < 0.05 were considered statistically significant, except where otherwise noted. All data are shown as means ± SEM.

The authors declare that they have no competing interests.