Marta Kulis

Angelika Merkel

Simon Heath

Ana C Queirós

Ronald P Schuyler

Giancarlo Castellano

Renée Beekman

Emanuele Raineri

Anna Esteve

Guillem Clot

Nuria Verdaguer-Dot

Martí Duran-Ferrer

Nuria Russiñol

Roser Vilarrasa- Blasi

Simone Ecker

Vera Pancaldi

Daniel Rico

Lidia Agueda

Julie Blanc

David Richardson

Laura Clarke

Avik Datta

Marien Pascual

Xabier Agirre

Felipe Prosper

Diego Alignani

Bruno Paiva

Gersende Caron

Thierry Fest

Marcus O Muench

Marina E Fomin

Seung-Tae Lee

Joseph L Wiemels

Alfonso Valencia

Marta Gut

Paul Flicek

Hendrik G Stunnenberg

Reiner Siebert

Ralf Küppers

Ivo G Gut

Elías Campo

José I Martín-Subero

José I Martín
email: imartins@clinic.ub.es

Whole-genome fingerprint of the DNA methylome during human B cell differentiation

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

The multitude of cell types and tissues of an organism can be defined by their particular epigenetic makeup [START_REF] Ziller | Charting a dynamic DNA methylation landscape of the human genome[END_REF][START_REF] Bernstein | The mammalian epigenome[END_REF] . DNA methylation is an important component of the epigenome which is extensively modulated during regulatory and developmental processes, both in the context of physiological and pathological conditions [START_REF] Smith | DNA methylation: roles in mammalian development[END_REF][START_REF] Bergman | DNA methylation dynamics in health and disease[END_REF][START_REF] Bird | DNA methylation patterns and epigenetic memory[END_REF] . Although recent reports have analyzed the DNA methylation profiles of various cell types at the whole-genome scale 1,6- 16 , the DNA methylome of a single human cell type during its complete differentiation process has not been defined so far. The B-cell lineage represents a paradigmatic cellular model to study the dynamic epigenome during cell development and specification because major B-cell maturation stages have distinct phenotypic and gene expression features and can be isolated in sufficient numbers from hematopoietic tissues [START_REF] Lee | A global DNA methylation and gene expression analysis of early human B-cell development reveals a demethylation signature and transcription factor network[END_REF][START_REF] Lai | DNA methylation profiling in human B cells reveals immune regulatory elements and epigenetic plasticity at Alu elements during B cell activation[END_REF][START_REF] Shaknovich | DNA methyltransferase 1 and DNA methylation patterning contribute to germinal center B-cell differentiation[END_REF] .

B-cell lymphopoiesis is a complex and tightly coordinated process guided by a hierarchical expression of different stage-specific transcription factors and microenvironmental influences [START_REF] Matthias | Transcriptional networks in developing and mature B cells[END_REF][START_REF] Kurosaki | B cell signaling and fate decision[END_REF] . The process starts in the bone marrow, where hematopoietic stem cells differentiate into multipotent progenitors and common lymphoid progenitors, which then commit to the B-cell lineage and give rise to precursor B cells. These precursors gradually rearrange their immunoglobulin genes and differentiate into mature naive B cells, which leave the bone marrow to enter the blood stream. Resting naive B cells transit through lymph nodes and, eventually, they are activated by specific antigens via activation of the B-cell receptor, which induces the germinal center reaction. Germinal center B cells further rearrange and mutate their immunoglobulin genes, rapidly proliferate and differentiate. Finally, the germinal center reaction gives rise to plasma cells producing large amounts of high-affinity antibodies and memory B cells. Plasma cells exiting the lymph nodes migrate to the bone marrow where they can reside for extended periods of time, and longlived memory B cells recirculate through the blood and lymphoid organs, providing the basis for enduring humoral immunity [START_REF] Manz | Lifetime of plasma cells in the bone marrow[END_REF][START_REF] Slifka | Humoral immunity due to long-lived plasma cells[END_REF] . Hence, an interesting feature of the B-cell maturation process is that it entails a variety of cells with different functional features, proliferation abilities, microenvironmental influences and life spans, providing an exceptional opportunity to study the epigenome in the context of different biological processes, and to provide insights into the fields of cell differentiation, B-cell biology, cancer and aging.

RESULTS

Whole-genome DNA methylation maps of B-cell subpopulations

We generated unbiased DNA methylation maps of uncommitted hematopoietic progenitor cells (HPCs) and five B-cell lineage subpopulations, including pre-B-II cells (preB2Cs), naive B cells from peripheral blood (naiBCs), germinal center B cells (gcBCs), memory B cells from peripheral blood (memBCs) and plasma cells from bone marrow (bm-PCs), by wholegenome bisulfite sequencing (WGBS) (Fig. 1a and Supplementary Table 1). We sequenced two biological replicates of each subpopulation and a total of 2,217 billion base pairs (bp) of which 85-95% could be mapped (mean depth of 54-fold per sample) (Supplementary Table 2). On average, we measured methylation levels of 22.7 million CpGs per sample (ranging from 21 to 25 million). Unsupervised principal component analysis (PCA) of CpG methylation levels showed that B-cell subpopulations segregate according to their developmental stage (Fig. 1b).

Globally, B-cell differentiation is accompanied by a gradual widespread demethylation of the genome, which was more pronounced at late differentiation stages such as memBC and bm-PC (Fig. 1c-e). The global methylation status of CpGs was largely bimodal in all sorted cell populations and the level of partially methylated regions increased to 19-24% in advanced maturation stages (Fig. 1e). This result contrasts to other WGBS studies using whole tissues, in which the proportion of partially methylated regions is usually high [START_REF] Hon | Epigenetic memory at embryonic enhancers identified in DNA methylation maps from adult mouse tissues[END_REF] , and highlights the importance of using purified cell subpopulations for DNA methylation studies.

The results obtained by WGBS were complemented with the analysis of 3 to 9 replicates of 10 different B-cell subpopulations by high-density DNA methylation microarrays [START_REF] Bibikova | High density DNA methylation array with single CpG site resolution[END_REF] (Fig. 1a). These subpopulations included those analyzed by WGBS as well as pre-B-I cells (preB1Cs), immature B cells (iBC), tonsillar naive B cells (t-naiBCs) and tonsillar plasma cells (t-PCs) (Supplementary Table 1). The biological replicates of each subpopulation analyzed by WGBS or microarrays showed high reproducibility (correlation coefficient > 0.95) (Supplementary Fig. 1), and both high-throughput techniques were further validated by bisulfite pyrosequencing (Supplementary Fig. 2). In line with WGBS, PCA of microarray data separated B-cell subpopulations mostly according to their developmental stage (Fig. 1f) and we observed a gradual global methylation loss throughout B-cell maturation (Fig. 1g).

Interestingly, samples clustered into "antigen-inexperienced" and "antigen-experienced" cells (Supplementary Fig. 3). This finding indicates that proliferative gcBCs start a massive reconfiguration of the DNA methylome [START_REF] Lai | DNA methylation profiling in human B cells reveals immune regulatory elements and epigenetic plasticity at Alu elements during B cell activation[END_REF][START_REF] Shaknovich | DNA methyltransferase 1 and DNA methylation patterning contribute to germinal center B-cell differentiation[END_REF] which continues in cell subpopulations with extended life span such as memBCs and bm-PCs.

Demethylation of non-CpG sites upon B-cell commitment

Cytosine methylation in mammals can occur outside CpGs, a phenomenon commonly observed in embryonic stem cells (ESCs) and neurons [START_REF] Lister | Global epigenomic reconfiguration during mammalian brain development[END_REF][START_REF] Lister | Human DNA methylomes at base resolution show widespread epigenomic differences[END_REF][START_REF] Varley | Dynamic DNA methylation across diverse human cell lines and tissues[END_REF] . In our study, considering non-CpG sites with high-confidence methylation estimates and having ruled out suboptimal bisulfite conversion and sequence variants, we did observe non-CpG methylation, which was primarily confined to the most undifferentiated cells (i.e. HPCs) and occurred mainly in a CpApC context (Fig. 2a-c). Significant non-CpG methylation in HPCs was detected by WGBS at 25,763 sites in replicate 1 and 16,838 in replicate 2, with a mean methylation levels of 25.1% and 24.7%, respectively. Non-CpG methylation in HPCs frequently targeted the same sites in the two biological replicates, fact that was also confirmed by bisulfite pyrosequencing in independent samples (Fig. 2d and Supplementary Fig. 4). Non-CpG methylation in HPCs was preferentially located in gene bodies (both introns and exons) and depleted in lamina-associated domains (Supplementary Fig. 5). Although the methylation microarrays only measured 3,091 non-CpG sites, we could confirm the WGBS results (Fig. 2e,f). Similarly to reports in ESCs [START_REF] Ziller | Genomic distribution and inter-sample variation of non-CpG methylation across human cell types[END_REF] , methylated non-CpGs were flanked by methylated CpGs. However, demethylation of non-CpGs and CpGs was not simultaneous. We observed a dramatic demethylation of non-CpGs upon B-cell commitment in preB2Cs, but 97% of these demethylated non-CpGs remained flanked by methylated CpGs (Supplementary Fig. 6). These data were confirmed by bisulfite pyrosequencing (Supplementary Fig. 4) and indicate that non-CpG methylation is passively erased from HPCs to preB2Cs without simultaneous demethylation of flanking CpGs (Fig. 2g).

Although previous studies have reported that high expression of DNMT3A and DNMT3B is associated with non-CpG methylation [START_REF] Ziller | Genomic distribution and inter-sample variation of non-CpG methylation across human cell types[END_REF][START_REF] Arand | In vivo control of CpG and non-CpG DNA methylation by DNA methyltransferases[END_REF] , we did not identify consistent differences in the expression levels of DNMTs between HPCs and preBCs (Supplementary Fig. 7).

Identification of dynamic DNA methylation patterns

Next, we focused our analysis on the modulation of CpG methylation throughout the complete B-cell maturation program. We identified dynamic methylation levels in 4.93 million CpGs (> 0.25 methylation change in one set of samples and > 0.1 in the other), which represent 30.6 % of the 16.1 million CpGs with methylation estimates in all 12 samples analyzed. The cell subtypes showing the most pronounced methylation changes compared to the preceding stage were gcBCs, memBCs and bm-PCs (Fig. 3a and Supplementary Fig. 8). Interestingly, multiple genes directly involved in B-cell differentiation (e.g. ARID3A, BCL2, BLK, EBF1 and IRF4) show a complex modulation of their DNA methylation profile across the gene length, with different regulatory elements losing methylation at distinct maturation stages (Fig. 3b and Supplementary Fig. 9). In a previous study, 5.6 million dynamic CpGs were identified to be differentially methylated in a wide range of human cell types and tissues [START_REF] Ziller | Charting a dynamic DNA methylation landscape of the human genome[END_REF] . Only one third of the dynamic CpGs in B-cell differentiation defined in our study overlapped with these CpGs (Supplementary Fig. 10). Although this result may in part be caused by differences in coverage and bioinformatic pipelines used in the two studies, it suggests that the majority of dynamic methylation in B cells may entail CpGs not previously detected to be differentially methylated in other cell types and tissues.

Similarly to WGBS, although at lower scale, microarray data revealed that 22.4% of the measured CpGs showed variable methylation levels during B-cell differentiation (n = 106,562, see online methods for details). Due to the more comprehensive set of samples, we used these data to define sets of CpGs sharing similar DNA methylation dynamics during B-cell differentiation. To that end, we based our next analysis on a linear model of B-cell maturation, from HPC to terminal differentiation of gcBC into t-PC and long-lived bm-PC. We defined 20 major modules containing at least 500 CpGs that could be classified according to four general DNA methylation modulation patterns during differentiation: i) methylation levels decrease (9 modules), ii) methylation levels first decrease and then increase (3 modules), iii) methylation levels first increase and then decrease (2 modules) and iv) methylation levels increase (6 modules) (Supplementary Data 1). Three of these 20 modules covered 57.8% of all dynamic CpGs: module 8 defined by CpG demethylation starting in gcBCs and continuing in bm-PCs (n = 34,604), module 9 showing demethylation only in bm-PCs (n = 13,044) and module 20 characterized by hypermethylation in bm-PCs only (n = 13,949) (Fig. 3c). These data confirmed that DNA hypomethylation mainly occurs in the gcBCs and bm-PCs. Furthermore, it could be appreciated that hypermethylation is a late event, mainly occurring in bm-PCs. Overall, we may conclude that 84.5% of the dynamic CpGs either gain or lose methylation (DNA methylation modulation patterns i and iv, n = 90,070) as B-cell differentiation progresses.

Hence, each B-cell differentiation stage has its specific DNA methylation pattern, but furthermore retains an epigenetic memory of the previous stages. Interestingly, although the methylation levels of the 20 modules are able to separate B-cell maturation stages, we performed a complexity reduction step and selected 5 CpGs in genes important for B-cell differentiation such as BLK, SEMA4B, ARID3A, AICDA and PRDM1 whose methylation levels could accurately classify each maturation stage (Supplementary Fig. 11).

Functional analysis of dynamically-methylated regions

In general, CpGs losing methylation at any B-cell maturation stage were preferentially located in introns, intergenic regions and repetitive elements (e.g. LINEs, SINEs and LTRs), and were enriched for genomic areas lacking CpG islands (CGIs) (Supplementary Figs. 12 and13).

In contrast, CpGs gaining methylation were enriched for CGIs and promoter regions (Supplementary Fig. 12). Additionally, we classified differentially methylated sites using a categorization of the genome into different chromatin states in immortalized mature B cells [START_REF] Ernst | Mapping and analysis of chromatin state dynamics in nine human cell types[END_REF] (Supplementary Fig. 14). Both WGBS and microarrays revealed that the majority of dynamic CpGs during B-cell differentiation were enriched for enhancer regions (mainly intragenic), polycomb-repressed regions or heterochromatin (Fig. 3d,e and Supplementary Figs. 15 and16). Demethylation in precursor B cells was mostly related to enhancer elements whereas that occurring exclusively from gcBCs onwards was preferentially located in heterochromatic regions. Gain of CpG methylation was a rare event in early B-cell differentiation but rather frequent in mature B cells, especially in bm-PCs. Such CpG hypermethylation targeted preferentially polycomb-repressed regions (Fig. 3d,e).

Next, we studied the mechanisms underlying enhancer demethylation in the B-cell differentiation process. We identified a significant enrichment (FC > 2, P < 0.01) in transcription factor binding sites (TFBSs) of key B-cell transcription factors such as BCL11A, EBF1, IRF4, MEF2A, MEF2C, PAX5 or TCF3 (E2A) (Fig. 4a, Supplementary Fig. 17 and Supplementary Data 2 and 3). As B-cell commitment is associated with the expression of lineage-specific TFs [START_REF] Matthias | Transcriptional networks in developing and mature B cells[END_REF] , we analyzed the transition from HPC to preB1C (module 1) in detail and, globally, we observed an inverse correlation between TF expression and the methylation levels of their binding sites (Fig. 4b,c). We further investigated the association between TF expression and the methylation status of their binding sites during the entire differentiation program, and in general, once a TFBS becomes demethylated at any B-cell differentiation stage, it remains unmethylated in subsequent stages, suggesting an epigenetic memory of TF binding 24 (Supplementary Fig. 18).

At the functional level, genes within microarray-based methylation modules enriched for enhancer elements were involved in multiple immune system-related functions (Supplementary Fig. 19 and Supplementary Data 4). In contrast, CpGs in modules enriched for heterochromatin or polycomb-repressed regions, targeted genes not involved in the immune system but rather affected terms such as development, locomotion or behavior (Supplementary Fig. 19 and Supplementary Data 4). Based on this observation, we hypothesized that differential methylation in enhancer elements may be globally associated with gene expression, while that affecting inactive elements (i.e. heterochromatin or polycomb-repressed regions) may not. We initially explored the transcriptome of B-cell subpopulations and observed that they clustered separately using an unsupervised approach (Fig. 4d). We then calculated the number of genes differentially expressed and differentially methylated by comparing adjacent cell subpopulations. Globally, there was a poor association, and large transcriptional changes could be related to a minor modulation of the DNA methylome and vice versa (Fig. 4e). We further explored the association between DNA methylation patterns and gene expression by exploring genes with dynamic methylation in enhancers, heterochromatin and polycomb-repressed regions. We observed that both variability and mean expression levels of genes containing dynamic CpGs in enhancer elements were much higher in comparison with genes showing modulation of CpG methylation in nonfunctional chromatin states such as heterochromatin and polycomb-repressed regions (P < 0.001; Fig. 4f-h). However, regardless of the DNA methylation pattern throughout B-cell maturation, dynamic CpGs target genes with higher expression levels and variation throughout the whole B-cell differentiation. Thus, these data suggest that there is no direct correlation between DNA methylation and gene expression, as previously shown elsewhere [START_REF] Hovestadt | Decoding the regulatory landscape of medulloblastoma using DNA methylation sequencing[END_REF][START_REF] Kulis | Epigenomic analysis detects widespread gene-body DNA hypomethylation in chronic lymphocytic leukemia[END_REF] , but rather that dynamic CpGs affecting functional elements target immune system-related genes whose expression is modulated during B-cell differentiation.

Next, we sought to analyze the functions of genes with dynamic enhancer methylation in more detail. As B-cell receptor signaling is a key element specific for B-cell differentiation, we studied the methylation pattern of 41 genes involved in this function. Out of those, 38 (93%) had a total of 234 CpGs with dynamic methylation and were preferentially located in enhancers (Supplementary Fig. 20). Furthermore, we observed that 38% of all enhancers (P < 0.001) belonged to a B-cell specific functional gene network 29 (Supplementary Fig. 21).

As our study comprises cell subpopulations isolated from three different compartments (bone marrow, peripheral blood and tonsil), we aimed at detecting particular DNA methylation imprints related to these different locations but we did not identify any consistent pattern. However, we did observe that naiBCs isolated from different compartments (i.e. peripheral blood or tonsils) showed a massive modulation of their transcriptome while their methylomes remained virtually identical (Supplementary Figs. 22 and23). Based on this observation, we may hypothesize that the gene expression changes in t-naiBCs are essential to optimize antigen recognition in the tonsil, followed by successful B-cell activation. However, when no antigen is found, the t-naiBCs will re-enter the blood stream.

Hence, the t-naiBCs state has to be reversible, which may be the reason why no changes were observed at the level of the DNA methylome.

Epigenetic link among long-lived B cells, cancer and aging

We observed that long-lived B cells such as memBCs and bm-PCs [START_REF] Slifka | Humoral immunity due to long-lived plasma cells[END_REF] show an extensive perturbation of their DNA methylome as compared to gcBCs (Supplementary Fig. 8).

Remarkably, although bm-PCs and newly generated t-PCs share a similar transcriptome, they have different methylomes (Supplementary Fig. 22). The major DNA methylation changes in memBCs and bm-PCs affect presumably non-functional elements such as heterochromatin and polycomb-repressed regions. Thus, our findings suggest that part of the epigenetic makeup of memBCs and bm-PCs is caused by an epigenetic drift associated with their long life span [START_REF] Slifka | Humoral immunity due to long-lived plasma cells[END_REF] . To determine whether such drift may be related to the expression of DNMTs, we analyzed DNMT1, DNMT3A and DNMT3B by qRT-PCR in sorted gcBCs, memBCs, t-PCs and bm-PCs.

DNMT1 showed low levels in memBCs, t-PCs and bm-PCs (Supplementary Fig. 7). As these cell types are considered to be non-proliferative, our finding implies that demethylation occurs either through an active mechanism or passively if these cell types still proliferate at a low rate [START_REF] Tooze | A Replicative Self-Renewal Model for Long-Lived Plasma Cells: Questioning Irreversible Cell Cycle Exit[END_REF] . Gain of methylation of polycomb-repressed regions has been linked to the activity of DNMT3A and DNMT3B [START_REF] Vire | The Polycomb group protein EZH2 directly controls DNA methylation[END_REF][START_REF] Cedar | Linking DNA methylation and histone modification: patterns and paradigms[END_REF] . We detected that as compared to gcBCs, DNMT3A was upregulated in t-PCs and bm-PCs, and at a lower extent also in memBCs. Thus, these results may suggest that minor increases of this enzyme in the context of long-lived cells could result in hypermethylation of polycomb-repressed areas. To obtain further insights into the mechanisms of chromatin repression by hypermethylation at polycomb-repressed regions, we performed bisulfite sequencing of chromatin immunoprecipitated DNA (ChIP-BS) with a H3K27me3 antibody [START_REF] Brinkman | Sequential ChIP-bisulfite sequencing enables direct genome-scale investigation of chromatin and DNA methylation cross-talk[END_REF][START_REF] Statham | Bisulfite sequencing of chromatin immunoprecipitated DNA (BisChIP-seq) directly informs methylation status of histone-modified DNA[END_REF] . This experiment suggests that in memBCs, H3K27me3 and DNA methylation coexist both in CpG-rich and CpG-poor regions, and that the DNA methylation levels within nucleosomes containing H3K27me3 were lower than regions outside such nucleosomes (Supplementary Fig. 24).

Remarkably, demethylation of heterochromatin (in part bound to the nuclear lamina) and hypermethylation of polycomb-repressed regions are among the most frequent epigenetic changes in solid and hematological tumors [START_REF] Berman | Regions of focal DNA hypermethylation and long-range hypomethylation in colorectal cancer coincide with nuclear lamina-associated domains[END_REF][START_REF] Martin-Subero | New insights into the biology and origin of mature aggressive B-cell lymphomas by combined epigenomic, genomic, and transcriptional profiling[END_REF][START_REF] Ohm | A stem cell-like chromatin pattern may predispose tumor suppressor genes to DNA hypermethylation and heritable silencing[END_REF][START_REF] Schlesinger | Polycomb-mediated methylation on Lys27 of histone H3 premarks genes for de novo methylation in cancer[END_REF][START_REF] Widschwendter | Epigenetic stem cell signature in cancer[END_REF] . To analyze whether those CpGs showing differential methylation in long-lived B cells overlap with those becoming differentially methylated in cancer, we used methylation data from various lymphoid neoplasms including acute lymphoblastic leukemia [START_REF] Busche | Integration of high-resolution methylome and transcriptome analyses to dissect epigenomic changes in childhood acute lymphoblastic leukemia[END_REF] (ALL, n = 46), chronic lymphocytic leukemia [START_REF] Kulis | Epigenomic analysis detects widespread gene-body DNA hypomethylation in chronic lymphocytic leukemia[END_REF] (CLL, n = 139), diffuse large B-cell lymphoma (DLBCL, n = 40) and multiple myeloma [START_REF] Agirre | Whole-epigenome analysis in multiple myeloma reveals DNA hypermethylation of B cell-specific enhancers[END_REF] (MM, n = 104). We analyzed the methylation levels of CpGs losing methylation (modules 8 and 9) and those gaining methylation in memBCs and bm-PCs (modules 19 and 20). Indeed, we observed that, considering these CpGs, neoplastic cells show a DNA methylation profile similar to memBCs and bm-PCs (Fig. 5a-d and Supplementary Fig. 25). To further evaluate the epigenetic link between normal B-cell differentiation and neoplastic transformation, we compared the DNA methylome of B-cell neoplasms with their normal cell counterparts, e.g. ALL [START_REF] Busche | Integration of high-resolution methylome and transcriptome analyses to dissect epigenomic changes in childhood acute lymphoblastic leukemia[END_REF] vs. preBCs, DLBCL-GCB type vs. gcBCs, and MM 40 vs. PCs (Fig. 6a and Supplementary Data 5). The results indicate that a large fraction of the CpGs differentially methylated in tumors are dynamically methylated during normal B-cell differentiation, ranging from 53% to 82% for hypermethylated and from 29% to 84% for hypomethylated sites (Fig. 6b). Interestingly, we identified that hypomethylation in ALL was enriched in enhancers whereas hypomethylation in DLBCL and MM predominantly affected heterochromatin (Fig. 6c). Additionally, although ALL cells are arrested at the preBCs stage, they acquire hypermethylation in polycomb-repressed regions, which is characteristic of more mature differentiation stages (Fig. 6d). MM cells, in contrast, do not acquire hypermethylation of polycomb-repressed regions because their cellular origin already shows this feature, but, as they downregulate the B-cell program, they acquire hypermethylation of B cell-specific enhancers [START_REF] Agirre | Whole-epigenome analysis in multiple myeloma reveals DNA hypermethylation of B cell-specific enhancers[END_REF] .

Finally, global hypomethylation accompanied by local hypermethylation of polycomb targets is also a molecular hallmark of aging [START_REF] Lopez-Otin | The hallmarks of aging[END_REF][START_REF] Horvath | Aging effects on DNA methylation modules in human brain and blood tissue[END_REF][START_REF] Hannum | Genome-wide methylation profiles reveal quantitative views of human aging rates[END_REF][START_REF] Heyn | Distinct DNA methylomes of newborns and centenarians[END_REF][START_REF] Maegawa | Widespread and tissue specific age-related DNA methylation changes in mice[END_REF] . We analyzed DNA methylation values of CpGs within modules 8, 9, 19 and 20 using data from 694 peripheral blood samples from donors ranging from 0 to 101 years [START_REF] Hannum | Genome-wide methylation profiles reveal quantitative views of human aging rates[END_REF][START_REF] Heyn | Distinct DNA methylomes of newborns and centenarians[END_REF] . Indeed, we detected a significant correlation with age (P < 0.01) following the same tendency as in long-lived cells and cancer cells, i.e. methylation of heterochromatin diminished with age whereas methylation of polycomb targets increased (Fig. 5e-f; Supplementary Fig. 25). We next compared B cells with short and long lifespan (naiBCs and memBCs, respectively) isolated simultaneously from individuals of different age and we detected that memBCs acquire hypomethylation of heterochromatin and hypermethylation of polycomb targets regardless of the chronological age of the donor (Supplementary Fig. 26). As the cellular composition of blood changes with age, our results imply that a relative increase of long-lived cells in older individuals may represent a confounding variable in age-related methylation studies [START_REF] Jaffe | Accounting for cellular heterogeneity is critical in epigenomewide association studies[END_REF] .

Discussion

The B-cell maturation process is an orchestrated program integrating internal and environmental signals to finally give rise to plasma cells and memory B cells that play an essential role in adaptive immunity. Although previous reports have studied epigenetic changes in the context of B-cell differentiation, they only studied partial DNA methylomes either of precursor or mature B cells [START_REF] Lee | A global DNA methylation and gene expression analysis of early human B-cell development reveals a demethylation signature and transcription factor network[END_REF][START_REF] Lai | DNA methylation profiling in human B cells reveals immune regulatory elements and epigenetic plasticity at Alu elements during B cell activation[END_REF][START_REF] Shaknovich | DNA methyltransferase 1 and DNA methylation patterning contribute to germinal center B-cell differentiation[END_REF][START_REF] Deaton | Cell type-specific DNA methylation at intragenic CpG islands in the immune system[END_REF] . With the exception of few cell subpopulations such as transitional B cells, CD5+ B cells and splenic marginal zone B cells, our study comprises all major B-cell differentiation stages and represents the first whole-genome epigenetic characterization of a complete human cell lineage from progenitor to terminally-differentiated cells. The comprehensive nature of our study has allowed us to provide epigenetic insights into different scientific fields and offers a resource for researchers working in different areas of cell differentiation, B-cell biology and related diseases, cancer and aging, both at single gene and genome-wide levels.

Our study points to a massive perturbation of the DNA methylome during B-cell differentiation, affecting 30% of all autosomal CpG sites. These changes follow an accumulative pattern in which each B-cell maturation stage, although characterized by a particular signature, keeps an epigenetic memory of past differentiation stages. In contrast to other reports in hematopoietic precursors [START_REF] Hodges | Directional DNA methylation changes and complex intermediate states accompany lineage specificity in the adult hematopoietic compartment[END_REF][START_REF] Jeong | Large conserved domains of low DNA methylation maintained by Dnmt3a[END_REF] , we did observe non-CpG methylation in HPCs, which virtually disappeared upon B-cell commitment in regions lacking simultaneous CpG demethylation.

Precursor B cells showed relatively small losses of CpG methylation, which mostly affected enhancers containing binding sites of B cell-specific TFs. The functional link between TF binding, CpG demethylation and enhancer activation has been recently analyzed during stem cell differentiation [START_REF] Tsankov | Transcription factor binding dynamics during human ES cell differentiation[END_REF][START_REF] Stadler | DNA-binding factors shape the mouse methylome at distal regulatory regions[END_REF] as well as in hematopoietic cells [START_REF] Schmidl | Lineage-specific DNA methylation in T cells correlates with histone methylation and enhancer activity[END_REF][START_REF] Tagoh | Dynamic reorganization of chromatin structure and selective DNA demethylation prior to stable enhancer complex formation during differentiation of primary hematopoietic cells in vitro[END_REF] , non-hematopoietic cells [START_REF] Wiench | DNA methylation status predicts cell type-specific enhancer activity[END_REF] and cancer [START_REF] Agirre | Whole-epigenome analysis in multiple myeloma reveals DNA hypermethylation of B cell-specific enhancers[END_REF][START_REF] Taberlay | Reconfiguration of nucleosome-depleted regions at distal regulatory elements accompanies DNA methylation of enhancers and insulators in cancer[END_REF] .

Interestingly, more than half of all enhancers defined in immortalized mature B cells show dynamic DNA methylation throughout the B-cell differentiation process, and 38% of all genes with dynamically-methylated enhancers are included in a regulatory network associated with human B cells [START_REF] Lefebvre | A human B-cell interactome identifies MYB and FOXM1 as master regulators of proliferation in germinal centers[END_REF] . Although mainstream research on DNA methylation still remains centered on promoter regions, our results imply that DNA methylation changes in enhancers seems to be more closely related to cell specification and maturation [START_REF] Ziller | Charting a dynamic DNA methylation landscape of the human genome[END_REF][START_REF] Gifford | Transcriptional and epigenetic dynamics during specification of human embryonic stem cells[END_REF][START_REF] Schlesinger | De novo DNA demethylation and noncoding transcription define active intergenic regulatory elements[END_REF] . However, similar to other recent studies 6,10,57 , we rarely observed a direct correlation between gene expression and DNA methylation, even in regulatory elements. Our study also suggests that at later stages of B-cell differentiation (from naiBC onwards), DNA methylation changes are more guided by other mechanisms than the intrinsic program of B-cell TFs. Upon antigen encounter, the germinal center reaction is induced and, at this stage, gcBCs start experiencing a wave of global demethylation, mostly affecting late-replicating regions such as heterochromatin and DNA repeats, and local hypermethylation of polycomb-repressed regions. This finding can be partially explained by the high proliferation rate, as normal proliferative tissues tend to lose methylation at late-replicating regions [START_REF] Aran | Replication timing-related and gene body-specific methylation of active human genes[END_REF] . However, downstream B-cell subpopulations derived from gcBCs, such as non-proliferative memBCs and PCs, which recirculate through the body and reside in bone marrow, respectively, acquire additional epigenetic changes in heterochromatin and polycomb-repressed regions. We postulate that these additional changes may be related to a potential epigenetic drift in the context of longevity that may be mediated by downregulation of DNMT1 and slight upregulation of DNMT3A.

Hypomethylation of heterochromatin and hypermethylation of polycomb-repressed regions have been previously described as an epigenetic hallmark of organismal aging, cellular senescence and cancer [START_REF] Bergman | DNA methylation dynamics in health and disease[END_REF][START_REF] Lopez-Otin | The hallmarks of aging[END_REF][START_REF] Horvath | Aging effects on DNA methylation modules in human brain and blood tissue[END_REF][START_REF] Hannum | Genome-wide methylation profiles reveal quantitative views of human aging rates[END_REF][START_REF] Heyn | Distinct DNA methylomes of newborns and centenarians[END_REF][START_REF] Teschendorff | Age-dependent DNA methylation of genes that are suppressed in stem cells is a hallmark of cancer[END_REF][START_REF] Cruickshanks | Senescent cells harbour features of the cancer epigenome[END_REF][START_REF] Rakyan | Human aging-associated DNA hypermethylation occurs preferentially at bivalent chromatin domains[END_REF] . Here, we observe that this signature starts in proliferating gcBCs and becomes particularly enhanced in non-proliferative long-lived B cells. Based on our results, we hypothesize that not all the cells of the organism are subjected to epigenetic drift as an influence of time but only those with long life span. At last, one of the most relevant implications of our study is related to the field of cancer. We demonstrate that B-cell tumors and long-lived cells share in part similar DNA methylation signatures. Furthermore, comparing various B-cell tumors with their normal cellular counterparts, we observed that a large proportion of the differentially methylated sites overlap with those undergoing dynamic methylation during normal differentiation, especially with those changing in memBCs and bm-PCs. Interestingly, similarly to preBCs, hypomethylation in ALL is enriched for enhancer elements whereas in line with gcBCs and bm-PCs, hypomethylation in DLBCL and MM is mostly enriched for heterochromatin. In general, this finding suggests that the epigenetic configuration of a particular maturation stage influences the DNA methylation changes acquired during its clonal expansion and neoplastic transformation. This novel strategy of analyzing the DNA methylome of B-cell tumors in the context of the entire differentiation program may allow us to provide new insights into the role of DNA methylation in cancer. We postulate that changes shared during neoplastic transformation and normal differentiation may represent epigenetic passengers whereas those exclusively taking place in tumor cells should contain epigenetic drivers with a potential functional impact in the disease. RNA was extracted with TRIzol (Invitrogen) following the manufacturer's recommendations. RNA quality was assessed with the Agilent 2100 Bioanalyzer.

URLs. European Genome

Whole-genome bisulfite sequencing

We performed the WGBS on two independent sets of biological replicates from six B-cell differentiation stages. Briefly, genomic DNA (1-2 μg) was spiked with unmethylated λ DNA (5 ng of λ DNA per μg of genomic DNA) (Promega). The DNA was shared by sonication to 50-500 bp using a Covaris E220 and fragments of size 150-300 bp were selected using AMPure XP beads (Agencourt Bioscience Corp.). Genomic DNA libraries were constructed using the Illumina TruSeq Sample Preparation kit (Illumina Inc.) following the lllumina standard protocol: end repair was performed on the DNA fragments, an adenine was added to the 3' extremities of the fragments and Illumina TruSeq adapters were ligated at each extremity. After adaptor ligation, the DNA was treated with sodium bisulfite using the EpiTexy Bisulfite kit (Qiagen) following the manufacturer's instructions for formalin-fixed and paraffin-embedded (FFPE) tissue samples. Two rounds of bisulfite conversion were performed to assure a conversion rate of over 99%. Enrichment for adaptor-ligated DNA was carried out through 7 PCR cycles using the PfuTurboCx Hotstart DNA polymerase (Stratagene). Library quality was monitored using the Agilent 2100 BioAnalyzer (Agilent), and the concentration of viable sequencing fragments (molecules carrying adaptors at both extremities) estimated using quantitative PCR with the library quantification kit from KAPA Biosystem. Paired-end DNA sequencing (2x100bp) was then performed using the Illumina Hi-Seq 2000. Amounts of sequence reads and the proportion of aligned reads are shown in Supplementary Table 2.

Read mapping and estimation of cytosine methylation levels

Read mapping was carried out using the GEM aligner (v1.242) [START_REF] Marco-Sola | The GEM mapper: fast, accurate and versatile alignment by filtration[END_REF] against a composite reference containing two copies of the human GRCh37 reference and two copies of the NCBI viral genome database (v35). For both the human and viral references, one copy had all C bases replaced by T and the other had all G bases replaced by A. The names of the contigs in the combined reference FASTA file were modified by adding #C2T or #G2A to the end of the contig names depending on the conversion performed. Before mapping was performed the original sequence of each read was stored. The first read of each pair then had all C bases replaced by T, and the second read had all G bases replaced by A. Read mapping with GEM was performed allowing up to 4 mismatches per read from the reference. After read mapping the original sequence of each read was restored. Estimation of cytosine levels was carried out on read pairs where both members of the read mapped to the same contig with consistent orientation, and there was no other such configuration at the same or less edit distance from the reference. After mapping, we restored the original read data in preparation for the inference of genotype and methylation status. We estimated genotype and DNA methylation status simultaneously using software developed at the Centro Nacional de Análisis Genómico, taking into account the observed bases, base quality scores and the strand origin of each read pair. For each genome position, we produced estimates of the most likely genotype and the methylation proportion (for genotypes containing a C on either strand). A phred scaled likelihood ratio for the confidence in the genotype call was estimated for the called genotype at each position. For each sample, CpG sites were selected where both bases were called as homozygous CC followed by GG with a Phred score of at least 20, corresponding to an estimated genotype error level of <= 1%. Sites with >500x coverage depth were excluded to avoid centromeric/telomeric repetitive regions. A common set of called CpG sites for all analyzed samples was generated, and all subsequent analyses used this common set.

Microarray-based DNA methylation analysis with 450k arrays

We used the EZ DNA Methylation Kit (Zymo Research) for bisulfite conversion of 500 ng genomic DNA. Bisulfite-converted DNA was hybridized onto the HumanMethylation 450K BeadChip kit (Illumina) which covers 99% of RefSeq genes and 96% of CpG islands. The Infinium methylation assay was carried out as described previously [START_REF] Bibikova | High density DNA methylation array with single CpG site resolution[END_REF][START_REF] Bibikova | Genome-wide DNA methylation profiling using Infinium(R) assay[END_REF] . Data from the 450k Human Methylation Array were analyzed in R using the minfi package [START_REF] Aryee | Minfi: a flexible and comprehensive Bioconductor package for the analysis of Infinium DNA methylation microarrays[END_REF] (version: 1.6.0), available through the Bioconductor open source software. To exclude technical and biological biases that might produce false results in further analyses, we developed and optimized an analysis pipeline with several filters (i.e. CpGs with low detection p values, sex-specific and individual-specific methylation or overlapping with SNPs). Taking into account the different performance of Infinium I and Infinium II assays we used the subset-quantile within array normalization (SWAN) [START_REF] Maksimovic | SWAN: Subset-quantile within array normalization for illumina infinium HumanMethylation450 BeadChips[END_REF] that corrects for the technical differences between the Infinium I and II assay designs and produces a smoother overall beta value distribution.

Detection of non-CpG methylation and differential methylation analysis

Cytosines in non-CpG context were identified as two adjacent nucleotides where the genotype of the first nucleotide was called with high confidence as homozygous C and the second was called with high confidence as a genotype not containing a G. Non-CpG cytosines were called methylated if they had at least two non-converted reads, at least six reads informative for methylation status, and a methylation probability greater than twice its standard deviation. The significance of the change in methylation levels between samples was assessed using the numbers of converted and non-converted reads in both samples with a chi square test, or Fisher's exact test when the chi square approximation was not appropriate. The difference in methylation levels between different stages of B-cell differentiation was calculated using 16.1 million CpGs with methylation estimates in all 12 samples analyzed by WGBS. The normal approximation to the binomial was used to test for significant differences of individual CpGs between samples. As we sequenced two biological replicates per cell subpopulation, we defined consistent DNA methylation changes between two differentiation stages if in one set of samples the methylation difference was above 0.25 and in the second set was at least a 0.1 difference in the same direction (hyper or hypomethylation). CpGs with dynamic methylation were defined as with differential methylation in comparisons of adjacent stages and in the comparison between HPCs and bm-PCs. To calculate the overlap between dynamically methylated CpGs in B cells and by Ziller et al [START_REF] Ziller | Charting a dynamic DNA methylation landscape of the human genome[END_REF] , we downloaded all the DMRs with dynamic CpGs from that study (GEO number GSE46644) and determined how many of the 4.93 million dynamic CpGs in B cells are located within DMRs identified by Ziller et al. We also defined CpGs that show variable methylation levels throughout B-cell development using 450k microarray data for HPC, preB1C, preB2C, iBC, naiBCs, t-naiBCs, gcBC, t-PC and bm-PCs. From gcBCs on, B-cell differentiation is branched into memBCs or PCs and therefore memBCs were not included in this linear analysis. We performed pair-wise comparisons between all these subsets of B-cell differentiation. We defined as "dynamic" those CpGs that presented a mean DNA methylation difference above 0.25 in at least one comparison (FDR<0.1, Wilcoxon test). Furthermore, we detected CpGs whose methylation shows a similar modulation pattern during the entire B-cell differentiation process and we group them into distinct modules. To allow for the identification of gradual changes (but consistent in different replicates) throughout the differentiation program, we applied a mean methylation difference between adjacent subpopulations of 0.1. The modules were grouped into four groups depending on their methylation tendency during differentiation (decrease, increase, decreaseincrease and increase-decrease). Differentially methylated sites between specific B-cell tumor entities and their normal cellular counterpart were identified by a mean DNA methylation difference above 0.25 and an FDR < 0.05.

Bisulfite pyrosequencing studies

Validation of DNA methylation levels generated by WGBS and 450k microarrays was validated by bisulfite pyrosequencing (BPS). Briefly, 500 ng of genomic DNA was bisulfite converted using the EpiTect 96 Bisulfite kit or the EpiTect Plus Bisulfite Conversion Kit (Qiagen) according to manufacturer's instructions. PCR amplification of the bisulfite treated DNA was done using the specific primers for each of the selected CpGs and non-CpGs (Supplementary Table 3). These primers were determined with the PyroMark Assay Design software (Qiagen). BPS and DNA methylation data analysis were performed with the pyrosequencer PyroMark Q96 ID and PyroMark CpG software (Qiagen).

Genomic and functional annotation of CpG sites

Both WGBS and 450k microarray data were annotated using the UCSC Genome Browser database (hg19). For the location relative to a gene, we used these categories: TSS 1500 (from 201 to 1,500 bp upstream of the transcriptional start site (TSS)), TSS 200 (from 1 to 200 bp upstream of the TSS), 5' UTR, first exon, exon (all exons excluding exon 1), intron, 3' UTR and intergenic regions. Owing to the presence of alternative transcription start sites and regions containing more than one gene, some of the CpGs were assigned multiple annotations. For the location relative to a CpG island (CGI), we used these groups: within CGI, in CGI shore (0-2 kb from the CGI edge), in CGI shelf (>2 kb to 4 kb from the CGI edge) and outside CGI. We also annotated all CpG probes using a recent categorization of chromatin and transcriptional states from the lymphoblastoid B cell line GM12878 [START_REF] Ernst | Mapping and analysis of chromatin state dynamics in nine human cell types[END_REF] (ChromHMM track of the UCSC Genome Browser), which has a DNA methylome similar to memory B cells and plasma cells (Supplementary Fig. 14). Regions with chromatin states 1-3 (active, weak and poised promoter) were considered as "Promoter regions", states 4-7 (strong and weak enhancer) as "Enhancer regions", state 8 as "Insulator", states 9 "Transcriptional transition", state 10 "Transcriptional elongation", state 11 as "Weak transcription", state 12 as "Polycombrepressed regions" and state 13 as "Heterochromatin (nuclear lamina)". Replication timing in GM12878 data was obtained from the UW Repli-seq track of the UCSC Genome Browser. Replication timing values for all sites from the background was divided into 3 bins: early-, mid-and late-replicating regions. Only early-and late-replicating regions were used for the analysis. The annotation of repeat elements was done based on the RepeatMasker Annotation, available at the UCSC Genome Browser.

  M.K., A.C.Q., G.Cas., R.B. and G.Clo. analyzed DNA methylation and gene expression arrays. L.A., J.B. and M.G. performed WGBS library preparation and sequencing. A.M., S.H., R.P.S., E.R., A.E. and M.D.-F. processed and analyzed WGBS data. M.K., N.V.-D. and R.V.-B. performed validation experiments. M.K., G.Cas., S.E., V.P., D.Rico. and A.V functionally characterized dynamically methylated genes. D.Rich., L.C., A.D. and P.F. were in charge of data management. I.G.G. and H.G.S. coordinated sequencing efforts and performed primary data analysis. H.G.S., R.S., R.K. and E.C participated in the study design and data interpretation. J.I.M.-S. conceived the study, J.I.M.-S. led the experiments and wrote the paper with predominant assistance of M.K. and R.B.

FIGURE LEGENDS Figure 1 .

 LEGENDS1 FIGURE LEGENDS

Figure 2 .

 2 Figure 2. Non-CpG methylation detected during B-cell differentiation. (a) Browser

Figure 3 .

 3 Figure 3. Dynamic DNA methylation during B-cell differentiation. (a) Differentially methylated

Figure 4 .

 4 Figure 4. Association between DNA methylation and gene expression in different chromatin

Figure 5 .

 5 Figure 5. DNA methylation changes during B-cell differentiation in the context of cancer and

Figure 6 .

 6 Figure 6. DNA methylation changes in various B-cell neoplasms as compared to their normal

  -Phenome Archive, http://www.ebi.ac.uk/ega/studies/; Blueprint project, http://www.blueprint-epigenome.eu; The R Project for Statistical Computing, http://www.r-project.org/; UCSC Genome Browser (http://genome.ucsc.edu/).

	Accession codes. WGBS data have been deposited in the European Genome-Phenome Archive
	(EGA) under the accession numbers EGAD00001001304 and EGAS00001000272. DNA
	methylation and gene expression microarray data are available from the EGA under accession
	numbers EGAS00001001196 and EGAS00001001197, respectively. Additional published data
	from precursor B cells are available through Gene Expression Omnibus under accession
	number GSE45461.

ACKNOWLEDGEMENTS

We thank Carlos López-Otín for critical reading of this manuscript, Marc Dabad Castellà for his assistance with the WGBS data analysis and Miguel A. Peinado for providing RNA from HCT116 DKO cells. This work was funded by the European Union's Seventh Framework Programme through the Blueprint Consortium (grant agreement 282510), the Spanish Ministry of Economy and Competitivity (MINECO) (project SAF2009-08663). Methylation microarrays were outsourced to the Spanish "Centro Nacional de Genotipado" (CEGEN-ISCIII). M.K. is supported by the AGAUR (Generalitat de Catalunya), EC is an Academia Researcher of the Institució Catalana de Recerca i Estudis Avançats, and J.I.M.-S. is a Ramon y Cajal researcher of the MINECO.

AUTHOR CONTRIBUTIONS

M.K., A.C.Q., N.R., M.P., X.A., F.P., D.A., B.P., G.Car., T.F., M.O.M., M.E.F., S.-T.L., J.L.W., provided samples from healthy donors and/or purified B-cell subpopulations.

COMPETING FINANCIAL INTEREST

The authors declare no competing financial interests.

ONLINE METHODS

Isolation of B-cell subpopulations

Precursor B cells were isolated from fetal bone marrow (22 week fetuses) using flow cytometry sorting. Early progenitors were isolated based on high levels of CD34 protein expression (CD34hi) and a lack of expression of the B-cell marker CD19. This population, designated as uncommitted hematopoietic progenitors (HPC), contains predominantly multipotent progenitors before lineage commitment and also common lymphoid progenitors and hematopoietic stem cells. B-cell-committed progenitors were isolated based on their expression of CD19 and CD34 (CD19+/CD34+), which were predominantly pre-B-I cells (preB1Cs). Two immature B-cell populations expressing CD19 and lacking CD34 were isolated and differentiated based on surface IgM (sIgM) expression: pre-B-II cells (preB2Cs) that express sIgM-/CD19+ and immature B cells that express sIgM+/CD19+ (iBCs). DNA methylation and gene expression data from these 4 subpopulations have been previously published [START_REF] Lee | A global DNA methylation and gene expression analysis of early human B-cell development reveals a demethylation signature and transcription factor network[END_REF] . Peripheral blood B-cell subpopulations, i.e. naive B cell (naiBCs) and memory B cells (memBCs) were obtained from buffy coats from healthy adult donors of age ranging between 28 and 66 years. After Ficoll-Isopaque density centrifugation CD19+ B-cells were isolated by positive magnetic cell separation by using AutoMACS system (Milteny Biotec, Auburn, CA). CD19+ cells were labeled with anti-CD27, anti-IgD, anti-IgM, anti-IgG and anti-IgA for 15 min at room temperature in staining buffer (PBS with 0.5% BSA). naiBCs (CD19+/CD27-/IgD+) and memBCs (CD19+/CD27+/IgA+ or IgG+) were obtained by FACS sorting on FACSAriaII (BD Biosciences). Plasma cells (t-PCs), germinal centre B cells (gcBCs) and naive B cells (t-naiBCs) were separated from tonsils of children undergoing tonsillectomies (ranging in age between 2 and 13), obtained from the Clinica Universidad de Navarra (Pamplona, Spain) or Children's Clinique La Sagesse (Rennes, France). Tonsils were minced extensively and after Ficoll-Isopaque density centrifugation, enrichment of B cells was performed with the AutoMACS system either by positive selection of CD19+ cells or by the B Cell Isolation Kit II (Milteny Biotec, Auburn, CA). t-PCs (CD20med/CD38high), gcBCs (CD20high/CD38med) and t-naiBCs (CD20+/CD23+) were separated by FACS sorting. In part, t-naiBCs were also selected using a slightly different marker combination (i.e. CD19+/CD27-/IgD+ or IgD+/CD38low/CD27-). gcBCs were also selected by this marker combination: IgD-/CD38hi/CD10+/CXCR4+. t-naiBCs and gcBCs isolated with different markers identified the same cell subpopulations, as evidenced by the fact that each subpopulation shows homogeneous DNA methylation and transcriptional profiles (Supplementary Fig. 22). Plasma cells from bone marrow (bm-PCs) were selected from healthy donors ranging from 20 to 30 years. After density gradient centrifugation, we performed a selective depletion of CD3+, CD14+ and CD15+ cells by immunomagnetic selection (Miltenyi Biotec, Germany), followed by flow cytometry cell sorting of CD45+/CD138+/CD38+ cells using a FACSAriaII (BD Biosciences) device. The purity of each of the isolated B-cell subpopulations exceeded 90% in all samples. DNA was extracted from purified samples by using a Qiagen kit (QIAamp DNA Mini Kit), following manufacturer's instructions, and was quantified using a Nanodrop ND-100 spectrophotometer. DNA samples for WGBS, 450k array and BPS experiments were derived from individual donors with the exception of those from bm-PCs, which were pooled from 4 different donors. Total

B-cell network analysis

We used the B-cell specific functional interaction network of Lefebvre et al., [START_REF] Lefebvre | A human B-cell interactome identifies MYB and FOXM1 as master regulators of proliferation in germinal centers[END_REF] containing 5,748 nodes (genes) and 64,600 unique edges (interactions) based on Entrez gene identifiers. We selected the 5,668 dynamic enhancer genes and mapped them to Entrez gene identifiers resulting in 5,658 unique Entrez gene identifiers. 2,154 of these Entrez genes are contained in the B-cell network and 1,993 of them are directly connected in the network by 11,741 edges. This subnetwork of 1,993 nodes and 11,741 edges was investigated further. We identified 9 communities in the subnetwork by using Gephi 66 and Louvain's method [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF] .

Chromatin immunoprecipitation coupled with bisulfite sequencing B cells were cross-linked with 1% formaldehyde during 8 min (RT) before FACS separation of memBCs. Chromatin preparation and ChIP were performed according to the BLUEPRINT Histone ChIP protocol with a anti-H3K27me3 antibody (C15410195, Diagenode). For wholegenome bisulfite library construction and sequencing the immunoprecipitated DNA (50 ng) was sheared on a Covaris™ E220 in order to reach the fragment size of 50-500 bp and size selected 150 to 500 bp fragments using AMPure XP beads (Agencourt Bioscience Corp.). 500 ng of unmethylated λ DNA (Promega) was treated in parallel on Covaris™ E220 and also size selected by AMPure XP beads to reach the same fragment sizes as the DNA sample. The unmethylated λ DNA was spiked into fragmented and size selected immunoprecipitated DNA (5 ng of λ DNA per 1 μg of DNA) and TruSeq Sample Preparation kit (Illumina Inc.) was used to prepare the Illumina library in order to add platform specific adaptors. After adaptor ligation 450 ng of fragmented and size selected unmethylated λ DNA was added to the library. Two rounds of bisulfite conversion was performed to reach > 99% conversion following the manufacturer's instructions for formalin-fixed and paraffin-embedded (FFPE) tissue samples (EpiTect Bisulfite kit; Qiagen). Adaptor-ligated DNA was enriched through ten cycles of PCR using the KAPA HiFi Uracil+ polymerase (Kapa Biosystems). The library was run in a fraction of a lane of HiSeq2000 flow cell (to reach 35 million paired end reads) with read length of 2x100 bp, according to standard Illumina operation procedures. Primary data analysis was carried out with the standard Illumina pipeline. The sequence reads were passed through the same read mapping and genotype/methylation calling pipeline as the conventional WGBS samples. In addition, the aligned reads were analyzed with the NucHunter package [START_REF] Mammana | Inferring nucleosome positions with their histone mark annotation from ChIP data[END_REF] , to provide predictions of the positions of H3K27me3containing nucleosomes. The average methylation of cytosines at all distances from 0-500 bp of each peak was calculated, taking into account the strand of the cytosine. The same analysis was performed on WGBS data from the same cell type (memBC), to allow comparison of the enriched (ChIP-BS) and non-enriched (WGBS) results. This analysis was repeated using only predicted nucleosome peaks that fell within predicted polycomb repressed regions from the Broad ChromHMM analysis of the ENCODE cell line GM12878.

Gene ontology analysis

The GOstat package [START_REF] Falcon | Using GOstats to test gene lists for GO term association[END_REF] available through Bioconductor was used to determine the enrichment of individual ontology terms in the different methylation modules as compared to all the genes analyzed in the 450k array. In Supplementary Data 4, the top 20 significant terms for each module (P < 0.001) are shown.

Analysis of transcription factor binding sites

Transcription factor binding site (TFBS) information was obtained by ChIP-seq data in the ENCODE project and available at the UCSC Genome Browser. A total of 79 TFBS were used for the analysis. The relative enrichment of each TFBS was calculated as compared to the background. A fisher exact test was used and an Odds Ratio and a P value were assigned to each comparison.

Gene expression analyses

RNA samples of HPCs, preB1Cs, preB2Cs, naiBCs, t-naiBCs, gcBCs, memBCs, t-PCs and bm-PCs were hybridized to Affymetrix Human Genome U219 Array Plates according to Affymetrix standard protocols. The analysis of scanned images for each probe set of the array was obtained with GeneChip Operating Software (GCOS, Affymetrix). Raw CEL files were processed and normalized with the robust multichip average (RMA) algorithm using R statistical software in conjunction with the Affy library 70 available through Bioconductor. The GeneChip Human Gene 1.0 ST Array data of progenitor B cells were downloaded through Gene Expression Omnibus under accession number GSE45461 and normalized using RMA (these data were only used for Fig. 4b,c). To evaluate the gene expression variability among B-cell subpopulations, we calculated the inter-cell subtype standard deviation per each of the Affymetrix tags. A global measure of the variability for a particular set of genes (within modules and chromatin states) was then calculated as the average of all the standard deviations. Differential expression between t-naiBCs and naiBCs was calculated using Limma library available through Bioconductor, using the criteria of FC > 1 (log2) between two groups and adjusted P < 0.05. Reverse transcription and quantitative PCR for DNMTs expression was done as follows: cDNAs from 100 ng RNA samples were synthesized in 20 µl reaction mix using oligo dT primers and SuperScript III enzyme according to the manufacturer recommendations (Invitrogen). Primers for DNMTs were taken from Fang et al. [START_REF] Fang | Epigenetic changes mediated by microRNA miR29 activate cyclooxygenase 2 and lambda-1 interferon production during viral infection[END_REF] whereas those for the housekeeping EEF2 gene are provided in Supplementary Table 3. PCR amplification was carried out with 1 µl of the 1:2 diluted reverse transcription sample with Power SYBR green PCR master mix according to manufacturer recommendations (Applied Biosystems). PCR reactions were run in triplicates on a StepOne System (Applied Biosystems).