B1a and B2 cells are characterized by distinct CpG modification states at DNMT3A-maintained enhancers

The B1 and B2 lineages of B cells contribute to protection from pathogens in distinct ways. The role of the DNA CpG methylome in specifying these two B-cell fates is still unclear. Here we profile the CpG modifications and transcriptomes of peritoneal B1a and follicular B2 cells, as well as their respective proB cell precursors in the fetal liver and adult bone marrow from wild-type and CD19-Cre Dnmt3a floxed mice lacking DNMT3A in the B lineage. We show that an underlying foundational CpG methylome is stably established during B lineage commitment and is overlaid with a DNMT3A-maintained dynamic methylome that is sculpted in distinct ways in B1a and B2 cells. This dynamic DNMT3A-maintained methylome is composed of novel enhancers that are closely linked to lineage-specific genes. While DNMT3A maintains the methylation state of these enhancers in both B1a and B2 cells, the dynamic methylome undergoes a prominent programmed demethylation event during B1a but not B2 cell development. We propose that the methylation pattern of DNMT3A-maintained enhancers is determined by the coincident recruitment of DNMT3A and TET enzymes, which regulate the developmental expression of B1a and B2 lineage-specific genes.

B cells are divided into two major sub-lineages. B2 cells, which constitute the majority of all B cells, are short-lived recirculating cells that are metabolically quiescent and can be activated in secondary lymphoid organs where they collaborate with helper T cells to generate specific antibodies. B1a cells, found predominantly in the peritoneal and pleural cavities, are longlived self-renewing cells that produce natural antibodies as part of innate-like immunity 1,2 . Given the crucial role of DNA methylation in controlling gene expression, lineage specification, and the maintenance of cellular states, we explored the B1 and B2 sublineages through the lens of whole-genome cytosine modification dynamics 3,4 . Intergenic CpG island methylation has been linked to lineage specification in the immune system 5 , and it has been suggested that DNA methylation might protect specific cell types from the aberrant activation of transcription factors in a lineagespecific manner 6 . Genome-wide studies of cytosine methylation have been undertaken in human B cells, and it has been shown that early human B-cell development is characterized by the loss of CpG methylation at enhancer sites; additionally, the differentiation of B cells into memory and plasma cells is characterized by the loss of CpG methylation in regions of heterochromatin and a gain of methylation in polycomb-repressed regions 7 . However, the identification of human B1 cells is still controversial and prior studies on human B-cell methylomes did not address the B1 vs. B2 lineage decision.
The enzymes that add methyl marks to cytosines include DNA methyltransferase 1 (DNMT1) and the DNA methyltransferases 3a (DNMT3A) and 3b (DNMT3B). DNMT1 is termed a maintenance methyltransferase because it recognizes hemi-methylated sites in newly replicated DNA and restores symmetric CpG methylation. DNMT3A and DNMT3B are de novo methyltransferases that are active during early embryogenesis and modify unmethylated CpGs symmetrically 8 . Although 5methylcytosine (5mC) is a stable modification, it can be converted to 5-hydroxymethylcytosine (5hmC) by methylcytosine dioxygenases of the Ten Eleven Translocation Family (TET1, TET2, and TET3) [9][10][11][12] . 5hmC itself is also a substrate for TET enzymes and can be progressively oxidized to 5-formyl (5fC) and 5-carboxy (5caC) cytosine, which are recognized and excised by thymine-DNA glycosylase, such that an unmethylated cytosine is restored by the base-excision repair pathway. Unlike 5mC, the 5hmC mark is not maintained by DNMT1 during cell division and appears to undergo DNA replication dependent loss, thereby providing an additional mechanism for DNA demethylation 13 . Thus, the TET enzymes play a critical role in active demethylation of mCpGs.
Hypomethylated CpGs mark gene regulatory regions such as promoters and some enhancers. The 5hmC mark may serve as an indicator of TET activity, but given that this mark is especially enriched in enhancer elements, it has been speculated that it may play an active, functional role beyond being just an intermediate in demethylation [14][15][16] . Several readers that specifically recognize methylated or unmethylated CpGs have been identified, but a specific reader for 5hmC has not yet been identified. Unlike DNMT1, which is associated with the DNA replication machinery, DNMT3A and the TET enzymes are recruited to specific genomic sites such as gene regulatory regions through their association with transcription factors. Engineered loss of Dnmt3a in mouse hematopoietic stem cells has been reported to result in malignant transformation in many hematopoietic lineages including the B lineage, resembling human B-CLL [17][18][19][20][21] . Leukemic transformation in the Eµ-Tcl1-transgenic model of CLL is also thought to be driven by the ability of TCL1 to inhibit DNMT3A 22 . Once a CpG methylation pattern is established across the genome by DNMT3A and TET enzymes during Blineage commitment, it is subsequently maintained across cell divisions by DNMT1. Mice lacking Tet2 and Tet3 in the B lineage retain normal numbers of follicular B2 cells but experience a partial block in early B-cell development and completely lack B1a and marginal zone B cells 23 , suggesting a continued role for DNA demethylation in establishing the B1a cell program after B lineage commitment. In this study, we elucidate the role of DNMT3A during and after B lineage commitment by studying Dnmt3afloxed mice in the CD19-Cre +/background at 4 weeks of age, well before the establishment of leukemia and at a time when the B1a cells are still polyclonal 20,21,24 , in order to examine the role of DNMT3A in the absence of the confounding effects of clonal somatic aberrations and widespread changes in CpG methylation that are known to accompany malignant transformation into CLL.
We profile cytosine modifications in B1a and B2 cells and their precursors and show that both B1a and B2 cells exhibit a dynamic layer of lineage-specific and DNMT3A-maintained CpG methylation that is superimposed on a shared foundational methylome. We find that continued action of DNMT3A is required to counteract TET activity in this dynamic methylome layer, and the sites whose CpG methylation is maintained by DNMT3A after B lineage commitment are highly enriched for enhancers, which we have termed DNMT3A-maintained enhancers (DMEs). DME methylation is linked to the regulated expression of lineagespecific genes and reveals a novel role for DNMT3A in both marking and maintaining DNA methylation at specific enhancers that modulate the fate of developing lymphocytes. Our results also provide an epigenetic framework for understanding the reported BCR-dependent plasticity in B-cell sub-lineages 25 .

Results
B1a cell development is characterized by a pronounced loss of CpG modification. We surveyed the genome-wide CpG modification profiles in murine B1a and B2 cell lineages, by performing whole-genome bisulfite sequencing (WGBS) of proB2 cells (Hardy fractions B and C from bone marrow), splenic follicular B2 cells, proB1 cells (Hardy fractions B and C from fetal liver) and peritoneal B1a cells; these populations will hereon be referred to as proB2, B2, proB1 and B1a cells, respectively. We identified about 40,000-50,000 hypomethylated regions (HMRs) that were dispersed across 1.7% to 1.9% of the genome, largely in intergenic and intronic regions in all four cell types 26 . Interestingly, the HMRs in wild-type B1a cells were found to be significantly longer than those in the B2 cells (median lengths of 830 bp and 668 bp, respectively; K-S test p < 1e-10), indicating that the B1a lineage may favor CpG demethylation at HMRs across the genome (Fig. 1a). Large HMRs over 3.5 kb have been previously referred to as canyons in studies on murine hematopoietic stem cells 27 . Using a similar cutoff of 3.5 kb, we identified methylation canyons in all four B-cell populations. Mature B1a cells exhibited 30-35% more canyons than all other B-cell stages examined (Fig. 1b). This developmental increase in the number of canyons in B1 cells was reminiscent of the canyon expansion observed in Dnmt3a-deficient murine HSCs 27 and suggested that changes in DNMT3A activity may contribute to the programmed demethylation observed during B1a cell differentiation.
We next identified differentially methylated regions (DMRs) that arise during the development of the B1a and B2 cell lineages. DMRs are regions with a high probability of differential methylation where an HMR exists in one methylome but not at the same location in the other. DMRs were particularly abundant between the proB1 and B1a stages. They encompassed~6 million bp across 5,392 genomic intervals in the B1a lineage (0.42% of the genome), but only about 0.5 million bp across 639 intervals in the B2 lineage (0.04% of the genome), suggesting that the regulation of CpG methylation is a far more prominent feature of B1a cell development, compared to B2 cell development (Fig. 1c-e). The majority of the B1a DMRs (~75%) exhibited a loss of CpG modification during B1a cell development (Fig. 1c, d). The median size of the observed B1a DMRs was 857 bp, and~12% of the hypomethylated B1a DMRs fell within methylation canyons in B1a cells (Fig. 1f) Programmed demethylation in B1a cells occurs at enhancers that are demethylated and re-methylated in B2 cells. The majority of B1a and B2 DMRs were present at a distance from the transcription start site (TSS) in intergenic or intragenic spaces, compatible with a possible enhancer function (Fig. 2a, b). Indeed, many prior studies of CpG methylomes in adult somatic tissues have found that variably methylated regions are strongly enriched at enhancer sites. Analysis of publicly available ChIP-Seq data for H3K4me1 and H3K4me3 modifications in mature B2 and B1a cells 28 showed that the DMRs that underwent developmental demethylation in both B1 and B2 cells co-localized with H3K4me1 peaks in mature B1a and B2 cells. This suggested that these DMRs might represent developmentally demethylated enhancers that modulate lineage-specific gene expression (Fig. 2c). To assess the significance of these enhancers in the context of hematopoietic development, we also examined H3K4me1 and H3K4me3 marks in HSCs using published ChIP-Seq data. The DMRs that were developmentally hypermethylated in mature B2 or B1a cells exhibited little H3K4me1 signal in B cells, but did bear an H3K4me1 signature in HSCs, suggesting  that these hypermethylated DMRs may represent enhancers that were developmentally silenced by CpG methylation (Fig. 2c). The B1a DMRs exhibited low levels of H3K4me1 in HSCs, suggesting that they may represent primed enhancers in HSCs that subsequently become fully active in the B lineage. In B2 cells, only a small subset of developmentally demethylated DMRs that were present in the close proximity of annotated transcription start sites (TSS) exhibited an H3K4me3 signal, consistent with a promoter function (Fig. 2a). However, in B1a cells, most DMRs that were developmentally demethylated acquired some H3K4me3 trimethylation marks (Fig. 2c). Although H3K4me3 is generally associated with promoters, it can be found in enhancers in some contexts 29,30  and prevent CpG methylated enhancers from acquiring ectopic H3K4me3 marks 31 . This CpG methylation barrier to H3K4me3 acquisition appears to break down at many enhancer sites in B1a cells, and the developmental acquisition of ectopic H3K4me3 marks at demethylated enhancers is unique to B1a cells.
Interestingly the H3K4me3 peaks were wider in B1a cells than B2 cells across the genome, consistent with a global increase in HMR lengths in B1a cells (Fig. 2d). WGBS revealed that 90% of the CpGs in DMRs that were developmentally demethylated in B1a cells were modified in B2 cells (Fig. 2c). As WGBS does not distinguish between cytosine methylation and hydroxymethylation, we also analyzed 5hydroxymethylation (5hmC) levels in mature B1a and B2 cells using TET-assisted bisulfite sequencing (TAB-Seq) 32 (Fig. 2c, far right). Given that the bulk of the bisulfite protection arose from 5mCpG rather than 5hmCpG, we focused on 5mCpG in our analyses. The majority of CpGs in B1 DMRs were modified in B2 cells. About 10 percent of these CpG were hydroxymethylated, indicating that TET enzymes are active at these sites even in B2 cells (Fig. 2e). In contrast, our analyses showed no appreciable level of 5hmC at the B1 DMRs in B1 cells. As previously described in HSCs, 5hmCpGs, a mark of TET activity, were found to be enriched at HMR boundaries in both B1a and B2 cells (Fig. 2f). TET activity at HMR boundaries plays a critical role in maintaining the HMRs by opposing the action of DNMT3A 27,33 .
B1 lineage-specific gene expression is linked to enhancer demethylation. Examination of gene expression data in B-cell subsets from the ImmGen Consortium indicated that the genes in the proximity of B1a cell-specific DMRs that lose methylation during B1a cell development are upregulated during differentiation from proB1 to mature B1a cells (either fetal liver Fraction E or peritoneal B1a cells) (Fig. 3a) 34 . These genes are enriched for pathways involved in B-cell receptor signaling and activation (Fig. 3b). In contrast, B1a cell-specific DMRs that gain methylation during development were not associated with statistically significant changes in the expression of neighboring genes (Fig. 3a). Unlike the hypomethylated DMR-proximal genes, the genes in the vicinity of hypermethylated or decommissioned enhancers are not linked to typical lymphocyte functions (Supplementary Table 1). B2 DMRs were less numerous and exhibited little correlation with gene expression changes during development ( Supplementary Fig. 1). Although some genes exhibited multiple DMRs in their vicinity, there was no correlation between the number of DMRs and gene expression ( Supplementary  Fig. 2). Genes whose expression was increased during B1a cell development exhibited a significant overlap with genes that were closest to DMRs that lost CpG methylation during B1a cell development (p < 0.01; Fisher's exact test) (Fig. 3c). Of these genes, 40% encode B1a lineage-specific genes, that include transcription factors such as Arid3a, Bhlhe41 and Zbtb32, and signaling regulators such as Siglecg and Rasgrp1, that are known to be required for B1a cell development and function (Fig. 3d) [35][36][37][38] .
Deficiency of Dnmt3a in the murine B lineage results in the selective expansion of B1a cells. In order to evaluate the possible contribution of DNMT3A to the B1 vs. B2 cell fate decision, we used Cd19-Cre to conditionally delete a floxed Dnmt3a allele (Dnmt3a 2lox ) in the B lineage, specifically from the proB1 and proB2 stages onwards 39,40 . The direct or indirect loss of DNMT3A in the B lineage has been shown to result in the leukemic transformation of B1a cells [20][21][22] . Similarly, Cd19-Cre +/-Dnmt3a 2lox/2lox (Dnmt3a -/-) mice also exhibited a marked expansion of B1a B cells that progressed to monoclonal leukemia by 9 months of age with a 100% penetrance (Fig. 4a). The increase in B1a cells was initially detectable in the peritoneal cavity as early as 4 weeks of age ( Fig. 4a) and became evident in the spleen and bone marrow after 10 weeks of age. The peritoneal B1b B-cell population was progressively displaced by the expanding B1a population with age. We observed no increase in the proportion of CD19 + B1a lineage progenitors or proB1 cells in the fetal liver 41 , suggesting that the expansion of B1a cells was driven by increased self-renewal or proliferation in the mature B1a cell compartment ( Supplementary Fig. 3a). There was a slight decrease in the number of follicular B cells in the spleen and no change in the number of developing B2 cells in the bone marrows of Dnmt3a -/mice ( Fig. 4b and Supplementary Fig. 3b). The number of marginal zone (MZ) B cells was also reduced at 8 weeks of age ( Fig. 4b and Supplementary Fig. 3c). Changes in the MZ B-cell compartment could not be assessed beyond 12 weeks of age, as the splenic marginal zone niche was obliterated by infiltrating B1a cells and MZ B cells were no longer detectable.
At 4 weeks of age, B1a cells lacking DNMT3A remained largely polyclonal, similar to wild-type B1a cells ( Fig. 4c and Supplementary Fig. 4). Staining of Igκ and Igλ light chains showed that the expanded B1a cell pool in the peritoneal cavity underwent a gradual restriction in repertoire from a polyclonal to monoclonal state. Of the eight mice analyzed beyond 30 weeks of age, seven showed only Igκ staining while one showed only Igλ staining. Dnmt3a -/-B1a cells exhibited higher Ki-67 staining, consistent with increased self-renewal in this compartment (Fig. 4d). IgH rearrangements from the clonal leukemias were amplified and sequenced, and like unmutated human CLL, were found to be free of somatic hypermutation. To avoid potential confounding effects arising from oncogenic somatic mutations that may have been selected in a leukemic Dnmt3a -/clone at a later time point, we used peritoneal Dnmt3a -/-B1a cells from 4-week old mice, when the cells were still in a polyclonal state, for all subsequent analyses.
DNMT3A loss reveals a B lineage-specific foundational methylome shared in B1a and B2 cells. We identified a total of 9,401 and 9,765 genomic intervals in the B2 and B1a cell lineages, respectively, that exhibited differential CpG modification across proB, mature B (wild-type) and Dnmt3a-deficient mature B cells, and henceforth refer to them as B2 and B1a DNMT3A-dependent DMRs (DDMRs), respectively (Fig. 5a). Remarkably, the Fig. 2 Extensive loss of enhancer methylation during B1a cell development. a Distribution of distance to nearest TSS for B1a (red) and B2 (blue) DMRs. B2 DMRs that overlap with H3K4me3 ChIP-Seq peaks in B2 cells (black) are separately depicted. b DMR genomic features. Asterisks depict genomic features that are significantly enriched with respect to HMRs. c The heatmap depicting the CpG modification ratio in B1a and B2 DMRs is reproduced from Fig. 1c to provide a frame of reference on the far left. H3K4me1 and H3K4me3 histone marks are correspondingly plotted + /-5 kb upstream or downstream of the B1a and B2 DMRs in HSCs 27, B1a and B2 cells 28. An aggregate signal derived from the 3 published replicates per condition is plotted. 5hmCpG levels + /-2 kb of the DMRs are shown on the far right as measured by TAB-Seq in B1a and B2 cells (one replicate per condition, using cells pooled from three mice each). d Distribution of H3K4me3 peak lengths in B1a (red) and B2 (blue) cells. e Average profile of 5hmCpG modification (% of total CpGs) in the B1a DMRs in B1a (red) and B2 (blue) cells. The DMRs are scaled to the same length and 5 kb flanking regions are depicted. f Average profile of 5hmCpG modification (% of total CpGs) in the HMRs in B1 and B2 cells and 5 kb flanking regions. Source data are provided as a Source Data file.
NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22458-9 ARTICLE NATURE COMMUNICATIONS | (2021) 12:2208 | https://doi.org/10.1038/s41467-021-22458-9 | www.nature.com/naturecommunications lineage-specific methylation patterns in mature B2 and B1a cells were almost completely erased in the absence of DNMT3A, unmasking a foundational methylome that was virtually identical in B1a and B2 cells (Fig. 5b). Indeed, there were only 34 differentially methylated intervals between Dnmt3a -/-B1a and B2 cells. Interestingly, the CpG methylation profile of wild-type B1a cells resembled Dnmt3a -/-B1a and B2 cells, suggesting that DNMT3A activity may be physiologically downregulated during B1a cell development (Fig. 5b). We found that the number of canyons increased in Dnmt3a-deficient B cells compared to their wild-type counterparts (Fig. 5c) and appeared to arise from the expansion and amalgamation of smaller HMRs present at the respective sites in wild-type B cells. The canyons in Dnmt3a-deficient B2 and B1a cells exhibited a very high degree of overlap (98%) and were strongly enriched for genes encoding transcription factors, such that for nearly one-third of these canyons, the closest gene encoded a transcription factor.
Comparison with publicly available WGBS data from murine wild-type and Dnmt3a-deficient bone marrow hematopoietic stem cells (HSCs) suggested that genomic loci exhibiting Dnmt3a-dependent CpG methylation showed very little overlap between HSCs and B cells (Fig. 5d). 27 . Of note, given the use of Cd19-Cre for deleting the Dnmt3a-floxed allele, the loss of DNMT3A occurs after B lineage commitment in the mice studied. Interestingly, the B lineage DDMRs were predominantly methylated in HSCs even in the absence of Dnmt3a (Fig. 5e). This suggests that the foundational and dynamic methylomes are recast following differentiation into the B lineage whereupon the CpG methylation at the B lineage DDMRs becomes dependent upon DNMT3A, perhaps to counteract the activity of TET2 as suggested by the enrichment of 5hmC at these sites ( Fig. 5f, g, h).
The DDMRs were stratified into 12 groups based on their patterns of CpG modification; six patterns accounted for >97% of all DDMRs in B2 and B1a cells ( Supplementary Fig. 5a, f). As previously observed, the majority (~75%) of the B1a DDMRs exhibit a loss of CpG modification during the differentiation of proB1 cells into B1a cells (Fig. 5f). Notably, the DDMRs that undergo partial demethylation from the proB1 to the mature B1a stage exhibited almost complete demethylation in Dnmt3a -/-B cells (Pattern 7 in Fig. 5f), suggesting that even the reduced level of CpG methylation observed in the DDMRs in wild-type B1a cells is maintained by DNMT3A. In B2 cells, which exhibited >65% CpG modification in the DDMRs, about 10% of the CpGs in the B2 DDMRs were 5-hydroxymethylated indicative of ongoing TET activity at these sites (Fig. 5g). In contrast, our analysis showed no appreciable level of 5 hmCpG at the B1a DDMRs in B1a cells (Fig. 5h). This suggests that while TET can convert a portion of mCpG generated by DNMT3A at B2 DDMRs to 5hmCpG in B2 cells, the developmentally dampened activity of DNMT3A at B1a DDMRs in B1a cells may lead to a  paucity of the mCpG substrate for TET in B1a cells. In a majority of the B lineage DDMRs, DNMT3A appears to function as a maintenance methyltransferase at the DDMRs even in nondividing B2 cells, in a manner that is distinct from the replicationcoupled maintenance methylation function performed by DNMT1. This set of DDMRs may, therefore, be viewed as the DNMT3A-maintained dynamic methylome in B1a and B2 cells.

TRANSCRIPTION FACTORS
DNMT3A-maintained enhancers (DMEs) represent a novel methylation-sensitive subset of enhancers. The DDMRs were predominantly in intergenic and intronic regions at a distance from the TSS, consistent with an enhancer function (Fig. 5i, j). To assess this, we examined chromatin features at the DDMRs using CUT&RUN with antibodies against H3K4me1, H3K4me3 and H3K27ac 42 and ATAC-Seq 43    DDMRs were enriched for the H3K4me1 histone modification, a widely studied mark of enhancers (Fig. 6a, b). Indeed, the DDMRs accounted for approximately 11% of all enhancers defined by non-TSS H3K4me1 peaks (7,990 of 74,344) in B cells.

B1a cells
Developmentally hypomethylated DDMRs (B1a DDMR patterns 1 and 7, and B2 DDMR pattern 7), acquired H3K27ac marks, which is indicative of an active enhancer state, during development ( Fig. 6a, b). B2 and B1a DDMRs overlapped with sites of enhanced chromatin accessibility. In both sets of DDMRs, chromatin accessibility was inversely correlated with CpG methylation (Spearman's correlation coefficient = −0.26 and −0.57, for B2 and B1a DDMRs, respectively). Given that the majority of the DDMRs were partially hypomethylated in wildtype B1a cells, the loss of DNMT3A had a more profound impact on the CpG methylation of DDMRs in B2 cells, concomitant with a gain in chromatin accessibility (Fig. 6b, c). This suggests that CpG methylation may constrain chromatin accessibility at these DNMT3A-maintained enhancers that we refer to as DMEs.
In B2 cells, as is typical of active enhancers, DMEs exhibited both 5mCpG and 5hmCpG, and were marked by H3K4me1 but not H3K4me3 (Fig. 6a). The few B2 DMEs (pattern 1) that were enriched in H3K4me3 were also enriched in TSS, suggesting that they are promoters ( Supplementary Fig. 6). In contrast, in B1a cells, the developmentally hypomethylated DMEs (especially pattern 1) exhibited both H3K4me1 and H3K4me3 marks, although they exhibited a minimal overlap with promoters of annotated transcripts (Fig. 6b, d, and Supplementary Fig. 6). As illustrated in Fig. 2c using publicly available ChIP-seq data, and in Fig. 6b with a different antibody in a CUT&RUN assay, we show that acquisition of H3K4me3 marks at DMEs is a unique feature of B1a cell development. H3K4me3 is prominent at promoters, but detectable levels of this modification are also observed at active enhancers bound by RNA Pol II 44 . Although CpG methylation in the DMEs was largely erased in the Dnmt3adeficient B1a and B2 cells, the ablation of Dnmt3a had little impact on H3K4me1 or H3K4me3 levels in the DMEs in either B1a or B2 cells and did not result in the acquisition of H3K4me3 marks at DMEs in Dnmt3a-deficient B2 cells (Fig. 6a, b). Complete loss of CpG methylation in the DMEs in Dnmt3adeficient B1a cells was also associated with an increase in H3K27ac marks in the B1a DMEs bearing ectopic H3K4me3 enhancer marks (Pattern 1), suggesting that CpG methylation may constrain the activity of DMEs with H3K4me3 marks (Fig. 6b, d). In contrast, in B2 cells, the levels of H3K27ac, H3K4me1 or H3K4me3 in the DMEs were largely unchanged upon DNMT3A loss.
A small number of DDMRs (patterns 10 and 11) exhibited a paradoxical gain in methylation in Dnmt3a-knockout B2 cells compared to wild-type. This may be analogous to the aberrant hypermethylation driven by DNMT3B in the absence of DNMT3A in hematopoietic stem cells 45 . Interestingly, these sites were devoid of enhancer marks in our data (Fig. 6a), indicating that they are not functional enhancers in B cells. This is corroborated by publicly available ChromHMM data on immune cells types (such as CH12 or mouse spleen) where these intervals are largely unmarked, making it harder to ascertain a specific function for these intervals (Supplementary Fig. 7) 46 . We also observed that developmentally hypermethylated B1 DDMRs (Patterns 4 and 9) bore the H3K4me1 mark at the proB1 stage but lost it in mature B1a cells, suggesting that these DDMRs represent enhancers that were developmentally decommissioned by DNMT3A-mediated CpG methylation. Failure to decommission this set of proB1 enhancers in Dnmt3a-deficient B1a cells may have contributed to the preleukemic phenotype observed in Dnmt3a-deficient B1a cells. The developmentally methylated B1a DDMRs also do not exhibit 5hmCpG marks, suggesting that DNMT3A may not be required for their continued methylation. Failure to methylate patterns 4 and 9 DDMRs is not accompanied by H3K4me3 acquisition. The impact of DNMT3A loss on the pattern of histone marks and cytosine modifications in the two largest categories of B1a and B2 DDMRs is depicted in Fig. 5D.
Dnmt3a-dependent CpG methylation of DMEs controls lineage-specific gene expression in B cells. The global loss of DME methylation in Dnmt3a-deficient B cells resulted in widespread perturbation in gene expression in both B1a (2,872 genes) and B2 (2,018 genes) B cells. The expression of 1,002 of these genes was concordantly altered in B1a and B2 cells, and the expression of 1,710 and 946 genes was altered in a lineage-specific manner in B1a and B2 cells, respectively (Fig. 7a, b). In both Dnmt3a-deficient B1a and B2 cells, there was little correlation between the differential expression of genes and the presence of a DME in their vicinity (Fig. 7c). We surmised that this was because DME-proximal genes that are differentially expressed in Dnmt3a-deficient B cells secondarily influence genes lacking a DME in their vicinity. In order to explore this notion, we repeated this analysis after stratifying the differentially expressed genes into sets, which may be enriched for master regulatory genes such as lineage-defining genes or those encoding transcription factors. We found that B1a and B2 lineage-specific genes as well as genes that are upregulated during B1a and B2 cell development were significantly enriched in the vicinity of DMEs whose methylation was maintained by DNMT3A (Fig. 7c). Interestingly, genes in the vicinity of DMEs that preserved their CpG modifications despite the loss of DNMT3A (B1a pattern 5, and B2 patterns 10 and 11) were not associated with lineage-specific, developmental or Dnmt3a-dependent changes in gene expression. We had previously noticed that the gene closest to nearly one-third of the hypomethylation canyons in Dnmt3a -/-B cells encoded a transcription factor (TF). This also suggested to us that although the presence of DMEs was not directly correlated with global changes in gene expression upon DNMT3A loss, enhancer demethylation at TF genes may result in widespread secondary changes in gene expression. Indeed, we found that the TF genes that were Fig. 4 Dnmt3a deficiency in the B lineage results in the selective expansion of B1a cells. a Expansion of B1a cells (CD19 + IgM + CD11b + CD5 + ) in the peritoneal cavity, spleen and bone marrow of Dnmt3a -/mice (mean + SD). The top panel shows the gating strategy used to distinguish B1a and B1b cells. Five wild-type and Dnmt3a -/mice were analyzed in each group. Statistically significant differences are indicated with asterisks. b Number of bone marrow B2 cell precursors (Fractions A-C′, D, E, F) and mature splenic B2 cells (newly formed (NF), follicular (FO-I and FO-II), marginal zone precursor (MZP), and marginal zone (MZ) B cells) from 8 week old Dnmt3a -/mice and their littermate controls. The experiment was performed two times and a total of eight mice in each group were analyzed. The box and whiskers mark the 25th to 75th percentiles and the minimum or maximum values. Statistically significant differences (two-sided unpaired t-tests, FDR > 1%) are indicated with asterisks (p = 0.001104 (FO-I), p = 0.000087 (FO-II), p = 0.000175 (MZ)). c IgH repertoire of peritoneal B1a cells in wild-type and Dnmt3a -/mice at varying ages (top row). Each clonotype is represented by a distinctly colored rectangle, with an area proportional to the fraction of total reads. Kappa and lambda Ig light chain staining in B1a cells from wild-type and Dnmt3a -/mice at varying ages is also shown (bottom row). d Ki-67 levels in B1a cells from the peritoneal cavity and B2 cells from the spleens of 4-week-old wild-type and Dnmt3a -/mice. Source data are provided as a Source Data file. upregulated in Dnmt3a-deficient B cells were selectively enriched in the vicinity of DMEs in both B1a and B2 cells (Fig. 7c).
During development, lineage-specific enhancers arise at sites of silent chromatin progressing through "primed", or "poised" states to the "active" state, each state being marked by specific combinations of chromatin and cytosine modifications 44 . These transitions require site-specific recruitment of chromatin and cytosine-modifying enzymes through the action of TFs. We analyzed the enrichment of TF-binding motifs in the DMEs and interpreted the data in the context of cell-stage-specific methylation states. Indeed, DMEs with similar patterns of methylation exhibited similar TF motif enrichment, suggesting that particular sets of TFs may be responsible for specific methylation patterns (Fig. 8a)   DDMRs are plotted in Fig Supplementary Fig. 5c. g and   changes in gene expression (Fig. 7c)  Interestingly, TFs not only recruit writers and erasers of cytosine modifications, but some can also function as readers of CpG methylation as their binding can be influenced by the degree of CpG methylation at their recognition sites. The impact of cytosine methylation on TF binding has been experimentally determined in vitro for over 500 TF's 50 . About half of the TF binding motifs enriched in DMEs had been reported to be sensitive to CpG methylation in the Yin et al. study 50 . Interestingly, the enriched TF-binding motifs for TFs with vital roles in B lineage specification, such EBF1, TCF3 (E2A), MEF2A, MEF2C, SP1 (PU.1), PAX5, and NFKB1 were either not sensitive to methylation or lacked a CpG site 51 . They were enriched in both methylated and unmethylated intervals (Fig. 8a). ARID3A, which lacks CpGs in its binding motif and can selectively drive B1a lineage development, is selectively enriched in B1a DMEs in the pattern 5 intervals that are methylated during B1a development 35 . In contrast, TFs that are activated in response to BCR signaling such as FOS, JUN, and NFATC1 are sensitive to CpG methylation in their binding sites. The binding of FOS and JUN to their targets is decreased by CpG methylation and the binding of NFATC1 is increased by CpG methylation. FOS and JUN motifs are enriched in hypomethylated intervals in both B1a and B2 cells. Indeed, the developmental demethylation of DMEs in the B1a lineage may render a large number of DMEs pattern (P1-12) P 1 P 7 P 9 P 1 0 P 1 1 P 1 2 P 1 P 4 P 5 P 7 P 9 P  Fig. 7 B-lineage DMEs are maintained in distinct chromatin states in B2 and B1a cells. a Heatmap of differentially expressed genes between wild-type and Dnmt3a -/-B1a and B2 cells (two replicates per condition). b Genes differentially expressed between wild-type and Dnmt3a -/-B1a and B2 cells are shown as a Venn diagram. Statistically significant overlaps (p < 2.2e-16, Fisher's exact test) are marked with an asterisk. c Significance of overlap between genes in the proximity of B1a and B2 DMEs exhibiting distinct patterns of methylation and the following differentially expressed sets of genes (i) all genes whose expression is altered by the loss of Dnmt3a in B1a and B2 cells, or (ii) B2 and B1a lineage-specific genes, or (iii) genes whose expression is altered during development (from proB to mature B), or (iv) genes encoding transcription factors whose expression is altered by the loss of Dnmt3a. The significance of the overlap was assessed using two-sided Fisher's exact test and significantly overlapping pairs of gene sets (p < 0.001) are colored in pink (B1) or blue (B2) based on the p-value. Source data are provided as a Source Data file.
responsive to FOS and JUN, which are downstream of BCR signaling.
These B lineage DMEs exhibit elevated sequence conservation across 60 vertebrate species relative to their flanking regions ( Supplementary Fig. 5b). SNPs within the corresponding intervals in the human genome have been linked to human disease in genome-wide association studies (GWAS) studies 52 . In particular, the DME syntenic regions were strongly enriched for SNPs associated with multiple human immune-mediated disorders (Fig. 8b) 53 , but a similar enrichment was not seen for most neurodevelopmental, neuropsychiatric or metabolic disorders. Common variants implicated by GWAS are known to be enriched for regulatory variants and it is indeed possible that the disease-associated SNPs in these DME syntenic regions directly contribute to disease susceptibility by impacting their enhancer function 54 . From a functional standpoint, the DMEs identified in this study can be considered to be a conserved set of methylation-sensitive enhancers in the B lineage, whose enhancer function is influenced and maintained by DNMT3A-mediated methylation.

Discussion
Our study has uncovered the underlying foundational methylome that is shared among B cells and maintained, likely by DNMT1, in both B1a and B2 sub-lineages independent of DNMT3A. We also describe a dynamic methylome largely composed of enhancers that are maintained by DNMT3A in a B1a and B2 lineage-specific manner. DNMT3A is required to counteract ongoing TET activity at these enhancers. We therefore refer to these enhancers as DNMT3A-maintained enhancers or DMEs. DMEs undergo widespread programmed CpG demethylation during B1a cell development, which coincides with the induction of the expression of lineage-specific genes in the vicinity of these enhancers. These very same enhancers remain methylated by DNMT3A in B2 cells. However, DMEs are a site of ongoing TET activity in B2 cells, as evidenced by 5hmCpG marks, and are completely demethylated in the absence of DNMT3A. Our results suggest that DNMT3A can function in a maintenance role at DME sites exhibiting ongoing cycles of TET-dependent demethylation and DNMT3A-dependent maintenance re-methylation (Supplementary Fig. 8). We also show that the developmentally demethylated DMEs acquire H3K4me3 modifications in B1a cells. Recent studies analyzing the kinetics cytosine modifications in early passage ES cells following the loss of DNMT3A showed that somatic enhancers in ES cells exhibit signs of coincident DNMT3A and TET activities 55,56 . Given the unique biology and chromatin features at these enhancers, and their role in and sub-lineage determination, we propose that DMEs may represent a specialized category of enhancers and their role in other cell types remains to be explored. Although the underlying foundational CpG methylomes that are unmasked in Dnmt3a-deficient B1a and B2 cells are virtually identical in the two sub-lineages, the cells do not lose their identity. Indeed, differences in other epigenetic marks including the ectopic H3K4me3 marks in DMEs are also preserved despite the loss of DNMT3A. This suggests that other epigenetic features including the H3K4me3 marks in DMEs are likely to be important in maintaining B1a or B2 sub-lineage identity. H3K4me3 marks are overwhelmingly associated with transcription, and their presence at DMEs may indicate the presence of enhancer RNAs. Different levels of H3K4me1 and H3K4me3 may determine whether unidirectional or bidirectional transcription occurs at enhancer sites 57,58 . In B cells, transcription of novel eRNAs has been observed in the context of changes in cytosine modification at CTCF anchors 59 . However, the B1 and B2 DMEs that we identified were not enriched in the CTCF motif, and only 10% overlapped with the set of 4,516 polyadenylated long non-coding RNAs identified in a recent study of eight B-cell types in mice 28 . We did not find evidence of transcripts specifically mapping to DMEs bearing H3K4me3 marks in our transcriptomic data. Given that we enriched polyadenylated RNAs for sequencing, non-polyadenylated enhancer RNAs may have been missed and further studies will be required to assess transcription initiation at DMEs. The role of eRNAs in regulating gene expression is not yet fully understood and the relationship between DME methylation, eRNAs transcription and gene expression will be explored in future studies. However, the loss of DME methylation in Dnmt3a-deficient cells is not without consequence and results in a perturbed transcriptional program manifesting as a preleukemic state in B1a cells.
In B1a cells, DNMT3A activity at DMEs may be opposed, at least in part, by the acquisition of H3K4me3 marks. While some studies have described the acquisition of H3K4me3 marks at enhancers in CD8 + T cells and in a tumor context, these previous studies did not evaluate them in the context of cytosine modifications 29,30 . The well-characterized mutually inhibitory relationship between CpG methylation and H3K4me3 in mammalian cells may underlie the biochemical basis of the switch in CpG methylation and chromatin modification states of DMEs between B1a and B2 cells. In mouse embryonic stem cells, the CFP1 subunit of the SET1 H3K4 trimethylase complex, which binds unmodified CpGs, constrains H3K4 trimethylase activity to hypomethylated promoters. In the absence of CFP1, ectopic H3K4me3 can accumulate in otherwise methylated enhancers that interact with promoters, resulting in the increased expression of their target genes 31 . Conversely, H3K4me3, unlike H3K4me0/ me1, is unable to free the catalytic site of DNMT3A from its autoinhibitory ATRX-DNMT3-DNMT3L (ADD) domain, causing local inhibition of DNMT activity in its vicinity 60,61 . Therefore, we hypothesize that the established recruitment of the CFP1 subunit of the SET1 complex, the primary H3K4 trimethylase, to unmethylated DNA at promoters, also applies in part to unmethylated DMEs 31,62 . Once the H3K4me3 mark is established in the DMEs perhaps in response to a transient opportunity presented by TET-mediated demethylation, H3K4me3 marks may also subsequently prevent further CpG methylation by DNMT3A. Indeed, 5hmCpG is enriched at the boundaries of the DMEs in B1a cells (Fig. 5h). As 5hmCpG is not maintained during replication, the baseline self-renewal of B1 cells in the context of reduced DNMT3A may also contribute to the complete demethylation at the DMEs. Furthermore, mice lacking TET2/3 in the B lineage completely lack B1a and marginal zone B cells 23 . It is likely that TET enzyme mediated demethylation is necessary for the H3K4 trimethylase activity at B1a DMEs.
The importance of DMEs in the B lineage is further highlighted by the fact that they are highly enriched for binding sites of B lineage-determining TFs. The striking enrichment of SNPs linked to human immune-mediated disorders within the DME syntenic regions further highlights the importance of these DMEs in maintaining cellular function, and we speculate that DMEs are likely to be of broader biological significance. Future studies using tissue-specific deletion of DNMT3A in various tissues are likely to reveal key elements of DME-mediated gene regulation in various cell types. In studies of organismal development and cell differentiation, enhancers have been identified as poised, primed, or active, depending on their histone modifications and the functional states of genes in their vicinity 44 . We propose that a subset of enhancers are DMEs and likely also undergo a similar progression from a poised or primed to an active state. However, unlike conventional enhancers, DMEs continue to exhibit localized TET activity even in the active state, which is balanced by DNMT3A. 5hmC is known to be enriched in enhancers in ES cells and is considered to mark enhancer priming in models of ES cell differentiation 63,64 . TET2 has previously been shown to drive the generation of 5hmC at enhancers in proB cells 49 . Maintenance of the appropriate level of each CpG modification at DMEs may require persistent localization of DNMT3A and TET2 at these sites through interactions with chromatin and transcription factors, as observed in epidermal stem cells 65 . CpG methylation of the majority of B-lineage DMEs is maintained by DNMT3A in B cells, but not in HSCs, suggesting that the dynamic methylome is specific to a particular cell type or lineage. Changes in the transcription factor activities that recruit DNMT3A and TET enzymes to the DMEs may underlie the loss of CpG methylation seen in B1a cells. Given that PU.1 binds DNMT3A 47 and its binding motif is enriched in DMEs, it is likely to be one of the transcription factors that can recruit DNMT3A to DMEs. Interestingly, the loss of PU.1 (SPI1) in B cells has been shown to result in a switch from a B2 to a B1a cell fate, and B1a cells exhibit a lower level of PU.1 66 . The acquisition of H3K4me3 marks at DMEs in B1a cells may also be a consequence of constitutive BCR signaling, which was recently demonstrated to be sufficient to reprogram B2 cells into a B1a cell fate 25 . Indeed, the H3K4me3 bearing B1a DMEs are enriched in FOS, JUN and NFATC1 motifs, supporting the notion that BCR signaling may play a critical role in the acquisition of ectopic H3K4me3 in the DMEs. On the contrary, DNMT3A loss alone is not sufficient for H3K4me3 acquisition in the DMEs in B2 cells or in the developmentally methylated and decommissioned enhancers in B1a cells, suggesting that the acquisition of H3K4me3 requires a B1specific factor, likely constitutive BCR signaling or a B1a cellspecific combination of lineage-determining factors such as LIN28, ARID3A, BHLHE41 or ZBTB32.
The local balance between DNMT3A and TET2 has also been suggested to maintain hypomethylation canyon boundaries in murine HSCs. The loss of DNMT3A in HSCs results in the expansion of canyons, as well as increased self-renewal 67 . B1a cell development is also accompanied by a genome-wide increase in the size of HMRs and in the number of canyons observed during B1a cell development, suggesting that the reduction in DNMT3A activity during B1a cell differentiation may not be restricted to DMEs alone but may involve a global reduction in DNMT3A activity. During development of the B1 lineage, B1a cells spontaneously acquire a CpG methylation profile that closely resembles that seen in Dnmt3a-deficient B1a and B2 cells. It has been proposed that HSCs modulate their self-renewal capacity by suppressing the level of DNMT3A via miR-29, and the increase in DNMT3A caused by the loss of miR-29a in HSCs is associated with a reduction in HSC self-renewal 68 . We speculate that analogous mechanisms may control the levels of DNMT3A in B1a cells, and consequently their self-renewal capacity. Indeed, the direct or indirect loss of DNMT3A in the B lineage has been shown to result in the leukemic transformation of B1a cells [20][21][22] . The reduced level of CpG methylation at DMEs in B1a cells may be sufficient to facilitate self-renewal; complete demethylation in the absence of DNMT3A in the knockout context may lead to uncontrolled self-renewal.
Interestingly, two recent studies describe evidence for coordinated action of DNMT3A and TET enzymes at somatic enhancers in ES cells that are highly enriched for transcription factor binding sites 55,56 . Our study suggests that a similar pattern of regulation may also operate at enhancers in fully differentiated somatic cells such as B cells. Having comprehensively analyzed the CpG methylome in an untransformed differentiated hematopoietic cell lineage, our findings may be applicable to a broader range of somatic cells. This study also illustrates how DNMT3A functions as a maintenance methyltransferase across a broad swath of B lineage-specific enhancers in a manner that is distinct from the DNA replication-coupled maintenance methyltransferase function of DNMT1. Our results also provide additional insights into how perturbation of the dynamic methylome promotes the preleukemic phase of IgHV-unmutated chronic lymphocytic leukemia, a leukemia suspected to be of B1a cell origin.
As shown here in B cells, we surmise that large subsets of enhancers, modulated by DNMT3A are likely present in many differentiated cell types and the contribution of DMEs in the dynamic methylome to the differentiation and biology of other lineages will undoubtedly be very revealing. DMEs may provide an additional mode of developmental control of enhancer function, wherein programmed CpG demethylation and acquisition of H3K4me3 marks facilitate additional modulation of enhancer activity, particularly in the context of lineage-specific gene expression. Our work also suggests how the cellular phenotypic consequences of DNMT3A-mediated epigenetic modifications can be constrained to a specific developmental lineage and how the loss of such epigenetic regulatory mechanisms may contribute to a cell type-specific malignancy.

Methods
Mice. Generation of knockout mice bearing the conditional Dnmt3a allele, in which the exon encoding the catalytic domain of DNMT3A (exon 19) is floxed (Dnmt3a tm3.1Enl , originally referred to as Dnmt3a 2lox ) has been described 39 . The mice were a gift from Drs. Taiping Chen and En Li and are presently available from the RIKEN BioResource Center (BRC strain: RBRC03731). Dnmt3a 2lox/2lox mice were backcrossed into the C57BL/6J background (Jax strain:000664) for over ten generations. They were crossed into Cd19-Cre mice (Jax strain:006785) to obtain Cd19-Cre +/-Dnmt3a 2lox/2lox mice with a conditional deletion of Dnmt3a in B cells 40 . These mice are referred to as Dnmt3a -/mice in the manuscript for simplicity. C57BL/6J mice were used as controls. All animals were housed at the Massachusetts General Hospital under specific pathogen-free conditions, with the ambient temperature set-point at 70°F, relative humidity between 30-70% and a 12 h light-dark cycle. All rodent experimental procedures were reviewed and approved by the Institutional Animal Care and Use Committee (IACUC).
Flow cytometry and isolation of B-cell populations. Single-cell suspensions from spleens, bone marrow, and peritoneal washes were analyzed using flow cytometry. The clones, suppliers and catalog numbers of the antibodies used in this study are listed in the Reporting Summary. All antibodies were used at a dilution of 1:100, except for anti-mouse CD19, B220 and CD3, which were used at a dilution of 1:200. Cells were sorted based on the gating strategy employed by the ImmGen consortium ( Supplementary Fig. 9). Peritoneal B1a B cells were identified by the following surface phenotype: CD19 + B220 lo IgM ++ CD23 -CD43 + CD5 + . Splenic B2 (follicular) B cells were identified by the following surface phenotype: CD19 + B220 + IgM + IgD + AA4.1 -CD23 + CD43 -CD5 -. For purification of proB1 cells (Fraction BC from E16 fetal liver), fetal liver cells from a total of 14 male and female E16 mouse fetuses were pooled and stained with a cocktail of PE-Cy7conjugated antibodies against Ly6c, Ter119, CD11b and CD3. Lineage-negative cells were enriched using anti-Cy7 magnetic beads (Miltenyi Biotech). The enriched lineage-negative cell fraction was then stained with a cocktail of antibodies following the gating strategy used by the ImmGen Consortium ( Supplementary  Fig. 9). Similarly, proB2 cells (Fr-BC) were FACS sorted from the bone marrow of 4-week-old mice using the gating strategy used by the ImmGen Consortium (Ly6c -Ter119 -Cd11b -CD3 -CD19 + IgM lo CD93 + B220 lo CD43 + CD24 lo ).
Single-cell suspensions were made from spleens, bone marrow, peritoneal washes or fetal liver. RBCs were lysed with 2 mL ACK lysis buffer (Lonza). Prior to surface staining, 1 ×10 6 cells were labeled with the LIVE/DEAD Fixable Blue Dead Cell Stain kit (Thermo Scientific) at 1:1000 dilution in 1 mL of HBSS for 30 min at room temperature. Fc-receptors were blocked with 2.5 μg of clone 2.4G2 (Biolegend) for 20 min on ice. Surface staining was performed using appropriate dilutions of antibodies in 100 μL of cell staining buffer (Biolegend, CA) for 30 min in the dark at 4°C or on ice. Flow cytometry and FACS sorting was performed on a 4-laser BD LSR-II, a 5-laser BD Fortessa, or a BD FACS Aria (BD Biosciences), as appropriate. For RNA-Seq studies, 100,000 live cells from each subset were directly sorted into RLT plus buffer (Qiagen) containing 1% 2-mercaptoethanol (Sigma) and stored at −80°C until the RNA extraction. For ATAC-Seq, WGBS, and TAB-Seq, 200,000 live cells from each subset were sorted into media containing 1% BSA, spun at 500 × g at 4°C for 10 min. Cells were then processed immediately for ATAC-Seq library preparation or resuspended in DNA/RNA shield (Zymo Research) and stored at −80°C prior to DNA extraction using the QiaAmp DNA Mini kit for WGBS or TAB-Seq.
RNA-Seq library preparation. CD19 + CD23 -IgM + CD5 + B1a cells (35,000-100,000 cells) were FACS sorted from the peritoneal washes of two Dnmt3a -/mice and two littermate controls at 4 weeks of age; CD19 + IgD hi IgM lo CD21 + follicular B cells (100,000 cells per sample) were sorted from spleens of two Dnmt3a -/mice and two littermate controls at 8 weeks of age. RNA was isolated from the sorted cells using an RNeasy-plus Micro kit (Qiagen). RNA-Seq libraries were prepared using the KAPA Stranded RNA-Seq Library Preparation kit (KAPA Biosystems). Libraries were assessed for quality using High Sensitivity DNA chips on the Agilent Bioanalyzer, quantified using Qubit fluorometer (Thermo Fisher Scientific), as well as KAPA Library Quantification kit (KAPA Biosystems), and sequenced on the Illumina NextSeq 550 platform.
WGBS library preparation. DNA was isolated from peritoneal cavity B1a cells, splenic follicular B cells (Dnmt3a -/vs. WT), fetal liver proB1 cells, and bone marrow proB2 cells. WGBS libraries were prepared using 100 ng of DNA using the Pico Methyl-Seq™ Library Prep Kit (Zymo Research) according to the manufacturer's protocol. The quality of the libraries was assessed using High Sensitivity D5000 Screen Tapes (Agilent 4200 Tapestation) and quantified with KAPA Library Quantification kit. Libraries were sequenced on an Illumina NextSeq 550. Fewer CpGs were modified in proB1 and B1a cells (74.1% and 73.1%, respectively) than in the proB2 and B2 cells (75.9% and 76%, respectively) (Supplementary Table 2). Since CHH and CHG methylation (where H represents A, T, or C) represented <2% of all cytosine methylation in proB and B cells (Supplementary Table 2), we focused our analyses on CpG modifications.
TAB-Seq library preparation. TAB-Seq libraries were constructed using the TAB-Seq library preparation kit from Wisegene (catalog #K001) 69 . Briefly, 5hydroxymethyl cytosines in 500 ng genomic DNA with 0.1% spike-in control DNA (synthesized with 5-hydroxymethyl cytosines or 5-methyl cytosines exclusively) were protected using T4 β-glucosyltransferase, and 5-methyl cytosines were oxidized using recombinant TET1. Following TET1 oxidation, the DNA was subjected to bisulfite conversion and PCR-based sequencing adapter tagging as described for WGBS above (Pico Methyl-Seq Library Prep Kit, Zymo Research). The library was quantified and sequenced on an Illumina Nextseq 550. Genome-wide, 2.6% of CpGs were hydroxymethylated in B1a cells and 3.8% of CpGs were hydroxymethylated in B2 cells (Supplementary Table 2).
ATAC-Seq library preparation. ATAC-Seq libraries were prepared as originally described 43 . Fifty thousand freshly sorted cells were pelleted and washed with 50 μL chilled 1× PBS, and with 50 μL lysis buffer (10 mM Tris-HCl pH 7.4, 10 mM NaCl, 3 mM MgCl2, 0.1% IGEPAL CA-630). Three replicates were analyzed per condition. The nuclei were pelleted in a cold micro-centrifuge at 550 × g for 10 min and resuspended in a 50 μL transposition reaction with 2.5 μL Tn5 transposase (FC-121-1030; Illumina) to tagment open chromatin. The reaction was incubated at 37°C in a Thermomixer (Eppendorf) at 300 rpm for 30 min. Tagmented DNA was purified using a QIAGEN MinElute kit and amplified with 7 or 11 cycles of PCR, based on the results of a test qPCR. ATAC-Seq libraries were then purified using a QIAGEN PCR cleanup kit and quantified using KAPA library quantification kit (KAPA Biosystems, Roche) and sequenced on the Nextseq 550 platform.
CUT&RUN library preparation. CUT&RUN libraries were prepared using rabbit monoclonal antibodies against H3K4me1, H3K4me3 and H3K27Ac (Cell Signaling Technology catalog numbers 5326, 9751, and 8173, respectively) as previously described 42 . The protein-A Mnase fusion protein from a gift from Steven Henikoff.
Preparation and analysis of Ig Heavy Chain-Repertoire libraries. RNA from B1a cells was reverse transcribed and amplified using a Qiagen OneStep RT-PCR kit (Qiagen Inc., Valencia, CA, USA) and iRepertoire® mouse BCR heavy chain (MBHI) primers (iRepertoire Inc., Huntsville, AL, USA), following the iRepertoire user manual for the Illumina sequencing platform. Following manufacturer's instructions in the iRepertoire manual, Illumina adapters were added using a Qiagen Multiplex PCR kit, and a 350-500 bp bands corresponding to the Ig heavy chain amplicons were gel purified using a Qiagen Gel Extraction kit. The libraries were quantified using Qubit dsDNA quant (Life Technologies) and sequenced on the Illumina MiSeq platform (paired end, 250 bp). For samples from aged mice bearing clonal leukemias, the libraries were Sanger sequenced using the Illumina universal primers as the sequencing primer. Repertoire data was analyzed using the iRweb online analysis platform provided by iRepertoire Inc.
Analysis of WGBS and TAB-Seq data. The first 5 and last 2 bases were trimmed using cutadapt v1.11 and the trimmed reads were aligned to the mm10 mouse reference genome using Bismark v0.14.3 70 with the parameters: --bowtie2 --non_directional --un --ambiguous --phred33-quals. PCR duplicates were removed using Picard (v2.5.0). Approximately 5× coverage was obtained after removal of duplicates (Supplementary Table 2). The Bismark package was used to extract the coverage and methylation status of individual cytosines. We then used the Meth-Pipe v3.4.3 26 to identify differentially methylated regions. We implemented MethPipe in the following steps. First, we used the program to-mr to convert the output BAM file from Bismark to mr format. Second, we removed duplicate reads using duplicate-remover to get rid of the reads resulting from PCR overamplification before calculating the methylation level. Third, the bisulfite conversion rate, defined as the rate at which unmethylated cytosines in the sample appear as thymidines in the sequenced reads, was measured by the program bsrate from MethPipe based on the non-CpG cytosines, which are believed to be almost completely unmethylated in most mammalian cells. In all, 98-99% bisulfite conversion efficiency was observed. Finally, the methcounts program was used to estimate the methylation level for individual cytosines, calculated as a probability based on the ratio of methylated to total mapped reads at the respective cytosines. Symmetric CpGs were merged using symmetric-cpgs. The bedgraph output from Bismark was converted into the bigwig format for visualization on IGV and for plotting heatmaps using the SeqPlots package 71,72 . WGBS data for HSC (GSE49714) and CD8 + T cells (GSE107150) were additionally downloaded and reprocessed using the same methods and parameters as used in this study 27,73 . A similar approach was used for the analysis of TAB-Seq data. Analysis of TAB-Seq spike-in controls showed an efficiency of~90% for TET oxidation and~98% for bisulfite conversion.
Detection of HMRs (hypomethylated regions) and canyons. HMRs were called using the program hmr in the MethPipe package, which utilizes a hidden Markov model (HMM)-based approach. A length cutoff of 3.5 kb or longer was used to identify canyons, as previously described 27 . Similarity in the HMRs among different cell types was assessed based on the Jaccard statistic using the BEDtools suite 74,75 . The Jaccard statistic between the two sets of genomic intervals ranges from 0 to 1 (no overlap to complete overlap, respectively).
Identification and annotation of DMRs. To compare two methylomes, we identified differentially methylated regions (DMRs) using the programs methdiff and dmr in the MethPipe package. The first step was to calculate the differential methylation score for each CpG site using methdiff. Then, the methdiff and hmr output files were used to compute DMRs using the program dmr. Finally, we obtained a filtered set of DMRs spanning at least 10 CpGs and having at least 5 significantly differentially methylated CpG as recommended by the MethPipe user manual. In addition, the roimethstat program was used to estimate the average methylation level in each DMR across all samples of interest.
To annotate the genomic features of DMRs (promoter, intergenic, intron, exon, 3′ UTR, 5′ UTR, TTS) and identify DMR-proximal genes, we used the annotatePeaks.pl script from HOMER (v4.9) 76 . We also used the HOMER scripts findGO.pl to perform functional enrichment tests. The findMotifsGenome.pl script was used to identify enriched vertebrate transcription factor binding motifs from JASPAR and HOMER databases in the DMR sets of interest (filtered using a pvalue cutoff <0.01 and a difference between test sets and background intervals >5%; the full set of DMRs was used as the background).
The DMRs were grouped into patterns based on the methylation dynamics across proB, mature B and Dnmt3a -/-B cells. We created a list of 12 pattern vectors of interest representing CpG methylation ratios (range [0,1]) across the three conditions (low, intermediate and high) ( Supplementary Fig. 5a). In each pattern vector, "low," "intermediate," and "high" values were set at CpG methylation ratios of 0.2, 0.5, and 0.8, respectively. Each vector represents the CpG methylation ratios across proB, wild-type B, and Dnmt3a -/-B conditions for a particular B-cell subset. We assigned each B1a and B2 DMR to the pattern vector that yielded the maximum Pearson correlation coefficient (PCC). No DMR was assigned to more than one pattern.
ChIP-Seq and ATAC-Seq analysis. We obtained H3K4me1, H3K4me3 ChIP-Seq data in B1a and B2 cells (GSE72017) and H3K4me1 ChIP-Seq data in mouse HSCs (GSE63276) from the NCBI GEO database 28,77 . The ATAC-Seq libraries were generated and sequenced in our laboratory as described above 43 . ChIP-Seq and ATAC-Seq reads were aligned to the mm10 mouse reference genome. We removed duplicated and mitochondrial reads, and generated bigwig files for plotting heatmaps using the SeqPlots package 72 . ATAC-Seq or ChIP-Seq peaks were called using MACS2 with a fold change >5 as a cutoff 78 .
Transcriptomic analysis. We aligned the RNA-Seq reads to the mm10 reference genome and performed transcript quantification using RSEM v1.25.0. Refseq gene annotations were used. We used EBSeq 79 to identify differentially expressed genes (DEGs) between two biological conditions; we used a cutoff of 0.95 for the posterior probability of differential expression reported by EBSeq to determine differentially expressed genes with a false discovery rate controlled at 5%. We used the HOMER findGO.pl script 76 to perform functional enrichment tests. Fisher's exact test was used to evaluate the significance of overlap between two gene sets, using all expressed genes as the background.
We also downloaded microarray data for developing and mature B1a and B2 cells (proB cells from bone marrow or fetal liver, follicular B cells from spleen, FrE from fetal liver or bone marrow, and B1a cells from peritoneal cavity) from the ImmGen consortium (GSE15907) and identified differentially expressed genes (p < 0.01) using the RankProd R package 34,80,81 .
SNPs enrichment test in DMRs. We used the probabilistic framework implemented by RiVIERA to infer the association of DMRs and SNPs of 27 human traits from four major disease categories (NeuroDegen, NeuroPsych, Immune, and Metabolic) 52 .
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
All sequencing datasets generated in the current study have been deposited in the NCBI GEO repository (GSE150497). Source data are provided with this paper.