Epigenetic Analysis of the Chromatin Landscape Identifies a Repertoire of Murine Eosinophil-Specific PU.1-Bound Enhancers

The long-held view of eosinophils as end-stage cells involved primarily in allergy and host protection against parasites has been transformed by demonstrating that eosinophils are pleiotropic, multifunctional leukocytes involved in a plethora of biological processes, including tissue homeostasis, mammary gland development, estrus cycling, Ag presentation, cytokine production, epithelial remodeling, fat metabolism, and cancer (1–8). Moreover, two subsets of regulatory eosinophils have recently been identified (2, 9), with transcriptomes and cellular phenotypes distinct from those of the conventional, inflammation-recruited eosinophil. The potency and resulting diversity fundamental to an eosinophil’s complex immunomodulatory functions and tissue specialization likely result from dynamic epigenetic regulation of the genome, a dynamic eosinophil regulome. However, the gene regulatory elements (REs) that orchestrate the eosinophil regulome remain largely unknown.

The importance of key transcription factor (TF) families in eosinophil lineage commitment and maturation has been well documented. Most notably, there are critical roles for GATA, C/EBPε, and PU.1 in eosinophil development, granule production, and function (10–12). However, the role of these TFs in regulating enhancers has not been deeply studied. The development of new technologies that enable the global identification of genomic elements and the TFs that regulate expression of associated genes has provided opportunities to discern how cis REs, such as enhancers, control gene expression programs under normal and disease states (13). Epigenetic features enable identification of enhancer-like regions within cell populations, facilitating systematic investigation of their importance in gene transcription and cellular phenotype. Promoters are often identified by the presence of open regions marked by histone H3 lysine 4 trimethylation (H3K4me3), surrounded by regions marked by histone H3 lysine 4 monomethylation (H3K4me1) and located proximal (within 1 kB) to the transcription start site (TSS). In contrast, enhancers are often identified by open regions marked by H3K4me1 and are distal (>1 kB) to the TSS (14–19). Subsequently, enhancers can be classified as “poised” (lacking enriched acetylation of histone H3 lysine 27 [H3K27ac]) or active (enriched for H3K27ac). Despite having a clear practical utility, global analysis and chromatin mapping of mature eosinophils has been limited, predominantly because of difficulties in isolating intact RNA and chromatin from these sensitive, lysis-prone, enzyme-rich granule-containing cells. At present, there is only one published study investigating H3K4me3 marks within promoters of eosinophils and eosinophil progenitors (20). Thus, the use of bioinformatic and epigenetic tools to investigate eosinophil cell populations is underrepresented in the literature.

Taking a global approach to identify the REs present in mature, bone marrow culture–derived murine eosinophils (bmEos), we performed chromatin immunoprecipitation (ChIP) coupled with massively parallel sequencing (ChIP-seq) for H3K27ac, H3K4me1, and H3K4me3; assay for transposase-accessible chromatin (ATAC) with massively parallel sequencing (ATAC-seq); and RNA sequencing (RNA-seq). Using these approaches, we defined a repertoire of eosinophil-specific enhancers that are located within 1–20 kB of genes expressed by mature eosinophils and then identified associated TFs. Throughout this manuscript, we use distance to classify promoter and enhancers, with promoters located proximally (within 1 kB) and enhancers located distally (1–20 kB) to the TSS of expressed genes. We focused our study on enhancers and excluded from analysis epigenetic marks detected within 1 kB of the TSS. Our analysis revealed the presence of numerous distal enhancers located predominantly within intronic regions. Gene ontology (GO) analysis revealed significant enrichment for genes associated with leukocyte migration, activation, chemotaxis, proliferation, degranulation, and cell death, highlighting that important functional genes may be dynamically regulated by distal enhancers. DNA binding site motifs for PU.1 were significantly enriched at enhancers that were predominantly near genes highly expressed in eosinophils, and PU.1 binding at these loci was further confirmed by ChIP-seq. Comparing our data with publicly available macrophage (21), neutrophil (22), and immature granulocyte (ImmGran) (23) datasets revealed that >25% of the PU.1-bound regions were unique to eosinophils. Furthermore, eosinophil-specific superenhancers [regions of the genome containing multiple strong enhancers that have been shown to correlate with cell identity (24)] were enriched for genes associated with genome-wide association study (GWAS)–derived risk loci for eosinophilia and allergic diseases.

Collectively, our data identify eosinophil-specific enhancers likely regulating key eosinophil genes through epigenetic mechanisms (open chromatin, H3K4me1, H3K4me3, and H3K27ac) and TF binding (PU.1). These findings are supportive of a dynamic, epigenetic eosinophil regulome underlying eosinophil diversity. Future comprehensive analyses of eosinophil epigenetic and transcriptional profiling from homeostatic, developmental, or disease states will further advance our understanding of allergic and eosinophil-related disorders and the potential development of therapeutic strategies.

Six- to eight-week-old BALB/c mice were used for all experiments. Animals were housed under standard specific pathogen-free conditions and handled under approved protocols of the Institutional Animal Care and Use Committee of Cincinnati Children’s Hospital Medical Center (IACUC2018-0028, IACUC2016-0083).

Murine bone marrow–derived eosinophils (bmEos) were generated from unselected bone marrow progenitor cells as first described by Dyer et al. (25, 26) and were used for all experiments. After 14 d, eosinophil maturity was assessed on the basis of cellular morphology and Siglec-F⁺CCR3⁺ positivity as determined by flow cytometry. Eosinophils were stained with Siglec-F–Alexa Fluor 647 (1:100, BD Biosciences) and CCR3-PE (1:100, R&D Systems) for 25 min and Zombie NIR Fixable Viability (1:1000, BioLegend) for an additional 5 min. Cells were washed prior to analysis on FACSCanto (BD Biosciences). Cell cultures that were >95% Siglec-F⁺CCR3⁺ were harvested and used for downstream experiments.

Eosinophil RNA was isolated from 2–5 × 10⁶ bmEos by double chloroform extraction prior to purification with the Direct-zol RNA MiniPrep (Zymo Research). RNA quality was assessed via Agilent RNA Nano Chip, and only samples with an RNA quality number >8 were sent for sequencing analysis. Libraries were constructed with the TruSeq Stranded mRNA kit (Illumina) and sequenced on the Illumina NovaSeq 6000 (S1 Flow Cell) or Illumina HiSeq 2500, paired-end, 100 bases for each end by the Cincinnati Children’s Hospital Medical Center DNA Sequencing and Genotyping Core (Cincinnati, OH), targeting ∼ 30 million reads per sample.

A total of 6.5 × 10⁴ viable bmEos was pelleted at 500 g for 5 min at 4°C prior to resuspension in 50 μL ATAC-Resuspension Buffer (10 mM Tris-HCl, 10 mM NaCl, and 3 mM MgCl₂ made up to 50 ml in ddH₂O) with added 0.1% NP-40, 0.1% Tween-20, and 0.01% digitonin. Samples were incubated on ice for 3 min prior to washing in 1 ml cold ATAC-Resuspension Buffer containing 0.1% Tween-20 only. Nuclei were pelleted at 500 g for 10 min at 4°C prior to the addition of 50 μL transposition mixture (2× Tagmentation DNA buffer, 100 nM transposase enzyme, 5 μL ddH2O, 1% digitonin, and 10% Tween-20 made up to 50 ml in 1× PBS) and incubated at 37°C for 30 min. After tagmentation, the transposased DNA was purified using the MiniElute Reaction Cleanup Kit (28204, QIAGEN). Libraries were constructed with the OMNI-ATAC protocol (27) and sequenced on the Illumina HiSeq 4000, paired-end, 150 bases each end (Novogene).

A total of 1–3 × 10⁷ viable bmEos was pelleted, and eosinophil chromatin was cross-linked by the addition of 1:10 of 10× cross-linking solution (8.8% formaldehyde, 0.1 M NaCl, 1 mM EDTA, 1 mM EGTA, and 50 mM HEPES made up to a total volume of 1 ml in ddH₂O) and incubated on ice for 5 min. The addition of 1:20 2.5 M glycine was used to stop the cross-linking reaction. The cells were centrifuged at 2000 g for 5 min at 4°C, after which the pellets were washed twice with PBS^−/− and frozen at −80°C until use. Following storage, the cell pellets were thawed on ice and resuspended in 1 ml L1 lysis buffer (50 mM HEPES, 140 mM NaCl, 1 mM EDTA, 10% glycerol, 0.5% NP-40, and 0.25% Triton X-100 made up to a volume of 50 ml in ddH₂O) with added protease inhibitor (1:100) and sodium butyrate (NaButyrate) (1:1000, stock 1 M) and rotated at 4°C for 10 min. Cells were centrifuged at 1900 g for 3 min at 4°C prior to the addition of 1 ml L2 buffer (200 mM NaCl, 1 mM EDTA, 0.5 mM EGTA, and 10 mM Tris made up to a volume of 50 ml in ddH₂O) with added protease inhibitor (1:100) and NaButyrate (1:1000, stock 1 M) and rotated at room temperature for 10 min. Cells were centrifuged at 2000 g for 5 min at 4°C to remove supernatant, followed by two additional washes with Tris-EDTA + 0.1% SDS to remove all excess salts, which inhibit the sonication process. The pellet was resuspended in 1 ml Tris-EDTA + 0.1% SDS and cooled on ice for 5 min prior to sonication using the Covaris S-series sonicator at 170.0 peak power, 10.0 duty factor, and 200 burst. A small amount (10 μL) of the chromatin suspension was saved for agarose gel analysis to confirm adequate sonication of the sample, and salts were readded (1% Triton X-100, 5% glycerol, 0.5M NaCl, 0.1% NaDOC, 1:100 protease inhibitors, and 1:1000 NaButyrate) to the remaining suspension. The chromatin was then precleared following incubation with Protein-A beads prior to performing the ChIP with Abs for H3K4me1 (C15410194, Diagenode), H3K4me3 (17-614, MilliporeSigma), H3K27ac (C15410196, Diagenode), and PU.1 (sc-352, Santa Cruz Biotechnology) overnight in the SX-8G IP-Star ChIP Robot (Diagenode). Libraries were constructed by ChIPmentation (28) and sequenced on the Illumina HiSeq 4000, single-end or paired-end, 150 bases per end (Novogene).

RNA-seq data were processed using the pipeline NextGenAligner (https://github.com/MarioPujato/NextGenAligner). Briefly, the raw sequence read quality was analyzed using FastQC/0.11.7 (https://github.com/s-andrews/FastQC). Adapter sequences were trimmed using Cutadapt/1.8.1 (29) and Trim Galore/0.4.2 (https://github/FelixKreuger/TrimGalore) prior to alignment of sequence reads to the mm9 genome and Reference Sequence database–based transcriptome using STAR/2.5 (30). Duplicates were removed using Picard/1.89 (http://broadinstitute.github.io/picard). FeatureCounts (31) and the calculateFPKM module were used to quantify gene transcript fragments per kilobase per million mapped reads (FPKM). To characterize gene expression by eosinophils, thresholds were determined using expression levels of known expressed eosinophil genes, with ≥5 FPKM in both replicates indicating genes that are expressed by eosinophils and <5 FPKM indicating genes that are silent, repressed, or minimally expressed. GO and GWAS analysis was performed using Enrichr (32, 33). ATAC and ChIP sequencing analysis was performed as above except for read alignment, which was performed with BOWTIE/2.3.4.1 (34). Peaks were called using MACS/2.1.0 (J.M. Gaspar, manuscript posted on bioRxiv) with a q value of 0.01. Read and peak information for ATAC-seq and ChIP-seq samples are in Supplemental Table I. The presence of ATAC-seq or ChIP-seq peaks within 20 kB of expressed eosinophil genes was determined using BEDtools/2.27.0 (35). ATAC-seq and ChIP-seq tracks were visualized using the Integrative Genomics Viewer (IGV) browser (36). Differential binding analysis of ChIP-seq and ATAC-seq peaks was performed with differential binding and occupancy analysis (DiffBind)/2.16.0 (37).

TF binding motif enrichment analyses of identified REs were performed using the HOMER software package (38). Briefly, HOMER uses a library of >7,000 TF binding models [in the form of position weight matrices, taken from the CisBP database (39)] to scan a set of input sequences for statistical enrichment of each position weight matrix. Calculations are performed using zero or one occurrence per sequence scoring coupled with hypergeometric enrichment analysis to determine motif enrichment. Input eosinophil enhancer sequences were also assessed for statistical enrichment of motifs for PU.1 binding sites using the findPeaks module and factor mode within HOMER. Significantly enriched TF binding site motifs are expressed as log p values. Superenhancer data were generated using the ROSE algorithm (24, 40) within Scientific Data Analysis Platform (https://scidap.com; Datirium) (41, 42).

Data were analyzed with GraphPad Prism software (v8; GraphPad Software), with flow cytometric data analyzed using FlowJo software (TreeStar). A p value < 0.05 was considered statistically significant.

The next-generation sequencing data are accessible through National Center for Biotechnology Information’s Gene Expression Omnibus (43) GSE144951 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE144951. The University of California, Santa Cruz Genome browser has been generated and is accessible at http://genome.ucsc.edu/s/feljj4/FELTON_mm9_bmEos.

To assess the chromatin landscape profile of mature, culture-derived, murine eosinophils (bmEos), open areas of chromatin and histone modifications for H3K4me1, H3K4me3, and H3K27ac were assessed across the genome. Prior to downstream analysis, eosinophils were evaluated for maturity, and >95% of cells were found to express CCR3 and Siglec-F (data not shown). Biological replicates for each of the four epigenetic marks showed a high correlation in peak location and strength as assessed by differential binding analysis of ATAC-seq and ChIP-seq peak data (DiffBind) output for correlation and replication of murine next-generation sequencing datasets (Supplemental Fig. 1A–B) and per-peak read count (Supplemental Fig. 1C). Biological replicates were subsequently pooled for downstream analysis. Next, we paired the chromatin accessibility and histone findings (ATAC-seq and ChIP-seq) with eosinophil transcriptomic data (RNA-seq).

Gene expression analysis revealed that 6288 genes were expressed (≥5 FPKM) by mature culture-derived bmEos. These genes included well-described marker genes for eosinophils, such as major basic protein (MBP, Prg2), eosinophil cationic protein (Ear1), eosinophil peroxidase (Epx), CCL11 (eotaxin) receptor (CCR3, Ccr3), IL-5Rα (Il5ra), Siglec-F (Siglec5), the TFs PU.1 (Spi1) and GATA-1 (Gata1), and the Ikaros TF family members Aiolos (Ikzf3) and Helios (Ikzf2). In contrast, >18,000 genes were considered to be silent, repressed, or minimally expressed (<5 FKPM, data not shown). Expressed eosinophil genes had higher chromatin accessibility and H3K4me1, H3K4me3, and H3K27ac levels in comparison with the silent, repressed, or minimally expressed genes (Supplemental Fig. 1D). Expressed genes were associated with the classical “see-saw” promoter profile (open chromatin, H3K27ac marks, and H3K4me3 marks surrounded by H3K4me1) within 1 kB of the TSS (44), whereas the silent, repressed, or minimally expressed genes did not present with this profile (Supplemental Fig. 1E).

Active promoters are traditionally identified by open regions of chromatin (ATAC) comarked with H3K4me3 and H3K27ac (Fig. 1A, top). Interestingly, when assessing the proportion of eosinophil genes found to possess active promoters within 1 kB of the TSS, we observed that there was an almost even split of active promoters lacking or possessing H3K4me1 at the H3K4me3 peak location (Fig. 1a, bottom). Of the 6288 expressed eosinophil genes, 2542 (40.43%) were found lacking H3K4me1 at the peak location (H3K4me1⁻), whereas 2776 (44.15%) were found to possess H3K4me1 (H3K4me1⁺) (Fig. 1B). When we compared the expression level of eosinophil genes associated with an H3K4me1⁺ or H3K4me1⁻ active promoter, it was observed that genes with an H3K4me1⁻ active promoter were expressed to a higher level than were those with an H3K4me1⁺ active promoter (p < 0.0001) (Fig. 1B). GO analysis of the 2542 genes with an H3K4me1⁻ active promoter revealed significant enrichment of genes associated with leukocyte-mediated immunity, migration, activation, and degranulation (Fig. 1C). In contrast, GO analysis of the 2776 genes with an H3K4me1⁺ active promoter revealed significant enrichment for genes associated with mRNA processing, gene expression, and translation (Fig. 1D). Interestingly, 970 (15.4%) of expressed eosinophil genes did not possess an H3K27ac-marked promoter within 1 kB of the TSS.

Next, we interrogated the regulation of key eosinophil genes by comparing the prevalence of active REs, identified by the presence of H3K27ac in open chromatin (REs lacking H3K27ac in open chromatin are classified as poised and were not included in this analysis) and TF binding motifs between the top 10% (653 genes, n = 2, average FPKM range 86.5–4870.8) and bottom 10% (653 genes, n = 2, average FPKM range 5–6.9) of expressed eosinophil genes. We hypothesized that highly expressed eosinophil genes would be enriched for H3K27ac relative to lowly expressed eosinophil genes. As expected, the top 10% of expressed eosinophil genes (Supplemental Fig. 2A, blue line) had an increase in H3K27ac signal compared with that of the bottom 10% of expressed eosinophil genes (Supplemental Fig. 2A, green line). Analysis of TF binding site motifs revealed preferential enrichment for PU.1 binding motifs in the top 10% of expressed eosinophil genes compared with the bottom 10% of expressed genes (Supplemental Fig. 2B, 2C). Highly expressed eosinophil genes were also preferentially enriched for YY1, FOS, JUN, and CEBP family members, albeit to a much lower extent than PU.1 (Supplemental Fig. 2B, 2C). Gene expression analysis revealed PU.1 to be the second most highly expressed TF in bmEos (Supplemental Fig. 2D), further supporting its importance in the eosinophil regulome.

To determine whether the proximal and/or distal REs were driving the PU.1 epigenetic signature, we examined DNA sequences within active (H3K37ac⁺) chromatin regions on the basis of genomic location; promoter regions (±1 kB from TSS) and enhancers (1–20 kB from TSS) were identified by H3K27ac density to assess TF binding site motif enrichment. Interestingly, HOMER analysis revealed that the DNA binding site motifs for PU.1 were the most enriched within the distal (1–20 kB from TSS) active chromatin regions identified as enhancers (Supplemental Fig. 2E) compared with the promoters (±1 kB of the TSS) (Supplemental Fig. 2F). Together, these data suggest that regulation of key eosinophil genes occurs by a select repertoire of TFs, including PU.1, at epigenetically regulated distal enhancers.

Enhancers are described as discrete regions of DNA that increase transcriptional activity of promoters from a distance, often many kilobases away, in part via the binding of lineage-specific TFs (45, 46). Because of 1) the key role of PU.1 in granulocyte maturation and function and 2) the observed enrichment of PU.1 in active enhancers distal (1–20 kB) to the TSS, we aimed to further characterize the chromatin modification profiles at these locations. Although the classical description of enhancers has an inverted profile of H3K4me1 and H3K4me3 (44), there are some noncanonical examples of H3K4me3 depositing at a subset of enhancers acting in lineage-specific functions (46). Considering this, we speculated that highly expressed genes associated with key eosinophil functional responses would have active enhancers with greater H3K4me3 deposition.

Notably, over 83% (5243/6288) of expressed eosinophil genes had at least one type of enhancer (active, poised, or inactive) within 1–20 kB of its TSS; genes associated with active enhancers were expressed at a higher level than those associated with poised enhancers (Fig. 2A). Furthermore, we observed that a greater proportion of active enhancers than poised enhancers were marked by H3K4me3 deposition (Fig. 2B). Poised enhancers were identified by regions of open chromatin comarked with H3K4me1 (with or without H3K4me3) and lacking H3K27ac (Fig. 2C, top). Active enhancers were identified by regions of open chromatin comarked with H3K4me1 (with or without H3K4me3) and possessing H3K27ac (Fig. 2C, bottom). Active enhancers were enriched for gene pathways associated with granulocyte-mediated immunity, migration, and phagocytosis (Fig. 2D) to a greater degree than were poised enhancers (Fig. 2E). Key functional genes with at least one active distal eosinophil enhancer located within 1–20 kB included Il5ra, Gata2, Epx, Ccr3, Pgr2, Mcl1, Ikzf1, Ikzf2, Ikzf3, Irf2, and Stat5b. These data identify epigenetically marked REs consistent with active distal enhancers in the eosinophil regulome and are suggestive of their involvement in key lineage-specific eosinophil functions.

To elucidate which TFs are involved in the regulation of the eosinophil enhancer epigenome, we examined DNA sequences within active enhancers for TF binding site motif enrichment. HOMER analysis revealed that the DNA binding site motifs for PU.1 were the most enriched (Fig. 2F, Supplemental Fig. 3A), with both the motif score (Fig. 2G, left; p = 0.0014) and proportion of predicted binding sites (Fig. 2G, right; p < 0.001, χ² test) being greater in active than poised enhancers. Other highly represented DNA binding site motifs included those for ELK:ELF family members, CEBP isoforms, IRFs, and GATA. De novo analysis revealed enrichment for DNA binding site motifs for KLF family members, AP1, ELK:ELF family members, and SMAD3 (Supplemental Fig. 3B).

We paired PU.1 binding sites (identified by ChIP-seq for PU.1) with enhancer data to determine the locations at which PU.1 is bound to enhancers across the eosinophil epigenome. The PU.1 motif was enriched in PU.1 ChIP-seq peaks compared with H3K27ac peaks alone or background (Fig. 3A). Over 95% (96.7%, 7970/8228) of active enhancers were bound by PU.1, in contrast to only 61.3% (3586/5848) of poised enhancers (Fig. 3B; p < 0.001, χ² test). To determine which TF may be facilitating PU.1 binding at the 7.2% of PU.1 peaks that did not possess a PU.1 or similar binding motifs, we performed de novo TF binding motif enrichment. First, all peaks with a GGAA within their genomic sequence were removed to ensure that ETS/ELK/ELF family TFs were not included because of the similarities in their core binding sequences. Interestingly, de novo motif enrichment revealed that FOXK1, members of the KLF family, STAT3, STAT5, SMAD3, and GATA1 were the most enriched de novo motifs (Fig. 3C), with a strong enrichment in the target sequences compared with background (Supplemental Fig. 3C). Of the total active enhancers with predicted PU.1 binding motifs identified, 83.2% (5374/6461) were confirmed by ChIP-seq (data not shown), including the active enhancer identified within the Il5ra gene (Fig. 3D).

To determine which of the identified PU.1-bound active enhancers were unique to murine culture-derived eosinophils, we compared our datasets with publicly available data (H3K27ac, PU.1 ChIP-seq and RNA-seq) from other myeloid cell types: bone marrow–derived macrophages [mBMDM: SRR502912, SRR502934, and RNA-seq SRR502948 (21)], neutrophils [SRR6751044, SRR6751056, and RNA-seq SRR6751062 (22)], and ImmGrans [SRR61173 and SRR611738 (23)]. No paired RNA-seq data were available for ImmGrans. Because of differences in mouse strain/background, we generated bmEos ATAC and H3K27ac datasets from mice on a BL6/C57 background to assess the correlation and overlap of called peaks between background strains. A strong correlation and overlap of peak locations was observed, with the most influential variable being epigenetic mark/assay type (PC1 = 61%; ATAC or H3K27ac) over mouse strain (PC2 = 14%; BL6/C57 or BALB/c) (Supplemental Fig. 3D–E), with 71.3–82.3% of ATAC and H3K27ac ChIP-seq peaks shared between BALB/c and BL6/C57 datasets (data not shown). Furthermore, the strongest correlation and clustering between PU.1 ChIP-seq datasets were within PU.1 myeloid datasets (Supplemental Fig. 3F).

DiffBind of active PU.1-bound bmEos enhancers compared with active PU.1 cobound regions in mBMDMs, neutrophils, and ImmGrans revealed 1691 PU.1-bound active enhancers that were unique to eosinophils (26%) (Fig. 4A). Eosinophils had ∼70% of their PU.1-bound enhancers in common with at least one of the other myeloid cell types. A total of 1445 PU.1-bound enhancers were common between all cells of the myeloid lineage assessed, describing a shared myeloid PU.1 regulome. The 1691 unique, active PU.1-bound, eosinophil-specific enhancers were associated with a total of 1300 genes with increased expression (FPKM) in bmEos (compared with mBMDMs [p = 0.0257] and neutrophils [p = 0.0219]) (Fig. 4B), and GO analysis revealed enrichment for gene pathways associated with myeloid differentiation, granulocyte activation, degranulation, and migration (Fig. 4C), including but not limited to Prg2, Epx (Fig. 4D), Ikzf1, Ikzf2, Ikf3 (Supplemental Fig. 4A), Cebpa, Cebpb, Cebpe (Supplemental Fig. 4B), Trib1 (Fig. 5D), Smad3, Cd44, Cd9, Ccl9, Ccl17, Il5ra (Supplemental Fig. 4C), Ear6, Grb7, Klf13, Nfkb1, Elmo1, and Gata1 (Supplemental Fig. 4D).

Together, these data demonstrate that relevant, functionally expressed genes are associated with PU.1-bound, lineage-specific, and cell-specific enhancers.

Last, we used the H3K27ac datasets to identify superenhancers, defined as regions of the genome containing multiple strong enhancers that have been shown to correlate with cell identity (24). To this end, we identified eosinophil superenhancers using the ROSE algorithm (24). Eosinophil superenhancers were associated with key genes involved in maturation and function, including Trib1, Ikzf2, Spi1 (PU.1), Il4, Ikzf1, Prg2 (proMBP), and Cebpe (Fig. 5A, top left). Superenhancers for mBMDMs, neutrophils, and ImmGrans were also identified using the ROSE algorithm (Fig. 5A). DiffBind analysis of superenhancers from each cell type (bmEos, neutrophils, mBMDMs, and ImmGrans) revealed a total of 643 eosinophil-specific superenhancers (Fig. 5B). Interestingly, eosinophils had the greatest proportion of cell type–specific superenhancers out of the four populations (Fig. 5C; p < 0.001, χ² test). Genes associated with eosinophil-specific superenhancers included but were not limited to Ccr3, Ikzf1, Ikzf2, Trib1 (Fig. 4E), Cebpe (Supplemental Fig. 4B), Smad3, and Prg2. Furthermore, murine, eosinophil-specific superenhancers were enriched in genes whose homologs are known to be associated with GWAS hits for eosinophil levels (eosinophilia, percentage of WBCs, or granulocytes) and allergic diseases in humans (Fig. 5D).

Together, these data demonstrate that relevant, functionally expressed genes are associated with lineage-specific and cell type–specific superenhancers, which are enriched in genes associated with allergic disease risk loci.

It is becoming increasingly appreciated that eosinophils are pleiotropic, multifunctional cells involved in a plethora of biological processes, suggesting the presence of complex cell autonomous pathways that determine eosinophil gene expression and function. As the genomic regulatory landscape of eosinophils was previously unknown, we took a global approach to analyze genome-wide transcriptome and epigenome of mature, murine, culture-derived bmEos to describe the cell’s enhancer profiles. Our aim was to identify and define the characteristics of unique, eosinophil-specific enhancers controlling expression of key genes, in addition to identifying TFs crucial for regulating eosinophil gene expression at the identified enhancers.

Active enhancers likely have roles in key eosinophil pathways and lineage determination. Highly expressed eosinophil genes were enriched for active (H3K27ac⁺) distal (1–20 kB of TSS) enhancers, particularly in genes associated with migration, activation, proliferation, and degranulation, suggesting that these key pathways may be regulated by enhancers. Furthermore, deposition of H3K4me3 was observed at a greater portion in active enhancers, suggesting that these could be involved in regulating lineage-defining functions (46). Russ et al. (46) observed a greater proportion of H3K4me3 deposition at predominantly poised enhancers (H3K27ac⁻) in naive T cells, yet we observed a greater proportion of H3K4me3 deposition in active (H3K27ac⁺) enhancers. When considering the terminally differentiated state of eosinophils in contrast to naive T cells, we speculate that the lineage-defining developmental processes regulated by these enhancer regions were activated at an earlier point in eosinophil development. Future studies will focus on delineating the epigenetic changes that occur during eosinophilopoiesis in addition to postmaturation activation using a range of cytokines and inflammatory mediators.

Correspondingly, the observed increase in the expression of genes possessing an H3K4me1⁻ active promoter, compared with those possessing an H3K4me1⁺ active promoter, further highlights the functional specificity of RE subtypes. These findings suggest that traditional (H3K4me1⁻) active promoters regulate highly expressed genes enriched for cell-specific functions (e.g., granulocyte activation and degranulation), whereas nontraditional (H3K4me1⁺) active promoters regulate lowly expressed genes involved in general cellular functions (e.g., regulation of translation and mRNA processing). Furthermore, the differential enrichment of TF binding motifs demonstrated between proximal promoters (±1 kB of the TSS) and distal enhancers (1–20 kB from TSS) demonstrates that select TF are involved in these discrete processes.

TFs have a major role in determining cell maturation and development in both cell type– and lineage-specific contexts. The most prominent of those involved in eosinophil maturation and lineage commitment are PU.1, GATA-1, and members of the CEBP family (10, 11, 20, 47, 48). Our data are consistent with these findings, showing high expression levels for all of these TFs, in addition to enrichment for their binding motifs within identified eosinophil enhancers. Interestingly, previous data published by our group showed enrichment for binding site motifs of members of the IKZF TF family (IKZF2 and IKZF3) in promoters of eosinophil progenitors and FACS-sorted, bone marrow–resident eosinophils (20). Despite the elevated transcript levels of these TFs and their association with both PU.1-bound eosinophil-specific enhancers and superenhancers, enrichment for their binding motifs in distal enhancers of cultured eosinophils was relatively low. This potentially suggests that 1) IKZF family members regulate eosinophil gene expression at an earlier time point in eosinophil development and maturation and may not have as prominent a role in mature, cultured cells under the direction of IL-5, or 2) IKZF family members have a more prominent role in the regulation of gene transcription in eosinophils at promoters (proximal REs). Similarly, the prevalence of PU.1 binding at enhancers in mature eosinophils supports the hypothesis that PU.1 is not essential for progenitor specification but instead has an important role in the terminal differentiation of mature eosinophils (49–53).

The roles of less prominent TFs in determining cell maturation and development in both cell type– and lineage-specific contexts are not to be discounted. Interestingly, de novo motif analysis of active enhancers revealed enrichment for TFs, including AP-1, KLF family members, STATs, SMAD3, and GATA. These same elements were also enriched in the subpopulation of PU.1-bound active enhancers without a canonical DNA binding motif for PU.1 within their sequence. Notably, AP-1–dependent elements have substantial overlap with risk loci for multiple immune diseases, including allergic disease (54). In addition, the activity of the IL5Ra P1 promoter has been shown to be mediated by AP-1 binding (55), with the consensus AP-1 binding motif present adjacent to the EOS1 enhancer-like element in the human IL5Ra promoter (56). KLF family members have been shown to interact with PU.1, CEBP, and GATA to facilitate myeloid and monocyte differentiation (57–59), in addition to regulating eosinophil function in adipose tissue (60). Furthermore, SMAD3 has known functional roles in the regulation of TGF-β signaling (61), and SMAD3-deficient mice exhibit attenuated pathology in a model of experimental eosinophilic esophagitis (62). The attenuated pathology comprises reduced levels of esophageal remodeling, fibrosis, and angiogenesis, despite no reduction in the eosinophil number (62), suggesting that SMAD3 may be more involved in eosinophil functional responses than eosinophilopoiesis. Our discovery of an eosinophil-specific superenhancer associated with SMAD3 substantiates a key potential role of SMAD3 in regulating gene expression in eosinophils; thus, modifying interactions between SMAD3 and active enhancers may be a novel strategy to reduce eosinophil-associated tissue pathology, including esophageal fibrosis and its associated complications, in experimental eosinophilic esophagitis (62). Collectively, these data highlight the regulation of cell type– and maturation-specific changes in the TF repertoire at eosinophil-specific REs.

In an effort to define the eosinophil-specific epigenome, we compared our findings with publicly available data describing PU.1-bound enhancers in macrophages (21), neutrophils (22), and ImmGrans (23). Differential binding analysis revealed that >25% of the identified PU.1-bound eosinophil active enhancers were unique to eosinophils and thus were described as PU.1-bound, eosinophil-specific enhancers; these enhancers were enriched for pathways involved in activation, degranulation, and granulocyte-mediated immunity. PU.1-bound, common myeloid lineage enhancers were also identified. Our discovery of both a sizeable common myeloid and a unique, eosinophil-specific PU.1-bound enhancer repertoire reinforces the importance of PU.1 in myeloid gene expression and development (10, 11, 20, 47, 50–53) and supports previous studies using PU.1-deficient mice; these studies showed an absence of Epx, Il5ra, and Prg2 (MBP) mRNA expression in fetal liver cells from embryonic lethal PU.1^−/− mice and a total loss of the eosinophil population within the bone marrow and spleen of viable PU.1^−/− mice (49). Further studies using gene-editing techniques, including CRISPR/Cas9, will allow investigation of the consequence of the loss of PU.1-bound, eosinophil-specific enhancers in eosinophil development and functional responses. We are currently in the process of generating an Il5ra enhancer deletion mouse to interrogate the impact of the loss of this PU.1-bound, eosinophil-specific enhancer on eosinophil function and development.

One limitation encountered routinely when using gene set enrichment analysis is the limited set of genes associated with GO terms outlining eosinophil functional responses; currently, there are only 38 genes associated with the eosinophil-related gene lists available in the GO database (63, 64), which may be due to the challenges in isolating RNA from these sensitive, lysis-prone cells. In contrast, there are 199 genes associated with neutrophil-related gene lists. Enrichment for gene pathways involved in neutrophil degranulation, activation, and immunity for genes associated with PU.1-bound eosinophil-specific enhancers emphasizes the gap in the literature describing and outlining the genes involved in eosinophil functional responses. Although these PU.1-bound, eosinophil-specific enhancers were absent in neutrophils, neutrophil-related terms were the most enriched because of limited GO terms outlining relevant degranulation, activation, and functional responses in eosinophils. Although we included these neutrophil-related terms in the analysis, we aim to use the current dataset to contribute to the GO consortium to develop terms relating specifically to eosinophil maturation and functional responses for more accurate distinction of these granulocyte populations within the literature and for the progression of the field of eosinophil-related research.

Finally, the description of the eosinophil-specific superenhancer landscape provides exciting opportunities for investigation into the role of regulator regions associated with driving cell identity (24). Parallel recent findings have identified the importance of Trib1 in the terminal differentiation of myeloid populations; Trib1 regulates both granulocyte precursor lineage commitment and mature eosinophil identity (65). From these findings, we speculate that ongoing studies investigating the tissue- and inflammatory-specific changes in the chromatin landscape will uncover novel key regulators of eosinophil tissue specificity and function, as well as elucidate potential approaches for treating allergic disorders.

Collectively, these novel data describe the epigenetic landscape of eosinophil-specific enhancers in mature, murine, cultured eosinophils, providing an initial step toward an understanding of the eosinophil chromatin landscape and its role in the eosinophil regulome. Future studies investigating the tissue- and inflammatory-specific changes in the chromatin landscape that drive eosinophil maturation, recruitment, tissue specificity, and function will elucidate potential opportunities and approaches for treating allergic disorders, such as manipulation of enhancer-regulated genes for therapeutic benefit.

We thank Shawna Hottinger for editorial assistance, Dr. Mark Rochman for thoughtful contributions, and Dr. Mario Pujato for continued assistance and development of bioinformatic pipelines. RNA-seq was performed in collaboration with the Cincinnati Children’s Hospital Medical Center DNA Sequencing and Genotyping Core (Cincinnati, OH).

M.E.R. is a consultant for Pulm One, Spoon Guru, ClostraBio, Celgene, Astra Zeneca, and Arena Pharmaceuticals and has an equity interest in the first three companies listed and royalties from reslizumab (Teva Pharmaceuticals), PEESSv2 (Mapi Research Trust), and UpToDate. M.E.R. is an inventor of patents owned by Cincinnati Children’s Hospital Medical Center. A.B. is a cofounder of Datirium. P.C.F. has received research funding from Knopp Biosciences. The other authors have no financial conflicts of interest.