Alternative Activation of Macrophages Is Accompanied by Chromatin Remodeling Associated with Lineage-Dependent DNA Shape Features Flanking PU.1 Motifs

Abstract

IL-4 activates macrophages to adopt distinct phenotypes associated with clearance of helminth infections and tissue repair, but the phenotype depends on the cellular lineage of these macrophages. The molecular basis of chromatin remodeling in response to IL-4 stimulation in tissue-resident and monocyte-derived macrophages is not understood. In this study, we find that IL-4 activation of different lineages of peritoneal macrophages in mice is accompanied by lineage-specific chromatin remodeling in regions enriched with binding motifs of the pioneer transcription factor PU.1. PU.1 motif is similarly associated with both tissue-resident and monocyte-derived IL-4-induced accessible regions but has different lineage-specific DNA shape features and predicted cofactors. Mutation studies based on natural genetic variation between C57BL/6 and BALB/c mouse strains indicate that accessibility of these IL-4-induced regions can be regulated through differences in DNA shape without direct disruption of PU.1 motifs. We propose a model whereby DNA shape features of stimulation-dependent genomic elements contribute to differences in the accessible chromatin landscape of alternatively activated macrophages on different genetic backgrounds that may contribute to phenotypic variations in immune responses.

FIGURE 1.

IL-4 stimulation leads to remodeling of open chromatin landscape in peritoneal macrophages. (A) PCA scores of individual ATAC-seq samples. PCA was performed using rlog-transformed ATAC-seq read counts of 30,856 regions with high variance (only regions with variance interquartile range >0.5 were retained) (n = 4–6 mice per macrophage population). Data points represent independent biological replicates. (B) Genome browser views of representative (boxed) constitutively accessible and IL-4–induced regions. (C) The contributions of individual accessible regions to principal components 1 and 2 are represented in the PCA loadings plot. Each data point is color-coded based on the direction of its IL-4 dependency. Hence, IL-4–induced regions (red) are highly associated with IL-4–stimulated macrophages, whereas IL-4–repressed regions (yellow) are highly associated with nonstimulated macrophages. (D) Venn diagrams indicating the number of constitutively and IL-4 induced accessible regions in (top) resident and (bottom) monocytic macrophages. We compared (E) enrichment levels for different types of genomic elements, (F) distance from a closest IL-4–induced gene, and (G) G/C content between constitutively accessible and IL-4–induced regions in AAM^mono and AAM^res, respectively. G/C content information is represented in two different ways: percentage of G/C bases in an accessible region (G, left panel) and CpG island enrichment for a given group of accessible regions (G, right panel). Number of IL-4–induced regions = 1572 in AAM^res and 1462 in AAM^mono; number of constitutively accessible regions = 8061 in AAM^res and 14,045 in AAM^mono. Enrichment p values are from binomial test, whereas two-class comparison p values are from two-sided Mann–Whitney test.

FIGURE 2.

IL-4–induced regions are associated with PU.1, KLF, and AP-1 motifs. (A) Heatmap visualizing the macrophage-specific IL-4–dependent regions. Each row represents one of the 2855 IL-4–dependent regions and each column a unique sample. Values are rlog-transformed, batch-subtracted read counts, scaled using a z-score transformation for each region. (B) Motifs discovered de novo from IL-4–induced regions in AAM^res and AAM^mono. (C) Frequency of IL-4–induced peaks delineated by the presence of de novo PU.1, KLF, and AP-1 motifs. (D) Twenty highly expressed TF genes related to the de novo motifs discovered from the IL-4–induced regions. Values are log₂ intensity values of microarrays (12). (E) Clustering analysis of de novo motifs and macrophage-specific motifs identified using an overrepresentation approach from the IL-4–induced regions of AAM^res (left) and AAM^mono (right). Asterisks indicate macrophage-specific motifs uniquely identified via the overrepresentation approach. Only macrophage-specific motifs with log₂ p value < −15 are included in this visualization.

FIGURE 3.

PU.1 motifs in AAM^res and AAM^mono are associated with macrophage-specific sequence features. (A) Comparison of PU.1 motif scores derived by FIMO in IL-4–induced regions of AAM^res versus AAM^mono, with horizontal lines in the violin plots representing values at 25th, 50th, and 75th percentiles. The p value is from a two-sided Mann–Whitney U test. Number of IL-4–induced regions = 1572 in AAM^res and 1462 in AAM^mono. (B) Macrophage-specific TF motifs found within IL-4–induced PU.1 motifs ±25 bp regions, represented using log₂ odds ratio values (two-sided Fisher test, adjusted p values < 0.1). Motifs are summarized as TF families and the specific TF with the maximum absolute log₂ odds ratio value is stated in parenthesis. When TF family annotation was not available, the log₂ odds ratio of the specific TF itself is used. (C) Average of DNA shape features at IL-4–induced PU.1 regions of AAM^mono and AAM^res. Scatter plots are centered on IL-4–induced PU.1 motifs (x-axes), with the solid lines representing Loess fit of predicted shape values at single-nucleotide resolution (dots). Boxplots are average predicted shape values over the PU.1 motif ±25 bp windows, and p values are from two-sided Kolmogorov–Smirnov test. (D) Predicted DNA shape at the eighth base pair of PU.1 motif of AAM^res and AAM^mono. Frequency distributions are represented by smoothed kernel density estimates. HelT, helical twist; MGW, minor groove width; ProT, propeller twist.

FIGURE 4.

Accessibility of IL-4–induced PU.1 regions can be altered through DNA shape change induced by sequence mutation without disruption of the PU.1 motif. Genomic regions containing PU.1 motif with flanking SNPs, selected from (A) AAM^res, (B) AAM^mono, and (C) nonstimulated BMDMs from Ref. . In all examples, sequences represent a unique PU.1 motif region, in which the PU.1 motif is underlined and shaded in gray, whereas SNPs are highlighted with red font. Line graphs of predicted shape values are centered on the PU.1 motif (sequences between the vertical lines). Genome browser tracks illustrate the strain-specific accessibility of the region containing the matching sequence (shaded in gray). ATAC-seq reads are visualized in (A) and (B), whereas ChIP-seq reads are visualized in (C). Boxplots in (A) and (B) represent size factor–normalized read counts, with p values are from DESeq2 (n = 4–8 mice per group) and adjusted by Benjamini–Hochberg procedure. RGMb expression in (C) is expressed as transcript per kilobase million, with values directly obtained from Ref. .

FIGURE 5.

Effects of SNP variants on chromatin accessibility can be predicted by the extent of DNA shape change and association with IL-4 inducibility. (A) Frequency of different types of sequence variants in strain-specific IL-4–induced PU.1 regions and (B) all IL-4–induced PU.1 regions (bottom). Number of strain-specific IL-4–induced PU.1 regions = 58 in AAM^mono and 37 in AAM^res; number of strain-common IL-4–induced PU.1 regions = 754 in AAM^mono and 828 in AAM^res. (C) Comparison of the amount of SNP-induced DNA shape change in 922 SNP(s)-containing strain-specific and 4776 SNP(s)-containing strain-common accessible regions, respectively. The shift in DNA shape for each sequence is quantitated using Euclidean distance, in which a larger Euclidean distance value indicates a larger shift in shape. The p values are from Mann–Whitney tests. (D) Coefficient values representing the relative contribution of each significant predictor in distinguishing between strain-specific and strain-common accessible regions. (E) Enrichment values of IL-4–induced and constitutively accessible regions in strain-specific regions for AAM^mono and AAM^res. Enrichment p value is based on binomial test. (F) Scatter plots illustrating the positive correlations between SNP frequency and amount of DNA shape change undergone by a specific region. Red, strain-specific regions; blue, strain-common regions. Euclidean distances in this study are scaled across each DNA shape feature for visualization purpose. HelT not shown as not correlated. Correlation coefficient and p values are based on the Spearman test using unscaled Euclidean distances.

FIGURE 6.

AAMs from C57BL/6 and BALB/c are functionally distinct. (A) PCA of 7431 genes with high variance (only genes with variance interquartile range >0.5 were retained). Data points represent independent biological replicates. (B) Venn diagrams indicating the number of genes that were commonly and uniquely identified as significantly differential (FDR < 0.1) in different comparisons—(left) macrophage-specific genes in C57BL/6 and BALB/c AAMs and (right) strain-specific genes in AAM^res and AAM^mono. (C) Enrichment values from Ingenuity Pathway Analysis visualized as—log₁₀ p value for the four different groups of genes—1) BALB/c specific in AAM^res, 2) C57BL/6 specific in AAM^res, 3) BALB/c specific in AAM^mono, and 4) C57BL/6 in AAM^mono. Only the top 10 pathways (as defined by enrichment p values) are included in this visualization. Specific pathways are highlighted for clarity. (D) Representative flow cytometric analysis of F4/80 and PD-L2 surface expressions in AAM^mono of C57BL/6 and BALB/c mice. Boxplots show frequency of CD11b⁺ F4/80⁺ PD-L2⁺ singlet, live cells. The p value is based on a two-sided unpaired t test. (E) Expression of the Pdcd1lg2 gene in AAM^mono of C57BL/6 versus BALB/c mice, represented by size-factor normalized read counts. The p value is from DESeq2 and adjusted by Benjamini–Hochberg procedure. (F) Representative flow cytometric analysis of F4/80 and MHCII surface expressions in AAM^mono of C57BL/6 and BALB/c mice. Boxplots show frequency of CD11b⁺ F4/80 ⁺ MHCII⁺ singlet, live cells. The p value is based on a two-sided unpaired t test. (G) Expression values of all MHCII genes in AAM^mono of C57BL/6 versus BALB/c mice, represented by size-factor normalized read counts. The p values are from DESeq2 and adjusted by Benjamini–Hochberg procedure. Hinges of all boxplots correspond to values of the 25th, 50th, and 75th percentiles, whereas boxplot whiskers extend to no more than 1.5 × interquartile range, beyond which the outlier data points will be plotted individually. Transcriptional profiling analysis: n = 8 AAM^res (C57BL/6), 4 AAM^res (BALB/c), 7 AAM^mono (C57BL/6), and 6 AAM^mono (BALB/c). Flow cytometric analysis: n = 6 AAM^mono (C57BL/6) and 6 AAM^mono (BALB/c).