Genome-wide profiling of transcribed enhancers during macrophage activation.

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.

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, 7–9]. 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, 24–26]. 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, 75–77] (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.

Additional file 1: Table S1. 184 macrophage samples used in this study. Table S2. 744 FANTOM5 non-macrophage mouse samples that were used as a background set for the calculation of macrophage-specific expression. Table S3. Established regulatory associations between transcribed enhancers, promoters and protein-coding genes. Table S4. Gene set enrichment analysis of 1306 genes associated with three or four enhancers. Table S5. Gene set enrichment analysis of 4149 genes not associated with any enhancer. Table S6. Gene set enrichment analysis of top 500 genes with the highest expression in macrophages among 1481 genes associated with macrophage-specific enhancers. Table S7. Gene set enrichment analysis of top 500 genes with the highest expression in macrophages among 1207 genes associated with non-macrophage-specific enhancers. Table S8. Regulatory associations between 115 M(IFN-γ)-responsive enhancers and 105 M(IFN-γ)-responsive genes. Table S9. Regulatory associations between 131 M(IL-4/IL-13)-responsive enhancers and 98 M(IL-4/IL-13)-responsive genes. Table S10. Over-representation analysis of binding sites in macrophage-specific and non-macrophage-specific enhancers. Table S11. Over-representation analysis of binding sites in E1 and E2 enhancers.

Additional file 2: Table S12. Numbers of overlapping regions between the CAGE-based enhancers and classes of ChIP-seq-based enhancers. Table S13. A matrix that shows whether each CAGE-based transcribed macrophage enhancer overlaps (value of 1) or does not overlap (value of 0) ChIP-seq-based enhancers of different classes.

Additional file 3: Figure S1. Comparison of 222,870 TAD-based E–P pairs to a subset of 64,891 correlation-based E–P pairs. Figure S2. 1844 macrophage-specific and 8923 non-macrophage-specific genes. Figure S3. Expression of macrophage-specific and non-macrophage-specific genes associated with different number of enhancers. Figure S4. KEGG pathway maps significantly enriched for G1 and G2 genes. Figure S5. Overlaps of M(IFN-γ)- and M(IL-4/IL-13)-responsive and macrophage-specific genes and enhancers. Figure S6. M(IFN-γ) marker enhancer associated with Cxcl9, Cxcl10, and Cxcl11 M(IFN-γ) marker genes. Figure S7. Time-course expression of Arg1 and associated M(IL-4/IL-13)-specific enhancer. Figure S8. Igf1 marker gene. Figure S9. M(IL-4/IL-13) marker enhancer associated with Igf1 M(IL-4/IL-13) marker gene. Figure S10. Macrophage-specific enhancer, associated with Spi1 gene.