Open Access

Genome-wide profiling of transcribed enhancers during macrophage activation

Epigenetics & Chromatin201710:50

https://doi.org/10.1186/s13072-017-0158-9

Received: 10 September 2017

Accepted: 13 October 2017

Published: 23 October 2017

Abstract

Background

Macrophages are sentinel cells essential for tissue homeostasis and host defence. Owing to their plasticity, macrophages acquire a range of functional phenotypes in response to microenvironmental stimuli, of which M(IFN-γ) and M(IL-4/IL-13) are well known for their opposing pro- and anti-inflammatory roles. Enhancers have emerged as regulatory DNA elements crucial for transcriptional activation of gene expression.

Results

Using cap analysis of gene expression and epigenetic data, we identify on large-scale transcribed enhancers in bone marrow-derived mouse macrophages, their time kinetics, and target protein-coding genes. We observe an increase in target gene expression, concomitant with increasing numbers of associated enhancers, and find that genes associated with many enhancers show a shift towards stronger enrichment for macrophage-specific biological processes. We infer enhancers that drive transcriptional responses of genes upon M(IFN-γ) and M(IL-4/IL-13) macrophage activation and demonstrate stimuli specificity of regulatory associations. Finally, we show that enhancer regions are enriched for binding sites of inflammation-related transcription factors, suggesting a link between stimuli response and enhancer transcriptional control.

Conclusions

Our study provides new insights into genome-wide enhancer-mediated transcriptional control of macrophage genes, including those implicated in macrophage activation, and offers a detailed genome-wide catalogue of transcribed enhancers in bone marrow-derived mouse macrophages.

Keywords

Transcriptional enhancerseRNATranscriptional regulationMacrophage activation

Background

Macrophages are innate immune system sentinel cells that mediate homeostatic and protective functions, including host defence against invading pathogens [1]. Macrophages respond to a wide range of external stimuli by acquiring heterogeneous activation states that exert functional programmes tailored for specific microenvironments [2]. A spectrum of macrophage phenotypes has been observed, with macrophages activated in response to interferon-γ, M(IFN-γ), and interleukin-4/interleukin-13, M(IL-4/IL-13), representing two extreme states [3, 4].

M(IFN-γ), often referred to as classically activated macrophages, is pro-inflammatory macrophages characterized by efficient antigen presentation, high bactericidal activity and production of pro-inflammatory cytokines, reactive oxygen and nitrogen intermediates [5, 6]. M(IL-4/IL-13), often classified as alternatively activated macrophages, is predominantly regulatory macrophages involved in homeostasis, angiogenesis, wound healing, tissue remodelling and parasitic and bacterial infection [1, 2, 79]. M(IL-4/IL-13) macrophages release anti-inflammatory cytokines and show less efficient antigen presentation and decreased production of pro-inflammatory cytokines [2, 7]. Macrophage activation is driven by specific transcriptional changes and is controlled by complex cellular mechanisms [6, 10].

Imbalance in populations of macrophages with opposing pro- and anti-inflammatory roles has been implicated in disease progression [1]. Intracellular pathogen Mycobacterium tuberculosis, the causative agent of tuberculosis, interferes with classical activation of macrophages to avoid its antibacterial action and promotes alternative activation state [11, 12]. Tumour microenvironments promote phenotypic switches from pro- to anti-inflammatory macrophages, which might contribute to the tumour progression by inhibiting immune responses to tumour antigens [1, 2]. Conversely, the phenotypic switch from anti- to pro-inflammatory population of macrophages might contribute to obesity and metabolic syndrome [1, 2, 13]. Therefore, the development of techniques for manipulation and specific targeting of macrophage populations could ultimately improve diagnosis and treatment of inflammatory diseases [2]. To advance this area of research, the cellular mechanisms responsible for macrophage activation need to be further deciphered.

Gene expression in eukaryotic cells is a complex process guided by a multitude of mechanisms [14]. Precise regulation is required to ensure dynamic control of tissue-specific gene expression and to fine tune the responses to external stimuli [15]. One such level of control is facilitated via regulation of RNA transcription. This process is mediated by a complex transcriptional machinery with its components recognizing specific regulatory regions of DNA. Promoters represent a better-characterized class of such regions from which RNA transcription is initiated [16, 17]. They act in concert with other cis-regulatory DNA elements, including enhancers, which are believed to play key roles in transcriptional regulation [18].

Enhancers are defined as regulatory DNA regions that activate transcription of target genes in a distance- and orientation-independent manner [18]. According to the dominant model, transcriptional regulation by enhancers is exerted via direct physical interaction between enhancer and target gene promoter mediated by DNA looping [18, 19]. Recent identification of distinct properties of enhancer regions enabled novel approaches to enhancer profiling [18]. Enhancer regions are often distinguished by a specific combination of chromatin marks present at these locations, such as H3K4me1 and H3K27ac [20, 21]. Enhancer sequences contain transcription factor binding sites (TFBSs) that recruit transcription factors (TFs) to regulate target genes [22, 23]. In addition, enhancers are frequently bound by proteins such as histone acetyltransferase p300 and insulator-binding protein CTCF [21, 2426]. Large-scale profiling of these enhancer-associated signatures by chromatin immunoprecipitation followed by sequencing (ChIP-seq) [26, 27] has greatly advanced enhancer identification and enabled systematic and genome-wide enhancer mapping [28, 29]. Another group of methods such as chromosome conformation capture (3C) [30] and its variant Hi-C [31] has been employed to profile physical DNA contacts, including those between promoters and enhancers [32, 33]. However, none of these methods has become a gold standard of enhancer detection, and the field is still actively developing.

Recent studies have led to the unexpected finding that most active enhancers recruit RNA polymerase II and are bidirectionally and divergently transcribed to produce RNA transcripts, referred to as eRNAs [34, 35]. While the functionality of eRNA remains controversial, a recent study by Hon et al. showed that many enhancers are transcribed into potentially functional long-noncoding RNAs (lncRNAs) playing a role in inflammation and immunity [36, 37]. Recently, quantification of eRNA transcription laid the foundation for a novel method of large-scale enhancer profiling [38]. In their study utilizing cap analysis of gene expression (CAGE) [39], Andersson et al. [38] performed genome-wide mapping of transcriptional events followed by identification of enhancers based on co-occurrence of closely located divergent transcripts representing eRNAs. The capacity of CAGE to simultaneously profile the expression of eRNAs and genes became an additional advantage, since eRNA production was shown to positively correlate with the production of mRNAs of target genes [34, 40].

These and other studies unravelled the fundamental importance of enhancer regions as DNA regulatory elements in multiple cell types, including macrophages [28, 29, 35, 40, 41]. Enhancers are extremely widespread, with an estimation of up to one million enhancers in mammalian genomes [20, 23, 24, 42]. They are major determinants of gene expression programmes required for establishing cell-type specificity and mediating response to extracellular signals [23, 43, 44]. Our current understanding of these elements, however, remains incomplete. High tissue specificity of enhancers is a major hurdle towards establishing a comprehensive catalogue of the full enhancer population [23, 43]. Moreover, emerging evidence indicates that enhancers selectively act in a stimuli- or condition-specific manner [45, 46]. A major challenge is, therefore, to catalogue enhancers active in different tissues and conditions and link them to target genes.

Recently, we investigated the transcriptional regulatory dynamics of protein-coding and lncRNA genes during M(IFN-γ) and M(IL-4/IL-13) macrophage activation using CAGE data [10]. We showed that particular TFs, such as Nfκb1, Rel, Rela, Irf1, and Irf2, drive macrophage activation and are commonly activated but have distinct dynamics in M(IFN-γ) and M(IL-4/IL-13) macrophages [10]. Here, we extended the former study to understand the regulatory influence of enhancers in the macrophage activation process. Our genome-wide in silico study aimed at characterizing the transcribed enhancer landscape in mouse bone marrow-derived macrophages (BMDM) and studying its dynamic changes during M(IFN-γ) and M(IL-4/IL-13) activation. We used CAGE data and enhancer-associated chromatin signature to identify transcribed enhancer regions. We inferred regulatory associations between enhancers and target protein-coding genes using their spatial organization in topologically associating domains (TADs) [47] and correlation of CAGE-derived expression in our time course. With these data, we established a catalogue of transcribed enhancer regions linked to their target genes. This catalogue provides insights into genome-wide enhancer-mediated regulation of transcription in mouse BMDM. Furthermore, we highlight the role of enhancers during macrophage activation and report enhancers driving expression dynamics of known macrophage activation marker genes.

Methods

CAGE data and processing

Mouse genomic coordinates (mm10) and tag counts of CAGE transcription start sites (TSSs) were obtained from the FANTOM5 project [17] data repository [48]. Data for 969 mouse samples classified as “primary_cell”, “time course”, “tissue”, and “cell_line” were used. The set included 184 BMDM samples profiled by us as described elsewhere [10], which we used here to construct a macrophage enhancer–promoter interactome (see Additional file 1: Table S1 for the list of macrophage samples).

The DPI programme [49] was used as described in Forrest et al. [17] to cluster CAGE TSSs into CAGE peaks. Briefly, the algorithm uses independent component analysis to decompose regions with continuous CAGE signals into separate peaks based on their profile across different samples and tissues. With the default parameters, similarly to Forrest et al. [17], we obtained a list of all CAGE peaks and a subset of CAGE peaks enriched for promoter-associated signals. The latter file represents a subset of peaks meeting the FANTOM5 “robust” criteria, with a single TSS supported by 11 or more observations and one or more tags per million (TPM) in at least one experiment [17]. These two peak sets were used for identification of enhancers and annotation of protein-coding gene promoters, respectively. Tag counts of all TSSs composing a CAGE peak were summed up to derive a total tag count for that CAGE peak.

Annotation of protein-coding gene promoters

The set of “robust” CAGE peaks derived by DPI (see above) was used to annotate promoters of protein-coding genes. Ensembl gene models [50] version 75 downloaded from the UCSC Table Browser [51] on 11 August 2016 was used to obtain coordinates of protein-coding transcripts and genes. We assigned a CAGE peak to an Ensembl protein-coding transcript if its 5′ end was mapped within 500 bp of the 5′ end of the transcript on the same strand. The transcript annotation was extended to gene annotation by combining the CAGE peaks associated with all of the gene’s transcripts.

Calculation of gene and promoter expression

TMM normalization [52] of promoter tag counts was performed to derive normalized expression values in a form of tags per million (TPM). We excluded lowly expressed promoters from the analysis and retained only those with expression of at least one TPM in 10% of the macrophage samples. Expression of each gene was derived as a sum of expression of the gene’s promoters. The resulting set included 24,449 promoters of 10,767 protein-coding genes.

Identification of mouse enhancers with CAGE data

The full set of 3,188,801 DPI-derived CAGE peaks was used for identification of mouse enhancers. CAGE peaks located within 500 bp of protein-coding transcript start sites or within 200 bp of exons were excluded based on the Ensembl gene models [50] version 75. This filtering resulted in 1,890,465 CAGE peaks. Next, we used the same strategy as Andersson et al. to infer enhancer regions as clusters of closely located bidirectional divergent CAGE peaks and to derive the corresponding tag counts [38]. The resulting 42,470 regions were deemed mouse enhancer regions. The counts were normalized to tags per million (TPM) using TMM normalization procedure [52]. Enhancers with nonzero expression in at least 10% of our macrophage samples were deemed transcribed in our macrophage samples.

Selection of enhancers regulating protein-coding genes in macrophages

Enhancer-specific chromatin signatures were based on ChIP-seq profiling of H3K4me1 and H3K27ac histone marks and were obtained from a study by Ostuni et al. [53]. Transcribed enhancers with at least 1 bp overlap with the regions inferred by Ostuni et al. [53] were retained. Genomic coordinates of TADs in mouse embryonic stem cells were obtained from a study by Dixon et al. [47]. We selected pairs of enhancers and promoters where both features were entirely located within the same TAD. For each of these pairs, we calculated Spearman’s correlation coefficient between expression levels of enhancer eRNA and promoter across our macrophage samples and selected only pairs with positive correlation coefficient and FDR < 10−4 (Benjamini–Hochberg [54] procedure). We considered an enhancer to regulate a gene if it was associated with at least one of the gene’s promoters. All mm9 genomic coordinates were converted to mm10 using the liftOver program [55].

Gene set enrichment analysis

KEGG pathway maps [56] or GO biological process ontology [57] was used as sets of biological terms for GSEA. GO terms and associated genes were retrieved using the R package GO.db [58]. We used hypergeometric distribution to calculate the probability of obtaining the same or larger overlap between a gene set of interest and each biological term [59]. Derived p values were corrected for multiple testing using Benjamini–Hochberg procedure [54]. As a background, a set of 22,543 Ensembl protein-coding genes (version 75) was used [50].

Identification of macrophage-specific features

Normalized TPM expression data were used to calculate a z-score for each of our 184 macrophage samples for each enhancer and gene by subtracting the mean and dividing by the standard deviation of expression values of the same feature in 744 FANTOM5 non-macrophage mouse samples (Additional file 1: Table S2), similarly to Yao et al. [60]. Enhancers and genes with z-score > 3 (i.e. expressed more than 3 standard deviations above the mean of the non-macrophage samples) in at least 10% of macrophage samples were deemed macrophage specific.

TFBS over-representation analysis

TFBS data for mouse were obtained from ENCODE [61] and HT-ChIP [62]. Raw sequencing data were mapped to the mm10 genome build for each tissue and cell type separately, and peaks were called using MACS2 [63]. TFBS summits with FDR < 10−4 were retained. We used three different background sets: the whole set of identified mouse enhancers, the subset of enhancers not expressed in macrophages, and a set of random genomic regions located within TADs excluding gaps, repeated sequences, Ensembl coding regions, and mouse enhancers identified here. Gap and repeated sequence regions were obtained from the UCSC Table Browser [51] on 1 August 2016 (“gap” and “rmsk” tables of mm10 database). Significantly over-represented TFBSs were selected based on empirical p value < 0.01 from a Monte Carlo analysis of 1000 trials [64]. We retained only TFBS which showed p value < 0.01 using all three background sets and nonzero expression of the corresponding TF in our macrophages samples.

Identification of stimuli-responsive features

We calculated a z-score for each of 16 M(IFN-γ) and 16 M(IL-4/IL-13) macrophage samples for each enhancer and gene by subtracting the mean and dividing by the standard deviation of expression values of the same feature in ten non-stimulated macrophage samples, similarly to the approach for identification of macrophage-specific features. Genes and enhancers with z-score > 3 in more than 25% of the corresponding samples were deemed stimuli responsive. Of our associations between stimuli-responsive enhancers and genes (significant with positive correlation coefficient and FDR < 10−4), we selected those with a positive Spearman’s correlation coefficient of expression in the subset of either M(IFN-γ) or M(IL-4/IL-13) samples.

Identification of marker enhancers

We aimed to infer potential activation marker enhancers that regulate marker genes specifically during either M(IFN-γ) or M(IL-4/IL-13) activation. To identify marker enhancers in M(IFN-γ), we first selected enhancers which were deemed M(IFN-γ) responsive, but not M(IL-4/IL-13) responsive and were associated with known M(IFN-γ) marker genes. Second, a z-score for each of 16 M(IFN-γ) samples was calculated using 16 M(IL-4/IL-13) samples as a background. Enhancers with z-score > 3 in more than 25% of M(IFN-γ) samples were deemed potential activation marker enhancers in M(IFN-γ). The same strategy was used to infer activation marker enhancers in M(IL-4/IL-13) macrophages.

Differential expression analysis

Differential gene expression analyses were performed using the exact test implemented in edgeR [52], and the p values were adjusted for multiple hypothesis testing using the Benjamini–Hochberg procedure [54].

All analyses made extensive use of the BEDTools utilities [65] and the R software [66].

Results

Identification of transcribed mouse macrophage enhancers

Active enhancers were shown to be bidirectionally transcribed in mammals [34, 35], and eRNAs profiled by CAGE technology [39] were used before to reliably infer enhancer regions in human [38]. To identify transcribed enhancers in mouse tissue, we used the FANTOM5 collection of CAGE mouse samples [17] and a strategy developed before [38] (see “Methods”). This approach yielded 42,470 mouse enhancers, with 17,752 enhancers deemed transcribed in our BMDM samples (Fig. 1a, “Methods”). Ostuni et al. [53] defined enhancers in mouse macrophages based on two histone marks (H3K4me1 and H3K27ac) profiled by ChIP-seq before and after different types of macrophage stimulation. To refine our set, we sub-selected 11,216 (63%) CAGE-based transcribed enhancers that overlapped the ChIP-seq-based enhancer regions (Fig. 1a, “Methods”, Additional file 2). Of these 11,216 enhancers, only 6.4% overlap ChIP-seq-based enhancers that carry histone marks of poised not activated enhancers (only H3K4me1, see Additional file 2: Table S13), whereas the rest also carry active enhancer histone mark H3K27ac either before or after macrophage stimulation. The remaining 6536 CAGE-based enhancers did not overlap regions inferred by Ostuni et al. [53] and were excluded from further analysis. Notably, of all CAGE-based mouse enhancers not transcribed in BMDM, only 19% carry macrophage enhancer chromatin signatures, highlighting the specificity of enhancers in mouse tissues.
Fig. 1

Macrophage enhancer–gene interactome. a Pipeline for identification of enhancers and enhancer–promoter associations. b Number and expression of genes associated with different number of enhancers. Dashed line shows median expression of genes not associated with any enhancer. c KEGG pathway maps significantly enriched for genes associated with a single enhancer, FDR < 0.05. d Top 15 KEGG pathway maps with the lowest FDR enriched for genes associated with more than 8 enhancers. In c and d next to the bars are the numbers of genes in the KEGG pathway covered by our gene list; dashed lines indicate FDR = 0.05

Macrophage enhancer–gene interactome

We aimed at studying enhancers that regulate expression of protein-coding genes in BMDM. We first identified pairs of enhancers and promoters located within TADs [47], since this regulation is thought to be exerted via direct enhancer–promoter contact [18, 19]. Thereafter, we refined these pairs using CAGE expression data based on the observation that eRNA and their target expression are positively correlated [34] (“Methods”). This yielded 222,870 TAD-based enhancer–promoter (E–P) pairs, with 64,891 pairs showing significant positive correlation of expression in macrophages (Fig. 1a). These correlation-based regulatory associations formed the basis for our further analyses and included 8667 enhancers deemed actively transcribed in mouse BMDM. Interestingly, most of the TAD-based E–P pairs show positive expression correlation (Additional file 3: Figure S1a), which supports the definition of a TAD as a structural unit favouring internal regulatory interactions [67]. Our filtering approach further selected regulatory associations with correlation coefficient above 0.3 (Additional file 3: Figure S1a), which we considered more reliable. The median distance between enhancers and promoters in the correlation-based E–P pairs was significantly smaller at 191,033nt as compared to 278,735nt for all TAD-based pairs (Additional file 3: Figure S1b).

We further investigated associations between enhancers and target protein-coding genes (Additional file 1: Table S3). Of all 10,767 protein-coding genes with CAGE expression (see “Methods”), 4149 genes (38.5%) were not associated with any enhancer in our settings (Fig. 1b, upper panel). Given previous evidence of additive action of enhancers [18, 68], we asked whether genes regulated by different numbers of enhancers have different gene expression levels. Genes without associated enhancers were overall lower expressed than genes associated with one (two-sided Wilcoxon signed-rank test p value < 2.2 × 10−16) or more enhancers. A steady increase in gene expression concomitant with higher numbers of associated enhancers (Fig. 1b, Kruskal–Wallis rank sum test p value < 2.2 × 10−16) was observed, supporting the model of additive enhancer action.

We further asked whether genes associated with different numbers of enhancers within the enhancer–gene interactome show functional differences. Gene set enrichment analysis (GSEA) was performed for gene sets of similar size to avoid a size-related bias (“Methods”). The 1416 genes associated with a single enhancer were enriched for housekeeping pathways including “Aminoacyl-tRNA biosynthesis” and “Ubiquitin mediated proteolysis”, as well as a few inflammation-related pathways (Fig. 1c). In contrast, the 1358 genes associated with more than eight enhancers showed stronger enrichment for signalling pathways important for macrophage immune function, such as “Jak-STAT signalling pathway” and “Chemokine signalling pathway” (Fig. 1d). GSEA for 1306 genes associated with three or four enhancers showed enrichment for a combination of housekeeping and immune pathways (Additional file 1: Table S4). Finally, the larger set of 4149 genes not associated with any enhancer showed the strongest enrichment for housekeeping pathways (Additional file 1: Table S5). Hence, a shift towards stronger enrichment for pathways important for macrophage immune function was a concomitant of higher numbers of associated enhancers.

Macrophage-specific expression

We opted for a similar strategy as Yao et al. [60] (“Methods”) to uncover eRNAs and genes with higher expression in macrophages as compared to other FANTOM5 mouse tissues (further referred to as macrophage specific). We identified 1844 macrophage-specific and 8923 non-macrophage-specific genes (Additional file 3: Figure S2). These two sets showed significant differences in numbers of associated enhancers, with 65.6% of macrophage-specific genes being associated with more than one enhancer, whereas this proportion dropped to 44.7% for non-macrophage-specific genes (odds ratio 1.99, Fisher’s exact test p value < 2.2 × 10−16) (Fig. 2a). These results were in agreement with our observation of stronger enrichment for macrophage-related functions in genes associated with many enhancers. Similar to the trend observed above, both macrophage-specific and non-macrophage-specific genes showed higher gene expression concomitant with higher numbers of associated enhancers, with non-macrophage-specific genes showing lower expression levels than macrophage-specific ones (Additional file 3: Figure S3).
Fig. 2

Macrophage-specific enhancer and gene expression. a Percentage of genes associated with different number of enhancers. b Expression of 4739 macrophage-specific enhancer eRNAs and 1481 associated genes. c Expression of 3928 non-macrophage-specific enhancer eRNAs and 1207 associated genes. In b and c, expression is shown in 184 macrophage and 744 non-macrophage samples; asterisks denote Wilcoxon rank sum test p value < 2.2 × 10−16. d Top 15 KEGG pathway maps significantly enriched for genes associated exclusively to macrophage-specific enhancers. e KEGG pathway maps enriched for genes associated exclusively to non-macrophage-specific enhancers with FDR < 0.05. In d and e next to the bars are the numbers of genes in the KEGG pathway covered by our gene list; dashed lines indicate FDR = 0.05

Among 8667 active enhancers, 54.7% were deemed macrophage specific (“Methods”), in agreement with known tissue specificity of enhancers [23, 42, 43]. Interestingly, non-macrophage-specific enhancers still showed higher eRNA expression in macrophages as compared to the non-macrophage samples (Fig. 2c, left panel). This may be explained by the fact that for this analysis we excluded all enhancers that showed zero eRNA expression in the majority of our macrophage samples (“Methods”).

Next, we asked whether these two enhancer sets could regulate genes with different functions. Genes associated exclusively with macrophage-specific enhancers, as well as genes associated exclusively with non-macrophage-specific enhancers, were sub-selected. As expected, genes in the former set showed overall higher expression in macrophage samples as compared to the non-macrophage samples (Fig. 2b, right panel). In contrast, expression of genes in the latter set was lower in macrophage samples (Fig. 2c, right panel). Interestingly, the opposite was observed for non-macrophage-specific enhancers (Fig. 2c, left panel). Genes associated with macrophage-specific enhancers were enriched for both housekeeping and immune pathways (Fig. 2d). This observation reflects the fact that production of macrophage-specific factors and activation of housekeeping processes that facilitate it might be both regulated by the same set of enhancers. Genes associated with non-macrophage-specific enhancers were enriched for only four KEGG pathway maps with FDR < 0.05 (Fig. 2e), none of which can be considered a typical macrophage pathway. We obtained consistent results when we repeated the analysis for a subset of 500 genes with the highest expression in macrophages (Additional file 1: Tables S6 and S7). Taken together, these findings demonstrate that most of the identified active enhancers in macrophages show macrophage-specific eRNA expression and regulate genes with macrophage specific as well as housekeeping functions.

Stimuli-induced transcriptional changes

We set out to determine transcriptional changes that were dynamically induced in M(IFN-γ) and M(IL-4/IL-13) mouse macrophages and to infer enhancers important in these processes (Fig. 3a). M(IFN-γ)- and M(IL-4/IL-13)-responsive enhancers and genes were identified as those up-regulated upon stimulation; regulatory associations were retained for pairs with a positive correlation of expression in the corresponding activation state (“Methods”). In this manner, we discovered 115 M(IFN-γ)-responsive enhancers regulating 105 M(IFN-γ)-responsive genes (further referred to as sets E1 and G1), as well as 131 M(IL-4/IL-13)-responsive enhancers regulating 98 M(IL-4/IL-13)-responsive genes (sets E2 and G2) (Fig. 3b and Additional file 1: Tables S8 and S9). GSEA of G1 and G2 gene sets showed significant enrichment for GO and KEGG terms relevant to immune system and macrophage functions (Fig. 3c and Additional file 3: Figure S4). These results highlight the importance of enhancer regulatory control during macrophage activation and suggest a striking influence of cytokine stimulation on activation of enhancers, which, in turn, drive some of the transcriptional responses seen during M(IFN-γ) and M(IL-4/IL-13) activation.
Fig. 3

Stimuli-responsive genes and enhancers. a Time-course data used in this study. b Enhancer and gene sets. E1 and E2: M(IFN-γ)- and M(IL-4/IL-13)-responsive enhancers regulating M(IFN-γ)- and M(IL-4/IL-13)-responsive genes (G1 and G2), respectively; E1′ and E2′: M(IFN-γ)- and M(IL-4/IL-13)-responsive enhancers regulating non-stimuli-responsive genes; G1′ and G2′: M(IFN-γ)- and M(IL-4/IL-13)-responsive genes not regulated by stimuli-responsive enhancers. Black arrows denote regulatory associations between stimuli-responsive enhancers and genes; of all previously defined significant associations, here we retained only those with positive correlation of expression in the corresponding subset of macrophages. c GO biological process terms enriched for G1 and G2 genes (all terms with FDR < 0.05 for G1; six terms with the lowest FDR for G2 are shown); dashed lines indicate FDR = 0.05. d Expression of stimuli-responsive enhancer eRNAs (upper panel) and genes (lower panel) unique to M(IFN-γ) or M(IL-4/IL-13). Statistical significance was determined by Wilcoxon signed-rank test, asterisks indicate p value < 10−5. e Correlation of time-course expression of M(IFN-γ)-responsive (upper panel) and M(IL-4/IL-13)-responsive (lower panel) enhancers and genes. Vertical dashed lines show median values

To gain an understanding of the distribution of stimuli-responsive elements in macrophage-specific and non-specific genes and enhancers, we calculated their overlaps (see Additional file 3: Figure S5). We found that 15.6% of macrophage-specific genes and 6.8% of macrophage-specific enhancers were also stimuli responsive. Of note, of all stimuli-responsive enhancers, 70.1% were also macrophage specific (77% of E1 and 71% of E2 enhancers).

M(IFN-γ) and M(IL-4/IL-13) macrophages are known to possess different phenotypes and functions [2]. As expected, G1 and G2 sets had only 19 genes in common. Similarly, a small overlap of only 14 enhancers was observed for E1 and E2 sets. Moreover, enhancers and genes selected as stimuli-responsive for a single activation state showed significant differences in time-course expression in M(IFN-γ) and M(IL-4/IL-13) macrophages (Fig. 3d). These data indicate that M(IFN-γ) and M(IL-4/IL-13) macrophages not only differ in their gene expression profiles, but also differ in their active enhancer repertoire that likely drives observed gene expression changes.

Previous studies reported and exploited positive expression correlation of eRNA and target genes [34, 38, 60]. Hence, we compared expression correlation of E1–G1 and E2–G2 pairs in M(IFN-γ) and M(IL-4/IL-13) macrophages (Fig. 3e) to determine how correlations differ between conditions. E1–G1 pairs showed higher correlation in M(IFN-γ) macrophages as compared to M(IL-4/IL-13) (two-sided Wilcoxon signed-rank test p value = 1.633 × 10−6). Similarly, correlation for E2–G2 pairs was higher in M(IL-4/IL-13) macrophages (two-sided Wilcoxon signed-rank test p value < 2.2 × 10−16). Such stimuli-specific expression correlation suggests stimuli specificity of enhancer–gene regulatory associations in macrophages.

Marker genes of macrophage activation are regulated by stimuli-responsive enhancers

We further asked which known marker genes of classical and alternative macrophage activation [1, 2, 6, 7, 69] were identified in M(IFN-γ) and M(IL-4/IL-13) in our setting (Table 1). Among 20 examined marker genes of classical macrophage activation, we found eight genes in the G1 set; similarly, eight of examined 26 marker genes of alternative activation were found in the G2 set (significant overlap with hypergeometric test p value < 10−10) (Table 1). The G1′ set contained an additional four classical macrophage activation marker genes (Gpr18, Il12b, Il6, Inhba) and the G2′ set an additional three alternative activation marker genes (Il27ra, Klf4, Myc), which, although stimuli-responsive themselves, were not associated with stimuli-responsive enhancers.
Table 1

M(IFN-γ) and M(IL-4/IL-13) macrophage activation markers

M(IFN-γ) markers

M(IL-4/IL-13) markers

Gene (G1)

# enh. in E1

M(IFN-γ) marker enhancers in E1

Gene (G2)

# enh. in E2

M(IL-4/IL-13) marker enhancers in E2

Cd38

1

 

Arg1

1

chr10:25119065…25119466

Cxcl9

3

chr5:92368373…92368774, chr5:92369052…92369453, chr5:92374704…92375105

Ccl24

1

 

Cxcl10

4

chr5:92353639…92354040, chr5:92368373…92368774, chr5:92369052…92369453, chr5:92374704…92375105

Egr2

9

chr10:67595184…67595585, chr10:67598488…67598889, chr10:67628888…67629289, chr10:67636538…67636939, chr10:67694800…67695201, chr10:67695848…67696249, chr10:67712611…67713012, chr10:67713071…67713472, chr10:67715029…67715430

Cxcl11

5

chr5:92353639…92354040, chr5:92368373…92368774, chr5:92369052…92369453, chr5:92374704…92375105, chr5:92375350…92375751

Fn1

3

chr1:71938511…71938912

Nos2

1

chr11:78916390…78916791

Igf1

7

chr10:87731929…87732330, chr10:87753519…87753920, chr10:87805812…87806213, chr10:87830718…87831119, chr10:87832100…87832501, chr10:87839444…87839845

Ptgs2

1

 

Irf4

2

chr13:30714614…30715015

Socs3

3

 

Mrc1

2

chr2:14185406…14185807, chr2:14206798…14207199

Tnf

1

 

Socs2

2

chr10:95232562…95232963, chr10:95236240…95236641

The columns list marker genes in G1 and G2, number of associated enhancers in the corresponding activation state, and potential activation marker enhancers

Next, we investigated the enhancer regulation of these marker genes. Given that different enhancers can modulate expression of the same gene in different conditions, we aimed to infer potential marker enhancers that regulate marker genes specifically during either M(IFN-γ) or M(IL-4/IL-13) activation. Each of the 16 marker genes in G1 and G2 was associated with a minimum of one and maximum of nine enhancers in the E1 and E2 stimuli-responsive sets, respectively (Table 1). Of those, we identified enhancers that were selectively responsive in a single activation state and showed higher eRNA expression in this state as compared to the other one (“Methods”). A total of 13 M(IFN-γ) and 22 M(IL-4/IL-13) enhancers were inferred as potential activation markers (Table 1).

Interestingly, three marker genes found in M(IFN-γ), Cxcl9, Cxcl10, and Cxcl11 are located within one TAD and are co-regulated by a group of three marker enhancers (Fig. 4a–c, see also Additional file 3: Figure S6 for a genome browser view of one of these enhancers). These enhancers, along with the two marker enhancers regulating Cxcl10 or Cxcl11 but not Cxcl9 (Table 1), are located in close proximity, in the intronic regions of the Art3 gene (Fig. 4c). These enhancer regions were previously reported to show induced RNA polymerase II binding in macrophages upon stimulation with LPS, one of the known classical macrophage activators [35]. In addition, these marker enhancer regions were shown to carry H3K4me1 enhancer histone marks in untreated macrophages [53]. Moreover, H3K27ac modification, associated with active enhancers, is stronger enriched in these regions in M(IFN-γ) as compared to M(IL-4) and untreated macrophages [53] (Fig. 4c), providing further evidence of their functionality in macrophage M(IFN-γ) activation.
Fig. 4

Examples of macrophage activation marker genes and enhancers. a Expression of classically activated macrophage marker genes Cxcl9, Cxcl10, and Cxcl11. b eRNA expression of three potential marker enhancers that co-regulate Cxcl9, Cxcl10, and Cxcl11. c Genomic region of a TAD containing Cxcl9, Cxcl10, Cxcl11, and associated enhancers. Black links connect the marker genes with the three potential marker enhancers. Grey links denote other enhancer–gene interactions that we identified in macrophages. d Expression of alternatively activated macrophage marker gene Egr2 and two of M(IL-4/IL-13) marker enhancers associated with Egr2. e Genomic region of a TAD containing Egr2 and associated enhancers. Black links connect Egr2 with the nine M(IL-4/IL-13) marker enhancers. Grey links denote other enhancer–gene interactions that we identified in macrophages. In a, b, and d, data were averaged over replicates and log-transformed. Error bars are the SEM. In c and e, genes are split into two tracks based on the strand, and wide orange marks denote gene promoters; histone mark tracks show ChIP-seq peaks with the height of − 10 × log10(p value), data from [53]

Among marker genes found in M(IL-4/IL-13), Arg1 as expected is substantially expressed in M(IL-4/IL-13) macrophages but has extremely low expression in M(IFN-γ) and untreated macrophages (Additional file 3: Figure S7a). We found a single M(IL-4/IL-13)-responsive enhancer that might drive expression of Arg1 in M(IL-4/IL-13) macrophages and might serve as a marker enhancer (Table 1 and Additional file 3: Figure S7b). On the contrary, marker gene Egr2 found in M(IL-4/IL-13), a TF that activates macrophage genes [70], is associated with as many as nine M(IL-4/IL-13) marker enhancers (Table 1). Egr2 showed immediate up-regulation in response to both IFN-γ and IL-4/IL-13 stimulation; however, in M(IL-4/IL-13) macrophages the up-regulation sustained for up to 24 h, whereas in M(IFN-γ) macrophages expression dropped rapidly after 2 h (Fig. 4d, upper panel). Time-course eRNA expression for two Egr2 marker enhancers with the highest expression at 2 and 4 h is shown in Fig. 4d. The distribution of all nine Egr2 marker enhancers within a TAD (Fig. 4e) may suggest that the regions identified as nine individual enhancers potentially demarcate fewer regions of stretch enhancers [71, 72]. We observed a similar distribution for enhancers of marker gene Igf1, which is known to shape the alternatively activated macrophage phenotype and regulate immune metabolism [73] (Additional file 3: Figure S8, see also Additional file 3: Figure S9 for a genome browser view of one of these enhancers with the highest expression at 2 h). Importantly, in both Egr2 and Igf1, marker enhancer regions carried H3K4me1 in untreated macrophages and showed the strongest enrichment with H3K27ac in M(IL-4) as compared to M(IFN-γ) and untreated macrophages [53] (Fig. 4e and Additional file 3: Figure S8c).

Transcription factor binding sites are enriched in enhancer regions

To investigate whether our enhancer sets are enriched for known transcription binding sites (TFBSs), we performed an over-representation analysis of experimentally determined protein DNA binding sites established through ChIP-seq [61, 62] (“Methods”). The sets of macrophage-specific and non-macrophage-specific enhancers are both enriched for binding sites of general factors (p300, Tbp), as well as a range of TFs with well-established roles in macrophages, such as macrophage lineage-determining factor Spi1 (PU.1) [41, 74], Cebpb, required for macrophage activation [75], and Rela, regulating inflammatory genes [76] (Additional file 1: Table S10). Interestingly, TFBSs for Spi1 overlap 54.1% of macrophage-specific enhancers, but only 38% of non-macrophage-specific enhancers (overlap ratio of 1.4 for macrophage-specific/non-macrophage-specific enhancers). We observed similar and higher overlap ratios for other functionally important TFs in macrophages, including Stat1, Rela, Irf1, Junb, and Cebpb [6, 7577] (Table 2). Six out of the seven TFs in Table 2, except for Junb, showed macrophage-specific expression (see above). Moreover, four TFs (Stat1, Irf1, Spi1, Cebpb) were significantly differentially expressed and up-regulated (FDR < 0.01, log2FC > 2) in our 184 macrophage samples when compared to the 744 FANTOM5 non-macrophage mouse samples (“Methods”). Furthermore, we found that macrophage-specific enhancers associated with Rela (chr19:5931597…5931998, chr19:6157210…6157611) and Spi1 (chr2:91086328…91086729, chr2:91204688…91205089) also carry TFBSs of the corresponding TFs, suggesting the formation of enhancer-mediated positive feedback loops, where a TF may induce its own enhancers. A genome browser view of one of these enhancers (chr2:91086328…91086729) shows a distribution of TFBSs and CAGE tags in the corresponding genomic region, which, notably, is located within an intron of the enhancer target gene Spi1 (Additional file 3: Figure S10). In addition, we show a genome browser view for M(IFN-γ) marker enhancer with non-macrophage-specific eRNA expression and M(IL-4/IL-13) marker enhancer with macrophage-specific eRNA expression (Additional file 3: Figures S6 and S9). Although both regions carry multiple TFBSs, we note that non-macrophage-specific enhancers can be characterized by a high occurrence of general enhancer proteins Ctcf and Ep300, as reflected in the TFBS over-representation analysis results (Additional file 1: Table S10).
Table 2

TFs regulating more macrophage-specific than non-macrophage-specific enhancers

TF

Mean TPM in macrophages

% Macrophage-specific enhancers

% Non-macrophage-specific enhancers

Ratio

Stat1

284.6

15.8

8.5

1.9

Rela

130.7

26.6

14.5

1.8

Atf4

329.4

25.5

14.3

1.8

Irf1

508.7

21.7

12.2

1.8

Junb

176.3

11.2

7.4

1.5

Spi1

599.9

54.1

38

1.4

Cebpb

647.9

48

34.8

1.4

Columns show TF name, average expression in macrophage samples (in TPM, see “Methods”), percentage of enhancers overlapping corresponding binding sites, and a macrophage-specific/non-macrophage-specific percentage ratio. TFBSs are statistically significantly enriched in both enhancer sets

The E1 and E2 stimuli-responsive enhancer sets are enriched for TFBSs of known macrophage TFs including Spi1, Cebpb, Rela, Irf, and Stat families [6, 78, 79] (Table 3, Additional file 1: Table S11). Interestingly, TFBSs of Stat1, Rela, and Irf1, involved in classical macrophage activation [6, 80], overlap a higher percentage of E1 enhancers as compared to E2 and macrophage-specific enhancers (Tables 2, 3, and Additional file 1: Tables S10 and S11). For instance, Irf1 TFBSs overlap 21.7% of macrophage-specific enhancers, 26.7% of E2 but 44.3% of E1 enhancers. In addition, Stat1 and Irf1 were deemed M(IFN-γ) responsive and were significantly differentially expressed and up-regulated (FDR < 0.01, log2FC > 1) in M(IFN-γ) when compared to either control or M(IL-4/IL-13) macrophages (“Methods”). Taken together, these results provide an additional layer of support for our regions as functionally important macrophage enhancers and implicate key macrophage TFs in modulating their activity. These findings further reflect that enhancers are selectively activated depending on the transcriptional machinery involved in the cellular response.
Table 3

TFs with binding sites enriched in both E1 and E2 enhancer sets

TF

Mean TPM in M(IFN-γ)

% Enhancers in M(IFN-γ)

Mean TPM in M(IL-4/IL-13)

% Enhancers in M(IL-4/IL-13)

Ratio

Stat1

726.3

40.9

218.6

22.9

1.8

Irf1

1345.4

44.3

426.7

26.7

1.7

Ets1

1.4

39.1

7

25.2

1.6

Jun

74.3

14.8

70

9.9

1.5

Rela

131.6

42.6

102.8

29.8

1.4

Atf4

277.3

33.9

208.3

27.5

1.2

Junb

122

13

99.7

13

1

Irf4

2.3

13

18.7

13.7

0.9

Spi1

766.3

52.2

681.7

56.5

0.9

Cebpb

424.6

51.3

262.9

59.5

0.9

Atf3

233.9

9.6

242.9

12.2

0.8

Columns show TF name, average expression in M(IFN-γ) and M(IL-4/IL-13) macrophages (in TPM, see “Methods”), percentage of enhancers overlapping corresponding binding sites, and a M(IFN-γ)/M(IL-4/IL-13) percentage ratio

Discussion

In this study, we investigated the transcribed enhancer landscape in mouse BMDM and its dynamic changes during M(IFN-γ) and M(IL-4/IL-13) activation. Using CAGE data combined with ChIP-seq, we identified 8667 actively transcribed enhancers forming 64,891 regulatory associations with protein-coding gene promoters in mouse BMDM. We highlighted tissue and stimuli specificity of both enhancers and their regulatory interactions. The enhancer–gene interactome established here supports a model of additive action of enhancers [18, 68], with higher gene expression concomitant with higher numbers of associated enhancers. Moreover, we observed a shift towards stronger enrichment for signalling pathways important for macrophage immune function in genes associated with many enhancers. Cytokine stimulation had a striking influence on enhancer transcription, which highlights the importance of enhancers in macrophage activation. In addition, we inferred potential macrophage activation marker enhancers. Finally, we find that binding sites of inflammatory TFs are enriched in enhancer regions, proposing a link between the response to stimuli and enhancer transcriptional activation.

Several studies have previously reported on enhancer landscape in mouse macrophages. Different populations of tissue macrophages were shown to be highly heterogeneous and to possess distinct sets of active enhancers, as defined by ChIP-seq profiling of histone modifications [28, 29]. Studies by Ostuni et al. [53] and Kaikkonen et al. [40] used ChIP-seq experiments to uncover enhancers that were established in macrophages de novo in response to a range of stimuli. In contrast to previous studies, we combined two complementary data types, transcriptomic (CAGE-derived) and epigenomic (ChIP-seq-derived, profiled by Ostuni et al. [53]) data, to infer more reliable transcribed active enhancers in mouse BMDM. Ostuni et al. separated macrophage enhancer regions into different enhancer classes based on the enhancer response to a range of stimuli. Our data show that in naïve macrophages 31% of active ChIP-seq-based enhancers show transcriptional activity. Of the poised not activated ChIP-seq enhancers, only 7.1% showed transcriptional activity in our set of 42,470 mouse enhancers. Both these observations support the idea of enhancer transcription being associated with histone-mark-based active states of enhancers.

Importantly, our analysis extended beyond identification of enhancers and characterization of their nearest genes. Here, instead of a widely used linear proximity-based approach [35, 38, 45], we employed TAD data to infer enhancer–gene associations. Accumulating evidence suggests that many enhancers regulate distal genes, bypassing the nearest promoter [81, 82]. At the same time, TADs have emerged as units of chromatin organization that favour internal DNA contacts [67], and the majority of characterized interactions between enhancers and target promoters occur within the same TAD [67, 82, 83]. Hence, our TAD-based approach enriched with correlation-based filtering enabled us to establish a more reliable mouse BMDM enhancer–gene interactome.

Our interactome covers 8667 actively transcribed enhancers. Of these, 70% overlap RNA polymerase II ChIP-seq peaks in untreated mouse macrophages [53]. Our enhancer regions show significant enrichment for binding sites of histone acetyltransferase p300, an enhancer-associated marker [26], and known inflammatory TFs. Hence, the regions identified here show a range of known enhancer properties, generally supporting our approach. Most of the active enhancers show macrophage-specific eRNA expression, in line with known tissue specificity of enhancers [23, 42, 43]. Kaikkonen et al. [40] identified mouse macrophage enhancers using ChIP-seq against H3K4me2. Our findings based on CAGE-seq extend their enhancer repertoire by additional 3974 transcribed enhancer regions. A comparison expanded to non-macrophage enhancers shows that 39.8% of our enhancers overlap a set of cis-regulatory elements from 19 non-macrophage mouse tissues identified by Shen et al. [42]. In another recent study, Schoenfelder et al. [33] employed a Capture Hi-C approach to identify enhancers and their target promoters in mouse foetal liver cells and embryonic stem cells. Even though enhancers and enhancer–promoter interactions are known to be highly tissue specific, their data still support 24.8% of our 64,891 E–P pairs.

Recent reports suggested that genes regulated by multiple enhancers were higher expressed than those regulated by a single enhancer, proposing that enhancers might contribute additively to the expression of their target genes [18, 68]. In support of this, we observed a steady increase in gene expression concomitant with increasing numbers of associated enhancers, with the genes not associated with any enhancers showing the lowest overall expression. A study of 12 mouse tissues has reported the enrichment for tissue-specific functions in genes associated with enhancers that transcribe eRNAs as compared to genes associated with non-transcribed enhancers [84]. Jin et al. recently showed that genes that did not interact with distal enhancers were enriched for housekeeping genes and also suggested that cell-specific genes were extensively controlled by cis-regulators [85]. We showed in this study that genes associated with many enhancers were more enriched for macrophage-related functions as compared to genes associated with only few or no enhancers. This finding might reflect a more fundamental principle of genome organization and evolution, such as the importance of multiple enhancers for fine-tuned and redundant control of cell specialization and cell-specific responses.

Studies by Ostuni et al. [53] and Kaikkonen et al. [40] revealed stimuli-specific epigenomic changes in enhancer regions in mouse macrophages and introduced a concept of stimuli-specific enhancer activation. In our study, we focused on enhancers and genes that responded to the stimuli with increased expression in order to further investigate this phenomenon. Notably, many stimuli-responsive genes were associated with stimuli-responsive enhancers, highlighting the importance of enhancer regulation in macrophage activation. As expected for such a cell-type-specific process as macrophage activation, most of the responsive enhancers showed macrophage-specific eRNA expression, and genes were enriched for macrophage-specific functions. In addition, our study suggests stimuli specificity of enhancer–gene regulatory associations in macrophages.

As an important example, we assessed 20 and 26 marker genes of classically and alternatively activated macrophages to characterize their enhancer regulation [1, 2, 6, 7, 69]. Of those, seven markers (Ccl20, Fpr2, Ido1; Chi3l3, Chi3l4, Alox12e, Chia) were not expressed in our data. For a total of 16 marker genes, we identified associated enhancers. Moreover, for 11 of them we found enhancers that might regulate these genes specifically in M(IFN-γ) or M(IL-4/IL-13) stimulation (Table 1). Hence, these enhancers present new potential markers for a particular macrophage activation status. Seven additional marker genes, identified as stimuli responsive, were not associated with any stimuli-responsive enhancer (Gpr18, Il12b, Il6, Inhba in the G1′ set; Il27ra, Klf4, Myc in the G2′ set). The remaining marker genes were not deemed stimuli-responsive. Of those, classically activated macrophage markers Il1b, Cd86, Marco, and Il23a and alternatively activated macrophage markers Mmp12, Tgm2, Clec4a2, Stab1, F13a1 were associated with at least one enhancer in macrophages. Ccr7, Retnla, Ccl17, Ccl22, Chi3l1, Cxcl13, and Ccl12 were not associated with any enhancers in macrophages.

We observed a particular genomic distribution of potential marker enhancers associated with Egr2 and Igf1 marker genes in M(IL-4/IL-13), which suggested that these regulatory DNA regions might represent stretch enhancers. Parker et al., in a recent study, investigated stretch enhancers in human cells and proposed that such extended regions could serve as molecular runways to attract tissue-specific TFs and focus their activity [72]. Similarly to Parker et al., potential stretch enhancer regions identified here were associated with cell-type specific genes and were demarcated by broad H3K27ac signals, specifically higher enriched in M(IL-4) as compared to M(IFN-γ) and untreated macrophages (Fig. 4e, Additional file 3: Figure S8c). Therefore, we propose that stretch enhancers might be involved in the regulation of macrophage activation. However, further studies are required to investigate this phenomenon in more detail.

Our approach inferred M(IFN-γ)- and M(IL-4/IL-13)-responsive enhancers that were strongly enriched for TFBS of known inflammatory TFs. These results are in line with previous reports in mouse macrophages. For example, Spi1 (PU.1) has been extensively studied as a crucial TF involved in macrophage differentiation and transcriptional regulation [41]. Moreover, Spi1 was deemed a pioneering or lineage-determining TF in macrophages, which defines enhancer regions and occupies many enhancers in macrophages [28, 35, 41, 74]. Furthermore, Heinz et al. suggested that collaborative action of Spi1 with Cebpb was required for the deposition of enhancer-associated chromatin marks [74]. Ghisletti et al. reported enrichment for NF-kB (Rel) and Irf TFs in enhancers induced by LPS in mouse macrophages [41]. Likewise, transcribed enhancers induced by LPS and IFN-γ stimulation showed enrichment for NF-kB/Rel, Irf, and Stat1 binding motifs [35]. In addition, we previously showed that TFs including Rela and Irf1 drive expression of protein-coding and lncRNA genes during macrophage activation [10]. Taken together, our results link enhancer activation to the transcriptional programme induced by IFN-γ and IL-4/IL-13 stimuli.

Conclusions

In this study, we have established a genome-wide catalogue of enhancers and enhancer–promoter regulatory interactions in mouse BMDM. In contrast to previous studies of enhancer landscape in mouse macrophages, we focused on transcribed enhancers and employed an improved method for identification of enhancer target genes, based on location within a TAD and correlation of expression. Hence, our study represents the most comprehensive analysis of transcribed enhancer activities in mouse macrophages to date and extends current knowledge of transcriptional regulation in macrophages in general and during activation in particular.

Abbreviations

CAGE: 

cap analysis of gene expression

ChIP-seq: 

chromatin immunoprecipitation followed by sequencing

E–P pairs: 

enhancer–promoter pairs

FDR: 

false discovery rate

GSEA: 

gene set enrichment analysis

TADs: 

topologically associating domains

TFBS: 

transcription factor binding site

TF: 

transcription factor

TPM: 

tags per million

Declarations

Authors’ contributions

ED performed computational analyses. SS designed the study. SS and ED analysed data, interpreted results, and wrote the manuscript with input from all authors. RG, MMM, HS, and FB helped interpret results and provided data. All authors read and approved the final manuscript.

Acknowledgements

None.

Competing interests

The authors declare that they have no competing interests.

Availability of data and materials

The dataset analysed in the study is available in the FANTOM5 repository, http://fantom.gsc.riken.jp/5/datafiles/reprocessed/mm10_v2/basic/. The datasets supporting the conclusions of this article are included within the article and its additional files.

Consent for publication

Not applicable.

Ethics approval and consent to participate

Not applicable.

Funding

This work was supported by grants from the South African National Research Foundation (NRF) and from the Department of Science and Technology, South African Research Chair Initiative (SARCHi) and South Africa Medical Research Council (SAMRC) to FB, grant from the Japan Society for the Promotion of Science (JSPS) and National Research Foundation of South Africa to FB and HS, and grants from the South African National Research Foundation (NRF) Competitive Programme for Unrated Researchers (CSUR) to RG.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Authors’ Affiliations

(1)
Institute of Natural and Mathematical Sciences, Massey University
(2)
Division of Immunology, Institute of Infectious Diseases and Molecular Medicine (IDM), South African Medical Research Council (SAMRC) Immunology of Infectious Diseases, Faculty of Health Sciences, University of Cape Town
(3)
International Centre for Genetic Engineering and Biotechnology (ICGEB), Cape Town Component
(4)
Gene Expression and Biophysics Group, CSIR Synthetic Biology ERA
(5)
Division of Chemical Systems and Synthetic Biology, Institute of Infectious Diseases and Molecular Medicine (IDM), University of Cape Town
(6)
Gene Expression and Biophysics Unit, Instituto de Medicina Molecular, Faculdade de Medicina da Universidade de Lisboa
(7)
Division of Genomic Technologies, RIKEN Center for Life Science Technologies

References

  1. Murray PJ, Wynn TA. Protective and pathogenic functions of macrophage subsets. Nat Rev Immunol. 2011;11:723–37.View ArticlePubMedPubMed CentralGoogle Scholar
  2. Mosser DM, Edwards JP. Exploring the full spectrum of macrophage activation. Nat Rev Immunol. 2008;8:958–69.View ArticlePubMedPubMed CentralGoogle Scholar
  3. Ginhoux F, Schultze JL, Murray PJ, Ochando J, Biswas SK. New insights into the multidimensional concept of macrophage ontogeny, activation and function. Nat Immunol. 2016;17:34–40.View ArticlePubMedGoogle Scholar
  4. Murray PJ, Allen JE, Biswas SK, Fisher EA, Gilroy DW, Goerdt S, et al. Macrophage activation and polarization: nomenclature and experimental guidelines. Immunity. 2014;41:14–20.View ArticlePubMedPubMed CentralGoogle Scholar
  5. Sica A, Mantovani A. Macrophage plasticity and polarization: in vivo veritas. J Clin Investig. 2012;122:787–95.View ArticlePubMedPubMed CentralGoogle Scholar
  6. Martinez FO, Gordon S. The M1 and M2 paradigm of macrophage activation: time for reassessment. F1000Prime Rep. 2014;6:13.View ArticlePubMedPubMed CentralGoogle Scholar
  7. Gordon S, Martinez FO. Alternative activation of macrophages: mechanism and functions. Immunity. 2010;32:593–604.View ArticlePubMedGoogle Scholar
  8. Herbert DR, Holscher C, Mohrs M, Arendse B, Schwegmann A, Radwanska M, et al. Alternative macrophage activation is essential for survival during schistosomiasis and downmodulates T helper 1 responses and immunopathology. Immunity. 2004;20:623–35.View ArticlePubMedGoogle Scholar
  9. Jenkins SJ, Ruckerl D, Thomas GD, Hewitson JP, Duncan S, Brombacher F, et al. IL-4 directly signals tissue-resident macrophages to proliferate beyond homeostatic levels controlled by CSF-1. J Exp Med. 2013;210:2477–91.View ArticlePubMedPubMed CentralGoogle Scholar
  10. Roy S, Schmeier S, Arner E, Alam T, Parihar SP, Ozturk M, et al. Redefining the transcriptional regulatory dynamics of classically and alternatively activated macrophages by deepCAGE transcriptomics. Nucleic Acids Res. 2015;43:6969–82.View ArticlePubMedPubMed CentralGoogle Scholar
  11. Lugo-Villarino G, Verollet C, Maridonneau-Parini I, Neyrolles O. Macrophage polarization: convergence point targeted by Mycobacterium tuberculosis and HIV. Front Immunol. 2011;2:43.View ArticlePubMedPubMed CentralGoogle Scholar
  12. Guler R, Parihar SP, Savvi S, Logan E, Schwegmann A, Roy S, et al. IL-4Ralpha-dependent alternative activation of macrophages is not decisive for Mycobacterium tuberculosis pathology and bacterial burden in mice. PLoS ONE. 2015;10:e0121070.View ArticlePubMedPubMed CentralGoogle Scholar
  13. Odegaard JI, Ricardo-Gonzalez RR, Goforth MH, Morel CR, Subramanian V, Mukundan L, et al. Macrophage-specific PPARgamma controls alternative activation and improves insulin resistance. Nature. 2007;447:1116–20.View ArticlePubMedPubMed CentralGoogle Scholar
  14. Komili S, Silver PA. Coupling and coordination in gene expression processes: a systems biology view. Nat Rev Genet. 2008;9:38–48.View ArticlePubMedGoogle Scholar
  15. Weake VM, Workman JL. Inducible gene expression: diverse regulatory mechanisms. Nat Rev Genet. 2010;11:426–37.View ArticlePubMedGoogle Scholar
  16. Lenhard B, Sandelin A, Carninci P. Metazoan promoters: emerging characteristics and insights into transcriptional regulation. Nat Rev Genet. 2012;13:233–45.PubMedGoogle Scholar
  17. Forrest ARR, Kawaji H, Rehli M, Baillie JK, de Hoon MJL, Haberle V, et al. A promoter-level mammalian expression atlas. Nature. 2014;507:462–70.View ArticlePubMedGoogle Scholar
  18. Shlyueva D, Stampfel G, Stark A. Transcriptional enhancers: from properties to genome-wide predictions. Nat Rev Genet. 2014;15:272–86.View ArticlePubMedGoogle Scholar
  19. Gorkin DU, Leung D, Ren B. The 3D genome in transcriptional regulation and pluripotency. Cell Stem Cell. 2014;14:762–75.View ArticlePubMedPubMed CentralGoogle Scholar
  20. Creyghton MP, Cheng AW, Welstead GG, Kooistra T, Carey BW, Steine EJ, et al. Histone H3K27ac separates active from poised enhancers and predicts developmental state. Proc Natl Acad Sci USA. 2010;107:21931–6.View ArticlePubMedPubMed CentralGoogle Scholar
  21. Heintzman ND, Stuart RK, Hon G, Fu Y, Ching CW, Hawkins RD, et al. Distinct and predictive chromatin signatures of transcriptional promoters and enhancers in the human genome. Nat Genet. 2007;39:311–8.View ArticlePubMedGoogle Scholar
  22. Spitz F, Furlong EE. Transcription factors: from enhancer binding to developmental control. Nat Rev Genet. 2012;13:613–26.View ArticlePubMedGoogle Scholar
  23. Romanoski CE, Link VM, Heinz S, Glass CK. Exploiting genomics and natural genetic variation to decode macrophage enhancers. Trends Immunol. 2015;36:507–18.View ArticlePubMedPubMed CentralGoogle Scholar
  24. Heintzman ND, Hon GC, Hawkins RD, Kheradpour P, Stark A, Harp LF, et al. Histone modifications at human enhancers reflect global cell-type-specific gene expression. Nature. 2009;459:108–12.View ArticlePubMedPubMed CentralGoogle Scholar
  25. Calo E, Wysocka J. Modification of enhancer chromatin: what, how, and why? Mol Cell. 2013;49:825–37.View ArticlePubMedGoogle Scholar
  26. Visel A, Blow MJ, Li Z, Zhang T, Akiyama JA, Holt A, et al. ChIP-seq accurately predicts tissue-specific activity of enhancers. Nature. 2009;457:854–8.View ArticlePubMedPubMed CentralGoogle Scholar
  27. Barski A, Cuddapah S, Cui K, Roh TY, Schones DE, Wang Z, et al. High-resolution profiling of histone methylations in the human genome. Cell. 2007;129:823–37.View ArticlePubMedGoogle Scholar
  28. Gosselin D, Link VM, Romanoski CE, Fonseca GJ, Eichenfield DZ, Spann NJ, et al. Environment drives selection and function of enhancers controlling tissue-specific macrophage identities. Cell. 2014;159:1327–40.View ArticlePubMedPubMed CentralGoogle Scholar
  29. Lavin Y, Winter D, Blecher-Gonen R, David E, Keren-Shaul H, Merad M, et al. Tissue-resident macrophage enhancer landscapes are shaped by the local microenvironment. Cell. 2014;159:1312–26.View ArticlePubMedPubMed CentralGoogle Scholar
  30. Dekker J, Rippe K, Dekker M, Kleckner N. Capturing chromosome conformation. Science. 2002;295:1306–11.View ArticlePubMedGoogle Scholar
  31. van Steensel B, Dekker J. Genomics tools for unraveling chromosome architecture. Nat Biotechnol. 2010;28:1089–95.View ArticlePubMedPubMed CentralGoogle Scholar
  32. de Laat W, Klous P, Kooren J, Noordermeer D, Palstra RJ, Simonis M, et al. Three-dimensional organization of gene expression in erythroid cells. Curr Top Dev Biol. 2008;82:117–39.View ArticlePubMedGoogle Scholar
  33. Schoenfelder S, Furlan-Magaril M, Mifsud B, Tavares-Cadete F, Sugar R, Javierre BM, et al. The pluripotent regulatory circuitry connecting promoters to their long-range interacting elements. Genome Res. 2015;25:582–97.View ArticlePubMedPubMed CentralGoogle Scholar
  34. Kim TK, Hemberg M, Gray JM, Costa AM, Bear DM, Wu J, et al. Widespread transcription at neuronal activity-regulated enhancers. Nature. 2010;465:182–7.View ArticlePubMedPubMed CentralGoogle Scholar
  35. De Santa F, Barozzi I, Mietton F, Ghisletti S, Polletti S, Tusi BK, et al. A large fraction of extragenic RNA Pol II transcription sites overlap enhancers. PLoS Biol. 2010;8:e1000384.View ArticlePubMedPubMed CentralGoogle Scholar
  36. Denisenko E, Ho D, Tamgue O, Ozturk M, Suzuki H, Brombacher F, et al. IRNdb: the database of immunologically relevant non-coding RNAs. Database. 2016. doi:10.1093/database/baw138.PubMed CentralGoogle Scholar
  37. Hon CC, Ramilowski JA, Harshbarger J, Bertin N, Rackham OJ, Gough J, et al. An atlas of human long non-coding RNAs with accurate 5′ ends. Nature. 2017;543:199–204.View ArticlePubMedGoogle Scholar
  38. Andersson R, Gebhard C, Miguel-Escalada I, Hoof I, Bornholdt J, Boyd M, et al. An atlas of active enhancers across human cell types and tissues. Nature. 2014;507:455–61.View ArticlePubMedPubMed CentralGoogle Scholar
  39. Kanamori-Katayama M, Itoh M, Kawaji H, Lassmann T, Katayama S, Kojima M, et al. Unamplified cap analysis of gene expression on a single-molecule sequencer. Genome Res. 2011;21:1150–9.View ArticlePubMedPubMed CentralGoogle Scholar
  40. Kaikkonen MU, Spann NJ, Heinz S, Romanoski CE, Allison KA, Stender JD, et al. Remodeling of the enhancer landscape during macrophage activation is coupled to enhancer transcription. Mol Cell. 2013;51:310–25.View ArticlePubMedPubMed CentralGoogle Scholar
  41. Ghisletti S, Barozzi I, Mietton F, Polletti S, De Santa F, Venturini E, et al. Identification and characterization of enhancers controlling the inflammatory gene expression program in macrophages. Immunity. 2010;32:317–28.View ArticlePubMedGoogle Scholar
  42. Shen Y, Yue F, McCleary DF, Ye Z, Edsall L, Kuan S, et al. A map of the cis-regulatory sequences in the mouse genome. Nature. 2012;488:116–20.View ArticlePubMedPubMed CentralGoogle Scholar
  43. Kieffer-Kwon KR, Tang Z, Mathe E, Qian J, Sung MH, Li G, et al. Interactome maps of mouse gene regulatory domains reveal basic principles of transcriptional regulation. Cell. 2013;155:1507–20.View ArticlePubMedGoogle Scholar
  44. Thomas GD, Hanna RN, Vasudevan NT, Hamers AA, Romanoski CE, McArdle S, et al. Deleting an Nr4a1 super-enhancer subdomain ablates Ly6Clow monocytes while preserving macrophage gene function. Immunity. 2016;45:975–87.View ArticlePubMedGoogle Scholar
  45. Arner E, Daub CO, Vitting-Seerup K, Andersson R, Lilje B, Drablos F, et al. Transcribed enhancers lead waves of coordinated transcription in transitioning mammalian cells. Science. 2015;347:1010–4.View ArticlePubMedPubMed CentralGoogle Scholar
  46. Mukhopadhyay S, Ramadass AS, Akoulitchev A, Gordon S. Formation of distinct chromatin conformation signatures epigenetically regulate macrophage activation. Int Immunopharmacol. 2014;18:7–11.View ArticlePubMedGoogle Scholar
  47. Dixon JR, Selvaraj S, Yue F, Kim A, Li Y, Shen Y, et al. Topological domains in mammalian genomes identified by analysis of chromatin interactions. Nature. 2012;485:376–80.View ArticlePubMedPubMed CentralGoogle Scholar
  48. FANTOM consortium: FANTOM5 Data. http://fantom.gsc.riken.jp/5/datafiles/reprocessed/mm10_v2/basic/ (2016). Accessed 07 Sept 2016.
  49. Kawaji H. Decomposition-based peak identification (DPI): source code. https://github.com/hkawaji/dpi1/ (2016). Accessed 04 Sept 2016.
  50. Flicek P, Amode MR, Barrell D, Beal K, Brent S, Chen Y, et al. Ensembl 2011. Nucleic Acids Res. 2011;39:D800–6.View ArticlePubMedGoogle Scholar
  51. Karolchik D, Hinrichs AS, Furey TS, Roskin KM, Sugnet CW, Haussler D, et al. The UCSC Table Browser data retrieval tool. Nucleic Acids Res. 2004;32:D493–6.View ArticlePubMedPubMed CentralGoogle Scholar
  52. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.View ArticlePubMedGoogle Scholar
  53. Ostuni R, Piccolo V, Barozzi I, Polletti S, Termanini A, Bonifacio S, et al. Latent enhancers activated by stimulation in differentiated cells. Cell. 2013;152:157–71.View ArticlePubMedGoogle Scholar
  54. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B (Methodol). 1995;57:289–300.Google Scholar
  55. UCSC Genome Browser Utilities: Batch Coordinate Conversion (liftOver). https://genome.ucsc.edu/cgi-bin/hgLiftOver. Accessed 15 Sept 2016.
  56. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28:27–30.View ArticlePubMedPubMed CentralGoogle Scholar
  57. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000;25:25–9.View ArticlePubMedPubMed CentralGoogle Scholar
  58. Carlson M. GO.db: a set of annotation maps describing the entire gene ontology. 2015.Google Scholar
  59. da Huang W, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009;37:1–13.View ArticleGoogle Scholar
  60. Yao P, Lin P, Gokoolparsadh A, Assareh A, Thang MW, Voineagu I. Coexpression networks identify brain region-specific enhancer RNAs in the human brain. Nat Neurosci. 2015;18:1168–74.View ArticlePubMedGoogle Scholar
  61. Yue F, Cheng Y, Breschi A, Vierstra J, Wu W, Ryba T, et al. A comparative encyclopedia of DNA elements in the mouse genome. Nature. 2014;515:355–64.View ArticlePubMedPubMed CentralGoogle Scholar
  62. Garber M, Yosef N, Goren A, Raychowdhury R, Thielke A, Guttman M, et al. A high-throughput chromatin immunoprecipitation approach reveals principles of dynamic gene regulation in mammals. Mol Cell. 2012;47:810–22.View ArticlePubMedGoogle Scholar
  63. Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 2008;9:R137.View ArticlePubMedPubMed CentralGoogle Scholar
  64. North BV, Curtis D, Sham PC. A note on the calculation of empirical P values from Monte Carlo procedures. Am J Hum Genet. 2002;71:439–41.View ArticlePubMedPubMed CentralGoogle Scholar
  65. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.View ArticlePubMedPubMed CentralGoogle Scholar
  66. R Core Team. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing. 2015. http://www.R-project.org/.
  67. Rocha PP, Raviram R, Bonneau R, Skok JA. Breaking TADs: insights into hierarchical genome organization. Epigenomics. 2015;7:523–6.View ArticlePubMedPubMed CentralGoogle Scholar
  68. Chepelev I, Wei G, Wangsa D, Tang Q, Zhao K. Characterization of genome-wide enhancer–promoter interactions reveals co-expression of interacting genes and modes of higher order chromatin organization. Cell Res. 2012;22:490–503.View ArticlePubMedPubMed CentralGoogle Scholar
  69. Jablonski KA, Amici SA, Webb LM, Ruiz-Rosado J, Popovich PG, Partida-Sanchez S, et al. Novel markers to delineate murine M1 and M2 macrophages. PLoS ONE. 2015;10:e0145342.View ArticlePubMedPubMed CentralGoogle Scholar
  70. Laslo P, Spooner CJ, Warmflash A, Lancki DW, Lee HJ, Sciammas R, et al. Multilineage transcriptional priming and determination of alternate hematopoietic cell fates. Cell. 2006;126:755–66.View ArticlePubMedGoogle Scholar
  71. Hnisz D, Abraham BJ, Lee TI, Lau A, Saint-André V, Sigova AA, et al. Super-enhancers in the control of cell identity and disease. Cell. 2013;155:934–47.View ArticlePubMedGoogle Scholar
  72. Parker SC, Stitzel ML, Taylor DL, Orozco JM, Erdos MR, Akiyama JA, et al. Chromatin stretch enhancer states drive cell-specific gene regulation and harbor human disease risk variants. Proc Natl Acad Sci USA. 2013;110:17921–6.View ArticlePubMedPubMed CentralGoogle Scholar
  73. Spadaro O, Camell CD, Bosurgi L, Nguyen KY, Youm YH, Rothlin CV, et al. IGF1 shapes macrophage activation in response to immunometabolic challenge. Cell Rep. 2017;19:225–34.View ArticlePubMedPubMed CentralGoogle Scholar
  74. Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol Cell. 2010;38:576–89.View ArticlePubMedPubMed CentralGoogle Scholar
  75. Ruffell D, Mourkioti F, Gambardella A, Kirstetter P, Lopez RG, Rosenthal N, et al. A CREB-C/EBPbeta cascade induces M2 macrophage-specific gene expression and promotes muscle injury repair. Proc Natl Acad Sci USA. 2009;106:17475–80.View ArticlePubMedPubMed CentralGoogle Scholar
  76. Saliba DG, Heger A, Eames HL, Oikonomopoulos S, Teixeira A, Blazek K, et al. IRF5:RelA interaction targets inflammatory genes in macrophages. Cell Rep. 2014;8:1308–17.View ArticlePubMedPubMed CentralGoogle Scholar
  77. Fontana MF, Baccarella A, Pancholi N, Pufall MA, Herbert DR, Kim CC. JUNB is a key transcriptional modulator of macrophage activation. J Immunol. 2015;194:177–86.View ArticlePubMedGoogle Scholar
  78. Tugal D, Liao X, Jain MK. Transcriptional control of macrophage polarization. Atertio Thromb Vasc Biol. 2013;33:1135–44.View ArticleGoogle Scholar
  79. Taniguchi T, Ogasawara K, Takaoka A, Tanaka N. IRF family of transcription factors as regulators of host defense. Annu Rev Immunol. 2001;19:623–55.View ArticlePubMedGoogle Scholar
  80. Roy S, Guler R, Parihar SP, Schmeier S, Kaczkowski B, Nishimura H, et al. Batf2/Irf1 induces inflammatory responses in classically activated macrophages, lipopolysaccharides, and mycobacterial infection. J Immunol. 2015;194:6035–44.View ArticlePubMedGoogle Scholar
  81. Li G, Ruan X, Auerbach RK, Sandhu KS, Zheng M, Wang P, et al. Extensive promoter-centered chromatin interactions provide a topological basis for transcription regulation. Cell. 2012;148:84–98.View ArticlePubMedPubMed CentralGoogle Scholar
  82. Lupianez DG, Spielmann M, Mundlos S. Breaking TADs: how alterations of chromatin domains result in disease. Trends Genet. 2016;32:225–37.View ArticlePubMedGoogle Scholar
  83. Symmons O, Uslu VV, Tsujimura T, Ruf S, Nassari S, Schwarzer W, et al. Functional and topological characteristics of mammalian regulatory domains. Genome Res. 2014;24:390–400.View ArticlePubMedPubMed CentralGoogle Scholar
  84. Cheng JH, Pan DZ, Tsai ZT, Tsai HK. Genome-wide analysis of enhancer RNA in gene regulation across 12 mouse tissues. Sci Rep. 2015;5:12648.View ArticlePubMedPubMed CentralGoogle Scholar
  85. Jin F, Li Y, Dixon JR, Selvaraj S, Ye Z, Lee AY, et al. A high-resolution map of the three-dimensional chromatin interactome in human cells. Nature. 2013;503:290–4.PubMedPubMed CentralGoogle Scholar

Copyright

© The Author(s) 2017

Advertisement