Agathe Duchateau

Névé Auperin

Federico Miozzo

Anne Le Mouel

Olivier Kirsh

Myriame Mohamed

Sascha Ott

Délara Saberan-Djoneidi

Valérie Mezger
email: valerie.mezger@u-paris.fr

Immediate perturbation of DNA methylation upon acute prenatal alcohol exposure in the mouse developing brain cortex

Keywords: Prenatal alcohol exposure, Fetal alcohol spectrum disorder (FASD), immediate alcohol-induced modifications in DNA methylation, differentially methylated regions, differentially expressed genes, imprinted genes, protocadherins, neurodevelopment, mice embryonic cortices, epigenetics, methylome capture

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

The prenatal developing brain is particularly vulnerable to detrimental events that generate neurodevelopmental defects and potentially have long-term consequences in the adulthood [START_REF] Bale | Early Life Programming and Neurodevelopmental Disorders[END_REF][START_REF] Schang | Epigenome and transcriptome landscapes highlight dual roles of proinflammatory players in a perinatal model of white matter injury[END_REF]. Among a diversity of adverse in utero stresses, prenatal alcohol exposure (PAE) is a leading cause of non-genetic mental retardation in the Western world [START_REF] Popova | Health Care Burden and Cost Associated with Fetal Alcohol Syndrome: Based on Official Canadian Data[END_REF][START_REF] Popova | Comorbidity of fetal alcohol spectrum disorder: a systematic review and meta-analysis[END_REF]. Depending on many parameters, such as the timing of exposure, the drinking pattern (chronic or acute) and the amount of alcohol consumed, PAE gives rise to a wide range of neurodevelopmental defects, referred to as part of the clinical picture of fetal alcohol spectrum disorders (FASD), whose prevalence is estimated around 9 for 1000 live births [START_REF] Burd | Recognition and management of fetal alcohol syndrome[END_REF][START_REF] Jones | Recognition of the fetal alcohol syndrome in early infancy[END_REF][START_REF] Kleiber | Neurodevelopmental alcohol exposure elicits long-term changes to gene expression that alter distinct molecular pathways dependent on timing of exposure[END_REF][START_REF] Lemoine | Les enfants de parents alcooliques[END_REF][START_REF] Mattson | Fetal Alcohol Spectrum Disorders: Neuropsychological and Behavioral Features[END_REF][START_REF] Popova | Health Care Burden and Cost Associated with Fetal Alcohol Syndrome: Based on Official Canadian Data[END_REF][START_REF] Popova | Comorbidity of fetal alcohol spectrum disorder: a systematic review and meta-analysis[END_REF]. These defects lead to impairment in cognition, behavior, executive function, attention (linked or not to hyperactivity), learning, judgment and social adaptation [START_REF] Gibbard | The neuropsychological implications of prenatal alcohol exposure[END_REF]. In addition to these primary defects, individuals diagnosed with FASD are at high risk for neuropsychiatric disorders, including anxiety disorders, depression, and addiction in their adulthood [START_REF] Gibbard | The neuropsychological implications of prenatal alcohol exposure[END_REF] reviewed in [START_REF] Kodituwakku | Defining the behavioral phenotype in children with fetal alcohol spectrum disorders: a review[END_REF]O'Connor and[START_REF] O'connor | Psychiatric conditions associated with prenatal alcohol exposure[END_REF], in line with the DOHaD concept (Developmental Origins of Health and Disease; Schang et al., 2018b). The most severe form of FASD, called fetal alcohol syndrome (FAS), is defined by characteristic facial dysmorphology, growth retardation and severe macroscopic structural abnormalities of the central nervous system (CNS; [START_REF] Jones | Recognition of the fetal alcohol syndrome in early infancy[END_REF]. Early diagnosis of PAE is necessary to enable efficient intervention [START_REF] Burd | Recognition and management of fetal alcohol syndrome[END_REF][START_REF] Paley | Behavioral interventions for children and adolescents with fetal alcohol spectrum disorders[END_REF]. Visible and well-characterized phenotypic defects observed for FAS newborns facilitated the early diagnosis of these individuals. In contrast, diagnosis of FASD children is more challenging, because the history of exposure is often unknown and clinical phenotypes are more cryptic, yet crucial since the functioning of the central nervous system damages is compromised [START_REF] Gibbard | The neuropsychological implications of prenatal alcohol exposure[END_REF]Popova et al., 2015). Indeed, PAE affects neurodevelopment at any stage, therefore impacting all neurodevelopmental processes, such as proliferation of neural progenitors, migration of young post-mitotic neurons and their differentiation, neuronal survival, synaptogenesis, neurotransmission and neuronal plasticity [START_REF] Kleiber | Neurodevelopmental alcohol exposure elicits long-term changes to gene expression that alter distinct molecular pathways dependent on timing of exposure[END_REF][START_REF] El Fatimy | Heat shock factor 2 is a stress-responsive mediator of neuronal migration defects in models of fetal alcohol syndrome[END_REF][START_REF] Guerri | Foetal Alcohol Spectrum Disorders and Alterations in Brain and Behaviour[END_REF][START_REF] Hashimoto-Torii | Roles of Heat Shock Factor 1 in Neuronal Response to Fetal Environmental Risks and Its Relevance to Brain Disorders[END_REF][START_REF] Ishii | Variations in brain defects result from cellular mosaicism in the activation of heat shock signalling[END_REF].

Although defects caused by prenatal alcohol exposure are well identified, the exact molecular mechanisms underlying these alterations and their persistence are still unclear. However, PAE is known to long-lasting impact on gene expression levels. These transcriptional disturbances, observed in the adult brain, occur in a manner dependent on the developmental stage at which the exposure has occurred. Indeed, transcriptional changes affect genes falling into distinct GO (Gene Ontology) categories, depending on the time-window of exposure in preclinical models, (equivalent to the first to third trimester of pregnancy in women) [START_REF] Kleiber | Neurodevelopmental alcohol exposure elicits long-term changes to gene expression that alter distinct molecular pathways dependent on timing of exposure[END_REF][START_REF] Kleiber | Long-term genomic and epigenomic dysregulation as a consequence of prenatal alcohol exposure: a model for fetal alcohol spectrum disorders[END_REF]. These long-lasting disturbances therefore affect processes governed by gene networks that are active during the gestation, at the time of exposure. This "memory" of the exposure could occur via Epigenetic mechanisms.

Remarkably, brain development is tightly controlled by epigenetic mechanisms, as underlined by the impact of mutations or variants in genes encoding epigenetic actors in the emergence of neurodevelopmental and neuropsychiatric disorders, such as Rett syndrome, Rubinstein-Taybi syndrome and autism spectrum disorders [START_REF] Bourgeron | From the genetic architecture to synaptic plasticity in autism spectrum disorder[END_REF][START_REF] Gräff | Epigenetic Regulation of Gene Expression in Physiological and Pathological Brain Processes[END_REF][START_REF] Lasalle | Epigenetic layers and players underlying neurodevelopment[END_REF].

The long-lasting and stage-specific transcriptomic alterations, described above and observed in the PAE-exposed brain, might be underlined by modifications in the epigenetic landscape, including perturbation of DNA methylation profile. Indeed, starting from seminal works on rodent models (e.g. [START_REF] Haycock | Exposure of Mouse Embryos to Ethanol During Preimplantation Development: Effect on DNA Methylation in the H19 Imprinting Control Region[END_REF][START_REF] Kaminen-Ahola | Maternal Ethanol Consumption Alters the Epigenotype and the Phenotype of Offspring in a Mouse Model[END_REF], a number of studies have identified disturbances in DNA methylation on candidate genes or loci, or in a genome-wide manner, in the mouse adult brain that was exposed to alcohol prenatally [START_REF] Kleiber | Neurodevelopmental alcohol exposure elicits long-term changes to gene expression that alter distinct molecular pathways dependent on timing of exposure[END_REF][START_REF] Laufer | Long-lasting alterations to DNA methylation and ncRNAs could underlie the effects of fetal alcohol exposure in mice[END_REF] reviewed in [START_REF] Lussier | Epigenetics studies of fetal alcohol spectrum disorder: where are we now?[END_REF]. Alterations in DNA methylation profile observed in mouse brain in response to PAE has been corroborated in peripheral tissues (cheek swabs) in cohorts of FASD children [START_REF] Laufer | Associative DNA methylation changes in children with prenatal alcohol exposure[END_REF][START_REF] Lussier | DNA methylation as a predictor of fetal alcohol spectrum disorder[END_REF]. PAE during the preimplantation period has also been associated to DNA methylation perturbations observed more than one week after and associated to neurodevelopmental defects [START_REF] Legault | Pre-implantation alcohol exposure induces lasting sex-specific DNA methylation programming errors in the developing forebrain[END_REF].

Strikingly, the short-term impacts of PAE on the genome-wide architecture of DNA methylation in the developing brain have been understudied and there is a need to investigate whether the aberrant DNA methylation events are the direct and immediate results of alcohol exposure. First, this is of importance for the field, because DNA methylation disturbances observed after birth or in the adulthood could (i) result from immediate PAE-dependent DNA methylation changes that would persist throughout life or (ii) be due to PAE-induced brain dysfunction. Indeed, in the last case, since neuronal activity can reshape DNA methylation throughout lifetime [START_REF] Guo | Neuronal activity modifies the DNA methylation landscape in the adult brain[END_REF], the remodeling of the DNA methylome could represent secondary effects of PAE, resulting from the PAE-induced alteration of neuronal activity in the infant or adult brain. Second, this is also a question of medical relevance. In many cases, the history of exposure to alcohol is unknown and FASD children are often diagnosed late, which compromised early intervention [START_REF] Gibbard | The neuropsychological implications of prenatal alcohol exposure[END_REF]. The search for accurate and relevant molecular biomarkers, especially for biomarkers of exposure, is thus necessary. Attempts to identify biomarkers of exposure linked to DNA methylation perturbations have been given special attention in the field. Ideally, such DNA methylation changes should occur quickly after exposure. This has remained to be determined.

We thus asked whether immediate changes in DNA methylation could be detected across the genome, in the developing mouse cerebral cortex, using a tailor-made methylome capture approach.

To capture immediate events of DNA methylation alteration and since binge drinking has become an increasing mode of alcohol consumption in young women (Dwyer-Lindgren et al., 2015;Tan et al., 2015;Popova et al., 2018), we chose a binge-drinking-like mouse model of PAE. We observe hundreds of differentially methylated regions (DMRs), as early as two hours after the last alcohol injection. We show that DMRs are significatively associated to genes or loci which, in physiological conditions, are dynamically regulated at the time of exposure, in terms of chromatin accessibility or expression, and which are linked to biological processes of importance for brain development and functions.

Moreover, DMRs are overrepresented in two families of mono-allelically expressed genes of importance for neurodevelopment -imprinted genes and genes of clustered protocadherins -and correlate with immediate modifications in their transcript levels. Our results therefore show that some genomic regions are the target of differential methylation events, very rapidly after exposure, that are associated with concomitant changes in gene expression. Interestingly, according to published results, DMRs associated with imprinted genes and protocadherins have been observed at temporal distance of PAE, in the mouse and/or human postnatal and adult brain, this suggest that these loci might represent immediate targets of PAE with long-lived DNA methylation perturbation and provide valuable and meaningful biomarkers of exposure in the future.

RESULTS

Working model: a methylome capture in a murine binge drinking model

To identify DNA methylation changes that would occur immediately after an acute prenatal alcohol exposure, we used a murine model of binge drinking at a developmental stage equivalent to the second trimester of pregnancy in human: mouse embryonic cortices were exposed to PAE, in utero, at embryonic days E15 and E16 (Figure 1A, see Materials and Methods). We chose this binge-drinking model to be as close as possible to a human mode of acute consumption. The quantity of injected ethanol corresponded to a dose known to induce brain defects that mimic FAS, in rodent pups [START_REF] Carloni | Caspase-3 and calpain activities after acute and repeated ethanol administration during the rat brain growth spurt[END_REF][START_REF] Ikonomidou | Ethanol-induced apoptotic neurodegeneration and fetal alcohol syndrome[END_REF][START_REF] Olney | Ethanol-induced apoptotic neurodegeneration in the developing C57BL/6 mouse brain[END_REF]. We searched for methylation modifications immediately occurring after PAE, using a capture of the DNA methylome (EpiCapture technology; see Materials and Methods), which allowed us to study, with a good resolution, a repertoire of 58,611 selected DNA regions corresponding to about 81.3 Mb. In particular, the capture design was composed of more than 75% of mouse promoters and all the enhancers that are active in the brain (Figure 1B, Supplementary Figure 1A; see Materials and Methods for details).

Quality control of the methylome capture

This methylome capture approach used sodium bisulfite (BS) conversion, to distinguish between methylated and unmethylated cytosines. First, we verified that the BS conversion had correctly occurred for all samples [START_REF] Clark | DNA methylation: bisulphite modification and analysis[END_REF][START_REF] Krueger | DNA methylome analysis using short bisulfite sequencing data[END_REF]. Abnormal per-base-sequence-content plot obtained with FASTQC was observed for each sample, as expected for DNA sequences converted with BS (Supplementary Figure 2A,B). Moreover, BS conversion rate of all samples was quantified, using spike-in DNA. Conversion rate was high (>97,8%) for all the samples, which means that no conversion biases were present in our data analyses (Supplementary Table 1).

To minimize inter-individual variability in the same group, we mixed half-cortices from 4 embryos from 4 different litters (2 left and 2 right half-cortices; Supplementary Figure 1B). In line with our sampling choice, we observed that intra-group variability was limited, according to hierarchical clustering analysis and Pearson correlation between samples (Supplementary Figure 2C,D). However, triplicates from the EtOH-treated group showed more differences between each other than did triplicates from control group (PCA analysis and hierarchical clustering analysis, Supplementary Figure 2C,E), suggesting that, even if samples were globally affected in the same manner by in utero alcohol exposure, some subtle differences in term of CpG methylation level could be observed. In addition, CpG base Pearson correlation values indicated that inter-group variability was also limited. This suggests that PAE does not majorly reshape the DNA methylome.

Using the bioinformatics workflow described in Supplementary Figure 2A, to compare control and alcohol-treated groups, we only identified three hyper-and one hypomethylated isolated (individual) CpGs, all having a high methylation differential percent (> 35%; Figure 2A; MethylKit analysis on normalised and filtered datasets, meth. diff < 5%, qval. < 0.05). Two of these differentially methylated cytosines (DMC) are located in intergenic regions (on chromosomes 5 and 10). One hypermethylated CpG is in the Kalrn gene, encoding the Kalirin RhoGEF kinase that is involved in various mechanisms, such as neuronal shape regulation, growth and plasticity. The third hypermethylated CpG is associated to Tiam2 gene, encoding a RAC1-specific guanine nucleotide exchange factor.

Definition of DNA methylated regions for capture analysis

Although the methylation of an individual CpG can affect gene expression [START_REF] Xu | Pioneer factor interactions and unmethylated CpG dinucleotides mark silent tissue-specific enhancers in embryonic stem cells[END_REF], DMRs are considered as more relevant in terms of impact on gene expression. In line with this consideration, studies of differentially methylated regions (DMRs) have been more frequently performed [START_REF] Bock | Analysing and interpreting DNA methylation data[END_REF]. We therefore pursued our analyses on the identification of DMRs. Bioinformatic analysis of this methylome capture required special attention, especially to define potential differentially methylated regions (DMRs) between control and alcohol-treated groups.

Indeed, to define such genomic regions, existing bioinformatics tools were often design for analyses considering the whole genome (or well-characterised arrays), as a reference. This kind of approach is not adapted to a customized capture based on a repertoire of selected regions. For this reason, we developed our own R function to define DMRs between control and alcohol-treated groups, which combined neighbouring CpGs that share similar differential methylation states, either all showing loss of methylation (hypo CpGs) or all showing gain of methylation (hyper CpGs). More precisely, we defined DMRs on the following basis (see Materials and Methods; Supplementary Figure 3): the association of at least 5 CpG having a same methylation state, in a region of maximum 2000 bases, with a maximum distance of 100 bases between two successive selected CpG, according to a p-value threshold (Supplementary Table 2). For this, a pre-selection of relevant CpGs was done, according to their p-value (CpGs with p-value < 0.07). This p-value threshold was determined using comparisons of real and random datasets (see below, Supplementary Figure 4, Supplementary 2B, Supplementary Dataset 1). Among these, DMRs showing gain of methylation were predominant, compared to DMRs corresponding to loss of methylation ( 257"hyper-DMRs"; 175 "hypo-DMRs"; Figure 2B). In either hyper-or hypo-DMRs, the median differential levels of DNA methylation in the majority of DMRs reached more than 10% (Figure 2B).

Even if there was a limited number of regions altered in their DNA-methylation levels, at least in the repertoire of the captured loci, our results suggest that DNA methylation is rapidly and locally redistributed in the cortices exposed to PAE. We performed hypergeometric tests either on all DMRs, or by separating hypo-DMRs and hyper-DMRs and observed that genomic regions were not randomly altered by PAE in their DNA methylation status. Indeed, brain active enhancers, characterized by the H3K27ac histone mark, were significantly over-represented among the DMRs identified in the capture, (345 on 432 DMRs i.e., 79.86% of DMRs; hypergeometric test, p-val < 0.05. Figure 2C,D; Figure 3). In contrast, promoters and other regions of interest, which represent respectively 134 and 35 DMRs, were not significantly affected by the binge drinking stress (Figure 2B,C).

To identify potential binding sites of transcription factors within the DNA sequences of enhancers located in DMRs, we explored the presence of DNA sequence motifs using findMotifsGenome (HOMER). We could not observe any over-representation of specific motifs. This result suggests that there are no sub-categories of enhancers (i.e. characterized by specific transcription factors), that would be more severely affected by PAE than others. Since these enhancers belong to an active repertoire in the mouse brain, we then performed a Gene Ontology (GO) analysis on all the enhancers containing DMRs, in order to unveil specific Biological Processes potentially rapidly altered by PAE. We observed that a significant number of genes associated to these enhancers played a role in neurogenesis, generation of neurons, cell development, nervous system development, and neuron differentiation (Figure 3A; Table 1). When considered separately, enhancers with a decrease in methylation levels upon PAE were not significantly associated to a specific Biological Process. In contrast, hypermethylated enhancers were predominantly associated with genes that are involved in the above-cited Biological Processes (Figure 3B,Table 2).

Regions showing chromatin remodelling during physiological brain development are significantly associated with PAE-induced DNA methylation events

The striking enrichment of DMRs in enhancer regions, which were specifically active in the brain, let us think that these DNA methylation perturbations could occur in regions whose chromatin accessibility was undergoing changes at cortical developmental stages around the PAE exposure. We thus both investigated the changes of global chromatin accessibility and gene expression dynamics throughout the physiological mouse forebrain development around E14.5 to E16.5. For this, we used a bioinformatics workflow, which allows for pairwise comparison of successive developmental stages, to analyse public ENCODE ATAC-seq and RNA-seq datasets (Supplementary Figure 2A, andMaterials and Methods, Supplementary Data notebook 2 and 3, Supplementary Tables 4 and5). Our results are described in Supplementary Dataset 1. The number of identified regions or genes in each dataset is indicated in Supplementary Figure 5A (for quality control of these data, see Supplementary Figure 5B,C).

We observed DMR enrichment in regions whose chromatin accessibility was changing around the PAE time-window, which concerned regions that gained chromatin accessibility between E14.5 and E15.5, in physiological conditions (Figure 4A,B; Supplementary Data notebook 4, Supplementary.

Dataset 2). This represented 5 regions out of 313 regions of the capture repertoire, among 1,442 total regions in the differentially opened regions (DOR) of the original ENCODE dataset (please see the legend of Figure 5 for normalization details taking into account the constraints of the capture approach). Interestingly, these 5 regions corresponded to active enhancers and were associated with with GO terms of importance for neurodevelopment and brain function: for example, cell motility and cell migration (Elmo1), synaptic signalling (Plcl2), regulation of multicellular organization process and cytoskeleton protein (Mapre2), and cell development, cell biogenesis, plasma membrane part, abnormal mouse morphology (Antxr). In addition, our analysis of the RNA-seq ENCODE dataset (Figure 4) revealed that DMRs were significantly associated with genes that were upregulated between E14.5 and E15.5, and between E15.5 and E16.6 under physiological conditions (Figure 4D; for quality control of these data, see Supplementary Figure 6A,B).

Immediate DNA-methylation and concomitant gene expression changes affect mono-allelically expressed gene upon PAE

First, when exploring CpG islands (CGis), which are relevant targets of DNA methylation and its potential perturbations [START_REF] Saxonov | A genome-wide analysis of CpG dinucleotides in the human genome distinguishes two distinct classes of promoters[END_REF], we found that only 36 out of the 432 identified DMRs were located into CGis. However, a quarter of these DMRs (9 out of 36; Figure 5A,B), significantly, corresponded to imprinted genes, which are tightly regulated by DNA methylation [START_REF] Perez | New Perspectives on Genomic Imprinting, an Essential and Multifaceted Mode of Epigenetic Control in the Developing and Adult Brain[END_REF].

As a whole, 13 DMRs were associated to 10 distinct imprinted genes in our analysis (9 DMRs in CGis, 4

DMRs not associated with CGis; Figure 5A,B). Moreover, these DMRs were located into the imprinting control regions (ICRs), which govern the expression of the imprinted locus (Table 3). This statistically significant over-representation DMR associated to imprinted genes among DMRs was confirmed by a hypergeometric test (see details in Materials and Methods). Interestingly, a majority of DMRs associated to imprinted genes were located within the promoter regions (8/13) and/or active enhancers in the developing brain (11/13, Figure 5A).

In addition, we found that a significative number of early DMRs were associated to clustered protocadherin genes (Pcdhs). Clustered Pcdhs are also mono-allelically expressed, but in a random manner, and are individually regulated by DNA methylation [START_REF] Phillips | Clustered protocadherin trafficking[END_REF]. Namely, we identified 10 DMRs in 9 distinct Pcdh, from both α-and γ-protocadherins clusters (Figure 5C). For 4 Pcdhs, DMRs observed immediately upon PAE are also located in a CpGi (Figure 5C).

The tight regulation of the monoallelic expression of imprinted genes and protocadherins by DNA methylation is crucial for neurodevelopment and normal brain functions [START_REF] Davies | What are imprinted genes doing in the brain?[END_REF]. For example, clustered Pcdhs play important roles in the modulation of dendrites arborisation and synaptogenesis or limitation of autapse formation through cell-cell recognition [START_REF] Light | δ-Protocadherins: Organizers of neural circuit assembly[END_REF][START_REF] Phillips | Clustered protocadherin trafficking[END_REF]reviewed in El Hajj et al., 2016;[START_REF] Matsunaga | Reelin transiently promotes N-cadherindependent neuronal adhesion during mouse cortical development[END_REF][START_REF] Molumby | γ-Protocadherins Interact with Neuroligin-1 and Negatively Regulate Dendritic Spine Morphogenesis[END_REF][START_REF] Yamagata | Cadherins Interact With Synaptic Organizers to Promote Synaptic Differentiation[END_REF]. The early alteration of their methylation levels by PAE could thus potentially affect their expression with potential defects typical of FASD. We thus investigated whether the establishment of DMRs in these genes were concomitant with changes in their expression levels. Using RT-qPCR, we examined the expression of four imprinted genes (Inpp5f_v2, Nap1l5, Peg13, and Zrsr1), three genes of protocadherin clusters (Pcdhα3, Pcdhα9, and Pcdhγa2). We also included three genes in our analysis, whose mutations in human lead to neurodevelopment disorders and which were also associated with DMRs with similar percentages of methylation differences (Supplementary Figure 7A):

Mid1, associated to the Opitz BBB/G syndrome reviewed in [START_REF] Winter | The MID1 protein is a central player during development and in disease[END_REF] andShank2 andAuts2 both linked to autism spectrum disorders (ASD; reviewed in [START_REF] Bourgeron | From the genetic architecture to synaptic plasticity in autism spectrum disorder[END_REF][START_REF] Pang | Untangle the Multi-Facet Functions of Auts2 as an Entry Point to Understand Neurodevelopmental Disorders[END_REF], with which FASD share overlapping traits [START_REF] Lussier | Prenatal Adversity Alters the Epigenetic Profile of the Prefrontal Cortex: Sexually Dimorphic Effects of Prenatal Alcohol Exposure and Food-Related Stress[END_REF][START_REF] Popova | Comorbidity of fetal alcohol spectrum disorder: a systematic review and meta-analysis[END_REF]. We showed that these genes which were associated with immediate DMRs exhibited concomitant disturbances of their expression in the prenatal cortex (Figure 6).

DISCUSSION

The impact of PAE on the DNA methylome in the brain has been almost exclusively interrogated in the postnatal and adult brain, at long temporal distance from the time-window of exposure [START_REF] Laufer | Long-term alterations to DNA methylation as a biomarker of prenatal alcohol exposure: From mouse models to human children with fetal alcohol spectrum disorders[END_REF] and the question whether perturbations of DNA methylation occur quickly after exposure has remained elusive. A recent study also investigated the effect of prenatal alcohol exposure during the preimplantation stages of the mouse embryo and, 10 days later, identified DMRs that interestingly affected genes important for neurodevelopment [START_REF] Legault | Pre-implantation alcohol exposure induces lasting sex-specific DNA methylation programming errors in the developing forebrain[END_REF]. To the best of our knowledge, our study represents the first analysis of the effect of PAE on the DNA methylome in the developing neocortex, addressing the question of the immediate deposition of aberrant DNA methylation marks, at very short temporal distance from the exposure.

Statistical intra-and inter-group variability

We observed that the intra-group variability was limited (Supplementary Figure 2D). This is likely due to our sampling mode which aimed at minimizing inter-individual variability in the same group, by mixing half-cortices from four embryos from four different litters in the same NGS sample (Supplementary Figure 1B). However, we also observe that the triplicates of the EtOH-treated group show more differences between each other than do triplicates from the control group. Interestingly, this is also observed in another paradigm prenatal stress due to neuroinflammation between samples within the exposed group [START_REF] Schang | Epigenome and transcriptome landscapes highlight dual roles of proinflammatory players in a perinatal model of white matter injury[END_REF], and more generally the exacerbated cell-to-cell variability provoked by PAE [START_REF] Ishii | Variations in brain defects result from cellular mosaicism in the activation of heat shock signalling[END_REF]. Indeed, this might reflect the immediate impact of cellular stress caused by ethanol exposure, and the disequilibrium, with some stochastic aspects, it imposes to cellular functions and molecular pathways, including the regulation of DNA methylation, before the system reaches a new equilibrium.

Alcohol exposure has multiple potential effects on DNA methylation in terms of availability of precursors of the methyl group [START_REF] Hutson | Folic acid transport to the human fetus is decreased in pregnancies with chronic alcohol exposure[END_REF][START_REF] Kleiber | Long-term genomic and epigenomic dysregulation as a consequence of prenatal alcohol exposure: a model for fetal alcohol spectrum disorders[END_REF][START_REF] Chater-Diehl | Changes to histone modifications following prenatal alcohol exposure: An emerging picture[END_REF], elevation of DNMT levels [START_REF] Miozzo | Alcohol exposure promotes DNA methyltransferase DNMT3A upregulation through reactive oxygen species-dependent mechanisms[END_REF], or potential redistribution of DNMTs in the genome, which would explain the fact that gain and loss of methylation can be observed depending on the genomic region considered. Our statistical analyses suggested that inter-group variability (control versus PAE-treated samples) is limited, and that PAE therefore might not majorly reshape the DNA methylome -at least in the repertoire of captured regions (Supplementary Figure 2D). Three observations could explain these results: i) it could result from the stochasticity of PAE effects on DNA methylation, and consequently, due to the fact that we capture only the most frequent and robust alterations in DNA methylation; this would not mean that the other perturbations of DNA methylation in one individual couldn't be meaningful in terms of clinical phenotypic consequences [START_REF] Cobben | DNA methylation abundantly associates with fetal alcohol spectrum disorder and its subphenotypes[END_REF]. ii) Another possibility is that, in spite of pleiotropic alcohol impacts, the DNA methylome is globally unchanged after such a binge-drinking-like mode of PAE due to rapid protective mechanisms, which would be interesting to unravel. iii) Alternatively, this apparent robustness of the DNA methylome might only reflect the fact that other perturbations of DNA methylation could occur later in life, for example due to abnormal neuronal plasticity due to PAE, that could reshape the DNA methylome profile in the postnatal and adult brain [START_REF] Guo | Neuronal activity modifies the DNA methylation landscape in the adult brain[END_REF]. Nevertheless, studies that explored the DNA methylation patterns at temporal distances in animal models or FASD cohorts, using capture-like or genome-wide approaches also identified a limited number of DMRs (a few hundred in most studies, around 6600 promoter regions in one study; [START_REF] Laufer | Long-lasting alterations to DNA methylation and ncRNAs could underlie the effects of fetal alcohol exposure in mice[END_REF]Laufer et al., , 2016;;[START_REF] Portales-Casamar | DNA methylation signature of human fetal alcohol spectrum disorder[END_REF][START_REF] Lussier | DNA methylation as a predictor of fetal alcohol spectrum disorder[END_REF][START_REF] Cobben | DNA methylation abundantly associates with fetal alcohol spectrum disorder and its subphenotypes[END_REF]. Whether these DMRs originate from early or late alterations of DNA methylation patterns, this means that the DNA methylome is altered in a limited number of regions, which does not mean that this could be detrimental for the brain.

Detection of individual differentially methylated cytosines, shortly after PAE

We detect few individual cytosines showing statistically differential methylation after PAE. Two of them, hypermethylated, reside into genes, KALRN and Tiam2, that are involved in a pathway important for neurodevelopment: guanine nucleotide exchange factors. These genes participate to neuronal shape and polarity, axon growth and/or neuronal plasticity [START_REF] Cahill | Kalirin regulates cortical spine morphogenesis and disease-related behavioral phenotypes[END_REF][START_REF] Honda | Rufy3 is an adapter protein for small GTPases that activates a Rac guanine nucleotide exchange factor to control neuronal polarity[END_REF]. KALRN misexpression has been linked to neuropsychiatric disorders like schizophrenia and addiction, and this gene has been found, in genome-wide association studies, to be related to ADHD, and schizophrenia, all pathologies being relevant for FASD (reviewed in [START_REF] Remmers | Abnormal kalirin signaling in neuropsychiatric disorders[END_REF]. Because the methylation of an individual CpG can affect gene expression [START_REF] Xu | Pioneer factor interactions and unmethylated CpG dinucleotides mark silent tissue-specific enhancers in embryonic stem cells[END_REF], these immediate DMCs could have deleterious effects on neuronal function.

Immediate DMRs are mostly found in brain active enhancers upon PAE

We have mainly focused on the identification and exploration of DMRs, which are, in general, more robustly associated to impacts on gene expression. We identified 432 DMRs, showing that DNA methylation is rapidly redistributed in the genome, shortly after PAE. We found that these immediate DMRs are not randomly distributed among the repertoire of capture regions: indeed, in contrast to promoters and other genomic regions, which are unaffected, enhancers that are known to be active in the adult brain are significantly over-represented, since more than 75 % of DMRs identified in our capture approach correspond to these genomic features. Interestingly, DMRs in these enhancers, that are particularly important for brain functions, if they would persist overtime in the adult, could disturb their activities, associated gene expression and explain part of the brain defects associated to FASD.

The fact that we could identified statistically robust immediate DMRs is suggestive that, in the perturbed, somewhat chaotic, environment, PAE triggers immediate DNA methylation alterations that are common to and could explained shared phenotypic traits among individuals diagnosed for FADS.

Conversely, other immediate changes that would occur more randomly amongst individuals might also have phenotypic consequences and be linked to the clinical variability observed in FASD patients [START_REF] Cobben | DNA methylation abundantly associates with fetal alcohol spectrum disorder and its subphenotypes[END_REF].

Immediate DMRs are linked to genes showing dynamic, physiological changes in chromatin accessibility and gene expression

We found that DNA methylation changes occur in specific regions that are, under no stress conditions, physiologically modified in terms of chromatin accessibility and gene expression, at the stage at which PAE was inflicted. First, DMRs are particularly enriched in regions whose chromatin gains accessibility between E14.5 and E15.5 in physiological conditions, and the corresponding enhancers are associated with genes whose GO terms underline importance for brain development.

Any perturbation of their methylation landscape by PAE might thus impair the physiological expression of these genes. Second, DMRs are also significantly associated with genes exhibiting upregulation of their expression, either between E14.5 and E15.5 or between E15.5 and E16.5. Notably, Satb2, one of the genes that undergoes upregulation during these two time-windows, in physiological conditions, is associated in our analysis with a PAE-induced hypermethylated DMR located in an enhancer (DMR in intron 2 on 10). This gene was previously identified by Hashimoto-Torii et al ( 2011) as the most downregulated gene in a similar paradigm of PAE (ethanol injections between E14 and E16). They also found that SATB2 expression was dysregulated in human cortices isolated from gestational week 15-18 fetuses, which were ex vivo exposed to ethanol. In their mouse model, the perturbation of Satb2 gene expression persists after birth (postnatal day 14, [START_REF] Hashimoto-Torii | Combined transcriptome analysis of fetal human and mouse cerebral cortex exposed to alcohol[END_REF], suggesting that the presence of such DMRs could potentially and persistently alter gene expression at temporal distance from the last ethanol injection. Similarly, the DMR-associated Mapre2 gene, which encodes a microtubule-associated protein, shows physiological increased chromatin accessibility between E14.5 and E15.5, and upregulated expression between E15.5 and E16.5. DNA methylation changes in Mapre2 DMR could therefore been followed by alteration in its expression. Interestingly, Tiam2, which contains one statistically significant DMC also belongs to the group of genes upregulated between E15.5 and 16.5.

In addition, we verified that three genes, Mid1, Auts2 and Shank2, whose mutations are associated with neurodevelopmental disorders including ASD for Auts2 and Shank2, are associated with DMRs and show correlative disturbances of expression, immediately after PAE.

Immediate DNA-methylation changes affect mono-allelically genes and correlates with misexpression upon PAE

Immediate DMRs are over-represented in imprinted genes in our PAE paradigm. They are mainly located within the promoter regions and active enhancers. We find that most of imprinted genes are ranked in the top of DMR-containing genes in terms of percentage of differential methylation and that these DMRs mostly overlap with the ICR of these genes, which is tightly controlled by DNA methylation and essential for parental monoallelic expression of these genes. The other family of genes prominently showing immediate alteration of DNA methylation patterns is the clustered protocadherin genes, which are individually regulated by DNA methylation. Interestingly, these two family of genes, which are key for neurodevelopment and neuronal function, have been already identified as associated to DMRs as in the adult brain in mice or FASD cohorts [START_REF] Laufer | Long-lasting alterations to DNA methylation and ncRNAs could underlie the effects of fetal alcohol exposure in mice[END_REF][START_REF] Lussier | Epigenetics studies of fetal alcohol spectrum disorder: where are we now?[END_REF] and in buccal swabs in cohorts of FASD children [START_REF] Cobben | DNA methylation abundantly associates with fetal alcohol spectrum disorder and its subphenotypes[END_REF][START_REF] Laufer | Associative DNA methylation changes in children with prenatal alcohol exposure[END_REF] therefore at temporal distance from the cessation of alcohol exposure.

Historically, alteration of DNA methylation patterns in imprinted genes has been also identified at the H19/Igf2 control region in mouse preimplantation embryos exposed in utero to alcohol at temporal distance of exposure, since material (whole embryos and placentas) was collected at E10.5 [START_REF] Haycock | Fetal Alcohol Spectrum Disorders: The Epigenetic Perspective[END_REF]. Early alteration in the DNA methylation status of imprinted genes has been suggested by [START_REF] Downing | Subtle decreases in DNA methylation and gene expression at the mouse Igf2 locus following prenatal alcohol exposure: effects of a methyl-supplemented diet[END_REF], who observed subtle decrease in DNA methylation and gene expression at the mouse Igf2 locus at E9, four hours after prenatal alcohol exposure in whole embryos and placentas, which were ameliorated by diet supplementation with methyl-group precursors [START_REF] Downing | Subtle decreases in DNA methylation and gene expression at the mouse Igf2 locus following prenatal alcohol exposure: effects of a methyl-supplemented diet[END_REF].

Other studies have also identified differential methylation status in imprinted genes after ex vivo ethanol exposure of cultured mouse embryos at the early neurulation stage (after 44 hours of exposure, [START_REF] Liu | Alcohol exposure alters DNA methylation profiles in mouse embryos at early neurulation[END_REF], as well as in human embryonic stem cells exposed to ethanol for 24 or 48 hours [START_REF] Khalid | Gene expression signatures affected by alcohol-induced DNA methylomic deregulation in human embryonic stem cells[END_REF]. These two ex vivo studies were suggestive that at least some imprinted genes could show early modifications in their DNA methylation profiles and thus, in line with our unbiased in vivo findings, using a methylome capture in the neocortex and two hours after PAE, only.

Among the genes whose expression is physiologically upregulated between E14.5 and E15.5 and that are associated with DMRs in our PAE paradigm, we identify Grb10, whose specific expression from the paternally inherited allele is involved in adult behaviour [START_REF] Dent | Brain-expressed imprinted genes and adult behaviour: the example of Nesp and Grb10[END_REF][START_REF] Perez | New Perspectives on Genomic Imprinting, an Essential and Multifaceted Mode of Epigenetic Control in the Developing and Adult Brain[END_REF].

Interestingly, Grb10 has already been identified as differentially methylated in adult mice, after PAE by Laufer et al. [START_REF] Laufer | Long-lasting alterations to DNA methylation and ncRNAs could underlie the effects of fetal alcohol exposure in mice[END_REF] and by [START_REF] Liu | Alcohol exposure alters DNA methylation profiles in mouse embryos at early neurulation[END_REF], as well as deregulated expression in response to very brief alcohol exposure in astrocytes [START_REF] Pignataro | The regulation of neuronal gene expression by alcohol[END_REF]. Interestingly, some DMRs are located in promoter region of several protocadherin genes, including Pcdhgb4 and Pcdhga12.

These observations suggest that expression of these particular genes could be affected both immediately and later after the PAE. PAE therefore quickly disturbs DNA methylation in these mono-allelically expressed genes, which are key for neurodevelopment. Moreover, we show that, at least for some of them, the presence of immediate DMRs in the ICR correlates with concomitant perturbation of their expression, which could potentially have rapid detrimental consequences on neurodevelopment. Note that the expression of the genes is disturbed in a relatively modest manner. However, even a 20% difference in gene expression upon PAE correlates with neurodevelopmental defects [START_REF] El Fatimy | Heat shock factor 2 is a stress-responsive mediator of neuronal migration defects in models of fetal alcohol syndrome[END_REF].

Moreover, as quoted above, at long temporal distances from PAE, DMRs were detected in several mono-allelically expressed genes of these two families, in the adult brain in mice [START_REF] Laufer | Long-lasting alterations to DNA methylation and ncRNAs could underlie the effects of fetal alcohol exposure in mice[END_REF][START_REF] Lussier | Epigenetics studies of fetal alcohol spectrum disorder: where are we now?[END_REF][START_REF] Legault | Pre-implantation alcohol exposure induces lasting sex-specific DNA methylation programming errors in the developing forebrain[END_REF] and in buccal swabs in cohorts of FASD children [START_REF] Cobben | DNA methylation abundantly associates with fetal alcohol spectrum disorder and its subphenotypes[END_REF][START_REF] Laufer | Associative DNA methylation changes in children with prenatal alcohol exposure[END_REF]. Our data thus suggest that these immediate alterations in DNA methylation may persist throughout life. In addition, the fact that we demonstrate that they occur so early upon PAE reinforces their potential as biomarkers of exposure. Notably, DMR persistence might not consist in a perfect conservation of the percentage of methylation for each CpG. Indeed, DNA methylation overtime might be modified by many factors, including brain activity itself [START_REF] Guo | Neuronal activity modifies the DNA methylation landscape in the adult brain[END_REF]. Rather, such persistence might be more qualitative than quantitative in the sense that a given DMR may globally exhibit differential methylation, both at short and long temporal distance from the end of exposure, but with different percentage of individual cytosine methylation at an early, compared to a late time-window.

CONCLUSION

To the best of our knowledge, our data demonstrate for the first time in an unbiased manner, that DNA methylation alterations in the developing cortex are detected immediately after the cessation of PAE across the genome. Although the number of DMRs that are formed right after PAE is limited, these alterations do not occur randomly, but target enhancer regions that are active at the time of alcohol exposure. Moreover, regions that undergo changes in chromatin accessibility and gene expression seem to be particularly vulnerable to perturbation of DNA methylation patterns and are associated with genes of importance for neurodevelopment. In addition, we show that genes of two gene families, whose monoallelic expression is tightly controlled by DNA methylation and crucial for brain development and function, carry immediate DMRs associated to synchronous dysregulation of their expression. These gene families have been previously identified as carrying or associated to DMRs in the adult brain, in human children or mouse models after PAE. Although they were suspected to be long-lasting epigenetic marks of PAE, there was no direct evidence that that they were established early in response to PAE. Our observations demonstrate that abnormalities in DNA methylation profiles are established immediately upon PAE and strongly suggest that they can persist at very long temporal distance from the end of exposure, and this could be considered as "scars" of exposure to prenatal stress. This makes these two gene families strong candidates as biomarkers of exposure, which are so urgently needed to better and earlier diagnose FASD children and accelerate and ameliorate follow-up and intervention.

MATERIALS AND METHODS

Mice mating and alcohol exposure

C57Bl/6N murine females of 2-to 4-months-old were time-mated and assessed for mating based on the presence of a vaginal plug. The noon of vaginal plug day was considered as embryonic day 0.5 (E0.5). For binge drinking stress, pregnant females received intra-peritoneal injections of ethanol (3g/kg, diluted in a final volume of 500 µL of PBS), at embryonic days E15, E15.5 and E16 (Figure 1A).

Control group received similarly PBS alone. Embryonic cortices were collected 2 hours after the last injection.

Tissue collection and genotyping

Embryonic cortices were harvested on ice, in cold L-15 medium (Leibovitz Gibco #11415-049).

Only male cortices were used for methylome analysis, whereas both male and female cortices were used for gene expression analysis. For methylome analysis, three replicates were generated per group.1M, 2M and 3M are ethanol-treated samples, while 4M, 5M and 6M design the PBS-treated ones (control). In order to reduce inter-litter variability, each replicate was composed of 2 right hemicortices and 2 left-hemi cortices, from embryos originated from 4 distinct litters, see Supplementary

Figure 1B.

Since sexual dimorphism was reported for behaviors of individuals exposed in utero to alcohol [START_REF] Hellemans | Prenatal alcohol exposure: Fetal programming and later life vulnerability to stress, depression and anxiety disorders[END_REF], we determined the sex of each embryo. Genomic DNA was extracted from animal tails incubation (95°C for 1h) in the extraction buffer (25 mM NaOH, 0.2 mM EDTA) then neutralized in Tris-HCl 40mM, pH5 (vol/vol). Sex of the embryos were identified by PCR genotyping using Ube1R (5'-CACCTGCACGTTGCCCTT-3') and Ube1F (5'-TGGATGGTGTGGCCAATG-3') primers that target Ube1X and Ube1Y genes [START_REF] Sugimoto | X Chromosome Reactivation Initiates in Nascent Primordial Germ Cells in Mice[END_REF]. These primers allow to amplify a 252bp amplicon from Ube1X and 334bp one from Ube1Y.

ATAC-seq and RNA-seq ENCODE data

Data mining from ATAC-seq and RNA-seq datasets was performed, using available public ENCODE data. These ATAC-seq and RNA-seq experiments were done using embryonic and newborn forebrains of C57BL/6N mice (experiments from Bing Ren laboratory, UCSD for the ATAC-seq and experiments of B. Wold, from Caltech laboratory for the RNA-seq).

ENCODE accession number of the data analysed

Methylome capture and high-throughput sequencing

To analyze the methylation of genomic regions that we selected (see below details about the design of regions, Figure 1B and Supplementary Figure 1A), a methylome capture was performed, using SeqCapEpi Developer Medium Enrichment kit (NimbleGen, Roche), following User's Guide and manufacturer's instructions. Briefly, gDNA samples were fragmented using Covaris technologies and non-directional libraries were prepared using KAPA Biosystems DNA Library Preparation Kits with

NimbleGen SeqCap Adapter Kits. Then, bisulfite conversion of these DNA samples libraries was performed using the Zymo EZ DNA Methylation-Lightning kit, before a DNA amplification using KAPA HiFi HotStart Uracil and ReadyMix. Capture of these bisulfite-converted sample libraries was done by hybridization, using SeqCap Developer M Enrichment, SeqCap Epi Accessory, and SeqCap HE-Oligo Kits.

In order to reduce variability between samples, the six libraries (one per sample) were multiplexed and captured all together on beads. This multiplexing during the 'capture step' has already been successfully tested with similar capture approaches [START_REF] Allum | Characterization of functional methylomes by next-generation capture sequencing identifies novel disease-associated variants[END_REF][START_REF] Van Der Werf | A robust protocol to increase NimbleGen SeqCap EZ multiplexing capacity to 96 samples[END_REF]. DNA was then amplified again using the KAPA HiFi HotStart ReadyMix.

Captured and amplified bisulfite-converted gDNA fragments were sequenced using an Illumina sequencing instrument (Paired-end, 150bp, NextSeq 500 High Throughput). All these steps were performed by Institut du Cerveau et de la Moelle épinière (ICM) sequencing platform.

Methylome capture composition

The capture is composed of 58,611 chosen genomic regions (Supplementary dataset 3), based on ENCODE available data or based on previous lab results (Figure 1B and Supplementary Figure 1A).

Majority of captured regions corresponds to active enhancers in adult (8 weeks-old) mouse cortex (characterized by the H3K27ac histone mark, regions of +/-500 bases from the middle of the peaks, ENCODE data). Capture also includes mouse promoter regions (+/-500 bases from the transcription start site, ENCODE data). Due to repetitive sequences that prevent the design of specific probes of some regions, and also because of annotation database information differences, only 75% of known promoter regions could be included in the capture. "Other regions" category is composed of genomic regions of interest based on previous lab results. It corresponds to regions potentially involved in stress response: (i) HSF2 binding sites found in unstressed mice cortices of E16.5 embryos (ChIP-seq data, VM lab) ; and to dynamic regions upon stress : (ii) differentially opened or closed regions (DOCR) identified in isolated oligodendrocyte precursors (O4+ cells) from cortices of 5 days-old mice, after an inflammatory stress from postnatal day P1 to P5, mimicking a neuroinflammation in the third trimester of pregnancy in human, (Schang et al., 2018, GSE197563); (iii) enhancers, characterized by the H3K4me1 histone mark (+/-100 bases from the middle of the peaks) of adult (8 weeks-old) mouse cortices among genes (and at +/-20 kb from genes ends) that were differentially expressed in response to another prenatal stress (neuroinflammation; Krishnan et al., 2017, GSE197563). Detail number on these sub-category "other regions" are indicated Supplementary Figure 1A. List of mm9 coordinates of all the genomic regions contained in the capture are available in Supplementary Dataset 3.

Bioinformatic analysis of methylome capture

Bioinformatic workflow used to analyze methylome capture is described in details in Supplementary data Notebook 1 and key figures of the analysis are shown in Supplementary Table 2.

Key steps of the analysis and associated tools are shown in Supplementary Figure 2A.

Trimming, mapping, methylation count table, and visualization. Briefly, reads quality was verified using FASTQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). To decrease methylation call errors from poor quality data, a trimming was carried out with Trimmomatic [START_REF] Bolger | Trimmomatic: a flexible trimmer for Illumina sequence data[END_REF] and Cutadapt [START_REF] Martin | Cutadapt removes adapter sequences from high-throughput sequencing reads[END_REF]. Trimmed reads were then paired-end mapped on the whole Mus musculus reference genome mm9, using bismark [START_REF] Krueger | Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications[END_REF] which is adapted for bisulfite converted reads. Libraries were considered as non-directional for the mapping. Deduplication of reads was performed using deduplicate_bismark function (from bismark suite)

Only cytosines on CpG context were studied. For each sample, number of methylated and unmethylated cytosines that covered a given cytosine site was obtained using bismark_methylation_extractor (BME) function (from bismark suite), with no_overlap option to prevent double counting of paired-end results. Then, to select only information that corresponds to capture regions, bedtools intersect [START_REF] Quinlan | BEDTools: a flexible suite of utilities for comparing genomic features[END_REF] was used to intersect capture regions with mapping files, or with CpG count tables. To avoid removal of potential interesting sequences, at the border of capture regions, each region was extended up and downstream with 150 bases, before the intersections, and were named enlarged capture). Mapping files restricted to the enlarged capture regions were sorted and indexed using samtools suite [START_REF] Li | The sequence alignment/map format and SAMtools[END_REF] before their visualization on IGV [START_REF] Robinson | Integrative genomics viewer[END_REF].

Detection of differentially methylated isolated CpG (DMCs).

Data from BME were formatted on R (The R Core Team, 2018) in order to use MethylKit R package [START_REF] Akalin | methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles[END_REF], for detection of significant differentially methylated isolated CpG between control and alcohol-treated samples. For this statistical analysis, formatted data were filtered: cytosine positions displaying abnormal high read coverage (> 99.9 th percentile) or low read coverage (number of reads < 10) were excluded from the analysis, to avoid bias [START_REF] Krueger | DNA methylome analysis using short bisulfite sequencing data[END_REF]. Read coverage was also normalized before statistical analysis, using median to calculate scaling factor. Only methylation state of CpG sites that were covered in all samples was investigated for the analysis (corresponding to 1,259,111 positions). DMCs were considered as significant when difference of methylation state between the two groups was higher than 5%, for a qvalue < 0.05.

Detection of differentially methylated regions (DMRs).

We generated an R function (called get_close_loci()) to define DMRs (Supplementary Figure 3,A-D). This function combines neighbouring CpGs that share similar differential methylation state, (all hypo-or all hypermethylated CpGs, Supplementary Figure 3A,C). Because MethylKit provides relevant statistical information about individual CpG, we used MethylKit statistical output as an input of get_close_loci() function (i.e.

position of the CpG site, percent of methylation difference between groups for a given site, p-value and q-value). Data from BME were thus formatted on R (The R Core Team, 2018) before using

MethylKit [START_REF] Akalin | methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles[END_REF]. Formatted data were filtered and normalized, with the same parameters than those used for DMC detection. Only methylation state of CpG sites that were covered in all samples was investigated for the analysis (corresponding to 1,259,111 positions). To define DMRs with get_close_loci(), only CpG sites having a reliable methylation state must be used, since one non relevant CpG is sufficient to affect DMR detection (Supplementary Figure 3C). Thus, a pre-selection of relevant CpGs was done, according to their p-value (CpGs with pvalue < 0.07), whose threshold was determined using comparisons of real and random datasets (see below, Supplementary Figure S4, Supplementary Table 2). To define DMRs, parameters were defined as follow: association of at least 5 CpG having a same methylation state, in a region of maximum 2000 bases, with a maximum distance of 100 bases between two successive selected CpG, according to p-value threshold (pval < 0.07,

Supplementary Table 2).

Randomization to define appropriate parameters for DMR detection. Random datasets (n=2) were generated and compared to real dataset, to define parameters of get_close_loci() function that were appropriated for DMRs detection in the real data. To obtain these random data, methylation states from real samples, were randomly switched, CpG position by CpG position, between samples and replicates using our own R script (Supplementary Figure S4, see details in Supplementary data notebook 1). Each methylation state is thus reassigned to a given sample, but methylation states between distinct CpG positions are not shuffled. This randomization assigns biologically plausible methylation values, actually observed at a given cytosine position (this is not the case if randomization was global, i.e redistribution of CpG methylation state, regardless of the CpG position, because there would be a risk of assigning, for example, high DNA methylation rates to a cytosine position that is always detected as non-methylated, regardless the conditions, PBS or EtOH-treatment).

Randomization was done on filtered but unnormalized real dataset. Read coverage normalization was performed after the randomization, using median to calculate scaling factor. To define appropriate parameters for DMR detection, different parameters were modulated and the number of DMRs found in the two datasets (real versus random) was compared (Supplementary Table 2). Similar results were obtained by using the two random datasets (Supplementary Table 2).

Annotation. To be able to compare DMRs with other bioinformatic analysis (ATAC-seq and RNAseq results), data obtained with mm9 coordinates were then converted with mm10 coordinates using LiftOver software (https://genome.ucsc.edu/cgi-bin/hgLiftOver), available on UCSC, with default parameters. Annotation file, based on mm10 coordinates was obtained using BiomaRt R package [START_REF] Durinck | Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt[END_REF][START_REF] Huang | Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources[END_REF]. DMRs were annotated in R according to this file information (see details in Supplementary data notebook 1). A file combining all information (DMRs coordinates, annotation and statistical information) were generated (Supplementary Dataset 1). To estimate if

DMRs are particularly located in a given capture categories (enhancers, promoters…), hypergeometric tests were done on R (using phyper() function), taking into account the number of each capture category that is represented in the methylome capture. Threshold chosen: pvalue < 0.05.

Motif enrichment. Transcription factors binding site enrichment within the DMRs located in active enhancer regions (H3K27ac mark), were explored using findMotifsGenome (HOMER suite, see details in Supplementary data notebook 1 ; [START_REF] Heinz | Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities[END_REF].

Over-representation of Imprinted genes among DMRs. To estimate whether DMRs are particularly located in imprinted genes (IG), hypergeometric tests were done on R (using phyper() function).

Threshold chosen: pvalue < 0.05. The exact number of capture regions that corresponds to imprinted genes are unknown, but we estimated that about 150 capture regions could be attributed to imprinted gene regions (there are almost 150 mouse IG, according to http://www.geneimprint.com/site/genesby-species.Mus+musculus website, and we potentially have, at least one promoter for each of this gene in the capture). Using this approximation, results of phyper test (phyper(11,150,58611-150,432))

is equal to 1.51e-09, which is significant. Even with an over-estimation of capture regions that correspond to imprinted gene regions, results of phyper test is significant (e.g. if we considered that 300 capture regions corresponds to imprinted genes regions (i.e. two regions per IG), phyher(11,300,58611-300,432) = 2.89e-06).

Estimation of sodium bisulfite conversion efficiency. To ensure that sodium bisulfite treatment correctly converted unmethylated cytosines, we evaluated the conversion efficiency by looking at conversion rate of a spike-in DNA (i.e. unmethylated known sequences of non-mammalian DNA), which was added in each sample preparation. This spike-in DNA, corresponding to sequences of lambda phage, was mapped on an appropriate reference genome, using bismark. The efficiency of sodium bisulfite conversion is obtained using deduplicated spike-in DNA mapped reads, by calculating the ratio of converted reads (containing thymines) to the total number of reads covering these cytosines (Supplementary Table S1).

Chromatin accessibility profile of the developing brain: ATAC-seq data time-course analysis

Bioinformatic workflow used to analyze ATAC-seq data is described in details in Supplementary data notebook 2. Key steps of the analysis and associated tools are showed in Supplementary Figure S2A and key features of the analysis are shown in Supplementary Table S4. Quality of sequenced reads was verified for each sample using FASTQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Reads having bad sequencing quality were trimmed using Trimmomatic [START_REF] Bolger | Trimmomatic: a flexible trimmer for Illumina sequence data[END_REF]. Trimmed reads were then paired-end mapped on Mus musculus reference genome mm9, using bowtie2 [START_REF] Langmead | Fast gapped-read alignment with Bowtie 2[END_REF]. Then, reads were deduplicated using samtools rmdup [START_REF] Li | The sequence alignment/map format and SAMtools[END_REF]. Reads that mapped to mitochondrial chromosome were removed using tools on Galaxy platform (Afgan et al., 2016, default parameters), by splitting the mapped reads per chromosome using bam-splitter tool [START_REF] Barnett | BamTools: a C++ API and toolkit for analyzing and managing BAM files[END_REF]. Then all the files were merged except the reads mapped on mitochondrial chromosome, using merge-bam tool.

MACS2 [START_REF] Zhang | Model-based analysis of ChIP-Seq (MACS)[END_REF] was then used to identify peaks corresponding to opened genomic regions, in each sample. For that, mapped reads from duplicates were merged, before this Peak calling. Then, table containing all opened regions, detected in at least one sample, was generated, using bedtools multiinter and bedtools merge [START_REF] Quinlan | BEDTools: a flexible suite of utilities for comparing genomic features[END_REF]. To obtain reads count at these regions, in each sample, htseq-count [START_REF] Anders | HTSeq--a Python framework to work with highthroughput sequencing data[END_REF] was runned. This tool prevents double counting of pairedend results since files containing mapped reads were sorted by read name. To identify regions where chromatin accessibility significantly changes during development, pairwise comparisons of successive developmental stages were performed using edgeR [START_REF] Robinson | edgeR: a Bioconductor package for differential expression analysis of digital gene expression data[END_REF], an R software (The R Core Team, 2018) and Bioconductor package [START_REF] Gentleman | Bioconductor: open software development for computational biology and bioinformatics[END_REF]. We adapted a script template called template_script_edgeR.r from Sartools R package [START_REF] Varet | SARTools: A DESeq2-and EdgeR-Based R Pipeline for Comprehensive Differential Analysis of RNA-Seq Data[END_REF], which implements some edgeR functions (see details in Supplementary Data notebook 2). To compare ATAC-seq results with other bioinformatic analyses (methylome capture and RNA-seq results), data obtained with mm9 coordinates were then converted with mm10 coordinates using LiftOver software (https://genome.ucsc.edu/cgi-bin/hgLiftOver), with default parameters. Annotation file, based on mm10 coordinates was obtained using BiomaRt R package [START_REF] Durinck | Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt[END_REF][START_REF] Huang | Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources[END_REF].

Differentially open and closed regions (DOCR) observed during brain development were annotated in R according to this file information (see details in Supplementary Data notebook 2). For each pairwise comparison, a file combining all information (DOCR coordinates, annotation and statistical information) were generated (Supplementary Dataset 1). To visualize data on IGV [START_REF] Robinson | Integrative genomics viewer[END_REF], files containing mapped reads were sorted and indexed using samtools suite [START_REF] Li | The sequence alignment/map format and SAMtools[END_REF].

Transcriptome profile of the developing brain: RNA-seq data time-course analysis

Bioinformatic workflow used to analyze RNA-seq data is described in detail in Supplementary data notebook 3. Key steps of the analysis and associated tools are showed in Supplementary Figure 2A and key features of the analysis are shown in Supplementary Table 5. For all samples, sequenced read quality was verified for each sample using FASTQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Adapter sequences and reads having bad quality were trimmed using TrimGalore (https://github.com/FelixKrueger/TrimGalore). Trimmed reads were then singled-end mapped using STAR on Mus musculus reference genome mm10, which was associated to its corresponding Gencode annotation during indexation, to improve accuracy of the mapping [START_REF] Dobin | Mapping RNA-seq Reads with STAR: Mapping RNAseq Reads with STAR[END_REF][START_REF] Dobin | STAR: ultrafast universal RNA-seq aligner[END_REF]. For each sample, htseq-count [START_REF] Anders | HTSeq--a Python framework to work with highthroughput sequencing data[END_REF] was used to obtain reads count of all genes contained in Gencode annotation file. To identify differentially expressed genes (DEG) during physiological brain development, pairwise comparisons of successive developmental stages were performed using edgeR [START_REF] Robinson | edgeR: a Bioconductor package for differential expression analysis of digital gene expression data[END_REF]. We used exactly the same script than those executed for ATAC-seq analysis (i.e. modified script from Sartools R package [START_REF] Varet | SARTools: A DESeq2-and EdgeR-Based R Pipeline for Comprehensive Differential Analysis of RNA-Seq Data[END_REF], see details in Supplementary data notebook 3). Annotation file, based on mm10 coordinates was obtained using BiomaRt R package [START_REF] Durinck | Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt[END_REF][START_REF] Huang | Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources[END_REF] and was used to annotate DEG observed during brain development with Unix commands (see details in Supplementary data notebook 3). For each pairwise comparison, a file combining all information (DEG coordinates, annotation, and statistical information) were generated (Supplementary Dataset 1). To visualize data on IGV [START_REF] Robinson | Integrative genomics viewer[END_REF], files containing mapped reads were sorted and indexed using samtools suite [START_REF] Li | The sequence alignment/map format and SAMtools[END_REF].

Gene ontology analyses

In order to perform Gene ontology analysis, gene lists (based on gene symbol names) were defined for each distinct datasets (i.e. genes associated to DMRs that are located in (or close to) captured active enhancers (H3K27ac) and genes associated to DOCR and DEG). For DOCR and DEG, only genes that are modulated between E14.5 and E15.5 and between E15.5 and E16.5 were studied. Gene ontology analyses were performed with WEB-based GEne SeT AnaLysis Toolkit (GESTALT, http://www.webgestalt.org/, using default parameters, except for the minimum number of genes for a category (set to one in our analysis). The reference gene list used for each analysis was either the genome for DOCR and DEG or all the genes of the studied capture category (i.e. in case of the "DMRs enhancer" list: all the genes associated to H3K27ac histone mark that are represented in the capture design (Supplementary Dataset 3).

RNA extraction and purification

Each E16.5 embryonic cortex in 500µL Trizol (Sigma) was lysed at 4°C using Precellys/Cryolys (Bertin) with 1.4mm ceramic beads (Mobio) during 10 seconds (5500 rpm). mRNAs were isolated using chloroform and isopropanol, followed by ethanol precipitation. Samples were resuspended into 50µL of ultrapure RNase-free Water (Invitrogen). Integrity of RNA was verified on an electrophoresis gel Reverse transcription was done using Biolabs kit (Lunascript RT SuperMix E010G), on less than 5µg of RNA, as recommended by the manufacter. cDNA concentration was normalized to 3ng/µL by diluting it into ultrapure water. To verify the purity of RNA extraction (absence of genomic DNA), a negative control without reverse transcriptase was also performed. Quality and Integrity of cDNA were verified by calculating 3'/5' ratio of Hprt1 gene, following MIQE recommendations [START_REF] Bustin | The MIQE Guidelines: Minimum Information for Publication of Quantitative Real-Time PCR Experiments[END_REF].

Only cDNA preparations with a ratio lower than 5 were used for RT-qPCR [START_REF] Bustin | The MIQE Guidelines: Minimum Information for Publication of Quantitative Real-Time PCR Experiments[END_REF].

Quantitative Real-Time RT-PCR

For real-time RT-qPCR, LunaGreen SYBR mix (Biolabs) was used, as recommended by the manufacturer.

Primers were purchased from Sigma-Aldrich. Probes sequences are listed in Supplementary Table 6.

384-wells PCR plates were filled using Tecan pipetting robot. RT-qPCR experiments were done using LightCycler LC480 (Roche). All reactions were made in quadruplicate for samples originated from at least three distinct litters.

Relative quantities of the target gene mRNAs were normalized to Hprt1 and Tpt1. Qbase+ software was used to define the values of calibrated and normalized relative gene expression (CNRQ).

Non parametric Mann-Whitney tests were performed using R software, to determine the pvalues when comparing mRNA levels between PBS and EtOH-treated samples. To take into account multiple tests, p-values were corrected with Benjamini-Hochberg correction (padj.), using p.adjust() R function.

For this analysis, a gene is considered as differentially expressed between PBS and EtOH-treated samples, when the adjusted p-value is below 0.05 (padj. < 0.05).

Sequences of primers used for RT-qPCR experiments. This screenshot was obtained using IGV. Scale is the same for all samples.

(C) Protocadherin genes that are associated to immediate DMRs upon PAE. As in (A). Results are generated using n=12 PBS samples and n=11 EtOH samples, from 3 distinct litters in each condition. Both males and females are considered together, since no sexual dimorphism is detected in terms of the expression of the studied genes (Supplementary Figure 7B,C).

Significance differences between PBS-and EtOH-treated samples were determined using nonparametric Mann-Whitney test. To take into account multiple tests, p-values were corrected into adjusted p-values (padj.), using Benjamini-Hochberg correction. Asterisks (*) denote significant differences between PBS-and EtOH-treated samples (when padj. < 0.05). ns, statistically nonsignificant results. Each dot represents calibrated normalized relative gene expression (CNRQ) value of each sample.

mRNA levels are normalized to the expression of Tpt1 and Hprt1 genes.

Results are generated using n = 6 females per condition (PBS or EtOH) and n = 5 or 6 males (6 for PBS, 5 for EtOH conditions) samples, from 3 distinct litters in each condition. Significant differences between PBS-and EtOH-treated samples were determined using non-parametric Mann-Whitney test.

To take into account multiple tests, p-values were corrected into adjusted p-values (padj.), using

Benjamini-Hochberg correction (padj. < 0.05 is considered as statistically significant result). For each gene, no statistically significant difference was observed between male and female samples (padj. always higher than 0.05). ,Pbx1,Sdccag8,Prox1,Ank3,Bcl11a,Kif3c,Eml1,Atxn1,Itga1,Farp1,Trappc9,App,Itsn1,Arid1b,Runx2,Efna5,Man2a1,Ptprm,Tcf4,Ablim1,Zeb2,Dclk1,Dclk2,Nfib,Ncdn,Spen,Cit,Cux2,Ctnna2,Mgll,Sox5,Zfp536,Inpp5f,Fgfr2,Shank2,Efnb2,Unc5d,Ephb1,Map4 GO:0022008 neurogenesis 1165 41 16.65 2.46 0.000137 Satb2,Pbx1,Sdccag8,Prox1,Ank3,Bcl11a,Kif3c,Eml1,Atxn1,Itga1,Farp1,Trappc9,App,Itsn1,Arid1b,Runx2,Efna5,Man2a1,Ptprm,Tcf4,Ablim1,Zeb2,Dclk1,Dclk2,Nfib,Ncdn,Spen,Cit,Cux2,Ctnna2,Mgll,Sox5,Zfp536,Inpp5f,Fgfr2,Shank2,Efnb2,Arhgef10,Unc5d,Ephb1,Map4 GO:0030182 neuron differentiation 990 35 14.15 2.47 0.001253 Satb2,Pbx1,Prox1,Ank3,Bcl11a,Kif3c,Itga1,Farp1,Trappc9,App,Itsn1,Arid1b,Runx2,Efna5,Ptprm,Tcf4,Ablim1,Zeb2,Dclk1,Dclk2,Nfib,Ncdn,Cit,Cux2,Ctnna2,Mgll,Sox5,Zfp536,Inpp5f,Fgfr2,Shank2,Efnb2,Unc5d,Ephb1,Map4 GO:0007399 nervous system development 1583 45 22.62 1.99 0.006082 Satb2,Pbx1,Sdccag8,Prox1,Ank3,Bcr,Apaf1,Bcl11a,Kif3c,Eml1,Atxn1,Itga1,Farp1,Trappc9,App,Itsn1,Arid1b,Runx2,Efna5,Man2a1,Ptprm,Tcf4,Ablim1,Zeb2,Dclk1,Dclk2,Nfib,Nfia,Ncdn,Spen,Mthfr,Cit,Cux2,Ctnna2,Mgll,Sox5,Zfp536,Inpp5f,Fgfr2,Shank2,Efnb2,Arhgef10,Unc5d,Ephb1,Map4 GO:0048468 cell development 1435 42 20.50 2.05 0.006082 Satb2,Pbx1,Prox1,Jmjd1c,Ank3,Bcl11a,Kif3c,Atxn1,Hrh2,Itga1,Farp1,App,Itsn1,Arid1b,Pacrg,Runx2,Efna5,Man2a1,Ptprm,Tcf4,Ablim1,Zeb2,Dclk1,Dclk2,Nfib,Ncdn,Spen,Cit,Cux2,Ctnna2,Mgll,Sox5,Sipa1l3,Zfp536,Inpp5f,Fgfr2,Shank2,Efnb2,Arhgef10,Unc5d,Ephb1,Map4 GO:0000904 cell morphogenesis involved in differentiation 549 23 7.84 2.93 0.006492 ,Jmjd1c,Ank3,Bcl11a,Hrh2,Farp1,App,Efna5,Ptprm,Ablim1,Zeb2,Dclk1,Nfib,Cit,Cux2,Ctnna2,Mgll,Sipa1l3,Fgfr2,Shank2,Efnb2,Unc5d,Ephb1 GO:0021953 central nervous system neuron differentiation 132 10 1.89 5.30 0.029620 Satb2, Prox1, Zeb2, Dclk1, Dclk2, Nfib, Sox5, Fgfr2, Unc5d, Ephb1 GO:0051239 regulation of multicellular organismal process 1836 47 26.23 1.79 0.029620 Cdc73,Pbx1,Prox1,Bcr,Grb10,Bcl11a,Sptbn1,Serpinf2,Kif3c,Ryr2,Atxn1,Hrh2,Nln,Farp1,Fam49b,Cblb,App,Itsn1,Runx2,Plcl2,Efna5,Man2a1,Ptprm,Tmem178,Prkce,Mapre2,Tcf4,Zeb2,Meis2,Nfib,Spen,Cmklr1,Cit,Cux2,Mgll,Iqsec1,Sox5,Zbtb32,Tshz3,Zfp536,Inpp5f,Fgfr2,Shank2,Efnb2,Unc5d,Gab1,Ephb1 GO:0048667 cell morphogenesis involved in neuron differentiation 448 19 6.40 2.97 0.029620 Ank3,Bcl11a,Farp1,App,Efna5,Ptprm,Ablim1,Zeb2,Dclk1,Nfib,Cit,Cux2,Ctnna2,Mgll,Fgfr2,Shank2,Efnb2,Unc5d,Ephb1 GO:0048812 neuron projection morphogenesis 494 20 7.06 2.83 0.030939 Ank3,Bcl11a,Itga1,Farp1,App,Efna5,Ptprm,Ablim1,Zeb2,Dclk1,Nfib,Cit,Cux2,Ctnna2,Mgll,Fgfr2,Shank2,Efnb2,Unc5d,Ephb1

Prox1

(

  detection of 18S and 28S ribosomal RNA). Extracts were DNAseI treated (Roche -ref 04 716 728 001) in presence of RNAse inhibitors (Promega -Ref N251A), before RNA purification using phenol chloroforme isoamyl-alcohol and chloroforme, followed by ethanol precipitation. Samples were resuspended into 20µL of ultrapure RNase-free Water. RNA was quantified with Denovix machine (DS-11 Denovix series).

Figure 4 .

 4 Figure 4. Immediate DMRs significantly fall into genomic regions that are remodeled in terms of chromatin accessibility or gene expression, during physiological brain development. (A) Number of capture regions that are associated to differentially opened or closed regions (DOCRs) or differentially expressed genes (DEGs) observed in the developing brain, under physiological conditions. DOCRs and DEGs are obtained in a pairwise-comparison manner of two successive developmental stages. Only DOCRs and promoters of DEGs strictly overlapping capture regions are taken into account (numbers indicated in black). Among these, DMRs strictly overlapping DOCRs/DEGs are considered for the analyses and their number of these regions is indicated in red. Asterisks indicated statistically significant results, i.e. when DMRs are found significantly over-represented (hypergeometric tests, see (B)). Total numbers of DOCRs and DEGs found during physiological brain development are indicated in Supplementary Figure 5A (i.e. independently of their presence in capture regions). Number of DOCR / DEG in Supplementary Figure 5A could be lower than the number indicated here, because a given DOCR or DEG could be represented by several capture regions. DOCR and DEG were identified using edgeR (using no logFC threshold but Benjamini-Hochberg p-value adjustment was performed and level of controlled false-positive rate was set to 0.05). (B) Hypergeometric test data showing the identification of over-represented DMRs among differentially opened or closed regions (DOCRs).

Figure 5 .

 5 Figure 5. Mono-allelically expressed genes presenting DMR immediately after PAE. (A) Imprinted genes associated to immediate DMRs upon PAE. Overlap between "PAE-affected" imprinted genes and capture regions corresponds to either strict overlap or to close vicinity (distance < 500bp). Negative values correspond to DMRs located close to capture regions with no overlap. Median pval.: median of p-values for each relevant CpG (with pval < 0.07) that composes the DMRs. (B) Illustration of a DMR observed in Peg10 imprinted gene. Visualization of deduplicated read coverage of EtOH-and PBS-treated samples. Genomic region containing a DMR with gain of methylation, located into the Peg10 gene is shown (thin blue bar at the bottom). Methylated cytosines in CpG context are in red, unmethylated ones in blue. Capture region is also indicated (thick blue bar).

Figure 6 .

 6 Figure 6. PAE-induced alterations of gene expression are concomitant with immediate perturbations of DNA methylation. Quantitative RT-PCR of 10 genes associated to immediate DMRs. Boxplots showing relative mRNA levels from E16.5 embryonic cortices upon PBS-or EtOH-treatment. Each dot represents calibrated and normalized relative gene expression (CNRQ) value of each sample. mRNA levels are normalized to the expression of Tpt1 and Hprt1 genes.

( B )

 B Comparison of mRNA expression does not reveal sexual dimorphism in the expression of these genes in fetal mouse cortices. (related to Figure 6) Quantitative RT-PCR of genes of interest, in (A) PBS-or (B) EtOH-treated samples. Boxplots showing Relative mRNA levels from E16.5 fetal male or female cortices, upon PBS-(A) or EtOH-treatment (B).

  Figure 1

Table 2). Immediate alterations in DNA methylation are detected in brain development upon PAE.

	Using this captured-specific bioinformatic workflow (

Supplementary Figure 2A; Supplementary Table 3, see Materials and Methods and Supplementary Data notebook 1) we

	identified 432 regions
	that were differentially methylated in embryonic cortices, immediately after PAE, among the 58,611
	regions included in the capture (Figure

gene (reference gene) Forward primer sequence Reverse primer sequence

	Legends of Figures
	Figure 1. Binge drinking model and DNA methylome capture design.
	(A) Binge drinking model. C57BL/6N pregnant mice received 1 intraperitoneal injection (IP)/day of
	ethanol (3g/kg) at embryonic days E15, E15.5 and E16. Control animals were similarly treated with IP
	of PBS. Fetal cortices were collected 2 hours after the last IP (at E16). For DNA methylome analysis,
	only male samples were studied. Three samples per group were generated (4M, 5M and 6M for the
	control group; 1M, 2M and 3M for the EtOH-treated one). To reduce variability between samples, each
	replicate was composed of 4 hemi-cortices, from 4 embryos of distinct litters (Supplementary Figure
	1B).
	Target gene (B) Design of the DNA methylome capture. The capture is composed of 58,611 chosen genomic Forward primer sequence Reverse primer sequence
	Auts2 regions, based on ENCODE available data (H3K27ac histone mark, a hallmark of enhancer regions and TCTGTCGGCTGAAATCCGAG CCTGGCTTACCGAGCTTCAA
	Mid1 promoter regions) and on our previous results (for other regions). For more details about capture CCATCTGGACTGTACGGAGC TTGGGGAAGATGTCACGACG
	Pcdhα3 composition, see Supplementary AGAAGACAAACTGGTTGGAGACAT CACTGGATACTGTTGGCCACT
	Pcdhγa2	CGTGCTTTCCTGCAGACCTA	TTGCTGAGGGGTTTCTTCTCTT
	Peg13	CCGGCCAGCAATCCTTATCT	TGAGGCACCCAAGTGGAATC
	Zrsr1	GAAAGCACTGCAACTTCCTTCA	TTTTACCAGAGGAGCCAGTCC
	Nap1l5	ATTCTGGCTAGCCCGTTTCT	GGATTTCAACGTGTGACGCA
	Inpp5f	TACACACACGGACTGGCTTC	GCGTCTGTCCCATTGGTTCT
	Pcdhα9	CAGCGAGTATGCTCAGGAGAG	CTGGGCAAGGTGTAAGACTGG
	Shank2	AGTCTATGCAGACGGGCAAG	GCGTGTCTTGGTAAAGTGCG
	standardizer
	Hprt1	ACAAAGCCTAAGATGAGCGCAAG	CTAGGCAGATGGCCACAGGAC
	Tpt1	GTTGCTCTCCTGGACTACCGTG	GCAGCCAATTATGGTGACAGGT

Figure 1A and Materials and methods. Figure 2. Immediate modifications of DNA methylation are observed in fetal cortices upon PAE (A) Differentially methylated CpG (DMC) identified upon PAE.

  Only cytosines in CpG context that are covered in all samples were investigated. DMC were identified using MethylKit (see details in Materials

	and Methods). Meth.diff: differential methylation rate (%) observed between control and EtOH-
	treated groups. A negative percentage indicates a loss of methylation, compared to control. q-val. <
	0.05; a minimum of 5% of Meth.diff is required for a CpG to be considered a DMC.
	(B)

Number of differentially methylated regions (DMRs) immediately observed upon PAE. Number

	of DMRs filtered by methylation differential rate (indicated as percent) observed between EtOH-
	treated and control groups are also indicated. Regions with gain (Hyper-DMRs) or loss (Hypo-DMRs) of
	methylation; (Hyper+Hypo), regions with either gain or loss of methylation. All: all DMRs.
	(C) Number of DMRs in each category of capture regions. A DMR is considered to be associated to a
	category of capture regions, if it overlaps to or if it is closed to a capture region (distance less than 500
	bp). Asterisks indicate categories that are significantly over-represented among DMRs, in a statistical
	manner, according to hypergeometric tests (see (D)).
	(D)

DMRs immediately detected upon PAE are predominantly located in brain active enhancers.

  They differ very slightly from the numbers indicated in (C), because some DMRs overlap several capture regions, and are, thus, counted more than one time.

	Figure 3.
	Results of hypergeometric tests realized to determine whether specific categories of capture regions

nare over-represented among DMRs. Significant values are highlighted in red. For this analysis, both DMRs that overlap or are located closed to capture region categories) were considered (i.e. distance less than 500bp between the regions. The number (Nb) of DMRs considered in each analysis is indicated.

Gene ontology (GO) analyses performed on genes associated to DMRs that are observed into (or close to) active enhancers Volcano plots showing statistically relevant GO biological processes

  for genes that are associated to (A) all DMRs (corresponding to either gain or loss of methylation) (B) to DMRs with gain of methylation, only, and observed into or close to active enhancers. The distance between DMRs and active enhancers regions is less than 500 bp. The size and share of color of each dot is proportional to the number of overlap with one given region category[START_REF] Liao | WebGestalt 2019: gene set analysis toolkit with revamped UIs and APIs[END_REF]. The genes associated to the ten most relevant GO Biological Processes of these analyses are described in Table1(for GO performed on all DMRs observed in active enhancers, related to Figure2D) and Table2(for GO performed only on regions that gain methylation upon PAE in active enhancers). These data were obtained using Web-

based GEne SeT AnaLysis Toolkit (WebGestalt ;

[START_REF] Liao | WebGestalt 2019: gene set analysis toolkit with revamped UIs and APIs[END_REF]

.

Table 2 . GO results for genes associated to Hypermethylated regions, observed upon EPA, that are localized in active enhancers (H3K27ac)

 2

	GO category	Description	Size	Over-lap	Expec-ted	Enrichment Ratio	FDR	User ID
	GO:0048699	neurons generation of	1088	40	15.55	2.57	0.000126	Satb2

Acknowledgements

We are grateful to Slimane AIT-SI-ALI, (Université de Paris, CNRS Epigenetics and Cell Fate, Paris, France), Maxim GREENBERG (Université de Paris, CNRS Institut Jacques Monod Paris) and Deborah BOURC'HIS (Institut Curie, Paris) for helpful discussions. We are grateful to the VM team members, Aurélie de THONEL and Véronique DUBREUIL (Université de Paris, CNRS Epigenetics and Cell Fate, Paris, France) for helpful discussions and comments on the manuscript. We thank Renaud MASSART, former postdoc in the VM Lab, for designing methylome capture. We are grateful to Emeline MUNDWILLER and Yann MARIE (Brain and Spine Institute, Paris, France) for performing the methylome capture libraries and high-throughput sequencing of these samples. We thank Isabelle Le PARCO and the staff from the Buffon animal housing facility at the Jacques Monod Institute (Université de Paris, Paris, France). We thank Magali HENNION, Head of the Bioinformatics and Biostatistics Core Facility at the CNRS Epigenetics and Cell Fate Center (Université de Paris) for bioinformatics support.

Funding information

VM was funded by the Agence Nationale de la Recherche « HSF-EPISAME », SAMENTA ANR-13-SAMA-0008-01) and CNRS. AD was supported by a Doctoral Fellowship from the French Ministère de l'Enseignement Supérieur, de la Recherche et de l'Innovation (MESRI) and FM from a Doctoral Fellowship from the CNRS. This study contributes to the Université de Paris IdEx #ANR-18-IDEX-0001 funded by the French Government through its "Investments for the Future" program. DSD benefited from a CNRS Délégation de Recherche (2018-2020). OTT benefited from travel grants from Université Paris Diderot.

.

Ethical issues

The breeding and treatments of wild type C57BL/6N mice, used for the experimental protocols described in this study have been approved by the Institutional Animal Care and Use Ethical Committee of the Paris University (registration number CEEA-40). The project has been recorded under the following reference by the Ministère de l'Enseignement Supérieur et de la Recherche (#2016040414515579). All efforts were made to reduce stress and pain to animals.

Legends of Supplementary Figures

Supplementary Figure 1 "Other regions" correspond to HSF2 genome-wide occupancy, a stress-responsive transcription factor involved in FASD phenotypes [START_REF] El Fatimy | Heat shock factor 2 is a stress-responsive mediator of neuronal migration defects in models of fetal alcohol syndrome[END_REF]ChIP-seq data), and to differentially accessible regions and H3K4me1-marked enhancer regions associated to gene differentially expressed in response to another prenatal stress, neuroinflammation [START_REF] Schang | Epigenome and transcriptome landscapes highlight dual roles of proinflammatory players in a perinatal model of white matter injury[END_REF]Krishnan et al., 2017, GSE197563; see and Methods). that categories are not exclusive (i.e. some regions are counted several times, e.g. 8,144 regions that are both included in H3K27ac enhancers and promoters categories, are also among the 34,279 regions that are included in H3K27ac enhancers). Black dots indicated dataset(s) that are used for the counting.

(B) Sampling protocol. To minimize intra-group and inter-individual variability, each NGS sample was composed of four hemi-cortices of four different male embryos from four different litters. ATAC-seq and RNA-seq data mining analyses were performed using publicly available ENCODE datasets, obtained from mice fetal forebrains at distinct developmental stages and newborn samples in physiological conditions. Key steps and tools used for each bioinformatic analysis are described here.

Supplementary

For more details, please see Materials and Methods and detailed command lines in Supplementary notebooks 1 to 4. QC, quality control; PE, paired-end analysis; SE, single-end analysis; Genc. annot., Gencode annotation; BME: bismark_methylation_extractor; DMC: differentially methylated isolated CpG; DOCR: differentially opened or closed regions upon physiological brain development stages; DEG, differentially expressed genes upon physiological brain development stages.

(B) Per-base-sequence-content plots obtained with FASTQC. Illustrated for 1M sample. An imbalance is observed between bases proportions at each read position as expected after DNA conversion with sodium bisulfite (BS). The conversion of unmethylated cytosines into thymines leads to an overrepresentation of thymines and an underrepresentation of cytosines at a given position for R1 reads (for R2 reads, which are complementary to R1 ones, an over-representation of adenine and an underrepresentation of guanine are logically observed). These plots were obtained after the trimming of the data. Analyses represented on figures (B,C,D) were done using MethylKit (default parameters).

CpGs were filtered (CpGs covered by less than 10 reads or having coverage > 99.9 th percentile were excluded) and coverage was normalized (median method). Only information from CpG positions that are in all samples are kept (1,259,111 CpGs). (C) Selection of relevant CpGs is essential for the detection of DMRs. DMR detection could be compromised if CpGs with no reliable methylation state are kept by the method. In this example, the conservation of the two unreliably defined hypermethylated CpGs during DMRs detection process, prevents the detection of hypomethylated DMR, although this DMR exists. Therefore, the p-value threshold for the methylation status of each CpG must be defined to identify statistically relevant methylation information. Filled circles represent hypermethylated cytosines observed between two datasets (e.g. control group versus test group), whereas empty circles represent hypomethylated cytosines. *** significant p-value pval < 0.07); small empty circle, non-significant p-value (pval ≥ 0.07). samples, whereas inter-group variability is greater between E14.5 and E15.5 samples (Supplementary Figure 5B and Supplementary Figure 6A). Many DOCRs and DEGs are identified between E16.5 and P0 stages, compared to other pairwise comparisons. This result is not surprising: since a longer span of time (a few days) separates these two stages (versus only one day for the other stage comparisons), it is plausible that more biological events have occurred that can imply chromatin accessibility. This is consistent with scatterplots and Simple Error Ratio Estimate (SERE) values observed between these samples (Supplementary Figure 5B and Supplementary Figure 6A). Few regions are both identified as DOCRs and DEGs between each stage, compared to the number of observed DOCRs and DEGs taken separately. DOCRs and DEGs are identified using edgeR, without logFC threshold but Benjamini Hochberg p-value adjustment was performed and level of controlled false positive rate was set to 0.05).

Supplementary

(B) Pairwise scatterplots and matrix of pairwise Simple Error Ratio Estimate (SERE) values for ATACseq samples. Comparisons were performed between two successive developmental stages. SERE statistic was used as a similarity index (the more dissimilar the samples are, the higher the SERE value is). Replicates for a given stage, are globally similar, except at E13.5. Indeed, replicates at E13.5 are more dissimilar between each other than to E14.5 samples. This observation could explain why we do not obtain the expected p-values distribution profile for sample comparison (see C) and why identification of DOCRs between these two developmental stages is compromised (see (A)). High interand intragroup variabilities are observed for E15.5 and E16.5 samples. This could explain the lower number of DOCRs identified for these sample comparisons, compared to other developmental stages comparisons (Figure 4A). Scatterplots and SERE values were obtained using SARTools R package. statistic is used as a similarity index (the more dissimilar the samples are, the higher the SERE value is).

For each pairwise comparison, variability between replicates is lower than those of samples from distinct stages. Nevertheless, variability between replicates is higher than those observed for ATACseq data comparisons (Supplementary Figure 5). Inter-and intragroup variabilities are nearly equivalent for E15.5 and E16.5 samples. These similarities can explain the lower number of DEGs identified for these sample comparisons, compared to other developmental stages comparisons (Supplementary Figure 5A). Scatterplots and SERE values were obtained using SARTools R package.

(B) Raw p-values distribution obtained for ATAC-seq samples pairwise comparisons allowing DEG

identification. For all comparisons, p-values distribution profiles are as expected (i.e. a globally uniform distribution with a peak around 0), meaning that DEGs can be detected with edgeR statistic tool These p-value histograms were obtained using ggplot2 R package, with data from SARTools R package.

Supplementary Figure 7. Mono-allelically expressed genes and important neurodevelopmental genes show concomitant perturbations of methylation patterns and alterations of associated gene expression upon PAE. (A) Methylation status of monoallelic expression genes (imprinted genes and clustered Pcdh genes)

and some relevant genes for neurodevelopment Suppl. Figure 1 Duchateau et al. ,Klf7,Pbx1,Sdccag8,Prox1,Ank3,Syt1,Dbnl,Bcl11a,Kif3c,Eml1,Atxn1,Itga1,Farp1,Asap1,Trappc9,App,Itsn1,Arid1b,Anks1,Runx2,Efna5,Man2a1,Ptprm,Tcf4,Ablim1,Dab2ip,Zeb2,Dclk1,Dclk2,Nfib,Ncdn,Spen,Cit,Cux2,Cux1,Dlx5,Creb3l2,Ctnna2,Mgll,Itpr1,Cd9,Sox5,Mboat7,Zfp536,Shank1,Inpp5f,Fgfr2,Shank2,Efnb2,Lig4,Arhgef10,Unc5d,Gnao1,Snx1,Ephb1,Map4 GO:0048699 generation of neurons 1088 54 23.54 2.293537 0.000017 Satb2,Klf7,Pbx1,Sdccag8,Prox1,Ank3,Syt1,Dbnl,Bcl11a,Kif3c,Eml1,Atxn1,Itga1,Farp1,Asap1,Trappc9,App,Itsn1,Arid1b,Anks1,Runx2,Efna5,Man2a1,Ptprm,Tcf4,Ablim1,Dab2ip,Zeb2,Dclk1,Dclk2,Nfib,Ncdn,Spen,Cit,Cux2,Cux1,Dlx5,Creb3l2,Ctnna2,Mgll,Itpr1,Sox5,Zfp536,Shank1,Inpp5f,Fgfr2,Shank2,Efnb2,Lig4,Unc5d,Gnao1,Snx1,Ephb1,Map4 GO:0048468 cell development 1435 64 31.05 2.060957 0.000018

Paqr8, Satb2,Klf7,Pbx1,Prox1,Map7,Jmjd1c,Ank3,Chst11,Syt1,Dbnl,Bcl11a,Kif3c,Atxn1,Hrh2,Itga1,Farp1,Asap1,Cldn5,App,Itsn1,Arid1b,Pacrg,Anks1,Runx2,Efna5,Man2a1,Ptprm,Tcf4,Ablim1,Dab2ip,Zeb2,Dclk1,Dclk2,Nfib,Ncdn,Spen,Cit,Cux2,Rilpl1,Cux1,Dlx5,Creb3l2,Ctnna2,Antxr1,Mgll,Itpr1,Cd9,Sox5,Sipa1l3,Zfp536,Shank1,Inpp5f,Fgfr2,Shank2,Efnb2,Lig4,Arhgef10,Unc5d,Gnao1,Snx1,Myo1e,Ephb1,Map4 GO:0007399 nervous system development 1583 67 34.26 1.955847 0.000043 Satb2,Klf7,Pbx1,Sdccag8,Trp53bp2,Prox1,Ank3,Bcr,Apaf1,Syt1,Dbnl,Bcl11a,Abr,Rbfox3,Kif3c,Eml1,Atxn1,Itga1,Farp1,Mal2,Asap1,Trappc9,Cldn5,App,Itsn1,Arid1b,Anks1,Runx2,Efna5,Man2a1,Ptprm,Tcf4,Ablim1,Dab2ip,Zeb2,Slc1a2,Dclk1,Dclk2,Nfib,Nfia,Ncdn,Spen,Mthfr,Cit,Cux2,Cux1,Dlx5,Creb3l2,Ctnna2,Mgll,Itpr1,Cd9,Sox5,Mboat7,Zfp536,Shank1,Inpp5f,Fgfr2,Shank2,Efnb2,Lig4,Arhgef10,Unc5d,Gnao1,Snx1,Ephb1 The distance between ICR and DNM (bp) was calculated using the ClosestBed tool from the Bedtools suite ("0 "indicates that these regions overlap). The size of the overlap between ICR and DMR was obtained using Interval IntersectBed. The percentage of the ICR corresponding to the DMR was calculated with the following formula: (size of the overlap ICR-DMR / size of the DMR) x 100. Some genes have two names (Zrsr1,Peg13,Mcts2,Nespas,Nap1l5 et Inpp5f_v2), because they integrated into other genes (respectively, Commd1, Trappc9, H13, Gnas, Herc3 et Inpp5f;McCole et al. 2011). ICR, imprinted control region ; chr, chromosome ; pb, base pair.