Genetic variation in IL-4 activated tissue resident macrophages alters the epigenetic state to determine strain specific synergistic responses to LPS

TL;DR

Genetic differences between mouse strains affect how macrophages respond to IL-4 and LPS, leading to distinct immune responses.

Contribution

The study reveals how genetic variation alters epigenetic reprogramming in macrophages, causing strain-specific immune responses.

Findings

C57BL/6 macrophages show greater transcriptional and epigenetic responses to IL-4 compared to BALB/c.

IL-4-activated C57BL/6 macrophages exhibit enhanced synergy with LPS compared to BALB/c.

Transcriptional differences and synergy are cell-intrinsic, as shown by chimeric mouse experiments.

Abstract

How macrophages in the tissue environment integrate multiple stimuli will depend on the genetic background of the host, but this is poorly understood. Here, we investigated C57BL/6 and BALB/c strain specific in vivo IL-4 activation of tissue-resident macrophages (TRMs) from the peritoneal cavity. C57BL/6 TRMs are more transcriptionally responsive to IL-4 stimulation, with a greater association of induced genes with super enhancers, induced enhancers, and topologically associating domains (TAD) boundaries. IL-4-directed epigenomic remodeling revealed BL/6 specific enrichment of NF-κB, IRF, and STAT motifs. Additionally, IL-4-activated BL/6 TRMs demonstrated an augmented synergistic response upon in vitro lipopolysaccharide (LPS) exposure compared to BALB/c TRMs, despite naïve BALB/c TRMs displaying a more robust transcriptional response to LPS than naïve BL/6 TRMs. Single-cell RNA sequencing (scRNA-seq) analysis of mixed bone marrow chimeric mice indicated that transcriptional differences between BL/6 and BALB/c TRMs, and synergy between IL-4 and LPS, are cell intrinsic within the same tissue environment. Hence, genetic variation alters IL-4-induced cell intrinsic epigenetic reprogramming resulting in strain specific synergistic responses to LPS exposure.

Figures

F](#F1)) are positioned closer to TAD boundaries compared to BALB/c-specific genes (n=1205) ([Fig. 5C](#F5), left). Conversely, non-IL-4-induced genes exhibiting strain-specific differences only at basal levels do not show significant distinctions in the distance to TAD boundaries between BL/6 and BALB/c naïve macrophages ([Fig. 5C](#F5), right). Earlier studies have demonstrated that the mean expression of genes within TAD boundaries is higher than that of genes located outside these boundaries. Additionally, genes within boundary regions exhibit lower variability in expression compared to th

D](#F2)). To identify strain-specific IL-4-induced enhancers, we combined H3K27ac ChIP-seq and ATAC-seq analysis. Initially, we overlapped the IL-4-induced ATAC peaks with H3K27ac ChIP-seq data in both strains, allowing us to identify strain-specific induced enhancers ([Fig. 2F](#F2)). We identified 2,718 enhancers induced in IL-4 stimulated BL/6 TRMs and 1,457 enhancers induced in IL-4 stimulated BALB/c TRMs. Notably, the majority of these IL-4 induced enhancers are strain specific, with only 561 shared between strains ([Fig. 2G](#F2)). PU.1 and IRF motifs (PU.1:IRF8) were more enriched in BL

C](#F3), C2 and C4 gene clusters), in comparison to BALB/c TRMs. We identified the promoter-distal regions of C2 and C4 genes through DNA interactions obtained from Hi-C analysis, and MEME analysis was employed to pinpoint promoter-distal regions containing the NF-κB motif. Our findings revealed that 310 (out of 1361) genes from C2 (~23%) and 209 genes (out of 1030) from C4 (~20%) are regulated by promoter-distal regions associated with the NF-κB motif. Motif enrichment analysis conducted with HOMER in these regions demonstrated a significant enrichment of the NF-κB motif, verifying the region

B](#F4), [Fig. 1I](#F1)). In the accessible promoter-distal regions for these genes, we only found NF-κB and IRF motifs in BL/6 TRMs but not in BALB/c TRMs, which may partially explain why these two genes are only upregulated in BL/6 TRMs. We compared the expression levels of genes (n=181) containing an NF-κB motif in either the promoter or promoter-distal regions among the 513 genes specifically upregulated by IL-4 exclusively in BL/6 TRMs ([Fig. 1G](#F1)) with genes lacking NF-κB motifs (n=332). Notably, the expression level for genes with NF-κB motifs is higher than those that are not assoc

A](#F5)–[B](#F5)). In line with previous investigations, the majority (n=1941) of TAD regions are shared between BL/6 and BALB/c macrophages. However, notable differences exist, such as at the Chil3 locus, where gene expression is higher in BL/6, both at the basal level and after IL-4 treatment ([Fig. 5B](#F5), Fig. S3D).

F](#F6)) and subsequently separated for analysis using Signac to identify strain-specific transcription factor (TF) motifs (Fig. S6A).We found that bZIP motifs (e.g. FOS, JUNB) are more enriched in BALB/c macrophages (Fig. S6A, right), consistent with ATAC-seq analysis of bulk cells, whereas KLF motifs are enriched in BL/6 macrophages (Fig. S6A, left). Integrating analyses of chromatin accessibility and gene expression of individual cells can be used to infer gene regulatory networks (GRNs) through SCENIC+^[35](#R35)^. We employed SCENIC+ to construct GRNs for both BL/6 and BALB/c macrophages,

F](#F7), left), those exhibiting a synergistic response ([Fig.7F](#F7), middle), and those displaying a more pronounced response in the LPS treatment group in BALB/c macrophages compared to BL/6 macrophage ([Fig.7F](#F7), right) are shown. Genes identified to be synergistically upregulated by IL-4+LPS, such as *Ccl2*, *Ccl12*, *Ccl7*, and *Edn1*, had comparable percentage of cells expressing these genes between BL/6 and BALB/c macrophages. However, the induced level of expression per cell is higher in BL/6 macrophages compared with Ctrl ([Fig. 7F](#F7), middle), consistent with a cell intrinsi

Full text

Introduction

Tissue-resident macrophages (TRMs) play crucial roles in maintaining tissue homeostasis and responding to inflammation^1, 2, 3^. Like many TRMs, large peritoneal cavity TRMs are predominantly derived from embryonic progenitors seeded into tissues before birth and maintained by local proliferation^4^. Dysfunction of TRMs can cause severe and fatal developmental disorders, as well as dysregulated inflammatory responses^5, 6^. The effects of genetic variation on signal-dependent gene expression in bone marrow derived macrophages (BMDMs)^7^ has revealed inbred mouse strain-specific responses to stimulus, but strain-specific differences in TRMs remains poorly characterized. Understanding the diversity of TRMs response to activation stimuli across different genetic backgrounds may offer new insights into determinants of immune variation.

During type 2 immune response, the cytokines IL-4 and IL-13 signal through IL-4Rα to trigger a macrophage activation state (M(IL-4) or alternatively activated M2 macrophages)^8^ that promotes control of helminth infection^9^ and tissue repair^10, 11^. Importantly, different tissue macrophages will respond differently to IL-4 activation in vivo and we previously found that tissue-resident macrophages and inflammatory macrophages derived from monocytes are phenotypically different following IL-4 activation^3, 12^. IL-4 drives the synthesis and activation of signal-dependent transcription factors (SDTFs)^13^, which are mainly IRF4 and STAT6, while the SDTFs induced after LPS stimuli include IRF1, IRF5,IRF8,STAT1,STAT2, and NF-κB^14^. However, the exact SDTF repertoire may vary depending on differences in the initial cellular state prior to stimulation. The epigenomic and transcriptional programs of TRM responses to inflammatory signals are tightly and dynamically regulated by a combination of transcription factors including lineage-determining transcription factors (LDTFs) and SDTFs. Sequence variants can disrupt the collaborative binding of LDTFs such as PU.1 with SDTFs such as STAT6^15^. Motif mutation analysis of IL-4 activated BMDMs highlights the critical role of LDTFs, including PU.1, C/EBPβ, and AP-1 family members, as well as SDTFs such as STAT6 and PPARγ, in IL-4-driven transcription^16^. Importantly, the IL-4 activation state of macrophages is not a fixed differentiation state and can be converted to a classical activation state in vitro by subsequent treatment with IFNγ/LPS^17^. Recently, IL-4 primed BMDMs have been shown to undergo “extended synergy” with LPS treatment^18^. However, how genetic variation influences the epigenomic changes resulting from IL-4 activation, and subsequent exposure to toll-like receptor (TLR) ligands such as LPS remains unexplored. Furthermore, such systematic examinations of transcriptional regulation have thus far focused on in vitro derived BMDMs, which may be distinct from in vivo derived and activated TRMs.

Here, we uncovered how IL-4 activation can remodel the epigenomic organization differently in TRMs from BL/6 and BALB/c mice. IL-4 activation in BL/6 TRMs results in enhanced accessibility for NF-κB and IRF activity, as defined by ATAC-seq and H3K27ac ChIP-seq (promoter-proximal) or Hi-C (promoter-distal) analysis. This BL/6 specific NF-κB enrichment resulted in an enhanced synergistic effect with subsequent LPS exposure on IL-4-activated BL/6 TRMs compared with BALB/c TRMs. Hence, the background genetic variation can affect the integration of different stimuli within TRMs due to different epigenomic reorganization.

RESULTS

To better understand how genome structure and enhancer-promoter interactions affect transcriptional responses to IL-4 activation in vivo in tissue resident macrophages, we integrated multi-modal genomic data on F4/80^Hi^ large peritoneal macrophages from BL/6 and BALB/c mice after treatment with IL-4-Fc intraperitoneally (i.p.). Analyses included RNA-seq, ATAC-seq, ChIP-seq for H3K27ac, genome-wide chromatin conformation capture (Hi-C), scRNA-seq and scATAC-seq (Fig. 1A). All data generated are available in the Gene Expression Omnibus (GEO) series GSE248038.

Transcriptome analysis reveals strain-specific differences in IL-4-induced gene expression.

To evaluate strain specific transcriptional responses to IL-4 stimulation in vivo, RNA-seq was performed on F4/80+ purified peritoneal TRMs from individual mice (n = 3/group) and principal components analysis (PCA) showed sample replicates clustering tightly together (Fig. 1B). Strain specific differences on the PC1 axis explains more variance (61%) than the effects of IL-4 stimulation on the PC2 axis (34%). Notably, the difference between control (Ctrl) and IL-4 treated TRMs along the PC2 axis is greater for BL/6 than the BALB/c strain (Fig. 1B). The top 10 genes that are driving PC1 for BL/6 TRMs (Fig. 1C) include Gvin1, Marco and Fabp5. The top 10 genes that are driving PC2 (Fig. 1C) include genes associated with alternative activation (e.g. Chil3, Retnla and Arg1), but also include other genes (e.g. Rnase2a, Col1a2 and Col3a1). We next identified genes upregulated by IL-4 treatment in each strain (FDR≤0.05; ≥2-fold) relative to TRMs from naive untreated mice (Ctrl) and found a total of 530 genes upregulated in BALB/c and 1,227 genes upregulated in BL/6. Of those genes 304 overlapped between strains (Fig. 1D), resulting in 226 BALB/c and 923 BL/6 genes uniquely responding to IL-4 in the different strains. This result indicates that BL/6 peritoneal TRMs are more transcriptionally responsive to IL-4 activation than BALB/c TRMs. Since BL/6 and BALB/c TRMs responses could be shaped by differences in the tissue environment, we also examined transcriptional profiles from BMDMs activated in vitro by IL-4^16^ using the same bioinformatics analysis parameters (Fig. 1E). We identified 171 overlapping genes upregulated by IL-4 treatment in both strains. Additionally, BL/6 macrophages exhibited an additional 280 strain-specific upregulated genes, while BALB/c macrophages had 108 strain-specific upregulated genes (Fig. 1E). These results indicate that at least some component of BL/6 macrophages being more responsive to IL-4 stimulation than BALB/c should be cell intrinsic and genetically determined.

BL/6 and BALB/c TRM transcripts can be expressed at different levels under baseline conditions, hence defining only upregulated genes in each strain may miss aspects of strain specific IL-4 mediated transcriptional regulation. We used a multifactor DESeq2 model (further details in the method) to identify strain specific differentially regulated genes (FDR≤0.05) after IL-4 activation and k-means clustering to divide 3 main patterns of transcriptional changes (C1, C2 and C3), of which we focused on IL-4 upregulated genes (C2+C3) (Fig. 1F). 2006 BL/6 specific genes are significantly upregulated, of which 513 have more than 2-fold higher expression compared to naive macrophages (Ctrl). For BALB/c, there are 1205 strain specific upregulated genes, of which only 237 genes have more than 2-fold higher expression compared to naive macrophages (Ctrl). Fig. 1G displays a heatmap illustrating the fold-change values compared with Ctrl for genes list with more than a 2-fold change from Fig. 1F (C2+C3). Additionally, Fig. S1A (left) and Fig. S1A (right) show heatmaps of normalized reads for these genes specific to BL/6 and BALB/c, respectively. The results indicate that BL/6-specific genes are expressed at higher levels than BALB/c-specific genes. Gene ontology (GO) analysis of BL/6 or BALB/c specific genes indicates some genotype specific functional differences in chemotaxis and cell migration (Fig S1B). Examples of IL-4 induced genes that are not strain specific include Chil3 (Fig. 1H), whereas Trem2, which is important in Alzheimer’s disease and cancer^19, 20^, is an IL-4 induced gene specific for BL/6 macrophages (Fig. 1I). These results illustrate the divergent IL-4 transcriptional responses between BL/6 and BALB/c macrophages, with BL/6 macrophages being more transcriptionally responsive to IL-4 activation than BALB/c macrophages.

Next, we performed motif enrichment analysis on the promoter regions (−400bp to +100bp of transcription start site (TSS)) of BL/6 and BALB/c strain specific genes (Fig. 1G). De novo motif analysis identified the RELB motif (Fig. 1J) in BL/6-specific IL-4 activated genes which were induced more than 2-fold compared to Ctrl (Fig. 1J, left). No motifs were enriched in promoters of BALB/c-specific IL-4 activated genes (Fig. 1J, right). RELB constitutes one subunit of NF-κB family^21^. This implies that the genes upregulated specifically in BL/6 in response to IL-4 are more regulated by NF-κB signaling.

Helminth infections induce IL-4 production and Heligomosomoides polygyrus (H. polygyrus) infection can also induce the accumulation of IL-4 activated TRMs in the peritoneal cavity^22^. Hence, we compared the natural accumulation of peritoneal TRMs from H. polygyrus infected BL/6 and BALB/c mice. By RNA-seq, we found that the expanded BL/6 TRMs from H. polygyrus infected mice also have a stronger transcriptional response compared to infected BALB/c mice. 469 genes were uniquely upregulated (FDR≤0.05; ≥2-fold) in BL/6 H. polygyrus TRMs relative to TRMs from untreated mice, compared to 269 in infected BALB/c mice (Fig. 1K). By k-means clustering on the differentially regulated genes after DESeq2 analysis, we found 963 (FDR≤0.05) strain specific genes upregulated (C2+C3) in BL/6 H. polygyrus induced TRMs and only 545 (FDR≤0.05) strain specific genes upregulated (C2+C3) in BALB/c H. polygyrus induced TRMs (Fig.1L). We analyzed the promotor regions of the whole genes list that are strains specific upregulated compared to naïve TRMs (FDR≤0.05, ≥2-fold) for motif enrichment. IRF motifs are enriched in genes specifically upregulated in BL/6 H. polygyrus infected TRMs, but there were no significant motif enrichments in genes specially upregulated in BALB/c H. polygyrus TRMs (Fig.1M).

In summary, these results indicate that peritoneal TRMs from BL/6 mice exhibit a heightened type 2 transcriptional response compared to BALB/c, regardless of whether driven by IL-4 treatment or H. polygyrus infection. This heightened response is likely cell intrinsic, as this observation was confirmed in BMDMs activated in vitro. Additionally, motif analysis of the promoter region of induced genes indicated that NF-κB and IRF transcription factors may play a more substantial role as SDTFs in BL/6 TRMs following IL-4 activation in vivo than in BALB/c TRMs.

IL-4 activation results in strain-specific epigenetic states with distinct transcription factor accessibility in the enhancer landscape.

To investigate the effect of IL-4 activation on the chromatin structure in TRMs, we next performed ATAC-seq to measure open chromatin and enhancers, and H3K27ac ChIP-seq to measure enhancer activation (Fig. 2). First, we identified constitutively accessible regions which are not altered in accessibility in response to IL-4 treatment of which 20,188 regions were shared between strains and 18,336 were uniquely accessible in BALB/c mice (3,731 in BL/6 mice, respectively) (Fig. 2A). We further identified accessible regions that are inducible by IL-4 activation of which only 1,407 regions were shared between strains and 4,801 regions were uniquely accessible in BL/6 (3,660 in BALB/c mice, respectively) (Fig. 2B). Notably, BALB/c macrophages have more strain specific constitutive regions (n = 18,336 in BALB/c versus n = 3,731 in BL/6), whereas BL/6 macrophages have more IL-4 inducible accessible regions (n = 4,801 in BL/6 versus n = 3,660 in BALB/c) (Fig. 2A and 2B). Additionally, induced accessible regions from BL/6 macrophages are enriched in intron and intergenic regions, whereas constitutive accessible regions are enriched in promoters (Fig. 2C). However, this contrast between induced and accessible regions is less pronounced in BALB/c macrophages (Fig. 2C). Hence, BL/6 macrophages are more responsive to IL-4 driven chromatin remodeling than BALB/c macrophages, whereas BALB/c macrophages have more constitutively accessible regions.

To identify transcription factors associated with IL-4 driven chromatin remodeling, we performed transcription factor motif enrichment analysis with HOMER on the IL-4 induced ATAC-seq peak regions (Fig. 2D). This revealed significant enrichment of LDTFs such as Sfpi1 (PU.1) and AP1(Cebpd, Fols2), as expected in both strains. Notably, IRF and NF-κB SDTF motifs are significantly enriched in BL/6 specific IL-4 induced peaks, but not in BALB/c specific induced regions. In contrast, the homeobox motif (HOXA13) is only enriched in BALB/c induced ATAC peaks (Fig. 2D). These results suggest that chromatin remodeling after IL-4 activation in BL/6 TRMs is driven by IRF and the NF-κB family of transcription factors (TFs), while in the BALB/c background the bZIP and homeobox family of TFs are the driving factors. We next used chromVAR^23^ to identify enriched motifs that are strain specific and IL-4 responsive and grouped by hierarchical clustering (Fig. S2A, left). This analysis suggests that the putative activity level of PU.1, NF-κB, IRFs and STATs are enriched in IL-4 activated BL/6 TRMs, whereas AP-1(bZIP motif) activity is enriched in BALB/c TRMs under both baseline as well as IL-4 activated conditions (Fig. 2E). Other bZIP motif clusters, including C/EBP motifs such as C/EBPA, C/EBPB and C/EBPD were also more enriched in BALB/c TRMs after IL-4 treatment (Fig. S2A, right). These findings indicate that IL-4-induced chromatin remodeling enhances transcription factor accessibility differently in BL/6 and BALB/c TRMs and for BL/6 TRMs these include LDTFs like PU.1, along with SDTFs such as NF-κB, IRFs, and STATs.

Next, we applied H3K27ac ChIP-seq to study change of activated enhancers following IL-4 stimulation. We observed 141 enhancers uniquely activated in BL/6 and 347 enhancers uniquely activated in BALB/c (Fig.S2B), which showed different functional enrichment by using the Genomic Regions Enrichment of Annotations Tool (GREAT) on BL/6 and BALB/c specific enhancers (Fig.S2C). De novo motif analysis on strain-specific enhancer regions revealed that IRFs were exclusively enriched in BL/6-specific enhancers, whereas the homeobox motif (Meis) showed enrichment in BALB/c (Fig. S2D), which is consistent with other analysis (Fig. 2D). To identify strain-specific IL-4-induced enhancers, we combined H3K27ac ChIP-seq and ATAC-seq analysis. Initially, we overlapped the IL-4-induced ATAC peaks with H3K27ac ChIP-seq data in both strains, allowing us to identify strain-specific induced enhancers (Fig. 2F). We identified 2,718 enhancers induced in IL-4 stimulated BL/6 TRMs and 1,457 enhancers induced in IL-4 stimulated BALB/c TRMs. Notably, the majority of these IL-4 induced enhancers are strain specific, with only 561 shared between strains (Fig. 2G). PU.1 and IRF motifs (PU.1:IRF8) were more enriched in BL/6, whereas bZIP motifs (cEBP-like) were more enriched in BALB/c (Fig. 2G). We associated enhancers with genes specifically upregulated in BL/6 and BALB/c TRMs (Fig. 1G) by applying the enhancers to the nearest promoters of genes using the HOMER and found that the enhancers near the BL/6 specific upregulated genes (promoter-proximal genomic regions) are more enriched for SDTFs motifs (e.g. IRFs) (Fig. 2H). When we mapped BL/6 and BALB/c specific IL-4 upregulated genes to the nearby strain specific induced enhancers (Fig. 2I), we found that out of 513 BL/6 specific genes, 135 are associated with induced enhancers, whereas only 16 out of 237 BALB/c specific genes are associated with induced enhancers. Additionally, when we defined super-enhancer regions by H3K27ac ChIP-seq using HOMER^24^, after annotating the super-enhancer with the nearest promoter of gene, we found that 184 out of 2,006 BL/6 specific IL-4 upregulated genes are associated with super enhancers nearby whereas only 38 BALB/c specific genes are associated with super enhancers. The number of total super-enhancer regions in BL/6 (n = 578) and BALB/c (n = 555) macrophages are quite similar. These results may partially explain the increased IL-4 transcriptional responsiveness of BL/6 TRMs, as IL-4 upregulated genes in BL/6 TRMs are more associated with induced enhancers and super enhancer regions.

IL-4 primes BL/6 TRMs to synergize with lipopolysaccharide (LPS) signaling more effectively than BALB/c TRMs.

Recently, IL-4 priming was shown to increase NF-κB-p65 binding resulting in “extended synergy” of IL-4-polarized BMDMs upon lipopolysaccharide (LPS) exposure^18^. Since, we observe that in vivo IL-4 activation in BL/6 TRMs enhances accessibility to motifs associated with SDTFs such as NF-κB and IRFs, we hypothesized that synergy between IL-4 priming and LPS exposure would also be more pronounced in BL/6 TRMs relative to BALB/c TRMs. Hence, we compared the in vitro (3hrs) LPS response of TRMs isolated from Ctrl and IL-4 in vivo treated BL/6 and BALB/c mice by RNA-seq (Fig. S3A). PCA shows that the transcriptional response to IL-4+LPS exposure is more distinct from LPS treatment in BL/6 TRMs compared to BALB/c TRMs (Fig. 3A). However, the distance on the PC1 axis between LPS exposed and untreated TRMs is greater on the BALB/c background (Fig. 3A), indicating that the transcriptional effects of LPS on BALB/c TRMs may be greater. We next identified differentially expressed genes (DEGs) (FDR≤0.05 and ≥2-fold ) between untreated vs LPS treated TRMs (LPS vs Ctrl), as well as between LPS treated TRMs and LPS treated IL-4 activated TRMs (IL-4+LPS vs LPS) on both backgrounds (Fig. 3B). The majority of DEGs (n=863) between Ctrl versus LPS exposure in BL/6 TRMs are also affected in in BALB/c TRMs. Consistent with the PCA plot, more genes are differently upregulated between Ctrl and LPS in BALB/c macrophages (n=831) than BL/6 (n=178) indicating a stronger transcriptional response to LPS exposure compared to naive peritoneal BALB/c TRMs (Fig. 3B, left). However, in line with our observation in the PCA plot (Fig. 3A), when we compared the response of LPS treated IL-4 activated TRMs with LPS treated naïve TRMs (Fig. 3B, right), we found that more genes are upregulated in BL/6 (n=909) compared to BALB/c TRMs (n=113) with an only small number of genes affected in both strains (n=161). GO analysis revealed that ‘metabolite process’ was associated with the BALB/c-specific upregulated genes (LPS vs Ctrl, n=831), while ‘immune system process’ was linked to the BL/6-specific upregulated genes (IL-4+LPS vs LPS, n=909) (Fig. S3B).

Next, we employed k-means clustering to group the differentially regulated genes between IL-4+LPS and LPS alone, as identified by DESeq2 using fold-change values across all samples relative to LPS treatment (Fig. 3C). Additionally, an alternative k-means analysis based on fold-change values relative to Ctrl is presented in Fig. S3C. Both analyses revealed the presence of six distinct clusters. Notably, cluster 2 (C2:1361) and cluster 4 (C4:1030) exhibited a marked increase in LPS inducibility in IL-4-primed BL/6 TRMs compared to LPS treatment only (Fig. 3C), and cluster 2 exhibited a synergistic response following IL-4+LPS treatment (Fig. 3D). From the 1361 genes in Cluster 2 (C2), 568 genes showed a greater than 2-fold change in BL/6, while in Cluster 4 (C4), 266 genes met the filtering criteria. In contrast, among the 754 genes in Cluster 6 (C6), only 37 exhibiting a 2-fold difference between IL-4+LPS compared to LPS treatment only in BALB/c macrophages. This demonstrated that extended synergy of IL-4-polarized TRMs upon LPS exposure is much more pronounced in the BL/6 genetic background. Six representative genes of cluster C2 are shown in Fig. 3E. Of these genes, Edn1 and Ccl2 have previously been reported upregulated in IL-4-primed and LPS-exposed BL/6 BMDMs^18^. Induction of Nfkbia by LPS (Fig. S3D, left) and Chil3 by IL-4 (Fig. S3D, right) in both BL/6 and BALB/c macrophages shows that IL-4 and LPS are functionally active in both strains. GO analysis on the C2 and C4 gene clusters showed that the majority of immune response genes were concentrated in the C2 cluster, whereas C4 genes were predominantly associated with metabolic processes (Fig. S3E).

To elucidate the epigenomic state of IL-4-primed TRMs at the regulatory elements governing genes in C2 and C4, we examined the ATAC-seq signal at these loci. For both C2 and C4 genes, the promoter-proximal genomic regions, as defined by adjacent ATAC-seq peaks, exhibited enhanced ATAC-seq signals in IL-4-primed BL/6 TRMs. Notably, this heightened accessibility was only observed in BL/6 IL-4-treated TRMs, in contrast with BALB/c TRMs where such induced accessibility was not evident. This highlights a greater accessibility in BL/6 IL-4-treated TRMs compared to BALB/c TRMs in the promoter-proximal genomic regions of the synergistically regulated genes (C2) (Fig 3F).The accessible regulatory elements at the promotor-proximal regions of C2 and C4 genes are enriched for the Sfpi1 (PU.1) and IRF motifs (Fig 3G), and C2 genes have a significant enrichment of the NF-κB motif. These results indicate that IL-4 treatment in BL/6 modified the epigenome to increase accessibility for the IRF and NF-κB family of TFs, which may mediate the extended synergy of IL-4 primed LPS transcriptional responses. However, this epigenomic conditioning for synergistic responses does not occur in the BALB/c background.

Hi-C analysis identifies NF-κB binding sites in the promoter-distal regions associated with BL/6 specific IL-4 induced genes.

Regulatory elements can govern gene expression across extensive genomic distances, exerting influence on genes situated far from promoters—termed here as promoter-distal regions^25^. To define promoter-distal regions and investigate 3D genome structure, we performed Hi-C on IL-4 activated BL/6 and BALB/c peritoneal TRMs. We examined the role of long-range interactions of regulatory elements in controlling strain-specific IL-4 transcriptional responses through 3D interactions. We first combined Hi-C analysis with ATAC-seq data to define the promoter based regulatory elements. We determined the promotor-distal regulatory elements for BL/6 specific and BALB/c specific IL-4 upregulated genes, and next identified TF motifs enriched in these accessible regulatory elements. We found 535 promoter-distal regions were related to the BL/6 specific upregulated genes (513) and 143 promoter-distal regions were related to BALB/c specific upregulated genes (237). Motif analysis by HOMER revealed enrichment of IRF motifs (de novo) and NF-κB motifs (known) in promoter-distal regions specific to BL/6 (Fig.4A, left), but not in BALB/c promoter-distal regions (Fig.4A, right). These findings indicate that NF-κB and IRFs may play a role in the regulation of IL-4 upregulated genes specific to BL/6 through long-range promoter-distal interactions.

Utilizing MEME motif analysis^26^, we found that 18 out of 237 BALB/c specific IL-4 upregulated genes have NF-κB motifs in the distal regulatory regions (~7%), whereas 81 out of 513 BL/6 specific IL-4 upregulated genes have NF-κB motifs on the proximal-distal regulatory regions (~15%). For example, Trem2 and Wdfy3 are only upregulated in BL/6 TRMs after IL-4 activation (Fig. 4B, Fig. 1I). In the accessible promoter-distal regions for these genes, we only found NF-κB and IRF motifs in BL/6 TRMs but not in BALB/c TRMs, which may partially explain why these two genes are only upregulated in BL/6 TRMs. We compared the expression levels of genes (n=181) containing an NF-κB motif in either the promoter or promoter-distal regions among the 513 genes specifically upregulated by IL-4 exclusively in BL/6 TRMs (Fig. 1G) with genes lacking NF-κB motifs (n=332). Notably, the expression level for genes with NF-κB motifs is higher than those that are not associated with NF-κB motifs (Fig. 4C). Also, the level of induction by IL-4 compared to Ctrl is higher for genes associated with NF-κB motifs in BL/6 TRMs (Fig. 4C). Hence, NF-κB may play a role in BL/6 TRMs specific promoter-distal transcriptional regulation during IL-4 activation in vivo.

Subsequently, we investigated the interplay between 3D genome organization and genes that undergo upregulation upon IL-4 priming followed by LPS exposure in BL/6 TRMs (Fig. 3C, C2 and C4 gene clusters), in comparison to BALB/c TRMs. We identified the promoter-distal regions of C2 and C4 genes through DNA interactions obtained from Hi-C analysis, and MEME analysis was employed to pinpoint promoter-distal regions containing the NF-κB motif. Our findings revealed that 310 (out of 1361) genes from C2 (~23%) and 209 genes (out of 1030) from C4 (~20%) are regulated by promoter-distal regions associated with the NF-κB motif. Motif enrichment analysis conducted with HOMER in these regions demonstrated a significant enrichment of the NF-κB motif, verifying the region selection made by MEME (Fig. S4A). The ATAC-seq data revealed that the accessible regulatory elements within the promoter-distal regions of genes belonging to both C2 and C4 from Fig. 3C are also the regions that experience increased accessibility following IL-4 activation in BL/6 TRMs (Fig. S4B). For instance, the promotor distal region of Hbegf, which is important in would healing^27^ from cluster C2, has a NF-κB motif under BL/6 specific loop that may partly explain the synergistic LPS response (Fig. S3C, D). Hence, the heightened transcriptional activity upon LPS activation appears to be associated with IL-4 priming-induced chromatin remodeling in BL/6 TRMs. These findings imply the presence of a preexisting and comparatively stable landscape of BL/6-specific enhancer-promoter interactions in macrophages, whose regulatory function becomes activated in response to IL-4.

Genes activated by IL-4 in BL/6 TRMs exhibit closer proximity to topologically associating domain (TAD) boundaries and display a more homogeneous expression pattern.

Hi-C data also provides information on TADs. TADs are genomic regions that form units of three-dimensional (3D) nuclear organization whereby DNA physically interacts with each other, which can regulate gene expression through enhancer-promoter interactions within TADs^28^. We used the Juicer pipeline and TAD caller arrowhead^29^ to annotate TADs in BL/6 (n = 2362) and BALB/c (n = 2227) macrophages (Fig. 5A–B). In line with previous investigations, the majority (n=1941) of TAD regions are shared between BL/6 and BALB/c macrophages. However, notable differences exist, such as at the Chil3 locus, where gene expression is higher in BL/6, both at the basal level and after IL-4 treatment (Fig. 5B, Fig. S3D).

TAD boundaries (see methods) delineate the regions that lie between TADs^30^. Upon mapping IL-4-induced genes to TAD boundary regions, we observed that the promoters of BL/6-specific upregulated genes (n=2006, Fig. 1F) are positioned closer to TAD boundaries compared to BALB/c-specific genes (n=1205) (Fig. 5C, left). Conversely, non-IL-4-induced genes exhibiting strain-specific differences only at basal levels do not show significant distinctions in the distance to TAD boundaries between BL/6 and BALB/c naïve macrophages (Fig. 5C, right). Earlier studies have demonstrated that the mean expression of genes within TAD boundaries is higher than that of genes located outside these boundaries. Additionally, genes within boundary regions exhibit lower variability in expression compared to those situated in non-boundary regions^31^. Here, we also observed that BL/6-specific IL-4 upregulated genes (n=2006) exhibit higher expression levels, as indicated by normalized read counts, compared to BALB/c-specific IL-4 upregulated genes (n=1205) (Fig. 5D). Furthermore, to analyze transcriptional variability between strains, we calculated the Fano factor (a measure of gene expression noise^32^), and observed a lower Fano factor for BL/6-specific genes upregulated after IL-4 activation compared to BALB/c-specific upregulated genes (Fig. 5E). This aligns with the notion that gene expression variability tends to be lower in genes situated near TAD boundary regions. Using Weighted Gene Co-expression Network Analysis (WGCNA)^33^to identify modules of co-expressed genes, we found that BALB/c macrophages exhibited a higher number of modules (Fig. 5F), contrasting with BL/6 macrophages where a singular large module prevailed. Hence, there is more variability in the expression of genes in the BALB/c macrophages. When we analyzed strain-specific IL-4 induced genes (Fig. 1G, BL/6: n=513 and BALB/c: n=237), we observed that the promoter regions of BL/6-specific IL-4 upregulated genes (n=513) were closer to TAD boundaries compared to those in BALB/c mice (Fig. S5A), and, after IL-4 activation, the BL/6-specific genes exhibited higher expression values and lower variability, as indicated by the Fano factor (Fig. S5B, C).

In summary, our findings suggest that IL-4-induced genes specific to BL/6 TRMs tend to be located closer to TAD boundaries, exhibiting higher expression levels and lower variability. Conversely, BALB/c-specific IL-4-activated genes are situated further away from TAD boundaries with lower expression levels and higher variability. These differences may contribute to the stronger transcriptional response to IL-4 activation in BL/6 TRMs.

scRNA-seq analysis of BL/6 and BALB/c macrophages in the same tissue environment of chimeric F1 mice identifies cell intrinsic strain specific regulons.

The epigenetic state and transcriptional responses of BL/6 and BALB/c TRMs to IL-4 stimulation in vivo could be influenced by strain specific homeostatic differences in the tissue environment prior to IL-4 exposure. To investigate whether cell intrinsic genetic differences between BL/6 and BALB/c TRMs mediate IL-4 responses when they are generated within the identical tissue environment, we generated mixed bone marrow chimeric mice using CB6F1/J (F1) mice as recipients, reconstituted with mixed bone marrow (BM) from BALB/c and BL/6 mice (Fig. 6A). These mixed BM chimeric mice were then treated with IL-4 in vivo, and the total peritoneal exudate cells (PECs) were subjected to scRNA-seq. Based on the genetic variation between BL/6 and BALB/c transcripts, we used Souporcell^34^ to separate BL/6 and BALB/c PECs from chimeric F1 mice (Fig. 6B, upper). We subsequently focused on the macrophage subset (Fig. 6B, lower). As anticipated, Chil3 and Arg1 are expressed in both BL/6 and BALB/c TRMs, with expression observed in the majority of cells (Fig.6C, left). In contrast, Fabp5 demonstrates higher expression in BL/6 TRMs at both the population and average per-cell expression levels. Conversely, Spp1 exhibits higher population-level expression in BALB/c macrophages (Fig.6C). We then validated that the expression of IL-4-induced BL/6-specific genes, as indicated previously by bulk RNA-seq data, is also notably higher in BL/6 TRMs within the F1 mice. Notably, these genes include Trem2, Sgmas1, Per3, Zmiz1, and Mertk, whereby NF-κB motifs were identified in the promoter-distal regions through MEME-FIMO analysis, defined by chromatin loops from Hi-C data (Fig. 6D). Subsequently, we employed single-cell regulatory network inference and clustering (SCENIC) analysis to identify transcription factor (TF) regulons in BL/6 and BALB/c macrophages. This analysis revealed notable activity of NF-κB (Rela and Relb) and STAT (Stat1) regulons in a distinct subset of BL/6 macrophages, which is not detected in BALB/c macrophages (Fig. 6E). Notably, only a distinct subset of BL/6 TRMs is regulated by the NF-κB. Hence, differences between strains could potentially be driven by just a small subset of the macrophage population.

Next, we conducted single-cell ATAC-seq (scATAC-seq) on PECs from BL/6 and BALB/c mice that have been treated with IL-4 in vivo to examine the accessible chromatin landscape of peritoneal TRMs in BALB/c and BL/6 mice at single cell level. The macrophage populations were identified by integrating the scATAC-seq data with the scRNA-seq data from the mixed BM chimeric mice (Fig. 6F) and subsequently separated for analysis using Signac to identify strain-specific transcription factor (TF) motifs (Fig. S6A).We found that bZIP motifs (e.g. FOS, JUNB) are more enriched in BALB/c macrophages (Fig. S6A, right), consistent with ATAC-seq analysis of bulk cells, whereas KLF motifs are enriched in BL/6 macrophages (Fig. S6A, left). Integrating analyses of chromatin accessibility and gene expression of individual cells can be used to infer gene regulatory networks (GRNs) through SCENIC+^35^. We employed SCENIC+ to construct GRNs for both BL/6 and BALB/c macrophages, after integrating the scATAC-seq data with scRNA-seq data and subsetting the macrophages from the total PEC population. The enhancer-driven Gene Regulatory Networks (eGRNs) constructed for IL-4-activated BL/6 and BALB/c TRMs reveal distinct transcription factor networks. SDTF motifs such as STATs and ETS motifs such as Spib (PU.1) and E4f1 exhibit higher activity in BL/6 TRMs, whereas homeobox TFs such as Meis1 are more active in BALB/c TRMs. These findings align with the data obtained from bulk ATAC-seq analysis (Fig. 6G, H).

In summary, these findings validate that the major genomic differences observed between bulk TRMs of BL/6 and BALB/c mice can be largely reproduced when analyzing macrophages isolated from the exact same tissue environment at a single-cell level. Therefore, cell-intrinsic genetic variation significantly influences the strain-specific responses of peritoneal TRMs to IL-4 activation.

Cell-intrinsic strain-specific responses to IL-4 and LPS synergy are observed in the same tissue environment of chimeric F1 mice.

To assess whether the extended synergy between IL-4 priming and LPS activation results in distinct functional consequences in BL/6 and BALB/c mice, we subjected IL-4-treated BL/6 and BALB/c mice to an in vivo challenge with a sub-lethal dose of LPS. High dose LPS treatment in vitro overwhelmed the transcriptional signature of IL-4 primed synergy in both strains (data not shown). We used an adapted murine sepsis score (MSS) to evaluate the effects of sub-lethal dose LPS challenge (Fig. 7A). Consistent with a stronger transcriptional response to LPS in BALB/c peritoneal TRMs, we found that BALB/c mice have a higher MSS score than BL/6 mice after LPS treatment (p = 0.0108). However, IL-4 pre-treatment does not affect the MSS score in BALB/c mice after LPS challenge, but for BL/6 mice IL-4 pre-treatment trended towards having a higher MSS score (p=0.0867) (Fig. 7A).

To assess the cell intrinsic and strain specific in vivo transcriptional response to IL-4 and LPS synergy in the same tissue environment, we performed scRNA-seq analysis of PECs from mixed BM-chimeric mice (Fig.7B, C). We compared naïve untreated mice (Ctrl), mice treated with LPS alone, and those pre-treated with IL-4 and subsequently challenged with LPS. The separation of BL/6 and BALB/c cells from the mixed BM-chimeric F1 mice was achieved through Souporcell^34^ and cell types were annotated by singleR^36^(Fig 7C). As expected, LPS treatment resulted in reduction of the macrophage population (Fig. 7D, Fig. S6B). In Fig. 7D, the reduction is visually evident through the diminished green color, which represents the macrophage population. This phenomenon is commonly referred to as the “macrophage disappearance reaction”^2^. LPS treatment also increased the neutrophil population in PECs (Fig. 7D, Fig. S6B). IL-4 pre-treatment slightly reduced the loss of macrophages. However, there retained sufficient number of macrophages in all groups to allow for meaningful downstream analysis. To identify differential peaks among groups, our focus was exclusively on highly induced genes (logFC threshold = 0.25). Consistent with bulk cell analysis, we observed a more pronounced upregulation of genes by LPS in BALB/c macrophages compared to control naïve (Ctrl) macrophages, as opposed to BL/6 macrophages (Fig. 7E). However, in BL/6 macrophages, we identified a higher number of genes upregulated in the IL-4+LPS treatment group compared to the LPS-only group, indicative of a more robust synergistic response on the BL/6 background. Although the number of significant upregulated genes identified in both comparisons through scRNA-seq is considerably lower than in the bulk RNA-seq experiments (Fig. 3B), the overall pattern remains consistent.

Representative genes with no synergistic response (Fig.7F, left), those exhibiting a synergistic response (Fig.7F, middle), and those displaying a more pronounced response in the LPS treatment group in BALB/c macrophages compared to BL/6 macrophage (Fig.7F, right) are shown. Genes identified to be synergistically upregulated by IL-4+LPS, such as Ccl2, Ccl12, Ccl7, and Edn1, had comparable percentage of cells expressing these genes between BL/6 and BALB/c macrophages. However, the induced level of expression per cell is higher in BL/6 macrophages compared with Ctrl (Fig. 7F, middle), consistent with a cell intrinsic response. Next, we re-clustered the macrophages into distinct phenotypes to identify greater granularity in macrophage heterogeneity (Fig. 7G) and applied SCENIC analysis to identify TF regulons associated with the distinct macrophage clusters (Fig. 7H). Cluster 2 macrophages experiences the most significant depletion after LPS exposure and is distinguished by the activity of Gata6, Klfs, and Fli1 regulons (Fig. 7H, Fig. S6C). Cluster 6 macrophages are more abundant in the IL-4+LPS treatment group compared to the LPS-only group, particularly from the BL/6 background (Fig. S6C). SCENIC analysis of Cluster 6 reveals an enrichment of the NF-κB (Rel) regulon. As expected for naïve peritoneal macrophages (Ctrl), Gata6 regulon activity is predominant in the majority of these cells. However, after LPS treatment (including in the context of IL-4 pretreatment), the activity of SDTF regulons such as Nkb1, Stats, and Irf becomes activated in the remaining TRMs (Fig. 7I). Overall, these data from BALB/c and BL/6 macrophages that are primed by IL-4 treatment prior to exposure to LPS in the same tissue environment of a chimeric F1 mice clearly indicates that the strain specific differences in IL-4 and LPS transcriptional synergy are driven by cell intrinsic genetic variation.

Discussion.

We found that natural genetic variation between BL/6 and BALB/c mice shapes the epigenomic reprogramming of peritoneal TRMs by the type 2 cytokine IL-4, resulting in subsequent differences in extended synergistic responses to LPS activation. Our findings illustrate how the integration of multiple stimuli by macrophages in vivo could be shaped by the natural genetic differences between individuals. Additional complexity was observed by scRNA-seq analysis to determine cell-intrinsic effects of genetic variation in mixed chimeric F1 hybrid mice. These findings illustrate the importance of studying transcriptional regulation in mice of different genetic backgrounds and not just with the BL/6 background.

Our results are strikingly similar to recent studies on macrophages in the liver called Kupffer cells, whereby there was a clear bias for AP-1 and MAF family activation in BALB/c Kupffer cells and IRF, NF-κB and LXR activity in C57BL/6J (BL/6) Kupffer cells^37^. Additionally, recent work suggested that AP-1 plays a lesser role in the LPS synergy observed in IL-4-primed macrophages^18^. Our results here tie together the findings that AP-1 TFs have more activity in BALB/c mice, whereas NF-κB TFs are more important in BL/6 mice, which may explain why IL-4/LPS synergy is more significant in macrophages from BL/6. Additionally, we extend upon the important discovery that IL-4 could remodel the epigenome of BMDMs from BL/6 mice to expose NF-κB motifs that then mediate an extended synergistic response to LPS exposure^18^. We find here that IL-4 could also remodel the epigenome of peritoneal TRMs to expose NF-κB motifs, specifically on the BL/6 background but not on the BALB/c background. This results in a more heightened synergistic response to LPS exposure in BL/6 TRMs compared to BALB/c TRMs. In addition to confirming the preference of the BL/6 background for NF-κB motif enrichment and AP-1 motif enrichment in BALB/c macrophage^37^, we provide a functional consequence that these TFs preference affect how genetic background can remodel the epigenome to integrate subsequent stimuli, such as LPS exposure.

The impact of natural genetic variations on type 2 immunity is important in uncovering the underlying mechanisms contributing to phenotypic diversity across a wide spectrum of inflammatory diseases such as asthma and atopic dermatitis. This line of investigation provides insight into how individuals may exhibit distinct phenotypes when integrating signals from multiple stimuli. Understanding how genetic variation affects priming and shaping complex features of the immune response not only enhances our comprehension of signal integration within cells, but how this may shape disease heterogeneity. This may lay the future groundwork for developing more targeted and personalized therapeutic strategies. Studies to understand transcriptional programs in macrophage activation have primarily been performed with BMDMs stimulated in vitro^16^. There is a lack of studies investigating tissue-resident macrophages. Our study is a first step in filling this gap by investigating peritoneal TRMs that play key roles in the control of infections and inflammatory pathologies^38, 39^, as well as in the maintenance of immune response robustness^40^. We specifically focused on the large peritoneal cavity macrophages (LPM) that are F4/80^hi^MHCII^low^ and depend on the expression of Gata6 for their differentiation^40, 41^. In response to inflammation, LPM will aggregate at specific sites in the peritoneal cavity to contain infection or mediate tissue repair^42^. IL-4 treatment of infection with nematode parasites can expand the LPM population through proliferation^43^, while activating the macrophages to adopt a tissue repair phenotype. Hence, this is a powerful in vivo model to study type 2 macrophage activation in different genetic backgrounds.

In our genomic studies, we focused on strain-specific differences in enrichment patterns of LDTFs and SDTFs. In line with previous studies investigating IL-4 activation in BMDMs^16^, our findings from peritoneal TRMs exhibit similar motif enrichment pattern. Notably, in the realm of LDTFs, we observed that ETS (PU.1 or ELFs) motifs are more prominent in the BL/6 background. Additionally, IRFs and NF-κB emerge as SDTFs crucial in the BL/6 context after IL-4 activation and were important in fostering an enhanced synergistic response within peritoneal TRMs from BL/6 mice. These results may arise from the collaborative efforts of SDTFs with LDTFs, including PU.1, AP-1, and C/EBP (CCAAT/enhancer binding protein)^16^. This intricate network underscores their collective significance in governing the complexities of IL-4-induced transcriptional regulation. Notably, BMDM experiments also identified EGR2 except NF-κB as a key transcription factor influencing the synergistic impact of LPS on IL-4-primed macrophages^18^. While our current studies did not observe enrichment of EGR motifs, some preliminary experiments in thioglycolate-induced monocyte-derived macrophages suggests that the activity of EGR2 is more significant in inflammatory cells (data not shown). Hence, transcriptional patterns of SDTFs may differ substantially between embryonically derived peritoneal TRMs when compared to BMDMs or thioglycolate-induced monocyte-derived macrophages of hematopoietic origin. In contrast to BL/6 TRMs, bZIP motifs, including AP.1, Maf and Jund, are notably more enriched in BALB/c TRMs, as evidenced by both bulk and single cell sequencing analyses.

Integrative analysis of transcriptional changes and chromatin structure indicated that IL-4 activation is accompanied by increased activity of NF-κB and IRFs in regulatory regions associated with genes specifically upregulated in BL/6 mice. These regions predominantly encompass promoter-proximal regions defined by ATAC-seq and H3K27ac ChIP-seq signal. Furthermore, when considering the three-dimensional (3D) chromatin structures, promoter-distal regions defined by loops from Hi-C also indicate increased binding of SDTFs, such as NF-κB or IRFs, in promoter-distal regions associated with genes upregulated in BL/6. For example, genes like Trem2 and Wyd3, which show upregulation in BL/6, exhibit motifs for both NF-κB and IRFs in the promoter-distal regions. Hence, IRFs and NF-κB may exert influence on gene expression through both 1D (proximal to the promoter) and 3D (distal to the promoter) interactions, which may heighten the synergistic response observed after IL-4+LPS treatment on the BL/6 background.

It was unexpected to observe a role for NF-κB transcription factors in IL-4 activation, as these TFs are typically involved in TLR ligand (LPS) activation. NF-κB complexes are formed by members of the Rel family of transcription factors^21^. While LPS signaling is primarily dominated by RelA, it remains unclear which subunit of NF-κB complex dominates in IL-4-primed macrophages, despite some studies focusing on the knockout of NF-κB subunits in mice to understand their role in type 2 cytokine responses^44, 45^. Perhaps these genes are regulated by different subunits of NF-κB as different NF-κB subunits may competitively bind the same DNA sites^46^. Importantly, we found that variation in NF-κB motif sequences between BL/6 and BALB/c mice does not appear to contribute significantly to the differences in accessibility and transcription. But the cooperation between LDTFs and SDTFs is complex. For example, binding of NF-κB can be dependent on LDTFs binding, as mutations in PU.1 or C/EBPβ motifs could abolish signal dependent binding of NF-κB^15^. Hence, the differences in NF-κB activity may not be driven by mutations in the cognate recognition motif but could be driven by genetic variation in the motifs of other LDTFs or SDTFs that collaborate to mediate the different response between BL/6 and BALB/c strains.

In addition to cell-intrinsic genetic factors such as chromatin accessibility, 3D chromatin structure, and/or natural genetic variations^47, 48^, the cytokine milieu of the tissue environment plays a crucial role in shaping the innate immune response. Chromatin modifications can be regulated by crosstalk between the environment and cell ontogeny through just a small number of TFs^49^. How cell intrinsic and tissue environment interactions determines cellular state remains poorly understood. Recent studies indicate that the tissue environment plays a major role in determining macrophage phenotypes by regulating the expression and function of transcription factors that activate cis-regulatory enhancer elements^50, 51, 52^. The genetic background of the individual will thus determine how macrophages that are exposed to similar environmental cues can display different patterns of gene expression^51^. We examined mixed bone marrow chimeric F1 mice to include both BL/6 and BALB/c tissue resident macrophages in the same peritoneal environment. Notably, after IL-4 treatment and analysis by single cell sequencing, we confirmed that STATs, IRFs and NF-κB regulons were enriched specifically in BL/6 peritoneal TRMs. Hence, the enhanced activity of these SDTFs, as well as IL-4 and LPS synergy, in IL-4 treated BL/6 peritoneal TRMs is likely driven by cell intrinsic genetic factors, and not by the tissue microenvironment.

Our work underscores the substantial impact of inherent genetic variations between BL/6 and BALB/c mice on the distinct responses to in vivo IL-4 stimulation in TRMs. Notably, the synergistic interplay between IL-4 stimulation and TLR engagement via LPS is markedly more pronounced in BL/6 TRMs compared to their BALB/c counterparts. Epigenomic modifications in the chromatin structure of BL/6 TRMs, as observed through bulk and single-cell level analyses following IL-4 activation, resulted in heightened accessibility of NF-κB and IRFs motifs for binding, a feature absents in BALB/c TRMs. These strain-specific epigenetic responses to IL-4 results in distinct transcriptional responses to subsequent TLR engagement, which may contribute towards genetic variation exerting in vivo effects in a complex inflammatory milieu. Hence, genetic variation may be a key influencer of how repeated activating signals affect epigenomic memory and synergistic responses. Additionally, we explored the diverse responses observed in TRMs, attributing some of these distinctions to the inherent variability in RNA expression, particularly prominent in BALB/c cells. This variability potentially contributes to the divergent responses observed in the two strains. Moreover, our findings indicate that the peritoneal environment does not play a pivotal role in shaping the divergent responses following IL-4 treatment. A more detailed analysis with a combination of scRNA-seq and scATAC-seq analysis of BL/6 and BALB/c macrophages in the same tissue microenvironment may reveal how interactions of multiple genes and cell heterogeneity can quantitatively contribute to synergistic signaling in macrophages. These insights may particularly elucidate how the genetic background of specific individuals influences inflammatory hyperresponsiveness to TLR (LPS) activation in TRMs during type 2 inflammatory responses. A deeper understanding of these mechanisms holds the potential to enhance our capacity to customize immune-mediated treatments for individuals, offering a personalized approach to therapy for inflammatory diseases driven by type-2 cytokines.

Methods and materials

Mice and IL-4-Fc and LPS treatment

Male BALB/cJ and C57BL/6J mice, typically 6 to 8 weeks old, were purchased from The Jackson Laboratory and/or bred in specific-pathogen-free facilities at the NIH. A fusion protein of mouse IL-4 with the Fc portion of IgG1 (custom order with Absolute Antibody) was generated to extend half-life with similar effects to IL-4–anti–IL-4 mAb complex. Mice were injected i.p. with either PBS or 10 μg IL-4-Fc in 100 μl PBS on days 0 and 2, and peritoneal exudate cells (PECs) were harvested on day 4. PECs were enriched for >90% F4/80+ macrophages cells by using EasySep^™^ Mouse F4/80 Positive Selection Kit (100–0659, Stemcell Technologies). In vivo LPS (500 μg/kg) was administered after 4 days with or without prior IL-4-Fc injection, and the MSS score was observed the following day around 9:00 am. The experiments were performed under protocol number LPD16E. p-values for the MSS score were calculated using a two-sided Mann-Whitney test.

Generation of mixed BL/6 and BALB/c bone marrow chimera F1 mice

Mixed chimeric mice were generated by either irradiation or busulfan treatment. In one set of experiments, F1 mice: CB6F1/J (Jax no. 100007) were lethally irradiated with 950 rad and subsequently reconstituted i.v. with 2 × 10^6^ bone marrow cell mixture 1:1 from BL/6 and BALB/c donor mice after depletion of T lymphocytes by negative selection on CD90.2 MACS columns (Miltenyi Biotec). Alternatively, F1 were treated i.p. with 25 mg/kg of Busulfan on day −6, −4, −2 before receiving on day 0 bone marrow cell mixture. All BM recipients were maintained on antibiotic TMS water containing Trimethoprim (0.13 mg/ml), Sulfadiazine (0.5 mg/ml), and 0.67 mg/ml Sulfamethoxazole (0.67 mg/ml) for 8 weeks. The chimerism of reconstituted mice was confirmed by FACS using APC anti-H-2Kb (clone AF6–885.5.3 Invitrogen) and Super bright 436 anti-H-2Kd (clone SF1–1.1.1 Invitrogen). Data were collected on BD Symphony and analyzed by FlowJo. After 8 weeks, chimeric mice received two intraperitoneal injections of 10 μg IL-4-Fc per mouse 2 days apart. After 4 days of IL-4 treatment in vivo peritoneal cells were harvested for scRNA-seq. In other experiments, just before harvesting the peritoneal cells, LPS was injected at a dose of 1 mg/kg in vivo for 3 hours. The peritoneal cells were then isolated from the chimeric mice and used to prepare the scRNA-seq library.

ChIP and ChIP-seq library construction

ChIP experiments were performed as described previously^53^. Briefly, isolated macrophages were washed twice with PBS, pin samples for 10min, 300g at 4° C, fixed at room temperature with 1% formaldehyde in PBS for 10min. Reactions were quenched by adding glycine to a final concentration of 0.125 M, and cells were washed twice with cold PBS. Cells were resuspended in cell lysis buffer (1% SDS, 10mM EDTA, 50mM Tris-HCl pH8.1), incubated in ice for 1 hour. The crosslinked chromatin was fragmented by sonication (Covaris, M220) to sheer size 200–2000bp. The final sonication product was diluted 5 times by dilution buffer (0.01% SDS, 1.1% Triton X-100, 1.2 mM EDTA, 16.7mM Tris-Hcl pH 8.0, 167mM NaCl). ChIP was performed with 3 μg anti-H3K27ac (39133, Active Motif) antibody-bound Protein A beads (10001D, Invitrogen) and incubated at 4°C overnight. Beads were collected by centrifugation, washed, and incubated at 65°C for 4h in elution buffer (50 mM Tris-HCl, pH 7.5; 10 mM EDTA; 1% SDS) to reverse cross-linking. ChIP DNA was purified by ChIP DNA Clean & Concentrator (Zymo Research: D5205). For sequencing, ChIP-seq samples were pooled and sequenced on NextSeq2000 using TruSeq ChIP Library Prep Kit (IP-202–1012) and single-end sequencing. Two biological replicate ChIP-seq experiments were carried out. ChIP-seq data are available on the Gene Expression Omnibus (GEO) website (http://www.ncbi.nlm.nih.gov/geo/).

Analysis of ChIP-seq data

The Real Time Analysis software (RTA 3.4.4) was used for processing raw data files, the Illumina bcl2fastq v2.20 was used to demultiplex and convert binary base calls and qualities to fastq format. Samples were trimmed for adapters using Cutadapt (1.18) before the alignment. The trimmed reads were aligned with mm10 reference using Bowtie2(2.2.6) alignment. Library complexity is measured by uniquely aligned reads using picard (2.18.26)’s markduplicate utility. Peaks were called by MACS2(2.2.6) with –broad option on. HOMER (v4.11) findPeaks -style super was used to find super enhancers in ChIP-seq. Enrichment analysis of strain-specific enhancer regions based on biological process Gene Ontology (GO) terms was conducted using the Genomic Regions Enrichment of Annotations Tool (GREAT) (http://great.stanford.edu/). Bigwig files are generated by bamCoverage(deeptools/3.5.0), with the parameter --binSize 10 –normalize using RPGC.

ATAC-seq library construction

The ATAC-seq library was created followed previously established procedures with some minor adjustments^54^. To start, 100,000 macrophage cells were subjected to lysis using 100 μl of Lysis buffer (containing 10 mM Tris-Cl at pH 7.4, 10 mM NaCl, 3 mM MgCl2, 0.1% NP-40, and 0.1% Tween-20). Following a centrifugation at 500g for 10 minutes at 4°C, the resulting nuclei were reconstituted with 50 μl of a transposition mixture (FC-121–1030, Illumina; comprising 1X Tagment DNA buffer, Tn5 Transposase, and nuclease-free H2O), and this mixture was then incubated for 30 minutes at 37°C in a thermomixer. The transposed DNA was subsequently purified using MinElute columns (28004, QIAGEN) and then amplified using Nextera sequencing primers and the NEB high-fidelity 2X PCR master mix for 11 cycles (M0541, New England Biolabs). PCR-amplified DNA was purified using MinElute columns (28004, QIAGEN) and ATAC-seq samples were pooled and sequenced on NextSeq 2000 P2 using paired-end sequencing.

Analysis of ATAC-seq data

The Real Time Analysis software (RTA 3.9.25) was used for processing raw data files, the Illumina bcl2fastq v2.17 was used to demultiplex and convert binary base calls and qualities to fastq format. Samples were trimmed for adapters using Cutadapt (1.18) before the alignment. The trimmed reads were aligned with mm10 reference using Bowtie2(2.2.6) alignment. Library complexity is measured by uniquely aligned reads using picard (2.18.26)’s markduplicate utility. The peaks were called using Genrich (0.6.1) (genrich -t sample_sortedByRead.bam -o sample.narroPeak -j -y -r -v -d 150 -m 5 -e chrM,chrY -E blacklist_regions.bed-f sample.bdg -b sample.bed). Differentially accessible peaks after IL-4 treatment were analyzed by DiffBind (2.10)(http://bioconductor.org/packages/release/bioc/vignettes/DiffBind/inst/doc/DiffBind.pdf),w ith FDR ≤0.001 define induced peaks at different time point compared with Ctrl. HOMER (v4.11) findMotifsGenome.pl was used to investigate the motif enrichment of peaks region selected. chromVAR^23^ was used to find the motif enriched score in all ATAC-seq samples. Bedtools (2.30.0) intersect option are used to find the overlap of the ATAC-seq and H3K27ac ChIP-seq. Co-expression networks were constructed using the WGCNA (Weighted gene coexpression network analysis) algorithm implemented in R WGCNA package^33^. Raw expression values of selected differential changed genes were normalized with variance StabilizingTransformation (vst) function in R software. To obtain co-expressed modules, the parameters were adjusted to minModuleSize = 30 to cut the tree. For genome tracks, bigwig files for heatmap production were created from bam files using deepTools with the parameter --binSize 10 –normalize using RPGC. Heatmaps displaying normalized read densities of ATAC-seq and ChIP-seq samples were generated with the computeMatrix and plotHeatmap modules of the deepTools package^55^. Genome tracks were explored and visualized using the IGV browser.

RNA isolation and RNA-seq library construction

Total RNA was isolated using the RNeasy Plus Mini Kit (74134, QIAGEN). mRNA-seq samples were pooled and sequenced on NextSeq 2000 P2 using Illumina^®^ Stranded mRNA Prep and paired-end sequencing. Three independent replicates were used for RNA-seq analysis.

RNA-seq analysis

Reads of the samples were trimmed for adapters and low-quality bases using Cutadapt (1.18) before alignment with the reference genome (mm10) and the annotated transcripts using STAR (2.6.0c). The mapping statistics are calculated using Picard software (2.18.26). Library complexity is measured in terms of unique fragments in the mapped reads using Picard (2.18.26)’s MarkDuplicate utility. In addition, the gene expression quantification analysis was performed for all samples using STAR/RSEM (1.3.2) tools. EBSeq (1.24.0) was used to produce normalized reads for all samples and to identify differential RNA expression in BL/6 and BALB/c samples compared with Ctrl or LPS. To see the significantly trend different after IL-4 treatment or LPS treatment or IL-4+LPS treatment compared with Ctrl between BALB/c and BL/6 mice, we used DESeq2^56^. Differentially expressed genes were identified using DESeq2 following the model:

$$[eqn]$$

All differentially changed genes were selected based on an FDR ≤0.05 after DESeq2 analysis; both genes were selected that were more or less than 2-fold up for further analysis. k-means analysis of RNA expression data was carried out in MATLAB using fold-change compared with Ctrl for IL-4 treat samples in Fig.1F, fold-change compared with LPS treatment in Fig.3C, and fold-change compared with Ctrl in Fig.S3C, with correlation as the distance metric, repeat clustering set to 5, and other parameters set to default. The programs findMotifs.pl were used to identify transcription factor binding motifs within selected promoter regions (−400bp to +100 bp from TSS). For visualization, bigwig files were generated by deepTools with the parameter --binSize 10 --normalizeUsing RPGC. Genome tracks were explored using the integrative genomics viewer (IGV) browser.

Hi-C libraries and sequencing.

Chromatin interaction (Hi-C) libraries preparations were performed by Arima Genomics (https://arimagenomics.com/) using the Arima-HiC kit that uses two enzymes (P/N: A510008). The resulting Arima-HiC proximally ligated DNA was then sheared, size-selected around 200–600bp using SPRI beads and enriched for biotin-labelled proximity-ligated DNA using streptavidin beads. From these fragments, Illumina-compatible libraries were generated using the KAPA Hyper Prep kit (P/N: KK8504). The resulting libraries were PCR amplified and purified with SPRI beads. The quality of the final libraries was checked with qPCR and TapeStation Systems. HiC-seq samples were pooled and sequenced on NovaSeq and paired-end sequencing.

Hi-C data analysis.

All Hi-C data were processed with Juicer (v.1.6^29^) pipelines with default settings. Sites for Phase Genomics cutting enzyme (GATC) were detected using the generate_site_positions.py, Sites for Arima Genomics cutting enzyme (ĜATC, GÂNTC) were obtained from [Arima restriction enzyme files. Accessed 9 April 2020.ftp://ftp-arimagenomics.sdsc.edu/pub/HiCPro_GENOME_FRAGMENT_FILES]. Juicebox (2.3.0) was used to visualize the Hi-C heatmap. We defined TAD boundaries as regions 100 kb upstream of the TAD start and 100 kb downstream of the TAD end. Strain-specific loops were identified using Bedtools (2.30.0). Subsequently, we located promoter-based loops and annotated the promoter-distal regions using the HOMER (v4.11). Next, we overlapped the promoter distal regions with strain-specific highly expressed genes, as defined by DESeq2 analysis, for motif analysis in Fig. 4A. We employed MEME-FIMO (5.4.1) to scan for the locations of strain-specific active transcription factors (TFs) within promoter-distal regions (enhancers) that are linked to strain-specific loops associated with highly expressed genes. The presence of motif sites within the enhancer region indicates a regulatory relationship between the TF and the gene.

scRNA-seq libraries and sequencing

Using a Chromium Next GEM Single Cell 3’ Kit v3.1, 4 rxns (1000269, 10X Genomics), Chromium Next GEM Single Cell 3ʹ Gel Bead Kit v3.1, 16 rxns (PN-1000122) and Chromium Next GEM Chip G Single Cell Kit (PN-1000120, 10X Genomics), the peritoneal cells suspensions in PBS, at 2000 cells per ul, were loaded onto a Chromium single cell controller (10X Genomics) to generate single-cell Gel beads-in-emulsion (GEMs) according to the manufacturer’s protocol. Briefly, approximately 10,000 cells per sample were added to a chip to create GEMs. Cells were lysed and the bead captured poly(A) RNA was barcoded during reverse transcription in Thermo Fisher Veriti 96-well thermal cycler at 53°C for 45 min, followed by 85°C for 5 min. cDNA was generated and amplified. Quality control and quantification of the cDNA was conducted using Agilent Genomic DNA ScreenTape Analysis kit (5067–5366 for Genomic DNA Reagents and 5067–5365 for Genomic DNA ScreenTape) in the TapeStation system. scRNA-seq libraries were constructed using a Chromium Next GEM Single Cell 3ʹ Library Kit v3.1 (PN-1000158, 10X Genomics) and Single Index Kit T Set A, 96 rxns (PN-1000213, 10X Genomics). 10x Genomics Single Cell Gene Expression libraries were sequenced on a NextSeq 2000 run. The sequencing run was setup as a 28 cycles + 90 cycles non-symmetric run. Demultiplexing was done allowing 1 mismatch in the barcodes.

scRNA-seq data processing

The analysis was performed with the Cell Ranger 7.1.0 software using the default parameters. The number of cells captured ranges from 4,862 to 9,060 and mean reads per cell ranges from 18,629 to 52,533. Cells with extremely low number of UMI counts were filtered out. The results from Cell Ranger were processed in R v3.6.3 with Seurat v3.2.3^57^ using default parameters, unless otherwise specified. Souporcell (https://github.com/wheaton5/souporcell) was used to separate the cells from CB6F1/J (Jax no. 100007) into BL/6 and BALB/c cells with singularity exec function. Quality control filtering was applied to each sample to eliminate downstream analysis of empty droplets, low-quality cells, and potential doublets. Then, we used the Seurat (v3.2.3) package to perform calculation of the number of unique genes detected in each cell (nFeature_RNA), the total number of molecules detected within a cell (nCount_RNA), and the percentage of reads that map to the mitochondrial genome. The differential peaks in Fig. 7E were identified using the FindMarkers function in Seurat (v3.2.3) with the following parameters: logfc.threshold = 0.25, only.pos = TRUE, test.use = ‘MAST’. To identify the potential key transcriptional regulators after IL-4 treatment, we used pySCENIC (v0.10.)^58^ to explore the potential transcriptional regulators.

Multiome-seq assay (snRNA-seq and snATAC-seq) from peritoneal cells

Peritoneal cells were isolated from both BL/6 and BALB/c mice after IL-4 treatment. Nuclei isolation for snRNA-seq and snATAC-seq was optimized and performed following the demonstrated protocol CG000365 from 10x Genomics. The libraries were prepared following the 10x Genomics user guide CG000338 for Single Cell Multiome ATAC + Gene Expression Reagent Bundle (PN-1000285,10xGenomics). Only the ATAC Library Construction is performed by Sample Index Plate N, Set A (3000427, 10xGenomics). snRNA-seq libraries failed QC and were not sequenced. 10x Genomics Single Cell ATAC libraries were sequenced on a NovaSeq 6000 run.

scATAC-seq data analysis

The analysis was performed with the Cell Ranger ATAC 2.0.0 software using the default parameters. Quality control of scATAC-seq data from BL/6 and BALB/c strains was performed in R (v.4.3.0) using Signac (v.1.1.0)^59^. Rliger (1.0.0)^60^ was used to integrate scATAC-seq and scRNA-seq data. Parameters included k=20 for the optimizeALS function, n_neighbors=10, min_dist=0.3, and cosine distance for the runUMAP function. A resolution of 0.3 was used for the Louvain community detection. Cell types were labeled by cross-referencing the integrated scATAC-seq dataset with scRNA-seq datasets. Cell labels were transferred using the Seurat (4.9.9) package, focusing on the macrophage subset for downstream analyses.

Enhancer gene regulatory networks (eGRNs) were identified using SCENIC+^35^, an algorithm utilizing paired single-cell transcriptomic and open chromatin data. In our analyses, we used a search space of 150 kb for inferring region-gene relationships and retained regulons with a minimum of 10 target genes. For pycisTopic (1.0.3), we executed the essential steps of the SCENIC+ workflow, including object creation, topic modeling, dimensionality reduction, dropout imputation, and DAR (Differential Accessibility Region) inference using default parameters. The topic modeling process employed serial LDA (Latent Dirichlet Allocation) with Collapsed Gibbs Sampler, with 500 iterations. The generated topics spanned from two to 48, and the final model selection yielded 10 topics for BALB/c and 6 topics for BL/6.PycisTarget (1.0.3) was run using default parameters, incorporating cisTarget and differential enrichment of motif (DEM). The analysis used the bulk consensus peaks motif databases. SCENIC+ was run with default parameters using http://sep2019.archive.ensembl.org/ as the BioMart host. High quality regulons were selected based on the correlation between gene-based regulon AUC and region-based regulon AUC (>0.7).

Statistical analysis

In Fig. 7A, the p-values were calculated using a two-sided Mann-Whitney test. In others figures, p-values were calculated using an unpaired t-test, two-sided. Data are presented as mean ± standard deviation. Gene expression comparisons were reported by adjusted p values (i.e., FDR) from DESeq2^56^.