Mechanisms underlying divergent responses of genetically distinct macrophages to IL-4

To investigate how natural genetic variation affects the macrophage response to IL-4, we began by performing RNA sequencing (RNA-seq) in BMDMs derived from female BALB/cJ (BALB), C57BL/6J (C57), NOD/ShiLtJ (NOD), PWK/PhJ (PWK), and SPRET/EiJ (SPRET) mice under basal conditions and following stimulation with IL-4. Time-course experiments in C57 BMDMs indicated a progressive increase in the number of differentially expressed genes from 1 to 24 hours (fig. S1, A and B, and table S1). We therefore focused our analysis on the response to IL-4 in BMDMs from the five strains at this time point. Weighted gene coexpression network analysis (WGCNA) identified numerous modules of highly correlated mRNAs, most of which were driven by strain differences (Fig. 1B). Genes that were positively regulated by IL-4 across strains (red module, bottom) were enriched for functional annotations related to negative regulation of defense responses. Conversely, the purple (top) module captured genes that were negatively regulated by IL-4 and were enriched for pathways associated with positive regulation of inflammation (Fig. 1B).

Of the 693 genes induced >2-fold in at least one strain, only 26 (3.75%) were induced at this threshold in all five strains (Fig. 1D, fig. S1, C and D, and table S2). Conversely, more than half of the IL-4–responsive genes identified were induced >2-fold in only a single strain. NOD BMDMs were notable for a generally attenuated response to IL-4 (Fig. 1B, red module, and Fig. 1C, second panel). A similar pattern was observed for down-regulated genes (Fig. 1D). Despite these differences at the level of individual genes, similar pathways/gene programs were enriched in all strains for both induced and repressed genes (Fig. 1E). Substantial differences in IL-4 target gene expression across strains are illustrated by Arg1, Slc7a2, and Msx3 (Fig. 1F and fig. S1E). BMDMs from all strains exhibit a significant induction of Arg1 expression, but the absolute basal levels and induction folds vary by more than an order of magnitude. Slc7a2 exhibits similar levels of expression in C57 and NOD BMDMs after IL-4 treatment, but its differences at the basal level result in an 8.4-fold and 1.2-fold change, respectively. We refer to the pattern of reduced responsiveness to IL-4 in this comparison of C57 and NOD as being associated with “high basal” activity in the less responsive strain. Conversely, NOD and PWK BMDMs exhibit similar levels of basal Slc7a2 expression, but IL-4 only increased Slc7a2 expression more than twofold in PWK. We refer to this pattern of reduced responsiveness to IL-4 in NOD compared to PWK as being associated with “equal basal” activity. A third category is exemplified by Msx3, which is induced in C57 but not in PWK and SPRET BMDMs. In this case, lack of responsiveness is associated with low expression of Msx3 under basal conditions. We refer to this pattern as “low basal” in the less responsive strain. Quantitative analyses of pairwise comparisons indicate that 29% of the genes with decreased IL-4–induced gene expression were due to low basal expression, 36% had no differences before IL-4 stimulation (equal basal), and 35% were the result of a high basal expression level in the less responsive strain (Fig. 1G).

To investigate local versus distant effects of genetic variation on the differential responses to IL-4, we crossed C57 mice with the most genetically distinct SPRET mice to generate F1 offspring containing each parental chromosome. A total of 91.4% of parental-specific RNA-seq reads in the F1 strain are within twofold of their values in C57 and SPRET (blue data points) and considered to be due to local (cis) effects of genetic variation (Fig. 1H and fig. S1F), while only 2.8% were divergent between the parental strains but not in F1 BMDMs (green data points), indicating trans regulation. As NOD macrophages exhibited a broadly attenuated response to IL-4 on the level of gene expression, we followed the same strategy using F1 C57 × NOD macrophages. RNA-seq on IL-4–stimulated macrophages of F1 C57 × NOD macrophages showed strong convergence of expression of genes that were differentially regulated in the parental strain (green data points in fig. S1G), consistently a major contribution of trans regulation. To investigate the point at which this regulation occurs, we performed chromatin immunoprecipitation sequencing (ChIP-seq) for RNA polymerase II (Pol II) under control and IL-4–stimulated conditions. In contrast to mRNA levels, examination of the IL-4–dependent changes in gene body RNA Pol II indicated similar magnitude changes in all strains, including NOD (fig. S1H). These results suggest the presence of a transacting factor in NOD that acts downstream of transcription to attenuate mRNA levels. Collectively, these studies uncovered notable variation in the cell autonomous responses of BMDMs to IL-4 across these five strains, providing a powerful experimental system for investigating mechanisms by which natural genetic variation affects signal-dependent gene expression.

To investigate the impact of cis variation on putative transcriptional regulatory elements, we defined high-confidence IL-4–activated enhancers as intronic or intergenic open chromatin regions [based on assay for transposase-accessible chromatin using sequencing (ATAC-seq)] with at least 2.5-fold increase in H3K27ac (23) and RNA Pol II (24) after IL-4 treatment (fig. S2, A to D). In 24-hour IL-4–stimulated C57 BMDMs, 1093 regions exhibited a >2.5-fold increase in H3K27ac, whereas 441 regions exhibited a >2.5-fold decrease, corresponding to putative IL-4–activated and IL-4–repressed enhancers, respectively (Fig. 2A). Comparison of C57 enhancers to those of other strains under IL-4 treatment conditions revealed marked differences that scaled with the degree of genetic variation (Fig. 2B and fig. S2, E and F). We further subdivided these regions into “conventional enhancers” (Fig. 2C, blue) and “super enhancers” (Fig. 2C, orange), on the basis of the density distribution of normalized H3K27ac tag counts (25). Super enhancers represent regions of the genome that are highly enriched for cell-specific combinations of transcription factors and coregulators and control the expression of genes required for cellular identity and critical functions. In comparison to conventional enhancers, super enhancers exhibited significantly less variation in H3K27ac in response to IL-4 (Fig. 2D and fig. S2G). For example, IL-4 induction of the Ak2 super enhancer (Fig. 2E) is highly conserved between the five strains. In contrast, a typical example of strain specificity is provided by the conventional enhancers associated with the Msx3 gene. These enhancers are IL-4 inducible only in BALB, C57, and NOD and absent in PWK and SPRET macrophages (Fig. 2F).

We next compared the fractions of enhancers containing variants in strain-similar enhancers (<1.5-fold differences in H3K27ac between strains) to strain-differential enhancers at increasing levels of difference (fold differences >1.5 to >4; Fig. 2G). The fraction of enhancers containing variants at strain-similar enhancers ranged from 17 to 19% in the strains most similar to C57 (BALB and NOD) to 69 to 93% in the most genetically divergent strains (PWK and SPRET). As expected, the fraction of enhancers containing variants increased with increasing levels of difference, except for SPRET that may have reached a saturation of variation capacity (Fig. 2G). These findings are not only consistent with local variants affecting enhancer activity but also indicate that a substantial fraction of even strongly strain-differential IL-4–induced enhancers lack these variants, consistent with previous findings for strain-specific enhancers overall (14).

In an effort to distinguish silent variants from those affecting enhancer activities, we trained a DeepSEA convolutional neural network to classify enhancers as active or inactive under the 24-hour IL-4 condition on the basis of local sequence context (26). The training data consisted of enhancers active under IL-4 conditions (positive data) and random background (negative data). The area under the receiver operating characteristic curve (auROC) was 0.894 on test data. We then used DeepLIFT (27) to compute the importance score of each nucleotide on the basis of the model’s classification decision. Variants at positions with top importance scores within surrounding 300–base pair (bp) enhancer regions are hypothesized to affect enhancer activity. We considered variants residing in the top 20% of importance scores for each region as predicted functional variants. The Msx3 enhancer in Fig. 2F illustrates 4 predicted functional variants of 14 variants in PWK and SPRET (red dashed lines). By focusing on top-scoring variants rather than all local variants, we saw an expected overall decreased percentage of enhancers with top-scoring variants (Fig. 2H and fig. S2H). On the other hand, enrichment of predicted functional variants increases as a function of importance score threshold and is strongest for enhancers that show the highest differences across strains (Fig. 2I). This is true when considering all strains, including SPRET. These results reveal a quantitative impact of variants affecting enhancer under IL-4 treatment conditions and suggest the extent to which a deep learning approach can distinguish potentially functional variants from the silent variants.

Interpretation of effects of genetic variation on distal regulatory elements is facilitated by knowledge of cell-specific enhancer-promoter interactions (28). To identify connections of IL-4–responsive enhancers to target promoters, we performed HiChIP using an antibody to H3K4me3 (29) in C57 BMDMs under basal conditions and after 24 hours of IL-4 treatment. HiChIP interactions are exemplified in Fig. 3A at the Slc7a2 locus, a gene that becomes maximally activated after 24 hours of IL-4 treatment and connects primarily to an enhancer-like region within the Mtmr7 gene, which itself is expressed at negligible levels (Fig. 3A). Although we observed instances of IL-4–specific interactions (e.g., yellow loops), a differential interaction analysis was unable to identify significantly different interactions between basal and IL-4 conditions, supported by the high correlation of interaction intensity between the two conditions (fig. S3A). Moreover, enhancer-promoter interaction intensity did not correlate with IL-4–induced gene activity or the level of H3K4me3 at promoters (fig. S3, B and C). However, IL-4–activated promoters mostly interact with IL-4–activated enhancers (Fisher’s exact test, P = 2.2 × 10⁻¹⁶), and IL-4–repressed promoters strongly interact with IL-4–repressed enhancers (P = 1.2 × 10⁻¹⁵; Fig. 3B). These results suggest a preexistent and relatively stable landscape of enhancer-promoter interactions in macrophages, whose regulatory function was activated in response to IL-4.

Although the HiChIP assay is designed to capture promoter-enhancer interactions based on preferential occurrence of H3K4me3 at promoters, we also recovered 145,907 pairs of interactive enhancers (Fig. 3B), consistent with more than one enhancer being in local proximity of a target promoter. The H3K27ac correlations between interactive enhancers were significantly stronger than those between noninteractive enhancers (Fig. 3C and fig. S3D), consistent with their being functionally related. Noticeably, the closer enhancers, despite being noninteractive on the basis of our data, still have much stronger correlation than completely random enhancers, which might be due to more frequent contacts of nearby regions within the same interactive domain that were not captured by H3K4me3 HiChIP. On the basis of the high correlation of enhancer activity between connected enhancers, we hypothesized that enhancer-enhancer interactions could explain strain-differential enhancer when local genetic variants were absent (Fig. 2H). Among 224 interactive enhancers exhibiting a >4-fold difference in H3K27ac signal between BALB and C57 under the IL-4 condition, the original ~50% of strain-differential enhancers with predicted functional variants was further split into 20.5% that had top-scoring variants on both ends and 33.6% that had only local top-scoring variants (Fig. 3D, upper left). Depending on the strain comparison, an additional 8.2 to 19.9% of differential enhancers could be explained by genetic variants in interacting enhancers, indicating that enhancers may be affected by functional variants in other connected enhancers. Reducing the fold change requirement to twofold yielded a smaller proportion of strain-differential enhancers containing local variants overall but significantly increased the proportion having top-scoring variants on the connected ends only (15.4 to 26.5%, Fisher’s exact test P = 0.002 for BALB, 2.8 × 10⁻⁵ for NOD, 8.3 × 10⁻⁷ for PWK, and 3.3 × 10⁻¹⁰ for SPRET), suggesting that local variants have a stronger effect on inducing differential activation than variants at connected enhancers (Fig. 3D, bottom, and fig. S3E). Figure 3E illustrates an enhancer affected by genetic variants at the connected enhancer. The enhancer highlighted on the left is notably more active in C57 than NOD. This region lacks local variants in NOD but is connected to another enhancer ~100 kb away (highlighted on the right) containing multiple variants that are predicted to affect activity by deep learning. These findings are consistent with genetic variants at an enhancer influencing the activity states of other enhancers that lack local functional variants within the same connected network (30, 31).

IL-4 rapidly activates a set of enhancers, most of which exhibit maximal H3K27ac at 1 or 6 hours and returns to (near) basal levels by 24 hours (Fig. 4A, top three clusters) when most gene expression changes were found (fig. S1B). Others are long-lasting or become activated at later time points (Fig. 4A, bottom three clusters). De novo motif enrichment analysis of enhancers exhibiting >2.5-fold increase in H3K27ac and RNA Pol II at 1, 6, and 24 hours (fig. S2A) recovered a STAT6 motif as the most enriched motif for all time points (Fig. 4B). Motifs for the lineage-determining factors PU.1 and AP-1 (activator protein 1) family members were also recovered in all three classes of enhancers. Notably, an EGR2 motif was significantly enriched among enhancers induced at 24 hours.

As a genetic approach to identify functional transcriptional factor binding motifs, we assessed the quantitative impact of the genetic variation provided by the five different strains of mice on the IL-4 response of enhancers using the motif mutation analysis tool MAGGIE. MAGGIE associates changes of epigenomic features at homologous sequences (e.g., enhancer activation or enhancer repression) with motif mutations caused by genetic variation so that it can prioritize motifs that likely contribute to the regulatory function (32). This analysis identified more than a dozen motif clusters in which motif mutations were significantly associated with strain-differential IL-4–activated or IL-4–repressed enhancers (Fig. 4C and fig. S4A). The EGR motif was found as the top motif associated with enhancer activation at the 24-hour treatment time, as well as motifs of known SDTFs STAT6 and PPARγ and macrophage LDTFs PU.1, AP-1, and C/EBP (CCAAT/enhancer binding protein) (Fig. 4C). We also found Kruppel-like factor (KLF) motifs associated with IL-4 enhancer activation, which fits with increased KLF4 expression by IL-4 (fig. S4B), and an nuclear factor erythroid 2-related factor (NRF) motif associated with both enhancer activation (Fig. 4C).

The identification of STAT6 and PPARγ motif mutations as being functionally associated with strain-differential IL-4 activation is consistent with substantial prior work demonstrating the importance of these factors in regulating IL-4–dependent gene expression (21, 33). Out of the EGR family members, only Egr2 is expressed in unstimulated BMDMs and rapidly induced after IL-4 stimulation (Fig. 4D and fig. S4C). Egr2 has also been associated with late IL-4 enhancer activation in a recent study (22). Examination of the Egr2 locus indicates IL-4–induced binding of STAT6 and PPARγ to a set of upstream super enhancers that gain H3K27ac and RNA Pol II signal after IL-4 stimulation (fig. S4D). These super enhancers were observed in BMDMs of all five different strains (fig. S4E) that are strongly connected to the Egr2 promoter in C57 BMDMs as indicated by H3K4me3 HiChIP interactions. Overall, these findings suggest a functionally important role of EGR2 in contributing to IL-4–induced enhancer activation in BMDMs.

To incorporate EGR2 into a comprehensive analysis of the impact of genetic variation on IL-4 responses, we next performed ChIP-seq for EGR2 under basal and 24-hour IL-4 treatment conditions. This confirmed the prediction that mutations in EGR binding sites contribute to strain-differential enhancer activation by altering the binding of EGR2. An example is provided by the Btbd11 enhancer, which is IL-4 inducible in C57, but not in SPRET BMDMs (Fig. 4E). Consistent with the loss of EGR2 binding in SPRET, a C-to-T variant in SPRET mutated an EGR2 motif.

To study EGR2 binding and its effects on gene regulation over time, we additionally measured EGR2 binding at 1 and 6 hours after IL-4. We saw a marked expansion of the EGR2 cistrome after stimulation with IL-4 (Fig. 5, A and B). Most EGR2 binding sites had an increasing binding intensity over time and reached maximum values at 24 hours (Fig. 5A and fig. S5A). In contrast, STAT6 binding was strong immediately after 1 hour and slowly decreased in its binding intensity (fig. S5A). At 24 hours, there were over 20,000 newly gained EGR2 binding sites (Fig. 5B), which were associated with an increase in H3K27ac and RNA Pol II signals over time, supporting a major role of EGR2 in late enhancer activation (Fig. 5C and fig. S5B).

To extend these analyses, we crossed Egr2^flfl (34) (Egr2^WT) with LyzM-Cre⁺ mice to obtain LyzM-Cre⁺ Egr2^flfl (Egr2^MKO) mice. This resulted in efficient deletion of Egr2 in BMDMs (fig. S5, C and D). RNA-seq data from Egr2^WT and Egr2^MKO BMDMs indicated that at the mRNA level, EGR2 is regulating ~40% of the 24-hour IL-4–induced genes (Fig. 5D, fig. S5E, and table S3), many of which correspond to gene ontology terms related to cytokine production and signaling (Fig. 5D). At the regulatory level, we found that ~40% of the IL-4–induced enhancers at 24 hours had notably decreased activity in Egr2^MKO BMDMs (Fig. 5E, blue cluster, and fig. S5, F to H). In concordance, the IL-4 induction in H3K27ac was found to be decreased at IL-4–induced EGR2 binding sites in Egr2^MKO BMDMs (Fig. 5F). On the basis of a motif enrichment analysis of the EGR2-dependent IL-4–activated enhancers (blue cluster “deactivated”) with the EGR2-independent enhancers as background (green cluster “activated”), we found EGR2 as the most significantly enriched motif and SMAD3 and C/EBP motifs as the second and the third hits (Fig. 5E). When comparing the EGR2-independent activated enhancers (activated) to the EGR2-dependent ones (deactivated), the STAT6 motif was most significantly enriched (Fig. 5E), suggesting that STAT6 could work independently of EGR2 to maintain the late activation for a subset of enhancers. Together, these findings establish an essential role of EGR2 in IL-4–dependent enhancer activation and gene expression and are in full agreement with the recent studies of Daniel et al. (22).

Analysis of the genome-wide binding patterns of EGR2, STAT6, and PPARγ indicated intensive cobinding at IL-4–activated enhancers (fig. S6A). To achieve highly strain-differential responses to IL-4 that are observed at the level of gene expression (Fig. 1) and enhancer activation (Fig. 2), we hypothesize these factors to exert their transcriptional effects via correspondingly divergent genomic binding patterns. About 4000 EGR2 binding sites exhibited more than fourfold differences in normalized tag counts between C57 and BALB and about 10,000 between C57 and SPRET (Fig. 6A). Similar relationships are observed for STAT6 (fig. S6B). Within the strain comparisons of C57 to BALB and SPRET, C57-specific EGR2 binding sites are associated with stronger H3K27ac signal in C57 (Fig. 6B, green boxes). At these C57-specific EGR2 binding sites, H3K27ac signal is more strongly down-regulated in Egr2^MKO macrophages (Fig. 6B, orange boxes). An example of this is demonstrated in fig. S6C in which C57-specific binding of EGR2 is associated with an IL-4–induced increase in H3K27ac, which is absent in SPRET (where the EGR2 motif is disrupted; Fig. 4E) and is decreased in Egr2^MKO BMDMs.

The strain-differential binding patterns of SDTFs and LDTFs enabled motif mutation analysis to study the importance of motifs for SDTF and LDTF binding (Fig. 6C). As a validation, LDTF and SDTF binding depended on their own motifs (e.g., PU.1 motif mutation was significantly associated with PU.1 binding, indicating that when PU.1 binding is lost in one strain, it is often found that the PU.1 motif score is reduced in that strain compared to the other). In addition, mutations in the motifs of LDTFs PU.1, C/EBP, and AP-1 influence the binding of all LDTFs and SDTFs, which fits with earlier observations (13). We found that the PPAR motif is only significant for PPARγ binding, and likewise, the STAT6 motif is not associated with binding of other SDTFs or LDTFs but STAT6. We found that mutations in EGR2 motifs are significantly associated with binding of SDTFs STAT6 and PPARγ and LDTFs PU.1 and C/EBPβ under IL-4 conditions, but not under basal conditions. These analyses also provided evidence for functional roles of several additional transcription factors. Mutations in NRF motifs were strongly associated with the IL-4–dependent binding of all SDTFs and LDTFs. Both NRF1 and NRF2 are expressed in BMDMs (fig. S4B) and are involved in lipid metabolism and stress responses (35–37). Mutations in KLF motifs were strongly associated with EGR2 binding under both basal and IL-4 conditions. KLF2, KLF4, and KLF6 are expressed in BMDMs (fig. S4B), and KLF4 has previously been associated with anti-inflammatory roles in macrophages (38). Mutations in IRF motifs were moderately associated with IL-4–dependent STAT6 binding. Multiple IRFs, including IRF4, are expressed in BMDMs (fig. S6D), and IRF4 has previously been linked to macrophage polarization by IL-4 (39, 40).

A prediction emerging from the analysis results above is that EGR2 should have a small effect on the cobinding of SDTFs and LDTFs under basal conditions and a substantial effect following 24 hours of IL-4 treatment. To examine this prediction, we performed ChIP-seq for STAT6, PPARγ, PU.1, and C/EBPβ in Egr2^WT and Egr2^MKO BMDMs and evaluated their binding in the vicinity of IL-4–induced EGR2 binding sites. Deletion of Egr2 had little effect on PPARγ and STAT6 binding under basal conditions and a much greater effect following 24 hours of IL-4 treatment (Fig. 6D). As an example, in Egr2^MKO BMDMs, PPARγ and STAT6 binding was found decreased at the Mmp12 enhancer at sites where EGR2 normally binds (Fig. 6E). Similarly, PU.1 and C/EBPβ binding was more significantly affected by Egr2 deletion under the IL-4 condition than the basal condition (Fig. 6F). In concert, these findings provide evidence for collaborative interactions between PU.1, C/EBPs, AP-1, STAT6, PPARγ, and EGR2 as major drivers of late enhancer activation in response to IL-4. EGR2 is a strong collaborative factor as it promotes binding of LDTFs PU.1 and C/EBPβ and SDTFs STAT6 and PPARγ after IL-4 stimulation.

We next investigated the possibility that the mutational status of the dominant motifs recovered by MAGGIE analysis was sufficient to predict qualitative patterns of strain-differential responses of IL-4–induced enhancers. Following the classification of strain-differential mRNA responses (Fig. 1), we used H3K27ac to define three different categories of strain-differential IL-4–induced enhancers (Fig. 7A, left column): enhancers exhibiting lower levels of basal activity in the lowly induced strain (low basal); enhancers with a similar level of basal activity (equal basal); and enhancers in which a lack of IL-4–induced activity was associated with relatively higher basal activity compared with the more responsive strain (high basal). Using these criteria, we identified 760 low basal, 2797 equal basal, and 2013 high basal enhancers from all pairwise comparisons of the five strains that exhibited >2-fold differences in H3K27ac induction (Fig. 7B). The closest genes for enhancers of these three categories follow similar trends as observed for enhancer activity (fig. S7A). Low basal, equal basal, and high basal enhancers are exemplified by enhancers associated with the Treml2, Ripk2, and Cd36 genes, respectively (Fig. 7, C to E, and fig. S7, B to D).

Consideration of chromatin accessibility as determined by ATAC-seq further uncovered potential mechanisms that distinguished the three enhancer categories (Fig. 7A, right column). The enhancers in the low basal category showed low to absent basal ATAC signal in noninduced strains, suggesting a lack of LDTFs under the basal condition to preoccupy chromatin required for subsequent recruitment of SDTFs after IL-4 stimulation. In contrast, high basal enhancers exhibited a higher basal level of ATAC in noninduced strains compared with the induced strains (Fig. 7A, right column), suggesting stronger LDTF binding in noninduced strains under the basal condition. Different from the other categories, equal basal enhancers exhibited similar levels of chromatin accessibility under both basal and IL-4 conditions between comparative strains, suggesting that the recruitment of SDTFs might be the key determinant for the strain difference instead of basal LDTF binding.

To test the hypotheses above regarding the different determinants for the three categories of enhancers, we performed MAGGIE motif mutation analysis on each category of enhancers that contain motif mutations (fig. S7E). We found that mutations in motifs of LDTFs PU.1/ETS and C/EBP were associated with low basal enhancers and resulted in better motifs in induced strains, while mutations in motifs of SDTFs EGR, STAT6, PPAR, and NRF/MAF were associated with the equal basal category leading to better motifs in induced strains (Fig. 7F and fig. S7F). Mutations in EGR motifs were also associated with the low basal category, suggesting another role of EGR2 as a strong collaborative factor under the IL-4 condition, which is supported by the notable decrease in open chromatin under IL-4 conditions after deletion of Egr2 (fig. S5F). Of particular interest, the high basal category of enhancers was most strongly associated with negative significance scores for LDTF PU.1, C/EBP, and AP-1 as well as NRF/MAF, meaning higher motif affinity in noninduced strains (Fig. 7F).

We validated these findings with our ChIP-seq data by examining the binding profiles of PU.1, C/EBPβ, STAT6, PPARγ, and EGR2 in three categories of enhancers. In low basal enhancers, we saw significantly reduced binding of PU.1 and C/EBPβ in noninducible strains under both basal and IL-4 conditions (Fig. 7G and fig. S7G). This pattern was accompanied by significantly weaker binding of SDTFs STAT6, EGR2, and PPARγ after IL-4 stimulation (Fig. 7H and fig. S7H). The example in Fig. 7C showed the absence of C/EBPβ binding in NOD under the basal condition likely due to two local variants at high-scored positions according to DeepLIFT that together mutated a C/EBP motif. Upon IL-4 stimulation, neither C/EBPβ nor EGR2 was further recruited. For equal basal enhancers, we found that PU.1 and C/EBPβ binding was similar under basal conditions in induced and noninduced strains (Fig. 7G and fig. S7G). Upon IL-4 stimulation, the induced strains displayed significantly stronger binding of SDTFs STAT6, EGR2, and PPARγ (Fig. 7H and fig. S7H). In the example in Fig. 7D, STAT6 binding was strongly induced by IL-4 at the Ripk2 enhancer in PWK but was absent in SPRET. Despite the clear difference in STAT6 binding, none of the local variants between the two strains were predicted functional when using a neural network model trained with random genomic backgrounds. To better capture the sequence patterns relevant for enhancer activation, we retrained neural networks using noninduced enhancers as the background, which emphasized a relatively divergent set of k-mers and focused less on those matched with LDTF motifs (fig. S7I). As a result, our retrained model assigned a high DeepLIFT score to one of the nucleotides in a STAT6 motif that was mutated by a variant in SPRET (Fig. 7D). For high basal enhancers, we found stronger binding of not only the LDTFs PU.1 and C/EBPβ (Fig. 7G and fig. S7G) but also the SDTFs STAT6 and PPARγ (Fig. 7H and fig. S7H) in noninduced strains under basal conditions. For example, high basal levels of C/EBPβ and STAT6 binding were observed at the Cd36 enhancer in NOD mice (Fig. 7E). The only local variant in PWK was at a predicted functional position and mutated a C/EBP motif likely causing the low basal C/EBPβ binding in PWK. In concert, these analyses validated the importance of LDTF motif mutations as primary determinants of differential enhancer activation in low basal and high basal enhancers, while also demonstrating the expected consequences of SDTF motif mutations in determining strain-differential activation of equal basal enhancers (Fig. 7I).

To investigate the influence of genetic variation on signal-dependent gene expression, enhancer activation, and transcription factor binding, we performed RNA-seq, ATAC, and ChIP-seq to study the responses of macrophages derived from five different inbred mouse (C57, BALB, NOD, PWK, and SPRET) strains to the anti-inflammatory cytokine IL-4 (Fig. 1A).

Female and male breeder mice for C57, BALB, NOD, PWK, and SPRET mice were purchased from the Jackson laboratory. F1 C57 × SPRET mice were crossed, and Egr2^fl/fl mice were generously donated by Drs. V. Lazarevic and L. Warren (National Institutes of Health) and crossed to LyzM-Cre mice (the Jackson laboratory) to achieve myeloid-specific targeted deletion of Egr2. Mice were housed at the University of California San Diego animal facility on a 12-hour/12-hour light/dark cycle with free access to normal chow food and water. All animal procedures were in accordance with University of California San Diego research guidelines for the care and use of laboratory animals. Eight- to 12-week-old healthy female mice were used for all our experiments.

Femur, tibia, and iliac bones from the different mouse strains were flushed with Dulbecco’s modified Eagle’s medium (DMEM) high glucose (Corning), and red blood cells were lysed using red blood cell lysis buffer (eBioscience). After counting, 20 million bone marrow cells were seeded per 15-cm nontissue culture plates in DMEM high glucose (50%) with 20% fetal bovine serum (FBS; Omega Biosciences), 30% L929 cell–conditioned laboratory-made media [as source of macrophage colony-stimulating factor (M-CSF), as described before (14)], penicillin/streptomycin + l-glutamine (100 U/ml; Gibco), and amphotericin B (2.5 μg/ml; HyClone). After 4 days of differentiation, mouse M-CSF (16.7 ng/ml; Shenandoah Biotechnology) was added to the media. After an additional 2 days of culture, adherent cells were scraped and subsequently seeded onto tissue culture–treated petri dishes in DMEM containing 10% FBS, penicillin/streptomycin + l-glutamine (100 U/ml), amphotericin B (2.5 μg/ml), and M-CSF (16.7 ng/ml). Macrophages were left untreated or treated with mouse recombinant IL-4 (20 ng/ml; PeproTech) for 1, 6, or 24 hours.

Cells were fixed with Cytofix/Cytoperm Buffer (BD Biosciences, BD554714) for 10 min at room temperature. Cytofix/Cytoperm buffer was removed, and cells were washed twice with Hanks’ balanced salt solution containing 2% bovine serum albumin (BSA) and 1 mM EDTA. Cells were kept in permeabilization/wash buffer (BD Biosciences, BD554714) for 1 hour at 4°C or until the experiment was performed. Fixed cells were blocked using 3% BSA, 0.1% Triton–phosphate-buffered saline (PBS) for 30 min at room temperature and then with 1/200 of the EGR2 antibody (Abcam) overnight at 4°C. The next day, cells were washed with 0.1% Triton-PBS and incubated with 1/200 donkey anti-rabbit 555 (Thermo Fisher Scientific, no. A31572) secondary antibody and phalloidin (Abcam, ab176759) for staining actin filaments, and nuclei were counterstained with 4′,6-diamidino-2-phenylindole. After washing with 0.1% Triton-PBS, slides were mounted with ProLong Gold Antifade Reagent (Life Technologies, no. 10144). Images were taken using a Leica SP8 with light deconvolution microscope.

Total RNA was isolated from cells and purified using RNA Direct-zol microprep columns and ribonuclease (RNase)–free deoxyribonuclease digestion according to the manufacturer’s instructions (Zymo Research). Sequencing libraries were prepared in biological replicates from polyadenylate [poly(A)]–enriched mRNA as previously described (14). Libraries were polymerase chain reaction (PCR)–amplified for 9 to 14 cycles, size-selected using tris-boric acid-EDTA (TBE) gels or one-sided 0.8× AMPure cleanup, quantified by the Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific), and 75-bp single-end–sequenced on a HiSeq 4000 or NextSeq 500 (Illumina).

For histone marks, PU.1, C/EBPβ, and RNA Pol II ChIP-seqs, culture media were removed and plates were washed once with PBS and then fixed for 10 min with 1% formaldehyde (Thermo Fisher Scientific) in PBS at room temperature; reaction was then quenched by adding glycine (Thermo Fisher Scientific) to 0.125 M. For STAT6, PPARγ, and EGR2 ChIP-seq, cells were cross-linked for 30 min with 2 mM DSG (Pierce) in PBS at room temperature. Subsequently, cells were fixed for 10 min with 1% formaldehyde at room temperature, and the reaction was quenched with 0.125 M glycine. After fixation, cells were washed once with cold PBS and then scraped into supernatant using a rubber policeman and pelleted for 5 min at 400g at 4°C. Cells were transferred to Eppendorf DNA LoBind tubes and pelleted at 700g for 5 min at 4°C, snap-frozen in liquid nitrogen, and stored at −80°C until ready for ChIP-seq protocol preparation.

ChIP was performed in biological replicates as described previously (44). Samples were sonicated using a probe sonicator in 500-μl lysis buffer [10 mM tris-HCl (pH 7.5), 100 mM NaCl, 1 mM EDTA, 0.5 mM EGTA, 0.1% deoxycholate, 0.5% sarkozyl, and 1× protease inhibitor cocktail]. After sonication, 10% Triton X-100 was added to 1% final concentration, and lysates were spun at full speed for 10 min. One percent was taken as input DNA, and immunoprecipitation was carried out overnight with 20-μl Protein A Dynabeads (Invitrogen) and 2-μg specific antibodies for PU.1 (Santa Cruz Biotechnology, sc-352X), H3K4me2 (Millipore, 07-030), H3K4me3 (Millipore, 04-745), H3K27ac (Active Motif, 39135), RNA Pol II (GeneTex, GTX102535), STAT6 (Santa Cruz Biotechnology, sc-374021), EGR2 (Abcam, ab43020), and C/EBPβ (Santa Cruz Biotechnology, sc-150). Beads were washed three times each with wash buffer I (20 mM tris-HCl, 150 mM NaCl, 0.1% SDS, 1% Triton X-100, and 2 mM EDTA), wash buffer II (10 mM tris-HCl, 250 mM LiCl, 1% IGEPAL CA-630, 0.7% Na-deoxycholate, and 1 mM EDTA), tris-EDTA (TE) 0.2% Triton X-100, and TE 50 mM NaCl and subsequently resuspended 25-μl 10 mM tris-HCl (pH 8.0) and 0.05% Tween 20, and sequencing libraries were prepared on the Dynabeads as described below.

For PPARγ ChIP-seq, fixed cells were lysed in 500-μl radioimmunoprecipitation assay (RIPA) lysis buffer [20 mM tris-HCl (pH 7.5), 1 mM EDTA, 0.5 mM EGTA, 0.1% SDS, 0.4% Na-deoxycholate, 1% NP-40 alternative, 0.5 mM dithiothreitol (DTT), and 1× protease inhibitor cocktail (Sigma-Aldrich)] and chromatin was sheared using a probe sonicator. One percent was taken as input DNA, and immunoprecipitation was carried out overnight with 20-μl Protein A Dynabeads (Invitrogen) and 2 μg of both PPARγ antibodies (Santa Cruz Biotechnology, sc-271392 and sc-7273). Beads were then collected using a magnet and washed with 175-μl ice cold buffer as indicated by incubating samples on ice for 3 min: three times with RIPA wash buffer [20 mM tris-HCl (pH 7.5), 1 mM EDTA, 0.5 mM EGTA, 0.1% SDS, 0.4% Na-deoxycholate, 1% NP-40 alternative, 0.5 mM DTT, and 1× protease inhibitor cocktail (Sigma-Aldrich)], six times with LiCl wash buffer [10 mM tris-HCl (pH 7.5), 250 mM LiCl, 1 mM EDTA, 0.7% Na-deoxycholate, 1% NP-40 alternative, and 1× protease inhibitor cocktail (Sigma-Aldrich)], twice with TET [10 mM tris-HCl (pH 8.0), 1 mM EDTA, 0.2% Tween 20, and 1× protease inhibitor cocktail (Sigma-Aldrich)], and once with TE-NaCl [10 mM tris-HCl (pH 8.0), 0.1 mM EDTA, 50 mM NaCl, and 1× protease inhibitor cocktail (Sigma-Aldrich)]. Bead complexes were resuspended in 25-μl TT [10 mM tris-HCl (pH 8.0) and 0.05% Tween 20], and sequencing libraries were prepared on the Dynabeads as described below.

ChIP libraries were prepared while bound to Dynabeads using a NEBNext Ultra II Library preparation kit (NEB) using half reactions. DNA was polished, poly(A)–tailed and ligated after which dual UDI (IDT) or single (Bioo Scientific) barcodes were ligated to it. Libraries were eluted and cross-links reversed by adding to the 46.5-μl NEB reaction, 16-μl water, 4-μl 10% SDS, 4.5-μl 5 M NaCl, 3-μl 0.5 M EDTA, 4-μl 0.2 M EGTA, 1-μl RNAse (10 mg/ml), and 1-μl proteinase K (20 mg/ml), followed by incubation at 55°C for 1 hour and 75°C for 30 min in a thermal cycler. Dynabeads were removed from the library using a magnet, and libraries were cleaned up by adding 2-μl SpeedBeads 3 EDAC (Thermo Fisher Scientific) in 124-μl 20% polyethylene glycol, molecular weight 8000/1.5 M NaCl, mixing well and then incubating at room temperature for 10 min. SpeedBeads were collected on a magnet and washed two times with 150-μl 80% ethanol for 30 s. Beads were collected, and ethanol was removed following each wash. After the second ethanol wash, beads were air-dried and DNA-eluted in 12.25-μl 10 mM tris-HCl (pH 8.0) and 0.05% Tween 20. DNA was amplified by PCR for 14 cycles in a 25-μl reaction volume using NEBNext Ultra II PCR master mix and 0.5 μM each Solexa 1GA and Solexa 1GB primers. Libraries were size-selected using TBE gels for 200 to 500 bp, and DNA was eluted using gel diffusion buffer [500 mM ammonium acetate (pH 8.0), 0.1% SDS, 1 mM EDTA, and 10 mM magnesium acetate] and purified using ChIP DNA Clean & Concentrator (Zymo Research). Sample concentrations were quantified by the Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific) and 75-bp single-end–sequenced on HiSeq 4000 or NextSeq 500 (Illumina).

Approximately 80K cells were lysed in 50-μl room temperature ATAC lysis buffer [10 mM tris-HCl, (pH 7.4), 10 mM NaCl, 3 mM MgCl₂, and 0.1% IGEPAL CA-630], 2.5-μl DNA Tagmentation Enzyme mix (Nextera DNA Library Preparation Kit, Illumina) was added. The mixture was incubated at 37°C for 30 min and subsequently purified using the ChIP DNA purification kit (Zymo Research) as described by the manufacturer. DNA was amplified using the Nextera Primer Ad1 and a unique Ad2.n barcoding primer using NEBNext High-Fidelity 2× PCR MM for 8 to 14 cycles. PCRs were size-selected using TBE gels for 175 to 350 bp, and DNA was eluted using gel diffusion buffer [500 mM ammonium acetate (pH 8.0), 0.1% SDS, 1 mM EDTA, and 10 mM magnesium acetate] and purified using ChIP DNA Clean & Concentrator (Zymo Research). Samples were quantified by the Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific) and 75-bp single-end–sequenced on HiSeq 4000 or NextSeq 500 (Illumina).

For H3K4me3 HiChIP, 10 million formaldehyde cross-linked cells per condition in biological replicates were used. HiChIP was performed as described before (29). In our experiments, 375 U of Mbo I (NEB, R0147M) restriction enzyme was used for chromatin digestion. Shearing was performed in three Covaris microtubes per sample and using the following parameters on a Covaris E220 (Fill Level = 6, Duty Cycle = 5, PIP = 140, Cycles/Burst = 200, Time = 200 s). H3K4me3 IP was performed using 7.5 μg of antibody (Millipore, 04-745). Final PCR was performed using NEBNext High-Fidelity PCR MM and Nextera general Primer Ad1 and specific Nextera Primer Ad2.n. PCR product was run on a TBE gel (Invitrogen), and libraries were size-selected from 250 to 700 bp and cleaned up using 150-μl gel diffusion buffer [500 mM ammonium acetate (pH 8.0), 0.1% SDS, 1 mM EDTA, and 10 mM magnesium acetate] and purified using ChIP DNA Clean & Concentrator (Zymo Research). Samples were quantified by the Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific) and 75-bp paired-end–sequenced on a NextSeq 500 (Illumina).

Custom genomes were generated for BALB, NOD, PWK, and SPRET mice from the C57 or mm10 genome as before (14) using MMARGE v1.0 (45) and the variant call format (VCF) files from the Mouse Genomes Project (46). Data generated from different mouse strains were first mapped to their respective genomes using STAR v2.5.3 (47) for RNA-seq data or bowtie2 v2.2.9 (48) for ATAC-seq, ChIP-seq, and HiChIP data. Then, the mapped data were shifted to the mm10 genome using the MMARGE v1.0 “shift” function (45) for downstream comparative analyses.

For each strain, a differential gene expression analysis was performed to compare IL-4 to basal with Limma Voom (51). A linear model was fit for all five differential comparisons at once, and 1912 genes that were significant with q value below 0.05 and an absolute fold change of 1.5 in any comparison were included in a WGCNA (52). WGCNA was performed with a softpower value of 20, and a signed network was generated. Modules were cut with minimum module size of 50 and a cut height of 0.999 including PAM-stage. Nine modules were detected of which two genes were part of the gray (nonconnected) module that was subsequently excluded. Module Eigengenes (ME) were calculated and visualized using the verbose-boxplots function that also performed a Kruskal-Wallis significance test to test whether all ME values belong to the same distribution, and all modules were significantly different between conditions (all P values below <0.0012). Two modules exhibited consistent differential expression between IL-4 and notx across strains, while the other six modules were most prominently influenced in a strain-specific manner. Modules were annotated with Metascape (50).

On the basis of the HOMER tag directories created from mapped sequencing data, the reproducible ATAC-seq and transcription factor ChIP-seq peaks were identified by using HOMER to call unfiltered 200-bp peaks (parameters -L 0 -C 0 -fdr 0.9 -size 200) and running IDR v2.0.3 on replicates of the same sample with the default parameters (53). The levels of histone modifications and RNA Pol II were quantified within ±500 bp around the centers of reproducible ATAC-seq peaks using HOMER annotatePeak.pl with parameters “-size -500,500 -norm 1e7.” The transcription factor binding intensities were quantified within ±300 bp around the identified ChIP-seq peaks using parameters “-size -150,150 -norm 1e7.” For comparisons across multiple samples (e.g., different time points, mouse strains, and transcription factors), we merged the set of peaks first using HOMER mergePeaks “-d given” before quantifying the features above. To visualize the average profile of a dataset around a certain set of peaks, we used HOMER annotatePeaks.pl with parameters “-norm 1e7 -size 4000 -hist 20” to help compute the histograms of 20-bp bins within ±2000 bp regions.

IL-4–responsive enhancers were identified by the strong fold changes of H3K27ac and RNA Pol II at intergenic or intronic open chromatin. Reproducible ATAC-seq peaks called from each mouse strain for the basal and IL-4 conditions were first merged and then annotated for genomic positions and the enrichment of H3K27ac and RNA Pol II within ±500 bp using HOMER v4.11.1. On the basis of the genomic annotations from HOMER annotatePeaks.pl, we classified regions at promoter-transcription start sites (TSS) as promoters and regions at intergenic or intronic positions as enhancers. Regions with less than 16 normalized tags of H3K27ac or less than 8 normalized tags of RNA Pol II were filtered out. For the remaining promoters and enhancers, we computed the fold changes of the normalized tags of H3K27ac and RNA Pol II between basal and IL-4 conditions for each mouse strain. Regions were called IL-4–induced or IL-4–repressed if there were at least 2.5-fold increases or decreases, respectively, from basal to IL-4 state for both histone markers. Regions with less than 1.4-fold changes were called neutral elements.

We used ROSE to call super enhancers for the five mouse strains (25). The active enhancers were first merged within each strain for both basal and IL-4 conditions to obtain a set of starting conventional enhancers. Then, the ROSE algorithm was run for each strain on the mapped H3K27ac ChIP-seq data with the parameter “-t 2500” to exclude TSS. The overall activity of a super enhancer was quantified by the H3K27ac ChIP-seq read counts within the entire identified super enhancer region.

Significant interactions captured by HiChIP-seq were overlapped with previously identified active promoters and enhancers for the five mouse strains using HOMER mergePeaks “-d 2500” to identify three categories of interactive pairs: enhancer-enhancer, enhancer-promoter, and promoter-promoter. Enhancer-promoter interactions have enhancers on one end and promoters on the other end, while enhancer-enhancer or promoter-promoter interactions are the linked pairs of enhancers or promoters, respectively. We ended up with 145,907 enhancer-enhancer interactions, 81,411 enhancer-promoter interactions, and 10,710 promoter-promoter interactions. To better understand the regulatory landscape associated with IL-4 stimulation, we subsequently focused on enhancer-promoter interactions that contained IL-4–induced, IL-4–repressed, and/or IL-4–neutral promoters on one end and IL-4–induced, IL-4–neutral, and/or IL-4–repressed enhancers on the other end, and quantified the number of interactions between these possible promoter-enhancer combinations in nine categories as a contingency table. Fisher’s exact test was applied to the contingency table to determine whether any of the categories were significantly different for three comparisons of interest: IL-4–induced enhancer/promoter interactions versus noninduced enhancer/promoters; IL-4–repressed enhancer/promoter interactions versus nonrepressed enhancer/promoters; and IL-4–induced enhancer/promoter interactions versus IL-4–repressed enhancer/promoter interactions. For enhancer-enhancer interactions, we preselected enhancers that have at least fourfold difference in H3K27ac ChIP-seq tags between any two strains under the 24-hour IL-4 condition to obtain a set of strongly strain-differential enhancers. We then computed the Pearson correlation of H3K27ac tags across the five strains for every pair of interactive enhancers among the preselected set. To obtain noninteractive enhancers, we either randomly paired preselected enhancers on the same chromosome (same-chromosome random enhancers) or looked for enhancers within certain distances but not connected based on our data (distance-matched random enhancers).

Genetic variation between C57 and the other four strains at strain-differential enhancers was extracted using MMARGE annotate_mutations (45), which was based on the VCF files from the Mouse Genomes Project (46). Variants were searched within ±150 bp around the centers of enhancers. At least one genetic variant from the comparative strain needs to be present within the search area for such enhancer to be counted as having variants.

Among the strain-differential IL-4–induced enhancers as described above, we further split them into three categories based on the level of H3K27ac under the basal condition in noninduced strains. High basal enhancers have more than twofold stronger H3K27ac in noninduced strains, while low basal enhancers have more than twofold stronger H3K27ac in induced strains (lower basal H3K27ac in noninduced strains). Equal basal enhancers are those in between.

All sequencing data have been made available by deposition in the Gene Expression Omnibus (GEO) database: GSE159630. The UCSC genome browser was used to visualize sequencing data. The codes for neural network model training and interpretation are available on our Github repository: https://github.com/zeyang-shen/macrophage_IL4Response.

Two independent groups were tested using Mann-Whitney U test for medians and using Levene’s test for variance. Gene expression comparisons were reported by adjusted P values (i.e., q values) from DESeq2 (49). Enrichment was computed by odds ratio and tested by Fisher’s exact test. Effect sizes were reported by Cohen’s d. All gene expression data are displayed as means with 95% confidence interval. All data distributions are shown with means, 25th percentiles, and 75th percentiles.