IRF8 configures enhancer landscape in postnatal microglia and directs microglia specific transcriptional programs

Microglia are innate immune cells in the brain. Transcription factor IRF8 is highly expressed in microglia. However, its role in postnatal microglia development is unknown. We demonstrate that IRF8 binds stepwise to enhancer regions of postnatal microglia along with Sall1 and PU.1, reaching a maximum after day 14. IRF8 binding correlated with a stepwise increase in chromatin accessibility, which preceded the initiation of microglia-specific transcriptome. Constitutive and postnatal Irf8 deletion led to a loss of microglia identity and gain of disease-associated microglia-like genes. Combined analysis of scRNA-seq and scATAC-seq revealed a correlation between chromatin accessibility and transcriptome at a single-cell level. IRF8 was also required for microglia-specific DNA methylation patterns. Lastly, in the 5xFAD model, constitutive and postnatal Irf8 deletion reduced the interaction of microglia with Aβ plaques and the size of plaques, lessening neuronal loss. Together, IRF8 sets the epigenetic landscape, which is required for postnatal microglia gene expression.


Introduction
Microglia elicit innate immunity in the brain, providing antimicrobial protection 1,2 . Microglia also promote brain development through synaptic pruning, supporting neuronal survival, and vascular formation 3,4,5,6 . Microglia continuously surveil the entire brain and respond to changes in their environment 7,8,9 . Dysregulation of microglia is associated with neurodegenerative diseases such as Alzheimer's disease (AD) 10,11 .
IRF8 is a DNA binding transcription factor that binds to ISRE, and composite ETS/ISRE, the latter by interacting with PU.1 18 . IRF8 is expressed in developing myeloid cells in bone marrow, Ly6C+ macrophages, and CD8+ or plasmacytoid dendritic cells in the periphery 19,20 . IRF8 is essential for innate resistance against various pathogens 21 . Constitutive Irf8 knockout (IRF8KO) mice are susceptible to infection while having excess neutrophils 22 . In the brain, IRF8 is expressed only in microglia, except for a few macrophages 23 . Besides participating in early microglia development, IRF8 takes part in responding to cytokines and pain injury 13,24,25 . However, the mechanism of postnatal microglia maturation is largely unknown 14, 17 .
We show that IRF8 plays a decisive role in this process. CUT&RUN assay revealed that IRF8 binds to intergenic and intronic regions of postnatal microglia gradually reaching a maximum after P14, which led to the formation of microglia specific enhancers. IRF8 occupancy was a prerequisite for setting chromatin accessible sites and DNA methylation patterns unique to microglia.

IRF8 binds to the microglia genome during postnatal development in a stepwise fashion
The IRF8's role in postnatal development and adult microglia remains unclear. To gain insight into this question, we sought to identify genome wide IRF8 binding sites in postnatally developing microglia. To first verify IRF8 expression in postnatal microglia, we examined Irf8-GFP knock-in mice that express endogenous IRF8 fused to GFP 26 . FACS analyses of microglia from P9, P14, and P60 (adult) showed that IRF8 is expressed at similar levels, higher than that in peritoneal macrophages (Fig. S1A, S1B).
CUT&RUN experiments were carried out for microglia isolated from Irf8-GFP mice using anti-GFP antibody. Peritoneal macrophages from the same mice were also examined as a reference. We identified 16,935 IRF8 peaks in adult microglia and 12,062 in peripheral peritoneal macrophages (Fig. 1A).
About half of the IRF8 peaks were in the intergenic regions with >5 kb distal from the genes, while about 20% were at the promoter regions in both cell types ( Fig. 1A). De novo motif analysis revealed enrichment of PU.1/IRF composite motif (TTCC…G/CTTT) and classical ISRE (GAAAG, Fig. 1B). DEseq2-based statistics showed a differential IRF8 binding profile between microglia and peritoneal macrophages (Fig. 1C). Further, GREAT GO annotation indicated microglia specific IRF8 peaks (n=9,679) to be related to cell/leukocyte activation and differentiation ( Fig. 1D; upper panel). In contrast, macrophage specific IRF8 peaks (n=2,251) were related to immune responses and inflammatory responses against pathogens ( Fig. 1D; lower panel). In addition, IRF8 binding was found near microglia identity genes, such as Sall1 and Cx3cr1 11, 27 , but not in peritoneal macrophages (Fig. 1G).
IRF8 binding peaks on P9, P14, and adult microglia are shown in Fig.   1E. The IRF8 binding peaks on P9 microglia were very sparse. The binding increased at P14, although still less than on adult microglia, indicating a stepwise increase in IRF8 binding (Fig. 1E). Furthermore, cumulative distribution frequency graphs plotting the distance between two neighboring IRF8 binding sites (Fig. 1F, right) showed that the distance in the adult is the shortest, while it is longer in P14 microglia and peritoneal macrophages, and Saeki et al. 5 the longest in P9 microglia, which indicates that IRF8 on microglia could form arrays as illustrated in Fig. 1F (left). IGV tracks corroborated the stepwise increase in IRF8 binding, observed over the enhancer region of Sall1 and Cx3cr1 genes (Fig. 1G, left). These regions were not occupied by IRF8 in peritoneal macrophages (Fig. 1G, left). Rather, IRF8 was bound to regulatory regions of Clec4n, Lyz1/2, active in peritoneal macrophages (Fig. 1G, right).
These results reveal that IRF8 binding to the microglial genome is developmentally controlled and reaches a maximum when microglia gain full functional activity. It is of note that despite the stepwise increase in binding activity, the amount of IRF8 was virtually the same during postnatal stages (Fig.   S1B), suggesting a mechanism that directs IRF8's ability to gain access to the microglial genome.

IRF8 directs transcription cascades to activate microglia identity genes
To study how IRF8 binding controls microglia's transcriptional programs, we carried out bulk RNA-seq for microglia isolated from wild-type (WT) and IRF8KO mice. FACS analyses of brain cells revealed two CD45 + populations in the IRF8KO brain, one CD11b high CD45 + corresponding to WT microglia, and another CD11b low CD45 + subset, termed P3, which was absent in WT brain (Fig.   S2A). The former subset expressed CX3CR1, a typical microglial marker, although its levels were lower than WT microglia. P3 subset present only in the IRF8KO brain did not express CX3CR1, instead expressed Ly6C, unraveling the presence of an extra cell population (Fig. S2B). We isolated microglia from WT and IRF8KO brains and processed them for microglial transcriptome analyses. Statistical analysis revealed a large number of genes that were up or downregulated in the absence of IRF8 ( Fig. 2A). A set of microglia identity genes, e.g., P2ry12, Siglech, and Ccr5, encoding surface proteins and a transcription factor, Sall1, were drastically reduced in IRF8KO microglia (Fig. 2B left). In addition, many lysosomal genes represented by CLEAR genes were also reduced in IRF8KO cells, e.g., Ctsb, Ctsc, and Lamp1 (Fig. 2B right) 28 .
Given that IRF8 activates downstream transcription factors to promote monocyte function, we examined levels of transcription factors expressed in postnatal microglia 14, 20, 30 . Among them, Batf3 and Sall1 were markedly downregulated in IRF8KO cells (Fig. 2C, blue square) 14 . BATF3, a basic leucine zipper motif (bZIP) family of a transcription factor, is known to be critical for generating a dendritic cell subset 31 . Because the role of BATF3 in microglial transcriptome remained unelucidated, we next analyzed RNA-Seq in Batf3KO and WT microglia (Fig. 2D). As anticipated, genes downregulated in Batf3KO microglia, such as P2ry13 and Ctsc were also downregulated in IRF8KO cells, indicating that IRF8 and BATF3 control an overlapping set of target genes.
Consistent with this, BATF3-dependent genes showed a significant correlation with IRF8-dependent genes (Fig. 2E). Sall1 is a transcription factor of the Spalt family, expressed selectively in microglia. Buttgereit et al. reported that SALL1 is required for the expression of microglia signature 27 . The correlation biplot also showed a significant overlap between SALL1 dependent genes and IRF8 dependent genes, including microglia identity and Interferon related gene sets ( Fig. 2F). Thus, IRF8 activates Sall1 expression, which in turn activates or represses genes that shape microglia signature. These data indicate that IRF8 creates transcription cascades and activates Sall1 and Batf3 expression, thereby inducing downstream target genes necessary for microglia identity and function (Fig. 2G).

IRF8 configures microglia specific enhancers
Above results established that IRF8 binding represents a critical basis of microglia's transcriptional programs. To further characterize the chromatin environment in which IRF8 binds, we examined histone marks linked to active transcription, particularly those denoted for enhancers, H3K4me1 and H3K27ac 32, 33 . Of 16,935 IRF8 binding sites, about 73% and 66% of IRF8 colocalized with H3K4me1 and H3K27ac, respectively, and 64% of IRF8 binding sites carried both enhancer marks (Fig. 3A). Conversely, H3K4me1 and H3K27ac marks near IRF8 peaks were much lower in IRF8KO cells than in WT cells, indicating that IRF8 helps set an active chromatin environment in Saeki et al. 7 microglia (Fig. 3B). Differential deposition analysis for H3K27ac and H3K4me1 marks revealed that both WT and IRF8KO microglia exhibited cell type specific histone modification patterns (Fig. S3A). Since the H3K27ac is a surrogate marker for super-enhancers and super-enhancers neighbor many genes critical for defining cell type specificity and function 33, 34 , we examined super-enhancers using the ROSE program 33 . As expected, many IRF8-dependent microglial identity genes were neighbored by some of 837 super-enhancers, including Sall1, Siglech, and Cst3 (Fig. 3C). Indeed 735 of 837 (87.8%) super-enhancers were bound by IRF8 (Fig. 3D). While most H3K27ac high regions were populated by IRF8 binding in adult microglia, this was not the case in early postnatal P9 microglia and peritoneal macrophages (Fig. 3E). In addition to super-enhancers, IRF8 was abundantly enriched in smaller enhancers found in specific intergenic and intronic regions overlapping with H3K4me1 ( Fig. S3C-G). IGV tracks for the Sall1 and Batf3 genes show that clusters of IRF8 binding coincide with those of H3K27ac and H3K4me1 (Fig. 3F). These results indicate that IRF8 binding guides the formation of enhancers unique to microglia during postnatal development.

IRF8 confers chromatin accessibility on postnatal microglia and sets DNA methylation patterns
To investigate genomic sites accessible for transcription and whether IRF8 influences the process, we performed an assay for transposaseaccessible chromatin (ATAC) sequencing 35 . Differential accessibility analysis revealed that 10,514 peaks were lost in IRF8KO microglia while gaining 5,926 new peaks (Fig. 4A). As anticipated, gained and lost peaks were mostly localized to intergenic regions and introns (Fig. 4B). De novo motif analysis showed that lost peaks were enriched with the PU.1/IRF composite motifs (Fig.   4C, left). Whereas, gained peaks were enriched with motifs for CEBP, AP-1, and RUNX binding sites (Fig. 4C, right). Further, ATAC peaks lost in IRF8KO cells were enriched with WT-specific H3K4me1 and H3K27ac marks (Fig. S3A, 4D).
On the other hand, IRF8KO-specific H3K4me1 marks were found in gained the Sall1 and Batf3 genes in WT cells but few in IRF8KO cells (Fig. 4I). Further, in WT cells, H3K27ac and H3K4me1 marks were abundantly present at these ATAC-seq sites but reduced in IRF8KO cells (Fig. 4I). PU.1 binding sites also revealed a correlation with ATAC-seq peaks (Fig. 4I). However, Apoe, 4I, bottom two lanes). Sall1 and Batf3 were enriched with demethylated sites in WT cells, but not IRF8KO cells. Conversely, Apoe exhibited many methylated sites in WT cells, but demethylated sites in IRF8KO cells. The above data provide the first evidence that IRF8 suppresses DNA methylation in gene specific manner in microglia. Together our results indicate that the regulation of DNA methylation is an integral mechanism for setting epigenomic landscape in microglia, which IRF8 dictates.

Continuous IRF8 expression is required for the transcriptome in adult microglia
To assess whether IRF8 is required for postnatal microglia, we constructed Irf8 flox/flox Cx3cr1cre ERT2 mice to allow for conditional Irf8 deletion.
The mice also carried a Rosa26-loxP-stop-loxP-EYFP reporter, with which we could identify and isolate Cre-expressing microglia 40 (Fig. 5A). To our surprise, the efficiency of Irf8 deletion was very poor (28% ± 18%SD, Fig. S5A) when Tamoxifen was administered on 4-week-old mice five consecutive times. To mitigate this problem, we resorted to deleting Irf8 in younger mice, between P12 and P16, and found a considerable improvement in Irf8 deletion efficiency (62% ± 5%SD, Fig. S5A). The basis of poor deletion efficiency is currently unclear, but it may be due to a shift in Cre accessibility.
Similarly, the expression of transcription factors, Sall1 and Batf3, IRF8's downstream targets, and CLEAR Lysosome genes were much lower in IRF8cKO cells (Fig. 5C, 5D). Moreover, aberrant induction of DAM/NDG genes and Interferon-related genes was also observed in IRF8cKO microglia, mimicking changes in IRF8KO microglia (Fig. 5B, S7B). Together, these results show that continuous IRF8 expression is required for maintaining the transcriptional programs in postnatal microglia.

Single-cell RNA-seq analysis of microglia reveals a progenitor-like population
Next, we performed single-cell RNA-Seq (scRNA-Seq) analysis to assess transcriptome diversity. After filtration, 3,504 WT and 5,303 IRF8KO cells expressing a total of 15,376 genes were selected. Leiden clustering denoted three WT clusters, seven IRF8KO clusters, and a "Common" cluster consisting of WT and IRF8KO cells ( Fig. 6A and B). The SigEMD pipeline identified 6,149 DEGs between WT and IRF8KO cells, which covered 62% (3016/4848) of DEGs found in bulk RNA-seq analysis, lending credence to our data (Fig. S6B). Microglia identity gene expression was higher in all WT clusters than in all IRF8KO clusters (Fig. 6C, left). Population heterogeneity was evident within IRF8KO cells, in that interferon-related genes were exclusively found in  (Fig. 6C, S6A). Thus, this cluster likely represents the P3 population in IRF8KO brain (Fig. S2B). The Common cluster exhibited unusual features, in that most microglial genes were absent in this population (Fig. 6D).
This raised the possibility that Common cluster represents a population prior to microglia maturation. To delineate a relationship between Common cluster and other clusters, we applied the PAGA program, which provided a trajectory network (Fig. 6E) 41 . This network illustrated a tight connection with the Common cluster to all WT and IRF8KO clusters except IRF8KO Cluster 6. Another unusual feature of Common cluster was the expression of cell cycle genes such as Cenpa, Ccnb2, and Mki67, a marker for cell proliferation (Fig. 6F). If Common cluster represents a developmentally immature population, it would predict that IRF8 dependent transcription cascade would be present only when it differentiates into WT clusters. To test this, we examined genes downstream of IRF8, BATF3, or SALL1 in WT and IRF8KO clusters relative to Common cluster. Correlation analyses revealed a significant correlation only in WT clusters, not any IRF8KO clusters, as predicted. Together, these results support the idea that Common cluster represents a predecessor of microglia, possessing the capability to develop into functioning microglia.

Irf8 deletion ameliorates AD pathology in the mouse model
IRF8 has been shown to promote neuroinflammation and implicated in AD pathogenesis in mice 21,25,42,43 . To further address the role of IRF8 in AD, we first compared microglia transcriptome profiles of WT and 5xFAD mice at 9 months of age 44 . We found 699 genes upregulated in 5xFAD microglia over WT cells, which we regarded as the 5xFAD associated gene set. This gene set contained many DAM/NDG genes described earlier, such as Apoe, Tyrobp, and Cst7 ( Fig. S7A) 10,11 . Interestingly, many genes in this gene set, mostly DAM/NDG genes in an early stage of human and mouse AD microglia 29, 43, 45 , were also upregulated in IRF8KO microglia without the 5xFAD background (Fig.   7A). Further, microglia from 5xFAD/IRF8KO mice showed essentially the same expression pattern as IRF8KO microglia (Fig. 7A). These DAM/NDG genes were upregulated in IRF8cKO microglia even from much younger mice (three months old, Fig. S7B, S7C). Thus, Irf8 deficiency leads to the induction of DAM/NDG genes. In addition, scRNA-seq analysis found that DAM/NDG gene expression was most prominent in Cluster 1 of IRF8KO cells (Fig. 7B).
Transcription factor Runx3 was also upregulated in this cluster (Fig. 7B).
To further examine the role of IRF8 in AD, histological analyses were performed with a year-old 5xFAD/WT and 5xFAD/IRF8KO brains using Thioflavin S (ThioS) and 6E10 antibody 46 . The former detects polymerized Aβ sheets, while the latter stains both polymerized and unpolymerized Aβ.
Confocal images revealed that the number of ThioS positive aggregates was similar between 5xFAD/WT and 5xFAD/IRF8KO brains. However, the amyloid load and plaque size were smaller in the 5xFAD/IRF8KO than in the 5xFAD/WT brain (Fig. 7C, quantification in 7D). Moreover, in the 5xFAD/IRF8KO brain, many of the Aβ signals detected by 6E10 (green) did not directly colocalize with Iba1, a marker for microglia (Fig. S7D, red), whereas extensive colocalization was evident in WT cells. These results suggest that 5xFAD/IRF8KO microglia are defective in the recognition and perhaps internalization of Aβ in line with the report 42 . We next sought to assess the effect of Irf8 deletion on overall neuronal health in the 5xFAD brain. We thus stained cortical regions of a year-old 5xFAD/WT and 5xFAD/IRF8KO brains with NeuN antibody to detect intact neurons 54 . NeuN immunoreactivity is lower than those from 5xFAD/IRF8KO brains than the 5xFAD/WT brains (Fig. 7E, quantification in 7F). These results indicate that deletion of IRF8 lessens neuronal damage in the 5xFAD brain.
To further assess the role of IRF8 in AD, we checked Axl expression, as AXL is reported to regulate microglia's ability to phagocytose Aβ 47 . Axl levels were lower in IRF8KO microglia than in WT microglia with and without 5xFAD background (Fig. 7A, S7F). These results suggest that IRF8 would affect AD pathology in a complex manner, inhibiting Aβ spread on one hand, and facilitating phagocytic clearance of Aβ on the other.

Discussion
This study aimed at elucidating the mechanism of postnatal microglia development and the role IRF8 plays in the process, the subject that has long remained elusive. It was critical first to determine genome-wide IRF8 binding sites and epigenome structures in the developing microglia. We show that IRF8 binding occurs in a characteristic stepwise fashion. IRF8 gained binding activity gradually, from P9 through P14. Binding was sparse on P9, increased on P14, and reached the maximum in adults. Microglia from P9 to P14 are still immature, gaining full identity and functionality later 17 . Thus, full IRF8 binding coincided with the acquisition of microglia identity and functionality. It was striking that a stepwise increase in IRF8 binding was not due to corresponding increase in IRF8 expression, since IRF8 levels were constant during this period and through adulthood, suggesting an epigenome-based mechanism, yet to be deciphered. In adult microglia, IRF8 peaks were often in dense repeats and within or near H3K27ac and H3K4me1 peaks, indicating that IRF8 binds mostly in the enhancers, both large and small. Furthermore, IRF8 binding was obligatory for installing chromatin accessibility sites, as revealed by ATAC-seq analysis. Many ATAC-seq peaks were in the IRF8 dependent enhancers, as exemplified by those for Sall1 and Batf3, IRF8 dependent transcription factors endowing microglia's transcriptome profiles and phenotypes. These ATAC-seq peaks were lost in IRF8KO microglia, which in turn, gained alternate ATAC-seq peaks neighboring genes upregulated in IRF8KO microglia, including DAM/NDG genes, such as Apoe. We further show, by genome-wide bisulfite sequencing, that IRF8 sets DNA methylation patterns specific for microglia. The methylation patterns correlated well with the ATAC-seq peak patterns. Thus, Saeki et al.

13
Sall1 and Batf3 enhancers had unmethylated DNA. Conversely, in IRF8KO microglia, these regions of DNA were methylated, and the Apoe region was unmethylated. Together, the formation of enhances, chromatin opening, and DNA methylation patterns unique to microglia depend on developmentally regulated IRF8 binding. Thus, we suggest a model in which IRF8 binding is the primary initiating event that confers microglia specific epigenomic landscape (Fig. 8). Our data on IRF8cKO microglia further support this model, in that deletion of Irf8 from adult microglia led to the loss of microglia identity and function closely resembling the phenotype of microglia from constitutive IRF8KO mice, demonstrating that continuous IRF8 occupancy is necessary to maintain microglia specific chromatin and DNA methylation landscape.
Aside from this, it is noteworthy that Sall1 and Batf3, transcription factors that provide actual functionalities to microglia, are downstream of IRF8. Clearly, IRF8 forges transcription factor cascades to achieve microglia specific transcriptome programs. IRF8 promotes the differentiation of monocytes and dendritic cell subsets by activating downstream transcription factors 18,19 . These Irf8 deletion led to a strong reduction of the microglia identity gene set but unleashed marked expression of DAM/NDG and Interferon-related gene sets. Both gene sets are associated with AD and other neurodegenerative diseases 10,11 . Nevertheless, AD pathology in 5xFAD brains without IRF8 was less severe than that with IRF8, as evidenced by reduced large Aβ plaques and enhanced NeuN immunoreactivity. Microglia Saeki et al. 14 are known to constantly survey the brain space, recognize Aβ and distribute them to other brain regions 42 . IRF8KO microglia are deficient in these activities, since they lack the ability to scan the brain to recognize Aβ, as verified by the absence of colocalization of Aβ with microglia. d 'Errico et al. further reported that Aβ plaques did not spread in the IRF8KO brain relative to WT counterparts 42 . Pertinent to this, it has been shown that depletion of microglia by CSF1R inhibitors reduces plaque formation and neuronal loss 49 . Notwithstanding these results, IRF8 likely exerts a complex role, positive and negative, in AD progression, given that Axl and CLEAR/lysosomal genes, crucial for controlling AD, require IRF8 for full expression.
In summary, this study demonstrates the essential requirement of IRF8 for the adult development of microglia and its complex role in modulating AD progression.

Microglia preparation
Mice were euthanized by asphyxiating with CO2 and transcardially perfused with 10 mL of ice-cold PBS. The brain was harvested, chopped with a single-edge razor blade, and ground in ice-cold Hanks-balanced salt solution using a 7ml Dounce homogenizer, then filtered through a 70µm cell strainer.
The cells were resuspended in 30% isotonic Percoll and spun down 800xg for

Single-cell RNA-seq
Ten thousand cells were sorted out from 3 mice pools and immediately emulsified using a Chromium kit and 10x Genomics as the manufacturer's protocol. The libraries were made using SMART-Seq v4 Ultra Low Input RNA Kit followed by Nextera XT library prep kit and sequenced sequenced 26bp for R1 and 98bp for R2 on HiSeq2000. The fastq files were first processed with Cellranger v3.1 pipeline to make gene expression matrices, and the following analyses were performed with Scanpy 52 . Cells with less than 500 or greater than 6,000 detected genes were filtered from analysis. Genes detected in less than three cells were also removed from the data. The feature counts were normalized to 10,000 reads per cell and logarithmized. The data were regressed out based on the highly variable genes marked with  with the previous dataset (GSM1533906) 16 . The peak co-occurrence, genomic position, and motif analysis were performed using HOMER 32 with -d given option. To determine highly concordant peaks across biological replicates, we employed "the Majority rule," which picks up only peaks found within more than 50% of replicates 55 . The differential peaks were obtained using Diffbind-DESeq2 56 with FDR<0.05 cut-off. Then all bam files were merged and used for IGV visualization and depicting average plots and heatmaps with Deeptools 57 .
H3K27ac high large regions were identified as "super-enhancer regions" of the ROSE's readouts with default parameters 33 .

Bulk ATAC-Seq
ATAC-Seq libraries were constructed as described before 35 . Fifty thousand isolated microglia were lysed in lysis buffer (10mM Tris-HCl pH7.4, 10mM NaCl, 3mM MgCl2, 0.1% NP-40). Following lysis, nuclei were tagmented using the Nextera DNA library prep kit (Illumina) for 30 minutes. Then DNA was purified and amplified for 12 cycles. The libraries were sequenced with a 50bp paired end read. After removing the adapter and low-quality sequences, the data were aligned to the mm10 genome with Bowtie2. Peaks were called on aligned sequences using MACS2 with -p 0.01 --nomodel --shift -37 --extsize 73 -B --SPMR --keep-dup all options. Peaks were then evaluated with less than 0.01 irreversible differential rates. Differentially accessible regions were identified using Diffbind-EdgeR with FDR<0.05 56 . Deeptools and HOMER were employed for further analyses, including motif enrichment analysis with the indicated options.

Whole-genome bisulfite sequence
Genomic DNA was extracted from 40,000 male microglial cells using Quick-DNA Microprep plus kit (Zymo Research) as the manufacturer's protocol.
Then 50ng DNA was bisulfite-converted and processed to the library using Pico Methyl-Seq Library Prep Kit (Zymo Research), with 8-cycle PCR amplification.
The libraries were sequenced with 75bp paired end read on HiSeq3000 in the NHLBI sequencing core (9.2x coverage on average). Based on the Bismark Mbias plot, the adapters, low-quality sequence, and an additional 10bp of 5'-end and 2bp of 3'-end were removed from the sequencing data using TrimGalore.
The data were first aligned using Bismark 38 with -q --score-min L,0,-0.4 --ignorequals --no-mixed --no-discordant --dovetail --maxins options, and then unmapped sequences were aligned additionally with single-end mode. After deduplication and methylation calls, both files were merged into a single coverage file. Differential analyses were performed using R/DSS package 39 using FDR<0.05 for the differentially methylated loci (DMLs) and p<0.001 for the regions (DMRs).
Saeki et al.

20
Data availability: All high-throughput sequence datasets generated in this paper are available in GSE231406. Some gene sets were obtained from publications as denoted individually.
Acknowledgments: We thank Steven Coon, NICHD DIR MGC, for his guidance for scRNAseq, Yuesheng Li in NHLBI sequencing core for whole-genome bisulfite sequencing, and Fayaz Seifuddin for data analysis. We also thank Avinash Bhandola for providing Rosa26-EYFP mice, Lisa William-Simons for cell sorting, and the NICHD animal facility staff for caring for our mice and for other technical assistance. We gratefully acknowledge colleagues in the NIH and elsewhere for valuable advice and discussions. This work was supported by the NICHD Intramural programs ZIA HD008015-13.
Disclosure: There is no conflict of interest.    Representative genes are shown with a label.  for the WT or IRF8KO specific peak identification.