A ABSTRACT

Previous epidemiological and experimental studies have strongly implicated estrogens in breast cancer risk and Estrogen Receptor (ER), the transcription factor to which estrogen binds, is considered as the major molecular driver of around 70% breast cancers. The importance of the deregulated estrogen signalling is further highlighted by increasing evidence that current chemopreventive and therapeutic strategies that target hormonally responsive breast cancers frequently result in the development of resistance to anti-estrogens and metastatic progression, highlighting the need for understanding the molecular underlying mechanisms. While until recently, ER was believed to act as a stand-alone transcription factor, which can directly bind its motifs in DNA, it is now accepted that ER activity is a complex and dynamic process that requires highly concerted actions of a dozen transcriptional cofactors and various chromatin regulators at DNA. Recent studies focused on characterising ER-associated cofactors and their role in opening the chromatin provided a remarkable insight into transcriptional regulation mediated by ER. However DNA methylation and histone acetylation are poorly understood in the context of ER's dynamic binding.

In this thesis, I combined a cell culture protocol adapted for studying estradiol (E2) deprivation and re-stimulation in stricto sensu in ER-positive breast cancer cells with the latest methylation array, that allowed a genome-wide interrogation of DNA methylation (including a comprehensive panel of enhancers). I further investigated histone acetylation (ChIP-seq) and transcriptome (RNA-seq) after E2 deprivation and re-stimulation to better characterise the ability of ER to coordinate gene regulation. I found that E2 deprivation and re-stimulation result in time-dependent DNA methylation changes and in histone acetylation across diverse genomic regions, many of which overlap with enhancers. Further enrichment analysis of transcription factor (TF) binding and motif occurrence highlights the importance of ER tethering mainly through two partner TF families, AP-1 and FOX, in the proximity of enhancers that are differentially methylated and acetylated. This is the first study that comprehensively characterized DNA methylation at enhancers in response to inhibition and activation of ER signalling. The transcriptome and genome occupancy data further reinforced the notion that ER activity may orchestrate a broad transcriptional programme through regulating a limited panel of critical enhancers.

Finally, the E2 re-stimulation experiments revealed that although the majority of the observed epigenetic changes induced by E2 deprivation could be largely reversed when the cells were re-stimulated we show that DNA hypermethylation and H3K27 acetylation at enhancers as well as several gene expression changes are selectively retained. The partial reversibility can be interpreted as a sign of treatment efficiency but also as a mechanism by which ER activity may contribute to endocrine resistance.

This study provides entirely new information that constitutes a major advance in our understanding of the events by which ER and its cofactors mediate changes in DNA methylation and chromatin states at enhancers. These findings should open new avenues for studying role of the deregulated estrogen signalling in the mechanism underlying the "roots" of endocrine resistance that commonly develops in response to anti-estrogen therapy.

I INTRODUCTION

Breast cancer epidemiology

Despite the continuous efforts for prevention and surveillance, breast cancer (BC) remains the most common cancer in women across the world with an incidence ranging from 26 to 95 cases per 100,000 individuals (Figure 1) [START_REF] Bray | Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries[END_REF]. A small proportion of this incidence is explained by genetic predisposition attributed to family history and ethnicity, but BC is mostly a cancer associated to socioeconomic status and lifestyle practices characteristic of high-income countries [START_REF] Choi | Socioeconomic Inequalities in Cervical and Breast Cancer Screening among Women in Korea, 2005-2015[END_REF][START_REF] Kuchenbaecker | Risks of Breast, Ovarian, and Contralateral Breast Cancer for BRCA1 and BRCA2 Mutation CarriersRisks of Breast, Ovarian, and Contralateral Breast Cancer Among BRCA Mutation CarriersRisks of Breast, Ovarian, and Contralateral Breast Cancer Among BRCA Mutat[END_REF][START_REF] Lundqvist | Socioeconomic inequalities in breast cancer incidence and mortality in Europe-a systematic review and meta-analysis[END_REF]. Indeed, alcohol intake, reproductive factors, high-calorie diet, weight and lack of physical activity, are all factors that increase the risk of BC [START_REF] Schnitt | Breast Cancer[END_REF]. In addition to these, intake of estrogen-progestogen contraceptives and menopausal therapies as well as involuntary exposure to endocrine disruptive chemicals in critical developmental windows can disturb hormonal signalling and contribute to oncogenesis [START_REF] Rodgers | Environmental chemicals and breast cancer: An updated review of epidemiological literature informed by biological mechanisms[END_REF]. Mortality rates are lesser compared to incidence rates and also more evenly spread worldwide while they have been decreasing in high-income countries reflecting a discrepancy in healthcare systems [START_REF] Lundqvist | Socioeconomic inequalities in breast cancer incidence and mortality in Europe-a systematic review and meta-analysis[END_REF].

Adult mammary gland and molecular BC subtypes

The functional unit of breast is the mammary gland, which is build of a network of ducts and alveoli engulfed in stromal cells that are mostly adipocytes (Figure 2) [START_REF] Visvader | Mammary stem cells and the differentiation hierarchy : current status and perspectives[END_REF]. The luminal fraction of the mammary gland is heterogeneous and dynamic as the function of epithelial cells evolves throughout a woman's lifetime [START_REF] Tornillo | ERrrr…where are the progenitors? Hormone receptors and mammary cell heterogeneity[END_REF]. This is reflected by the expression of different molecular receptors, namely estrogen receptor (ER), progesterone receptors (PR) and human epidermal growth factor receptor 2 (HER2), that enable cells to be stimulated by different hormones and growth factors.

In breast tumours, the expression of the above receptors and Ki-67 -a marker of proliferation -at the protein level is mostly measured with antibody-based methods.

Based on these markers, BCs are classified in four main molecular subtypes [START_REF] Koboldt | Comprehensive molecular portraits of human breast tumours[END_REF]: luminal A (ER[+] and/or PR[+], HER-2 [-], low Ki-67), luminal B (ER [+] and/or PR[+], HER-2[+/-], high Ki-67), HER2-enriched (ER [-], PR [-], HER-2 [+]) and basal-like (ER [-], PR [-], HER-2 [-]). More than 70% of BCs are detected with high levels of ER making luminal subtypes the most frequent BC subtype [START_REF] Koboldt | Comprehensive molecular portraits of human breast tumours[END_REF]. In contrast, up to 50% of normal epithelial cells can express ER with only a minority expresses ER at systematically high levels which suggests that BC is more prone to develop from cells expressing ER (Cagnet et al., 2018;Fridriksdottir et al., 2015). Therefore a large body of the breast cancer scientific community has been putting their efforts into elucidating the molecular mechanisms of ER signalling ever since the identification of the protein in the late 1950s (Jensen, 1962).

Therapy and endocrine resistance of hormonal cancers

Because the proliferation of ER-expressing tumour cells is estrogen-dependent, a common approach to treat ER-positive BC is to target ER signalling. In this scope, two main drug strategies were implemented: a) to inhibit ER protein itself and b) to reduce the levels of estradiol, therefore reducing ER activity. The first performing antiestrogenic drug, tamoxifen, was developed already in the late 1960s and was later followed by ICI 182,780 (aka fulvestrant), both of which are still widely used in premenopausal BC patients [START_REF] Cole | A new anti-oestrogenic agent in late breast cancer. An early clinical appraisal of ICI46474[END_REF][START_REF] Thompson | ICI 164,384, a Pure Antagonist of Estrogen-stimulated MCF-7 Cell Proliferation and Invasiveness[END_REF]. In the same period, the first aromatase inhibitors (AIs) appeared as researchers started seeking for drugs that inhibit aromatases, the enzymes that metabolise androgens into estrogens [START_REF] Brodie | Inhibition of Peripheral Aromatization by Aromatase Inhibitors, 4-Hydroxy-and 4-Acetoxy-Androstene-3,17 Dione*[END_REF]. These treatments can be used sequentially or combined with other drugs but also with removal surgery of the BC tumour or ovaries where estrogens are produced.

Unfortunately, regardless of treatment type, overall survival is relatively low as almost 50% of patients present ultimately an intrinsic (from the first treatment) or an acquired (after secondary treatments) endocrine resistance to therapy [START_REF] Clarke | Endocrine resistance in breast cancer -An overview and update[END_REF]. The causes of endocrine resistance are still unclear although many ERmediated mechanisms have been suggested. The loss of ER expression in the tumour was proposed to be among potential mechanisms but this event is observed in only about 10% of patients [START_REF] Ellis | Outcome Prediction for Estrogen Receptor-Positive Breast Cancer Based on Postneoadjuvant Endocrine Therapy Tumor Characteristics[END_REF]. As a matter of fact, ER has an extended network of cofactors with which it cooperates meaning that even if drugs target ER, a large number of other cofactors may be affected in more or less predictable ways.

Therefore identifying ER cofactors and understanding ER signalling is crucial to the understanding of endocrine resistance and to the improvement of BC patient outcome.

Estrogen receptor

ERα and ERβ are members of the nuclear receptor (NR) super-family which contains 48 proteins that act both as transcription factors and as ligand receptors [START_REF] Germain | Overview of Nomenclature of Nuclear Receptors[END_REF]. In this work, only ERα will be considered, that is the most widely expressed form and which we will be referred to as ER for the sake of simplicity.

When dimerised, ER forms a DNA-binding domain that enables to bind directly to DNA, and a ligand-binding domain that binds a large panel of molecules, usually with steroid-like structure [START_REF] Heldring | Estrogen Receptors: How Do They Signal and What Are Their Targets[END_REF]Schwabe et al., 1993). 17β-estradiol (E2) is the most potent female estrogenic hormone whose production starts in the ovaries during puberty with hormone level fluctuation from that point due to menstrual cycles until menopause, when E2 levels drop [START_REF] Kuiper | Comparison of the Ligand Binding Specificity and Transcript Tissue Distribution of Estrogen Receptors α and β[END_REF]. ER comprises also two other main domains, AF-1 and AF-2 (that stand for "activation function") through which ER interacts with other TFs as well as receives post-translational modifications.

ER has also been described to have a ligand-independent activity that is rather associated with developmental functions [START_REF] Caizzi | Genome-wide activity of unliganded estrogen receptor-α in breast cancer cells[END_REF]. Nevertheless, most of ER's activity is ligand mediated and occurs in the nucleus where ER can bind to a multitude of cofactors and regulate gene expression (Figure 3). Until the late 90s, ER was studied mostly in the context of promoters and of direct DNA-binding to its canonical binding motif known as Estrogen Responsive Element (ERE, 5′-GGTCAnnnTGACC-3′) [START_REF] Driscoll | Sequence Requirements for Estrogen Receptor Binding to Estrogen Response Elements[END_REF]. Like for many other TFs and NRs, the emergence of genome-wide sequencing technologies opened the access to the entire ER binding landscape and revealed that ER binds in majority to cis-regulatory elements away from promoters, also known as enhancers [START_REF] Carroll | Genome-wide analysis of estrogen receptor binding sites[END_REF]. (From Heldring et al., 2007)

i. ER cistrome

The term cistrome was first introduced by [START_REF] Lupien | FoxA1 Translates Epigenetic Signatures into Enhancer-Driven Lineage-Specific Transcription[END_REF] and was defined as the "repertoire" alias the ensemble of binding sites of a transcription factor (TF) [START_REF] Lupien | FoxA1 Translates Epigenetic Signatures into Enhancer-Driven Lineage-Specific Transcription[END_REF]. This term is derived from the composition of "cis", a Latin prefix already used in chemistry and meaning "on the same side", and the suffix "ome" that initially comes from Greek (-ωμα) and refers to the sum of elements defined by the first compound of the term. The term sounds particularly familiar as it was inspired by the notorious word "genome" and all its derivatives. In the particular case of ER, the cistrome defines the totality of ER's binding sites which by extension intersects importantly the cistrome of ER cofactors. Indeed, through its association with other cofactors, ER can regulate a number of genes without directly binding on DNA. This way of binding is commonly referred to as "tethered" (Figure 3).

ii. Transcription factors

ER can regulate gene expression positively or negatively through its association with different types of cofactors. The use of protein identification methods (mass spectrometry) and genome binding occupancy assays (ChIP-seq) brought to light new ER cofactors such as FOXA1, GATA3 or GRHL2 to the already known AP-1 or AP2γ [START_REF] Mcpherson | AP2α and AP2γ: a comparison of binding site specificity and trans-activation of the estrogen receptor promoter and single site promoter constructs[END_REF][START_REF] Mohammed | Endogenous Purification Reveals GREB1 as a Key Estrogen Receptor Regulatory Factor[END_REF][START_REF] Theodorou | GATA3 acts upstream of FOXA1 in mediating ESR1 binding by shaping enhancer accessibility[END_REF]. In cell differentiation models, FOXA1 was the first TF to be termed a "pioneer" factor because it was detected as the chronologically first factor to access closed regions of chromatin [START_REF] Boller | Pioneering Activity of the C-Terminal Domain of EBF1 Shapes the Chromatin Landscape for B Cell Programming[END_REF][START_REF] Cirillo | Opening of Compacted Chromatin by Early Developmental Transcription Factors HNF3 (FoxA) and GATA-4[END_REF]. FOXA1 has been extensively studied in the context of ER signalling where it was shown to occupy more than half of ER binding sites [START_REF] Hurtado | FOXA1 is a key determinant of estrogen receptor function and endocrine response[END_REF]. Subsequently, it was shown that FOXA1 resides stably on DNA in time and this regardless of E2-activated ER which suggests that FOXA1 is already settled in chromatin in the context of differentiated BC cells in comparison to developmental models [START_REF] Glont | Comprehensive Genomic Analysis Reveals that the Pioneering Function of FOXA1 Is Independent of Hormonal Signaling[END_REF][START_REF] Swinstead | Steroid Receptors Reprogram FoxA1 Occupancy through Dynamic Chromatin Transitions[END_REF]. Furthermore, DNA binding motif search revealed that FOXA1 as well as GATA3 motifs occur at regular distances of <100bp from ER binding sites, suggesting that the three TFs co-bind in close proximity (Serandour et al., 2013). GATA3 is a TF that also plays a prominent role in ER signalling as it was found to form complexes with ER and FOXA1. Interestingly, up to 16% and 18% of breast invasive carcinomas harbour mutations or copy number alterations in GATA3 and GRHL2, another ER cofactor, against only 6% in FOXA1 and 3% in ER genes [START_REF] Ciriello | Comprehensive Molecular Portraits of Invasive Lobular Breast Cancer[END_REF][START_REF] Jozwik | FOXA1 Directs H3K4 Monomethylation at Enhancers via Recruitment of the Methyltransferase MLL3[END_REF] 2 . The number of possible genetic alterations is at least as large as the number of ER cofactors which makes breast oncogenesis particularly complex to study. Similar to ER, the genes encoding the components of the AP-1 complex, the transcription factor extensively studied in the context of ER signalling and endocrine resistance, rarely carry genetic alterations. AP-1 complex is a heterodimer composed of Fos (FOS, FOSB, FOSL1, FOSL2), Jun (JUN, JUNB, JUND), Maf and ATF subfamilies and can bind to ER through its AF domains [START_REF] Kushner | Estrogen receptor pathways to AP-1[END_REF]. AP-1 was also described as a pioneer factor in the context of glucorticoid receptor (GR) binding in a mouse cell line model [START_REF] Biddie | Transcription Factor AP1 Potentiates Chromatin Accessibility and Glucocorticoid Receptor Binding[END_REF]. AP-1 binds mainly genomic regions away from promoters and is a prerequisite to more than half of GR binding sites. In BC models, it was shown that anti-estrogenic treatments can generate AP-1 binding sites that do not rely on ER binding to DNA [START_REF] He | c-Jun/AP-1 overexpression reprograms ERα signaling related to tamoxifen response in ERα-positive breast cancer[END_REF][START_REF] Webb | The Estrogen Receptor Enhances AP-1 Activity by Two Distinct Mechanisms with Different Requirements for Receptor Transactivation Functions[END_REF]. In addition to ER/AP-1 interaction, JNK proteins that phosphorylate and activate AP-1 have been shown to be up-regulated in endocrine resistance [START_REF] Malorni | Blockade of AP-1 Potentiates Endocrine Therapy and Overcomes Resistance[END_REF], leading to the proposal that AP-1 can act independently from ER and contribute to endocrine resistance mechanisms.

Chromatin and epigenetic marks

Regulation of gene expression is the result of a complex molecular gearing where TFs, epigenetic regulatory factors (ERFs) and epigenetic marks work tightly together on chromatin. Chromatin can be compared to magical yarn where about 2 meters of DNA are compacted in tiny balls to fit a cell nucleus and where TFs and RNA polymerase II (RNApolII) can access many parts without creating knots or unravelling the remaining part. Chromatin evolves dynamically in the short and longterm and this is achieved thanks to changes in DNA methylation and posttranslational modifications of histones that are commonly referred to as epigenetic marks. These dynamic changes are catalysed by ERFs that have either the ability to add ("writers"), to remove ("erasers" or "editors") or simply to recognise ("readers") specific epigenetic patterns on promoters but also on gene bodies and enhancers.

i. Enhancers

Besides revealing the entire binding profiles of TFs, genome-wide sequencing also revealed a large number of enhancers that are in fact much more numerous than promoters [START_REF] Andersson | An atlas of active enhancers across human cell types and tissues[END_REF]. Enhancers play the role of docking platforms for TFs as well as ERFs and enable chromatin to create loops and consequently to bring in spatial proximity regions that are otherwise far in chromosome. While the identification of individual promoters and enhancers can be done with histone and TF binding assays (ChIP-seq), the construction of promoter-enhancer interaction networks requires chromatin conformation assays (Hi-C) [START_REF] Lieberman-Aiden | Comprehensive mapping of long-range interactions reveals folding principles of the human genome[END_REF]. This is also the case for ER where association with a variety of cofactors allows establishing chromatin interactions both in proximal and distant cis-regulatory elements, mediating concerted gene expression in response to hormonal stimulation [START_REF] Dily | Hormonecontrol regions mediate steroid receptor-dependent genome organization[END_REF][START_REF] Fullwood | An oestrogen-receptor-α-bound human chromatin interactome[END_REF][START_REF] Hsu | Estrogenmediated epigenetic repression of large chromosomal regions through DNA looping[END_REF][START_REF] Rafique | Estrogen-induced chromatin decondensation and nuclear re-organization linked to regional epigenetic regulation in breast cancer[END_REF].

While the location of promoters is easily defined by their proximity to the transcription start site (TSS), enhancers can operate from hundreds of thousands of kilobases away, therefore distance cannot be used as a criterion. In addition, multiple enhancers can regulate a given gene and enhancers are cell-type specific which further complicates enhancer identification [START_REF] Heinz | The selection and function of cell type-specific enhancers[END_REF][START_REF] Ko | Functional Enhancers As Master Regulators of Tissue-Specific Gene Regulation and Cancer Development[END_REF]. The location of enhancers is traditionally predicted by the presence of H3K4me1 as well as DNaseI and TF footprints, while enhancer activity is defined through p300 binding, H3K27 acetylation and enhancer RNA (eRNA) transcription [START_REF] Charlet | Bivalent Regions of Cytosine Methylation and H3K27 Acetylation Suggest an Active Role for DNA Methylation at Enhancers[END_REF][START_REF] Franco | Enhancer transcription reveals subtype-specific gene expression programs controlling breast cancer pathogenesis[END_REF][START_REF] Zhu | Predicting enhancer transcription and activity from chromatin modifications[END_REF]. Recently DNA methylation hyper-variability in relatively small genomic regions have been proposed to predict enhancer activity with high efficiency and fidelity [START_REF] Libertini | DNA methylation oscillation defines classes of enhancers[END_REF].

ii. DNA methylation DNA methylation (DNAm) is the covalent addition of a methyl group (-CH 3 ) on the 5 th carbon of a cytosine (5mC) basis that occurs predominantly in the context of a CpG sequence in mammalian cells [START_REF] Smith | DNA methylation: roles in mammalian development[END_REF]. In most cases, DNAm is symmetric which means that a given CpG in an individual cell can be either methylated or unmethylated on either strands. Nevertheless, DNAm is measured most of the time in groups of cells that do not necessarily all have the same methylation status at a specific CpG. For this reason, DNAm -that is expressed as the percentage or the proportion (aka β-value) of the number of methylated sites over the total number of CpG sites -can take any value between 0 and 100% or 0 and 1.

DNAm is edited by two main protein families: DNA methyltransferase (DNMTs) that transfer a methyl group from a donor S-adenosyl methionine onto the cytosine, and Tet methylcytosine dioxygenases (TETs) that contribute to DNA demethylation through a series of oxidation steps [START_REF] Denis | Regulation of mammalian DNA methyltransferases: a route to new mechanisms[END_REF][START_REF] Ecsedi | 5-Hydroxymethylcytosine ( 5hmC ), or How to Identify Your Favorite Cell[END_REF]. In total four DNMTs are known in mammalian cells: DNMT1 ensuring the maintenance of DNAm during cell cycle, DNMT3a and 3B that deposit DNAm de novo, and DNMT3L that has no catalytic domain. DNMTs are often found in the polycomb repressive complexes (PRC) which is in line with their role in repression of gene expression [START_REF] Jin | DNMT1 and DNMT3B modulate distinct polycomb-mediated histone modifications in colon cancer[END_REF]. Less is known about TETs since they were discovered more recently but more and more evidence show they play a significant role in chromatin regulation. The three members TET1, TET2 and TET3 can all catalyse 5mC into derivative forms, namely 5-hydroxymethylcytosine (5hmC), 5carboxylcytosine and 5-formylcytosine, all of which have different chemical properties and therefore can affect chromatin states and TF binding at different degrees.

DNAm has been studied mainly in the context of CpG-dense regions, known as CpG islands (CGI), that frequently coincide with gene promoters [START_REF] Saxonov | A genome-wide analysis of CpG dinucleotides in the human genome distinguishes two distinct classes of promoters[END_REF].

The current widely accepted view is that DNAm in CGIs is negatively correlated to gene expression, meaning that high DNAm levels in promoter are correlated with low expression and vice versa. But DNAm occurs also in other genomic contexts and in CpG-sparse sequences (namely shores, shelves and open sea regions) where it may affect gene regulation in indirect ways. Fore instance, gene-body DNAm, in association to H3K36me3, is believed to prevent RNA polymerase II from binding to DNA, thereby limiting spurious transcription [START_REF] Neri | Intragenic DNA methylation prevents spurious transcription initiation[END_REF] while DNAm in 1 st intron is inversely correlated to gene expression [START_REF] Anastasiadi | Consistent inverse correlation between DNA methylation of the first intron and gene expression across tissues and species[END_REF]. An inverse correlation is also observed in alternative splicing mechanisms where low DNAm of introns is linked to higher frequency of intron inclusion at the mRNA transcription level [START_REF] Kim | Population-dependent Intron Retention and DNA Methylation in Breast Cancer[END_REF][START_REF] Shayevitch | The importance of DNA methylation of exons on alternative splicing[END_REF].

DNAm levels can be stratified in unmethylated regions (UMRs, β-value <0.1), in lowly methylated regions (LMRs, 0.1< β-value <0.5), in partially methylated regions (PMRs, 0.5< β-value <0.85) and in fully methylated regions (FMRs, β-value >0.85), [START_REF] Nothjunge | DNA methylation signatures follow preformed chromatin compartments in cardiac myocytes[END_REF][START_REF] Stadler | DNA-binding factors shape the mouse methylome at distal regulatory regions[END_REF]. Low to intermediate levels of DNAm in CpG-scarce regions have been studied in several different contexts, but globally they have been associated to dynamic chromatin regions and to TF binding [START_REF] Sharifi-Zarchi | DNA methylation regulates discrimination of enhancers from promoters through a H3K4me1-H3K4me3 seesaw mechanism[END_REF]. The methylation status of TF binding motifs that contain CpG site can affect the TF's binding affinity either by attracting or repulsing TF [START_REF] Klose | Genomic DNA methylation: the mark and its mediators[END_REF][START_REF] Yin | Impact of cytosine methylation on DNA binding specificities of human transcription factors[END_REF]. Inversely, TF binding itself can also passively preserve existing DNAm levels simply by blocking the access to DNMTs [START_REF] Hervouet | Specific or not specific recruitment of DNMTs for DNA methylation, an epigenetic dilemma[END_REF][START_REF] Lienert | Identification of genetic elements that autonomously determine DNA methylation states[END_REF]. Nevertheless, changes in DNAm alone cannot explain TF binding at all sites in vivo [START_REF] Maurano | Role of DNA Methylation in Modulating Transcription Factor Occupancy[END_REF], which indicates that the importance of DNAm likely depends on other factors as well, such as post-translational modifications of TFs and histones.

iii. Histone modifications

Nucleosomes are spherical octamers of histones (twice of each H3, H4, H2A and H2B) that are positively charged and enable the negatively charged DNA to coil around them, thereby allowing DNA compaction in the nucleus. Many residues in the N-and C-terminal tails of histones are subjected to post-translational modifications (PTMs) that play a crucial role in condensing and loosening up chromatin [START_REF] Vanzan | Mechanisms of Histone modifications[END_REF]. Furthermore, nucleosomes directly obstruct the access to RNApolII and need to be evicted or reorganized during transcription initiation and elongation [START_REF] Workman | Alteration of Nucleosome Structure as a Mechanism of Transcriptional Regulation[END_REF].

There are currently hundreds of possible PTMs of histone tails (including covalent addition of methyl, acetyl, citrullin and other moieties). Among the most studied PTMs are H3K4me1, H3K4me3, H3K27ac and H3K27me3.

H3K4me1 and H3K4me3

Several histone methyltransferases (HMTs, such as SET1A, SET1B and MLL1 to MLL4) exist that can transfer up to three methyl groups from S-adenosyl methionine donors onto lysine and arginine residues of H3 and H4 tails [START_REF] Clarke | Protein methylation[END_REF].

Trimethylation of H3K4 (H3K4me3) by MLL1/2-containing complexes is strongly associated to promoters and is necessary for their activation [START_REF] Wang | Global Analysis of H3K4 Methylation Defines MLL Family Member Targets and Points to a Role for MLL1-Mediated H3K4 Methylation in the Regulation of Transcriptional Initiation by RNA Polymerase II[END_REF]. The addition of one methyl group (H3K4me1) is mostly catalysed by MLL3 and MLL4

and their presence is associated to enhancers [START_REF] Yan | Histone H3 Lysine 4 methyltransferases MLL3 and MLL4 Modulate Long-range Chromatin Interactions at Enhancers[END_REF]. In ER-expressing cells H3K4me1 deposition at enhancers largely depends on FOXA1 binding and is necessary to ER binding at these locations [START_REF] Jozwik | FOXA1 Directs H3K4 Monomethylation at Enhancers via Recruitment of the Methyltransferase MLL3[END_REF]. Both H3K4me1 and H3K4me3 usually stretch along large genomic windows, but unlike H3K4me3, H3K4me1 is not sufficient to activate an enhancer.

H3K27ac

To activate further enhancers and promoters, the 27 th lysine residue of H3 needs to be acetylated (H3K27ac) [START_REF] Creyghton | Histone H3K27ac separates active from poised enhancers and predicts developmental state[END_REF]. Histone acetylation is mediated by the recruitment of histone acetyltransferases (HATs) and more specifically, in mammalian cells acetyl is deposited by the catalytic bromodomain of CBP and p300 [START_REF] Raisner | Enhancer Activity Requires CBP/P300 Bromodomain-Dependent Histone H3K27 Acetylation. Cell Rep[END_REF]. H3K27ac marks appear in narrow windows and most of the time defining the edges of TF-bound regions. Although H3K4me1 and MLL3/4 are not sufficient to activate enhancers, their presence is essential to P300 recruitment [START_REF] Wang | Enhancer priming by H3K4 methyltransferase MLL4 controls cell fate transition[END_REF]. In breast cancer models, ER activation can lead to the recruitment of both p300 and AP-1 and highly co-localises with FOXA1 and GATA3 binding [START_REF] Jeffy | An estrogen receptoralpha/p300 complex activates the BRCA-1 promoter at an AP-1 site that binds Jun/Fos transcription factors: repressive effects of p53 on BRCA-1 transcription[END_REF][START_REF] Murakami | Dynamic assembly and activation of estrogen receptor α enhancers through coregulator switching[END_REF][START_REF] Rosell | Complex formation and function of estrogen receptor α in transcription requires RIP140[END_REF]. Histone deacetylases (HDACs) "erase" the acetyl group from histones and consequently lead to a decrease of distance between nucleosomes and ultimately to a condensation of chromatin structure [START_REF] Zentner | Regulation of nucleosome dynamics by histone modifications[END_REF].

Repressive PTMs

EZH2, that is part of the PRC2 complex along with SUZ12 and EED, is the only histone methyltransferase known to methylate H3K27me3 [START_REF] Comet | Maintaining cell identity: PRC2-mediated regulation of transcription and cancer[END_REF]. It can be found on promoters and enhancers and is not mutually exclusive to H3K27ac.

PRC2 is a large complex that may also contain HDACs and DNMTs and can therefore lead to an important gene expression down-regulation [START_REF] Hernández-Muñoz | Association of BMI1 with polycomb bodies is dynamic and requires PRC2/EZH2 and the maintenance DNA methyltransferase DNMT1[END_REF][START_REF] Marchesi | Targeting Enhancer of Zeste Homolog 2 as a promising strategy for cancer treatment[END_REF]. H3K9me3 is another repressive mark that is partially redundant to H3K27me3 since they are often superimposed. Both H3K9me3 and H3K27me3 can be read by heterochromatin protein 1 (HP1), a protein that plays a role in the constitution of heterochromatin [START_REF] Jamieson | Loss of HP1 causes depletion of H3K27me3 from facultative heterochromatin and gain of H3K27me2 at constitutive heterochromatin[END_REF][START_REF] Van De Werken | Paternal heterochromatin formation in human embryos is H3K9/HP1 directed and primed by sperm-derived histone modifications[END_REF].

iv. Chromatin states

Dynamic chromatin states are needed to establish tissue-and temporal-specific transcriptional programs in response to environmental and endogenous stimuli.

Initially, chromatin was approximately distinguished into euchromatic and heterochromatic domains based on the intensity of staining observed in microscopy [START_REF] Elgin | Heterochromatin and gene regulation in Drosophila[END_REF]. Ever since, the integration of genome-wide epigenetic, transcription and chromatin-conformation data has revealed a complex and multilayer regulation of chromatin.

Conventionally, DNAm is associated to heterochromatin. Nevertheless, intermediate DNAm is both positively and negatively correlated to H3K27ac at enhancers when considered in large genomic bins [START_REF] Charlet | Bivalent Regions of Cytosine Methylation and H3K27 Acetylation Suggest an Active Role for DNA Methylation at Enhancers[END_REF]. This lack of correlation could be accounted to intermediate DNAm levels, which include a large spectrum of DNAm levels (LMRs + PMRs). DNA methylation levels should also be considered at the single CpG level rather than in regions, especially in the context of enhancers. For instance, hypomethylation has been associated with large nucleosome-depleted regions (NDRs) while punctual methylation of histone-bound fraction of DNA seems to contribute to a more open structure [START_REF] Jimenez-Useche | DNA methylation Regulated Nucleosome Dynamics[END_REF][START_REF] Taberlay | Reconfiguration of nucleosome-depleted regions at distal regulatory elements accompanies DNA methylation of enhancers and insulators in cancer[END_REF].

The cross-talk between DNAm and histone marks is believed to be mediated by complexes that contain ERFs that are sensitive to DNA and/or histone methylation.

For instance, MLL1/2, that deposit H3K4me3, contains a zinc finger CXXC domain that has high affinity for unmethylated DNA and is often found on promoters [START_REF] Birke | The MT domain of the proto-oncoprotein MLL binds to CpGcontaining DNA and discriminates against methylation[END_REF][START_REF] Wang | Global Analysis of H3K4 Methylation Defines MLL Family Member Targets and Points to a Role for MLL1-Mediated H3K4 Methylation in the Regulation of Transcriptional Initiation by RNA Polymerase II[END_REF]. Inversely, MLL3/4 lack the CXXC domain and this can be reflected by the heterogeneous levels of DNAm found at enhancers [START_REF] Van Nuland | Quantitative Dissection and Stoichiometry Determination of the Human SET1/MLL Histone Methyltransferase Complexes[END_REF][START_REF] Sze | MLL3/MLL4/COMPASS Family on Epigenetic Regulation of Enhancer Function and Cancer[END_REF]. CTCF, that contributes to the definition of topologically associated domains (TADs) and to insulation mechanisms, was believed for a long time to be mutually exclusive to DNAm but recent integration of 5hmC data shows that the affinity of CTCF binding is more subtle than initially thought [START_REF] Dixon | Topological domains in mammalian genomes identified by analysis of chromatin interactions[END_REF][START_REF] Wiehle | DNA (de)methylation in embryonic stem cells controls CTCF-dependent chromatin boundaries[END_REF].

The global picture of chromatin dynamics proves to be highly complex, especially regarding the role of DNAm. Nevertheless, DNAm, histone marks and TF binding in specific genomic contexts can be combined to define four chromatin states [START_REF] Caglio | Keep Them Close: PRC2 Poises Enhancer-Promoter Interactions at Anterior Neuronal Genes[END_REF][START_REF] Calo | Modification of enhancer chromatin: what, how, and why?[END_REF] applicable to both promoters and enhancers but are more described in the context of the latter (Figure 4):

• Active enhancers carrying H3K4me1 as well as H3K27ac marks and exhibiting binding of MLL3/4 and p300 as well as producing eRNA.

• Primed enhancers maintaining H3K4me1 and showing lower, albeit not null, levels of H3K27ac, MLL3/4 and p300 and do not produce eRNA.

• Poised enhancers are similar to primed but in addition carry the repressive mark H3K27me3.

• Repressed enhancers showing no methylation on H3K4 and having high levels of H3K27me3 and/or H3K9me3.

ER, DNA methylation and breast cancer

BC, similar to a vast majority of cancers, can occur following epigenetic alterations and is characterised by a global hypomethylation and a focal hypermethylation [START_REF] Feinberg | Hypomethylation distinguishes genes of some human cancers from their normal counterparts[END_REF][START_REF] Herceg | Genetic and epigenetic alterations as biomarkers for cancer detection, diagnosis and prognosis[END_REF][START_REF] Szyf | DNA methylation and breast cancer[END_REF]. The fruition of large-scale sequencing efforts, spearheaded by The Cancer Genome Atlas (TCGA), and the Encyclopedia of DNA Elements (ENCODE) project has substantially accelerated our understanding of cancer (epi)genome. By matching high throughput data from hundreds of human tumour and normal (non-tumour) samples with similar data from cell line models, researchers have been gradually bridging the gap between clinical and basic research. For BC, the differences in methylome were shown to be largely driven by the expression of hormonal receptors that are used to define the classical molecular subtypes of BC [START_REF] Fernandez-Jimenez | Lowly methylated region analysis identifies EBF1 as a potential epigenetic modifier in breast cancer[END_REF][START_REF] Holm | Molecular subtypes of breast cancer are associated with characteristic DNA methylation patterns[END_REF]. For instance, the basal-like (triple-negative) subtype was shown to be the most hypomethylated BC subtype [START_REF] Koboldt | Comprehensive molecular portraits of human breast tumours[END_REF], although the underlying mechanism remain unknown.

In 2014, Ung et al. overlapped ER binding sites from ER-positive cell line MCF-7 with BC methylome data (Illumina Infinium HumanMethylation450 BeadChip) from TCGA and found that while ER-positive cancers were globally hypermethylated in comparison to ER-negative subtypes, DNAm at ER binding sites was substantially lower in ER-positive tumours [START_REF] Ung | Effect of estrogen receptor ERalpha binding on functional DNA methylation in breast cancer[END_REF]. A similar observation was subsequently reported for DNAm at FOXA1 binding sites that was found considerably lower in tumours expressing both FOXA1 and ESR1 [START_REF] Ciriello | Comprehensive Molecular Portraits of Invasive Lobular Breast Cancer[END_REF]. This is strongly in line with in vitro results showing that knock-down of FOXA1 in MCF-7 cells leads to an increased methylation on CpGs included within ±500bp away from FOXA1 peaks [START_REF] Zhang | Nucleation of DNA repair factors by FOXA1 links DNA demethylation to transcriptional pioneering[END_REF]. Interestingly, it was earlier suggested in a different cell line model that FOXA1 actually triggers an active DNA demethylation [START_REF] Sérandour | Epigenetic switch involved in activation of pioneer factor FOXA1-dependent enhancers[END_REF]. Binding of ER, FOXA1, GATA3 and FOS was found significantly enriched at functional enhancers that had lower DNAm levels in ER-positive compared to ER-negative BC [START_REF] Fleischer | DNA methylation at enhancers identifies distinct breast cancer lineages[END_REF]. This raises the possibility that the specific mechanisms are conserved between cell line models and clinical samples. Finally, long term E2 deprivation in MCF-7 was also shown to lead to a global hypermethylation, a small portion of which overlapped with ER-bound enhancers, based on which a mechanism playing a role in endocrine resistance was proposed [START_REF] Stone | DNA methylation of oestrogen-regulated enhancers defines endocrine sensitivity in breast cancer[END_REF]. Together, these results strongly indicate that ER and its cofactors may play a role into the maintenance of lower DNAm levels at enhancers.

Hypothesis and aims

Based on the current knowledge, I hypothesised that ER positively regulates gene expression by shaping the chromatin and more specifically by keeping DNA methylation low and histone acetylation high at ER binding sites (Figure 5). In the presence of E2, ER and its cofactors would occupy DNA and passively maintain these epigenetic marks by preventing repressive complexes to access these sites through simple steric hindrance. In the absence of E2, the ER-dependent complexes would not form permitting repressive complexes containing DNMTs and/or HDACs to change epigenetic marks, therefore to switch to repressed chromatin states and to downregulate gene-expression. Finally, I also put forward the hypothesis that if ER is reactivated by its ligand, the changes in epigenetic marks may not be reversible. The removal of E2 may mimic the mechanism of aromatase inhibition and therefore this project could to contribute to disentangling further the mechanisms underlying endocrine resistance.

The specific aims are:

• To show that modulating estradiol-mediated ER activity has an impact on genome-wide DNAm and histone acetylation.

• To address the functional impact of these epigenetic changes.

• To assess the reversibility of these epigenetic changes.

M MATERIAL AND METHODS

Cell culture and treatments i. Rationale behind a new stricto sensu E2-dependent cell culture protocol

In order to model the effect of steroid hormones in cell culture, many previous studies used the E2-deprivation/re-stimulation approach [START_REF] Ryan | Moving Toward Integrating Gene Expression Profiling Into High-Throughput Testing: A Gene Expression Biomarker Accurately Predicts Estrogen Receptor α Modulation in a Microarray Compendium[END_REF]. Briefly, an ERpositive cell line, typically MCF-7, is deprived for several days from E2 after which the cells are re-stimulated with E2 and the endpoint of interest is measured. This system allows exploring two distinct phenomena: 1) the effect of E2 removal from the medium, alias deprivation and 2) the effect E2 re-introduction, alias re-stimulation, both of which are mediated by ER. Nevertheless, if one aims to study the potential reversibility of epigenetic changes resulting from E2 deprivation subsequently followed by E2 re-stimulation the existing protocols present certain limitations. In the following section, I address some of these limitations, and I present a methodological approach that circumvent these shortcomings.

The very extensive use of E2 deprivation/re-stimulation protocol has the consequence of increasing number of parameters that vary within and between laboratories. The origin of the MCF-7 cell line, the composition of growth medium, the quality of charcoal stripping of serum, the duration of E2-deprivation or the E2 concentration used for re-stimulation, are some examples that can affect the estrogenic growth doseresponse and hence contribute to reproducibility issues [START_REF] Kleensang | Genetic variability in a frozen batch of MCF-7 cells invisible in routine authentication affecting cell function[END_REF][START_REF] Sikora | Endocrine Response Phenotypes Are Altered by Charcoal-Stripped Serum Variability[END_REF]. Firstly, it is important to know that although there are multiple cell culture alternatives, the main E2 deprivation protocol consists into passing ERpositive cells from medium containing standard serum to medium containing charcoal-stripped serum (CSS). Charcoal stripping is a chemical process that removes lipophilic compounds such as steroid hormones (E2, progestogens, androgens), several growth factors and some cytokines. E2 re-stimulation is subsequently performed by adding E2 on previously deprived cells. This method has been applied in a multitude of studies, some of which have focused on the effect of long-term E2 deprivation in MCF-7 on multiple readouts with the aim to characterise the mechanisms of endocrine resistance [START_REF] Aguilar | Biological reprogramming in acquired resistance to endocrine therapy of breast cancer[END_REF][START_REF] Jeng | Estrogen receptor expression and function in long-term estrogen-deprived human breast cancer cells[END_REF][START_REF] Santen | Adaptation to estradiol deprivation causes upregulation of growth factor pathways and hypersensitivity to estradiol in breast cancer cells[END_REF][START_REF] Stone | DNA methylation of oestrogen-regulated enhancers defines endocrine sensitivity in breast cancer[END_REF]. Long-term deprivation mimics the decrease of E2 blood levels that occurs in women with hormono-dependent breast cancer that either underwent oophorectomy or were treated with aromatase inhibitors [START_REF] Santen | Adaptation to estradiol deprivation causes upregulation of growth factor pathways and hypersensitivity to estradiol in breast cancer cells[END_REF]. In vitro, the protocol of long-term E2-deprivation (LTED) consists into culturing MCF-7 in medium containing CSS for durations ranging from 4 to 24 months and to compare the readout to that of MCF-7 maintained in regular-serum medium. This protocol design is problematic, given that the control group is cultured in medium containing regular serum containing all components previously mentioned.

This means that the comparison of the outputs does not reveal effects of E2 deprivation alone, but cumulates indistinguishably the effect of other steroid hormones and growth factors that affect MCF-7 cell processes [START_REF] Azeez | Progesterone regulates the proliferation of breast cancer cells -in vitro evidence[END_REF][START_REF] Macedo | Role of Androgens on MCF-7 Breast Cancer Cell Growth and on the Inhibitory Effect of Letrozole[END_REF]. Although the removal of E2 could be the main cause of the observed changes, it is not possible to exclude that the final LTED phenotype is the result of synergistic deprivation. This observation is of minor impact in the context of studies using short E2 deprivation -usually 3 days alias one passage -that are mostly aimed at elucidating ER dynamics at a mechanistic level.

In addition to the above limitation, one has to consider the capacity of cancer cell lines, in spontaneously evolving into different subclones [START_REF] Hara | Changes of proliferative activity and phenotypes in spontaneous differentiation of a colon cancer cell line[END_REF][START_REF] Qu | Evaluation of MCF10A as a reliable model for normal human mammary epithelial cells[END_REF]. Particularly MCF-7 is known to be a heterogeneous cell line that contains endocrine resistant subpopulations [START_REF] Devarajan | Human breast cancer MCF-7 cell line contains inherently drug-resistant subclones with distinct genotypic and phenotypic features[END_REF][START_REF] Kleensang | Genetic variability in a frozen batch of MCF-7 cells invisible in routine authentication affecting cell function[END_REF].

These endocrine resistant clones often lose their dependence to ER and do not respond to E2 modulation. To avoid the issue of spontaneous specification, it has been suggested to perform assays within less than 20 to 25 cell passages [START_REF] Nelson | Detection of Endogenous Selective Estrogen Receptor Modulators such as 27-Hydroxycholesterol BT -Estrogen Receptors: Methods and Protocols[END_REF].

Keeping in mind the considerations described above, the main concerns for designing the experiments were:

1. Focusing strictly on the effect of E2 deprivation and not the effect of charcoal stripped serum 2. Having a control that can be commonly compared to both E2 deprived and E2 re-stimulated condition to study the putative reversibility of the effects due to deprivation 3. Culturing cells long enough to obtain observable effects, but not too long to avoid loss of estrogen response In order to respect the above constraints and since removing E2 alone from CSS is not chemically possible, I opted for an inversed approach: to deprive cells from all components with the use of CSS and to maintain only E2 (and its solvent) in the medium by adding it on a regularly basis. To implement this, as soon as MCF-7 cells were thawed, they were placed in culture in phenol-red-free medium (another weak ER agonist) supplemented with CSS as well as a daily input of E2. They were maintained and subcultured in this medium until the beginning of the experiments, which usually took approximately 2 weeks and 4 passages. In this way, I made the assumption that by the beginning of the experiments, most changes resulting from the absence of charcoal-stripped components would be evened out, or at least minimised.

To specifically address the deprivation of E2 stricto sensu, MCF-7 were cultured in parallel in two groups; a control group (CTR) cultured in the above mentioned medium with E2, and a E2-deprived group (E2D), cultured in the same conditions, only lacking E2 (Figure 6). As for assessing the reversibility of the effects deriving from the deprivation, a third group of cells was simply re-stimulated with E2 (ReSt) after being deprived from it for a certain period. The presence of E2 being the only difference between CTR, E2D and ReSt groups, the comparison of all groups becomes possible and we can study mechanisms that are strictly ER-mediated. Concerning the duration of the treatment, we know from experience that changes in expression can be measured in a matter of hours while there is less evidence about the time needed for DNAm changes to be detectable. As a starting point, I followed the most employed duration of E2 deprivation that is used for mRNA expression measuring, that is 3 days of deprivation [START_REF] Stanislawska-Sachadyn | Transcriptomic Effects of Estrogen Starvation and Induction in the MCF7 Cells[END_REF]. For the E2 re-stimulation, the durations in literature vary from 30 minutes to 48h and I chose 24h. With these conditions, we ran a pilot 850k array (see Material & Methods section later).

For the sake of readability, the details of the pilot analysis will not be shown but briefly, no significant differences in DNA methylation were found using the canonical linear regression model that is used in the context of methylation array analysis [START_REF] Aryee | Minfi: A flexible and comprehensive Bioconductor package for the analysis of Infinium DNA methylation microarrays[END_REF]. Nevertheless, with the use of the more sensitive beta regression [START_REF] Ferrari | Beta Regression for Modelling Rates and Proportions[END_REF], 1885 DMPs were detected of which only 3% had differential methylation >10% after 3 days of E2 deprivation. The beta regression though proven to be also less specific as only 6 out of 10 tested DMPs were technically validated. Based on this pilot assay, I concluded that more samples would increase slightly the robustness of statistics but especially that a longer E2 deprivation and re-stimulation would increase the amplitude of DNA methylation changes.

Therefore, I increased the total duration of the assay to 14 days and the E2 restimulation particularly lasted for 10 days after 4 days of deprivation (Figure 7). All further experiments were based on this core design. ii. E2 deprivation and re-stimulation of MCF-7 MCF-7 HTB-22® cells were obtained from the American Type Culture Collection and all assays in this thesis were performed within 12 passages. Fresh E2 stock solutions at 1µM in 100% DMSO were regularly prepared, they were stored at -20°C in aliquots <50µL and were single-used daily to ensure a continuous availability of E2 in the medium. MCF-7 cultured for the E2 deprivation/re-stimulation and the siRNA assays were expanded for two weeks in phenol-red-free DMEM/F-12 supplemented with 10% (v/v) charcoal-stripped fetal bovine serum (csFBS, Sigma-Aldrich), 1% nonessential amino acids, sodium pyruvate and penicillin/streptomycin (Life Technologies) and with a daily addition of 10nM E2 and 0.1% DMSO (Sigma-Aldrich) at 37°C, in a humidified and 5% CO2-enriched atmosphere (E2+CSS medium). For the assay itself, control cells (CTR) were cultured continuously for 14 days in the above-mentioned conditions, while E2-deprived cells (E2D) were cultured in the same conditions only lacking E2. The re-stimulated cells (ReSt) were deprived from E2 for 4 days and re-stimulated for the 10 following. Cells were cultured in triplicates and collected at day 0, 4 and 14 for downstream analysis.

iii. ER inhibition and re-activation by antagonists

Another way to validate the effect of sequential ER inactivation and re-activation, is to apply and later remove ER inhibitors from the culture conditions. In order to keep ER deprivation and ER inhibition protocols comparable, I wanted initially to keep the same basic medium composition (phenol-red-free medium, CSS and E2) and add the ER antagonist in higher concentration to E2. To test if cells were responsive to these inhibitors, MCF-7 were exposed for 4 days to two well known antagonists separately at 1µM, fulvestrant or ICI 182,780 (ICI), a selective estrogen receptor degrader (SERD), and (Z)-4-hydroxytamoxifen (shortly tamoxifen, TAM), a selective estrogen receptor modulator (SERM) [START_REF] Cole | A new anti-oestrogenic agent in late breast cancer. An early clinical appraisal of ICI46474[END_REF][START_REF] Wakeling | A Potent Specific Pure Antiestrogen with Clinical Potential[END_REF]. During this exposure, the drugs were renewed twice, along with E2. For both drugs, expression of ER downstream targets is expected to change. After 4 days, I measured the expression levels of GREB1, a gene positively regulated by ER, as a proxy of ER activity (see Material & Methods section later). Surprisingly, I discovered that under the E2+CSS culture conditions, the expression of GREB1 remained unchanged (Figure 8A). After several tests changing concentrations of components, I abandoned the idea of culturing cells in E2+CSS medium and I switched to using medium containing regular FBS. Interestingly, this led in the expected result, meaning the down-regulation of ER target GREB1 in both ICI and TAM treatments (Figure 8B). I therefore performed downstream inhibition assays using medium containing regular FBS as follows. MCF-7 meant for exposure to the inhibitors were maintained for two weeks in DMEM supplemented with 10% fetal bovine serum (Eurobio), 1% nonessential amino acids, sodium pyruvate, penicillin/streptomycin and 0.1% DMSO.

The control cells (CTR) were cultured in the condition above, inhibited cells (ICI) were cultured in continuous presence of 1µM ICI and re-activated cells (ReAc) were cultured in presence of the inhibitor for 4 days followed by control conditions for another 10 days (Figure 9). Cells were cultured in triplicates and collected at day 0, 4 and 14 for downstream analysis.

Molecular biology i. DNA and RNA extraction

For all assays involving DNA and RNA, total DNA and RNA were extracted using

AllPrep DNA/RNA Mini kit (Qiagen). DNA was quantified by Qubit dsDNA HS kit (Invitrogen) while quality and quantities of RNA were assessed on a ND-8000 spectrophotometer (Nanodrop). Eluted DNA was stored at -20°C and RNA at -80°C until further use.

ii. Bisulfite conversion

In this work, DNA was used for the only purpose of measuring CpG methylation. To quantify the percentage of methylated cytosines, we used bisulfite conversion that causes unmethylated cytosines to be converted into uracil bases while methylated cytosines remain as such [START_REF] Vaissière | Quantitative analysis of DNA methylation after whole bisulfitome amplification of a minute amount of DNA from body fluids[END_REF]. For samples processed on the Infinium MethylationEPIC array (Illumina), 600 ng of DNA were converted with the EZ DNA Methylation kit (Zymo Research) 250 ng of which were used for the hybridisation. Samples processed by pyrosequencing were converted with a EZ DNA Methylation Gold kit (Zymo Research). In both cases DNA was eluted in 20 µL and stored at -20°C until further use.

iii. Pyrosequencing

Pyrosequencing was performed as previously established in our team [START_REF] Vaissière | Quantitative analysis of DNA methylation after whole bisulfitome amplification of a minute amount of DNA from body fluids[END_REF].10 to 25 ng of modified DNA was amplified by PCR in a final volume of 50 μL. 8 μL of PCR reaction was analysed on a 2% agarose gel to confirm efficacy of the reaction. The remaining volume of modified DNA was used for the pyrosequencing assay itself. One of the primers used for each PCR reaction carried a 5'-biotinyl modification. Briefly, the PCR amplicons were bound onto streptavidin-coated sepharose beads (Amersham-GE Healthcare), which attract the biotinylated strands.

Pyrosequencing was carried out in accordance with the manufacturer's protocol using PyroGold Reagent kit (Qiagen) on the vacuum-based workstation (PSQ 96MA, Biotage). The methylation levels of the CpGs of interest were evaluated by converting the resulting pyrograms to numerical values based on their peak heights. Primers for PCR and sequencing primers are presented in Table 1. iv. Chromatin immunopricipatation MCF-7 were cultured in duplicates in CTR, E2D and ReSt conditions as mentioned above and were formaldehyde-fixed at 1% for 8 minutes. The equivalent of 2 million cells were sonicated and processed with the iDeal ChIP-seq kit for histones as per the manufacturer's protocol (Diagenode). Chromatin from 1 million cells was immunoprecipitated with 2.5 µg of anti-H3K27ac, anti-H3K4me1 and anti-H3K27me3 antibodies (ab4729, ab8895, ab6002, Abcam). ChIPed and input DNA was extracted using the iPure kit (Diagenode). The same sonication conditions were used as for ChIP-seq (see below for details).

v. mRNA retro-transcription and quantitative RT-and ChIP-PCR

Reverse transcription reactions were performed using MMLV-RT enzyme (Invitrogen) and random hexamers on 500 ng of total RNA per reaction as per the manufacturer's protocol. Quantitative RT-PCR was done in triplicate for each condition. The PCR conditions used were: 95°C 5min, [95°C 15 s, 61°C 30 s] × 40 cycles, 95°C 1 min. The assays were performed using SYBR Green qPCR mix and a CFX96 RealTime PCR Detection System (Bio-Rad Laboratories). Expression levels were calculated using the 2-ΔCt method [START_REF] Livak | Analysis of relative gene expression data using real-time quantitative PCR and the 2-ΔΔCT method[END_REF]chip. Results from ChIP were calculated as a percentage of input (previously diluted 1:10). The primers used for RT and ChIP-qPCR are reported in Table 2 andTable 3

OMICs and data analysis i. Methylation bead arrays

250 ng of bisulfite converted DNA were hybridised on an Infinium MethylationEPIC array, referred to as 850k from this point (Illumina). Processing and normalization of data followed by identification of differentially methylated positions (DMPs) was performed as described earlier using Bioconductor package minfi [START_REF] Degli Esposti | Unique DNA methylation signature in HPV-positive head and neck squamous cell carcinomas[END_REF][START_REF] Fortin | Preprocessing, normalization and integration of the Illumina HumanMethylationEPIC array with minfi[END_REF]. Annotation of probes was done with the use of IlluminaHumanMethylationEPICanno.ilm10b2.hg19 [START_REF] Hansen | IlluminaHumanMethylationEPICanno.ilm10b4.hg19: Annotation for Illumina's EPIC methylation arrays[END_REF]. To determine the effect of E2D deprivation on DNAm, we used a linear regression model (limma)

where we compared E2D (n=6) treatment to CTR (n=9) while adjusting for timepoint and applying a differential methylation cutoff of 10% (false discovery rate, FDR <0.05; Δβ ≥10%) [START_REF] Ritchie | limma powers differential expression analyses for RNA-sequencing and microarray studies[END_REF]. Then, to test the reversibility of this effect, we applied a linear regression model between CTR and ReSt samples at d14. Drawing of QQ-plots and calculation of lambda (inflation factor) were done as per [START_REF] Yang | Genomic inflation factors under polygenic inheritance[END_REF]. FDb.InfiniumMethylation.hg19 and ChIPseeker R/Bioconductor packages were used for annotating the DMPs [START_REF] Triche | InfiniumMethylation.hg19: Annotation package for Illumina Infinium DNA methylation probes[END_REF][START_REF] Yu | ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization[END_REF].

ii. ChIP-sequencing

The library preparation requires an initial DNA fragment size ranging 200-800bp.

Because large fragments persisted with the increase of sonication cycles while small fragments became further fragmented, I decided to sonicate chromatin 14x (30sec ON:OFF, Bioruptor Pico) and to perform a fragment size selection downstream, after the ChIP (Figure 10). The very first step of the TruSeq library preparation consists into performing a small fragment clean-up using AMPure beads. After various tests, I found that a ratio of 0.5x AMPure beads (volume beads:DNA) successfully got rid of large-sized fragments (>800bp). I incorporated this extra step at the very beginning of library preparation and proceeded as per the manufacturer protocol with 10 ng of initial cleaned-up ChIPed DNA. For input, equal amounts from duplicate CTR at d4 and d14 were pooled to generate the input (1.67µL from each, 10µL in total diluted 1:10), which was de-cross-linked and further treated like the remaining ChIPed samples. In between the red lines is indicated the fragment range that will be finally sequenced (add ± 120bp for adapters and end-repair)

H3K27ac ChIP and input libraries were prepared using the TruSeq ChIP Sample Prep kit and we performed paired-end sequencing with the NextSeq 500/550 High Output Kit v2.5 (150 Cycles) (Illumina). Fastq files were trimmed using TrimGalore for adapters and a minimum quality of Q >30 (Babraham Bioinformatics3 ), then mapped on hg19 using BWA and further processed as previously described [START_REF] Jozwik | FOXA1 Directs H3K4 Monomethylation at Enhancers via Recruitment of the Methyltransferase MLL3[END_REF][START_REF] Li | Fast and accurate short read alignment with Burrows-Wheeler transform[END_REF]. Mapped read numbers ranged between 125 and 174 million, which is higher than the amount of reads recommended by ENCODE.

Narrow peaks were then called with MACS2 (2.1.1.20160309) against input and were filtered ad hoc for <1.5kb width, q <0.01 and fold-change from input >5 [START_REF] Zhang | Model-based Analysis of ChIP-Seq (MACS)[END_REF]. For downstream analysis we worked with irreproducible-discovery-rate (IDR) thresholded peaks. Differentially bound (DB) regions were identified with DiffBind between different groups at d14 and only DB peaks with |log 2 FC| >1 were considered [START_REF] Stark | DiffBind: differential binding analysis of ChIP-Seq peak data[END_REF].

iii. RNA-sequencing

Total RNA was extracted along with DNA using the AllPrep DNA/RNA Mini kit (Qiagen). One µg of RNA was used for TruSeq RNA library preparation (Illumina)

following the manufacturer's protocol and single-end sequencing was performed on

HiSeq 4000 yielding at least 20M 50bp-reads for each sample. Following quality control, reads were trimmed at 30bp and mapped on hg19 with hisat2 with mapping efficiencies of ≥87% [START_REF] Kim | HISAT: a fast spliced aligner with low memory requirements[END_REF]. Read counts were generated from bam files using htseq-count function and genes with ≤10 counts-per-million in at least two samples were eliminated [START_REF] Anders | HTSeq--a Python framework to work with high-throughput sequencing data[END_REF]. Differentially expressed genes were identified between CTR and E2D adjusting for time and between CTR and ReSt at d14 using edgeR (FDR <0.05, |log 2 FC| >1) [START_REF] Robinson | edgeR: a Bioconductor package for differential expression analysis of digital gene expression data[END_REF].

Other data analysis i. De novo motif analysis

A de novo TF motif enrichment analysis was performed after enlarging the genomic window to 250 bp around DMPs and to 1 kb for DB H3K27ac, using the findMotifsGenome.pl function from HOMER v.4.8.3 [START_REF] Heinz | Simple Combinations of Lineage-Determining Transcription Factors Prime cis-Regulatory Elements Required for Macrophage and B Cell Identities[END_REF]. Whole genomic ranges were kept as such for de novo motif search in publically available TF ChIP peaks.

ii. TF and histone enrichment

The ChEA database, provided online by the EnrichR collection of tools, was used for the first evaluation of the enrichment of TF binding on DEGs and genes mapping close to DMPs and DB H3K27ac regions [START_REF] Kuleshov | Enrichr: a comprehensive gene set enrichment analysis web server 2016 update[END_REF]. Only TF enrichments with an adjusted p-value <0.05 were considered. For overlaps with specific TFs, histone marks and ChromHMM annotations, we downloaded published histone and TF ChIP-seq datasets (Table 4) and we overlapped them with the genomic ranges of DMPs, DB H3K27ac peaks and DEGs produced in this study. All public datasets used were based on MCF-7 cells cultured in E2 containing media. Whenever possible, optimally IDR-thresholded peaks mapped on hg19 were downloaded.

Histone ChIP-seq datasets contained a large amount of peaks, many of which carried a low signal. To be more stringent, we ranked peaks by signal and we filtered out the second half of the peaks. E2-specific peaks were derived from ER peaks that were present in E2-treated MCF-7 and absent from E2 deprived cells [START_REF] Swinstead | Steroid Receptors Reprogram FoxA1 Occupancy through Dynamic Chromatin Transitions[END_REF]. The enrichment was calculated as the ratio of genomic-range overlap of the differential hits over the overlap with the total probes, peaks or expressed genes respectively for each dataset (differential/expected). More precisely, for DNAm we used the genomic ranges of all probes on the 850k array, for H3K27ac marks we used the IDR-thresholded peaks of CTR at d14 and for transcription we used the genomic range of expressed genes. Fisher's exact test was used to evaluate enrichments with a nominal significance level of p<10 -5 . This approach should be taken with some caution in the context of DEGs because the size of genes is variable and very large and therefore this particular analysis resulted in large fractions of overlap with mostly null fold-enrichments.

iii

v. Visualisation

For the visualization of H3K27ac peaks on IGV (2.4.19) [START_REF] Thorvaldsdóttir | Integrative Genomics Viewer (IGV): high-performance genomics data visualization and exploration[END_REF], bam files were normalized using the reads per genome coverage (RPGC) normalization method (Ramírez et al., 2016). BigWig or bedgraph format was used for publicly available datasets. Graphpad 6.0 was used for bar and scatter plots, R 3.5.1 was used for volcano and box plots, NMF package (0.21.0) was used for heatmaps and Inkscape 0.92 was used to assemble layouts.

vi. Data availability

All raw and processed data generated during this thesis project are made available under the accession number GSE132514 on the Gene Expression Omnibus (GEO)

platform and the access token can be provided the reviewers upon request. The same applies to tables of DMPs, DB H3K27ac peaks and DEGs as well as R scripts.

R RESULTS

A. E2 deprivation leads to a global epigenetic-mediated downregulation of enhancers

E2 deprivation is a well-established approach to study genome-wide response to ligand-dependent activity of ER. It is also a widely exploited in vitro model that mimics the molecular mechanisms behind aromatase inhibitor treatment used in clinical settings. Here I assess the impact of strict E2 deprivation on genome-wide DNAm, histone acetylation and gene expression using an optimised culture protocol of ER positive cell line MCF-7 (considered as the "workhorse" for studying ERpositive breast cancer) for a total duration of 14 days.

E2 deprivation leads to hypermethylation of ER-dependent enhancers i. Genome-wide methylation

The following section was written in an effort to give a comprehensive view of the bioinformatic processing steps and differential analysis behind this work. Please refer to the first section of Material and Methods for a complete view on the experimental design (p.26).

Preprocessing and normalisation

To analyse the effect of E2 deprivation and re-stimulation, samples in biological triplicates (see experimental design, Figure 7) were processed on an 850k array. The experimental design was summarised in the form of a table necessary for the further statistical analysis to which I will refer to as pdata (Table 5). The variables of interest are 'treatment' (CTR, E2D and ReSt) and 'time' (d0, d4 and d14) as well as the concatenation of these two, designated as 'group' After linking the pdata table and the raw data (.idat), samples were assessed for their quality. Median signal intensities for both methylated and unmethylated probes clustered together and all samples were well beyond the accepted intensity threshold (Figure 12). Based on these, all samples met the quality requirements and therefore they were all kept in downstream analysis. I then applied to all samples a functional normalisation (funnorm) [START_REF] Fortin | Functional normalization of 450k methylation array data improves replication in large cancer studies[END_REF], a method that normalises the signal based on some control probes that are not meant to measure biological signal but rather batch effects related to inter-array variability and bisulphite conversion. As a result, the global DNAm distribution is more homogeneous across samples (Figure 14A). For downstream analysis, it was also important to log-transform betas (M-values), as this process centres data around 0 while rendering the distribution more normal-like and hence more usable by standard statistical methods (Figure 14B). I removed a list of probes that was predicted in silico to hybridise to more than one genomic location (n=43254) [START_REF] Mccartney | Identification of polymorphic and off-target probe binding sites on the Illumina Infinium MethylationEPIC BeadChip[END_REF] as well as probes that had no signal in at least one sample (n=207), ending with 822693 working probes. In order to check for potential batch effects affecting the variability between samples, a principal component analysis was performed (Figure 15). None of the known variables were contained in the first principal component but treatment and time variables were the only two variables significantly explaining the second principal component (Table 6). Although multiple batch correction functions exist, we decided to proceed without batch correction to avoid inflating differences due to our limited sample size [START_REF] Rødland | Methods that remove batch effects while retaining group differences may lead to exaggerated confidence in downstream analyses[END_REF].

E2-deprivation leads to time-increasing hypermethylation

To determine whether E2-deprivation has an impact on DNAm, we compared E2D to CTR groups, not considering ReSt group. We used two approaches to test this:

Model (a)
a model where the effect of treatment was contrasted (CTR vs E2D) and where the effect of time was taken into account (d0, d4, d14) The comparison of E2D (n=6) against CTR (n=9) revealed a total of 2978 DMPs, the majority of which were hypermethylated in E2D (Model (a)). The lambda value, which is a statistical measure of inflation 4 , as well as the shape of the qq-plot indicate that the p-values obtained for the E2D DMPs are relatively inflated (Figure 16). A second observation was that time in culture also had an effect on DNAm and yielded 400 DMPs (d0 vs d14). The differences between the first and the last time point ranged from -32% to 23% Δβ. However, out of these 400 DMPs, only 9 were overlapping with the DMPs associated to the E2D treatment, indicating that the effect of treatment and of time were barely confounded (Figure 17). pdata$group <-relevel(pdata$group, ref = "CTR_d4") design <-model.matrix(~as.factor(group), data=pdata) fit <-lmFit(mval, design) and pdata$group <-relevel(pdata$group, ref = "CTR_d14") design <-model.matrix(~as.factor(group), data=pdata) fit <-lmFit(mval, design)

#
For the second analysis, we performed the differential analysis between CTR and E2D within each timepoint based on the assumption that if there is an effect of time in culture on DNAm, then this effect should be equivalent within the same timepoint.

Therefore, two regression models were ran for each time point d4 and d14 comparing every time CTR and E2D in triplicates, omitting the CTR at d0. The comparison of CTR and E2D at d4 revealed only one DMP while at d14, the same comparison revealed 3850 DMPs. This suggests that the effect of deprivation is mostly visible on d14. This is also reflected in the qq plots corresponding to the comparisons that were stratified by day (Figure 18).

E2-deprivation leads to hypermethylation of enhancers

In the previous section we observed that E2-deprivation leads to a majority of hypermethylation events that became more prominent with time (Figure 19, Figure 20). We next sought to examine if these changes are related to any particular functional elements. For this, we analysed the genomic features that overlapped with the DMPs and we observed that, compared to the distribution of all probes on the 850k array, the hyperDMPs were mostly found in genomic regions distant from both promoters and CpG islands (these elements often coincide) and were greatly enriched in intronic and especially in open sea contexts (Figure 21A and Figure 21B.). In contrast, hypoDMPs, represented by low number of DMPs, were not significantly enriched in any particular feature and were not further considered in this study.

In order to further investigate potential function of these DNAm changes, we overlapped the positions of DMPs with publicly available histone ChIP-seq data.

Among best-characterised histone marks are H3K4me1 and H3K4me3, which are location marks for enhancers and promoters respectively, and H3K27ac and H3K27me3, that are more dynamic in response to cellular and environmental stimuli distinguishing active from repressed regions, respectively. This analysis showed that hyperDMPs were significantly enriched in both H3K4me1 and H3K27ac marks (79.0% and 27.9%, respectively), which suggests that hypermethylation occurs in putative active regulatory regions (Figure 21C). In contrast, we observed that hyperDMPs were depleted (14.2%) in H3K4me3 marks, while no hyperDMPs were found in H3K27me3 (Figure 21C). To further characterise the functional context of the hyperDMPs, we overlapped the DMPs with ChromHMM5 annotations of predicted chromatin states in MCF-7 [START_REF] Taberlay | Reconfiguration of nucleosome-depleted regions at distal regulatory elements accompanies DNA methylation of enhancers and insulators in cancer[END_REF]. In this way, we confirmed that the non-promoter regions in which hyperDMPs are enriched could be likely enhancers, as hyperDMPS were significantly represented in enhancer and transcribed regions (49.8% and 14.9%, respectively), whereas they were lowly represented in promoter, CTCF-bound or repressed regions (1.9%, 0.8% and 0.0%, respectively) (Figure 21D). These results show that the increase in DNAm in response to estrogen deprivation mostly occurs within regulatory regions in introns that contain enhancers.

Hypermethylation associated to estrogen deprivation occurs predominantly in lowly and partially methylated locations

We observed that initial DNAm levels at sites that become hypermethylated had already low to intermediate methylation levels prior to E2 deprivation. More precisely, DNAm levels measured at d0 of the experiment ranged 40% and 71% (1 st and 3 d quartile) with an average of 55% (Figure 22). This is in line with past studies

showing that intermediate levels of DNAm are also associated with regulatory elements [START_REF] Nothjunge | DNA methylation signatures follow preformed chromatin compartments in cardiac myocytes[END_REF][START_REF] Stadler | DNA-binding factors shape the mouse methylome at distal regulatory regions[END_REF].

ii. Validation

To further validate the above changes in DNAm, I replicated experiments by repeating E2 deprivation in two independent assays but also by using two alternative methods to inhibit ER signalling: by (i) silencing of ESR1 transcription and (ii)

antagonist-mediated ER inhibition (being used in latter section of the thesis).

Biological validation of E2 deprivation

The E2 deprivation assay was repeated in two independent experiments where MCF-7

were deprived for 7 and 9 days (experiment 1 and 2, respectively) after which the DNAm of several top DMPs was measured by pyrosequencing. There was a visible increase of DNAm in time among all tested CpG sites (Figure 23). Intermediate DNAm levels at d0 were also observed in these assays.

Figure 23. DNAm levels measured after 7 days (exp1) and 9 days (exp2) of E2 deprivation. DNAm was measured by pyrosequencing and is shown as average ± SD of triplicates. The effect of treatment was evaluated using a Mann-Whitney test.

Silencing ESR1 gene expression

Finally, siRNA targeting ER's transcript, namely ESR1, was used as an alternative to reduce the activity of ER. The levels of ESR1 transcripts were efficiently knockeddown by 92.6% 5 days after the transfection (Figure 24). The efficiency of knockdown was also likely translated at the level of ER protein, since the expression of GREB1, TFF1 and PGR -that are canonical ER targets -was also strongly downregulated in the same samples, even in presence of E2 in the medium (Figure 24). Once the efficiency of the ESR1 knock-down was confirmed, I analysed DNAm of several CpG sites that were found hypermethylated in response to E2 deprivation.

Five days after knock-down, we observed an increase of average DNAm ranging from 6.7% to 13.8% for three CpG sites matching probes from the 850k array (Figure 25).

The sequencing window of the probe cg13374172 allowed us to capture the DNAm of a CpG site 10bp downstream of the 850k CpG site and interestingly its DNAm was not affected as much as its neighbouring CpG (Figure 25, last panel). Although changes appeared statistically significant and are in line with the DNAm changes observed in the E2 deprivation assay, this result should be interpreted with caution due to the limited sample size.

Figure 25. DNAm levels measured after 5 days of siRNA targeting ESR1 transcripts. DNAm was measured by pyrosequencing and is shown as average ± SD of triplicates. The effect of treatment was evaluated using a Mann-Whitney test.

Histone deacetylation in enhancers in response to E2 deprivation

The analysis of genome-wide DNAm revealed that E2 deprivation lead to the hypermethylation of CpGs located in intragenic enhancers. The latter were identified using two approaches: their occurrence within regions having a high ratio between

H3K4me1 and H3K4me3 and their overlap with enhancers that were previously annotated with ChromHMM in MCF-7 [START_REF] Jozwik | FOXA1 Directs H3K4 Monomethylation at Enhancers via Recruitment of the Methyltransferase MLL3[END_REF][START_REF] Taberlay | Reconfiguration of nucleosome-depleted regions at distal regulatory elements accompanies DNA methylation of enhancers and insulators in cancer[END_REF]).

Yet, the physical overlap with enhancers does not provide us with sufficient information on the degree of activity of these regulatory elements. Multiple methods have been established with time to measure enhancer activity, such as sequencing of enhancer RNAs, RNA Pol II binding or H3K27ac marks [START_REF] Franco | Enhancer transcription reveals subtype-specific gene expression programs controlling breast cancer pathogenesis[END_REF][START_REF] De Santa | A Large Fraction of Extragenic RNA Pol II Transcription Sites Overlap Enhancers[END_REF]. Based on the above, we sought to find whether the E2D DMPs were enriched within H3K27ac marks but also if H3K27ac levels varied in response to E2 deprivation.

i. ChIP of histone marks in proximity to DMPs

For this, we conducted a pilot experiment where MCF-7 were deprived for 10 days followed by measuring three histone marks in proximity (<100bp) of top DMPs detected earlier. We chose H3K4me1 to verify if DMPs were in an enhancer context as well as H3K27ac and H3K27me3, two histone marks indicating enhancer activation and repression, respectively (Figure 26). H3K4me1 and H3K27ac marks were detected in the proximity of all four DMPs tested while H3K27me3 levels were low. In addition, H3K27ac signal was on average higher in CTR group in contrast to E2D. As a positive control of active enhancer, we used a GREB1 enhancer known to be active under regular E2-containing culture conditions for MCF-7 [START_REF] Jozwik | FOXA1 Directs H3K4 Monomethylation at Enhancers via Recruitment of the Methyltransferase MLL3[END_REF]. Consistent with the picked regions, H3K27ac signal was enriched and was on average higher in CTR in comparison to E2D in this control region, although not significantly, while H3K27me3 was barely detected. In contrast, a higher H3K27me3 signal was detected in the control that mapped within the promoter of a gene that is neither expressed nor contains any DMPs. Together, the absence of H3K27me3 and the presence of H3K27ac signal suggest that the panel of picked DMPs is likely located in proximity to enhancers that are seemingly more active in CTR than in E2D conditions.

Figure 26. ChIP of H3K4me1, H3K27ac, H3K27ac marks and IgG in CTR at d11 and E2D at d10 across four top DMPs (title: 850k probe name). qPCR data was normalized over input DNA and presented as an average ± SD of duplicates. Data was analysed using a Wilcoxon ranked test within each histone mark.

ii. ChIP-seq of active histone mark H3K27ac

In order to have a genome-wide view of specific histone mark specific for active enhancers, MCF-7 were cultured in presence and absence of E2 for 4 and 14 days, and H3K27ac profiles were measured with ChIP-seq. Narrow peaks were called against a CTR input and the IDR method was used to keep highly consistent peaks between the biological duplicates. This yielded 31,090 consistent peaks for CTR, 32,995 peaks for E2D and 36,214 peaks for ReSt. The higher number of IDR-filtered peaks in the ReSt and E2D group simply reflects the higher read homogeneity between duplicates of these groups (PC1, Figure 27A). The comparison of CTR vs E2D led to the identification of 3053 differentially bound (DB) regions (|log 2 FC| >1, FDR <0.05) (Figure 27B) the large majority of which showed a lower signal in E2D cells compared to CTR cells (Figure 28). Total H3K27ac peaks were evenly detected across promoter, intronic and intergenic regions, but the DB regions were located in a similar genomic context to DMPs, that is, in intronic and intergenic regions (42.9% and 41.1%, respectively, Figure 29A). An overlap with the ChromHMM annotation revealed that 44.0% of these peaks were found in enhancer regions (Figure 29B).

Furthermore, under E2-stimulated conditions, 25.9% of DB peaks co-localised with ER binding events and almost half of the them (11.2% of total DB peaks) overlapped with E2-dependent ER peaks measured by [START_REF] Swinstead | Steroid Receptors Reprogram FoxA1 Occupancy through Dynamic Chromatin Transitions[END_REF] (examples shown in Figure 30). Taken together these results point to an inactivation of enhancers that could be in part mediated by direct ER binding. A B

Gene expression mostly decreases with E2 deprivation

To further characterise the genome-wide impact of ER activity on gene regulation, we performed an analysis of the transcriptome (RNA-Seq) of CTR and E2D conditions at d4 and d14 in the same samples used for the methylome analysis. Resulting expression data served two main purposes at this stage of the project:

• Measuring the global effect of E2 deprivation in stricto sensu and assessing the concordance of our protocol with other studies that were focused on ERdependent transcription

• Finding specific targets that may explain the observed epigenetic changes Overall 14,617 genes were expressed at detectable levels across all samples and the comparison of treatments (CTR vs E2D), after adjusting for time, revealed 5466 differentially expressed genes (DEGs) with a FDR = 0.05. The addition of an absolute cut-off of log 2 FC >1 limited the number of DEGs to 547, 71.1% of which were downregulated in response to E2 deprivation (Figure 31A). In a meta-analysis of ERregulated genes in MCF-7 cell line, Stanislawska-Sachadyn et al. identified a set of 24 genes that are commonly de-regulated in response to different exposure times and E2 concentrations [START_REF] Stanislawska-Sachadyn | Transcriptomic Effects of Estrogen Starvation and Induction in the MCF7 Cells[END_REF]. Of these, 19 DEGs were also significantly deregulated in this study, and 17 changed in the expected direction. In addition, among the top DEGs, I identified genes known to be directly upregulated by ligand-activated ER through their promoter-enhancer interactions, such as GREB1, TFF1 or PGR, supporting that our cell culture model is in line with the previously reported impact of inactivation of ER (Figure 31B) [START_REF] Fullwood | An oestrogen-receptor-α-bound human chromatin interactome[END_REF]. The expression of these particular genes was also validated by RT-qPCR (Figure 32).

A B

=2.1) and decreased the expression of PADI4 and TET2 (log 2 FC =-1.5 and -0.8 respectively). HDAC9 is a histone deacetylase and is part of the family of HDACs that are commonly associated to gene repression, PADI4 is a transcriptional coactivator, and TET2 is an enzyme catalysing CpG demethylation and is also associated to gene activation (Figure 33A) [START_REF] Plass | Mutations in regulators of the epigenome and their connections to global chromatin patterns in cancer[END_REF][START_REF] Zhang | Genome-Wide Analysis Reveals PADI4 Cooperates with Elk-1 to Activate c-Fos Expression in Breast Cancer Cells[END_REF]. The same genes were also significantly differentially expressed in the same direction in MCF-7 that were treated with ICI 182,270 for 4 and 14 days, a drug that is known for inducing ER protein degradation and consequently down-regulation of ER signalling (Figure 33B). This further consolidates the possibility that the deregulation of these ERFs is mediated by ER binding. Taken altogether these results show that E2 deprivation leads to a global gene down-regulation and to the deregulation of specific ERFs that are in line with the loss of enhancer activity. It is noteworthy that TET2 absolute expression levels were much higher in contrast to HDAC9 or PADI4 that had a low baseline expression, suggesting that a change in expression of TET2 could have a more consequent impact at the protein level compared to the other two transcripts.

Another interesting observation across all tested ER targets was that the expression changes were not noticeable between d4 and d14 unlike DNAm changes that were clearly time-dependent.

B. Integration of ER-dependent epigenetic and transcription changes

At this stage of the study we have in our hands three datasets (methylome, H3K27ac

ChIP-seq and transcriptome) that were all produced in ER-expressing breast cancer

A B

cells that were deprived of E2 for 14 days. The changes in DNAm profile (850k data)

and the binding profile of the active chromatin mark H3K27ac (ChIP-seq data)

suggested that the activity of ER-dependent enhancers decreases in response to E2 deprivation. In parallel, the transcription-wide data obtained by RNA-sequencing showed that the expression levels of 71% of DEGs decreased. In this section, we aimed to integrate these datasets in an attempt to find common molecular denominators associated to ER that may drive these epigenetic and transcriptional changes.

Associations of epigenetic and transcriptional changes i. Co-occurrence of epigenetic changes

In the onset of this study, we observed that a part of DMPs obtained in response to E2 deprivation were significantly overlapping with publically available H3K27ac peaks in MCF-7 (27.9% of hyperDMPs; 1.5x compared to 850k overlap, see p.51), which suggests that differential DNAm occurs in active chromatin regions. I therefore wanted to see if this overlap was also significant in our experimental setting. For this, I first overlapped DMPs with all filtered and IDR-thresholded H3K27ac peaks obtained under CTR conditions. 10.9% of E2-deprivation DMPs overlapped regions carrying H3K27ac marks in CTR conditions (Table 10). This was lower to the expected overlap that I observed earlier where 27.9% DMPs overlapped with ENCODE H3K27ac ChIP-seq data (see p.51), that was obtained from MCF-7 in similar culture condition to CTR in this study. One plausible explanation for this observation is that this specific ENCODE peak file contained peaks reaching up to 42kb, while I applied a maximal peak size of 1.5kb. According to this, there is a possibility that the first overlap of DMPs with H3K27ac peaks was over-estimated and thus advocate that more filters to publically available datasets need to be applied.

I then focused specifically on the 3053 DB regions and found only 3 DMPs overlapping within 500bp away from the centre of these H3K27ac peaks (Table 10).

A technical consideration that potentially contributes to the limited overlaps could be that the ChIP-seq and 850k arrays were performed not only in two independent experiments but also in two different batches of MCF-7. Nevertheless, after expanding the genomic window of both CTR and DB H3K27ac peaks, the percentage of overlapping DMPs considerably increases in an absolute way but also in comparison to the percentage of overlap of total available CpG sites (Table 10). This led us to the hypothesis that DNAm and H3K27ac levels at defined locations may be regulated by common chromatin remodelling complexes, and that these complexes do not necessarily operate in short genomic windows but rather in distant cis locations through chromatin looping events.

ii. Correlation of DNAm and gene expression

Although the overlap between DEGs and genes and carrying DMPs was not significant, I performed a more in-depth correlation analysis matching all DMPs and DEGs, this time with no differential cut-offs. I first correlated DNAm and expression within CTR and E2D samples (n=12) and found 666 DMPs that both negatively or positively correlated with expression (n=342 and 324, respectively) (Spearman |ρ| > 0.67, FDR <0.05, Figure 34). After filtering for highly variable DMPs and DEGs (|log 2 FC|>1 and |Δβ|>10%), only 51 DMPs strongly correlated -mostly negativelyto differential expression. These 51 DMPs match altogether to 40 genes that were identified in the pairwise overlaps above. Among these genes we found GREB1 and PKIB, both of which are common ER targets (Figure 35). Nevertheless, in order to find any causal relation between DNAm and expression, one would need to specifically target the DNAm at one of these sites for instance with the help of dCas9-Dnmt3a systems [START_REF] Liu | Editing DNA Methylation in the Mammalian Genome[END_REF] and check if the expression changes. iii.

Overlaps of nearest gene

Because the correlation of epigenetic and transcription changes within the same genomic location was limited, I tested if E2 deprivation had an impact in genes that were deregulated in at least two out of the three measured readouts, independently of the direction of changes (DNAm, H3K27ac and gene expression). For this, genomic regions carrying DNAm and H3K27ac changes were firstly annotated with the nearest gene, which sometimes was more than 1Mb away from DMPs and DB H3K27ac.

However, picking the nearest proximal gene in terms of genomic distance is a safe approach and commonly used to match regulatory elements such as enhancers and gene expression [START_REF] Corces | The chromatin accessibility landscape of primary human cancers[END_REF]. Next, I performed pairwise overlaps between genes corresponding to each differential dataset and tested if the common genes were significantly enriched compared to the background overlap of all available genes.

Genes containing or mapping closely to a DMP were significantly likely to also map near to a DB H3K27ac peak (FC=2.5x, Fisher p= 8.7x10 -30 , Figure 36). This was also the case for DEGs that were also significantly enriched in DB H3K27ac by 2.5x

(Fisher p= 7.6x10 -21 , Figure 36). Nevertheless the number of genes being altered at all three levels was extremely low and this low overlap was mainly driven by the limited overlap of DEGs and genes containing a DMP.

At this point, we found that the global down-regulation in gene expression is in line with the observed general increase of DNAm and loss of histone acetylation marks.

Nevertheless, the overlap between DEGs and genes containing DMPs or DB H3K27ac peaks in cis was relatively low compared to the total number of shared genes. Considering that most of these epigenetic changes in enhancers have a functional impact, we reasoned that they may not necessarily regulate the expression of the nearest gene but rather resulting from dynamic binding of TFs associated to ER. Together these results suggest that under E2-stimulated conditions, ER shapes epigenetic marks both directly and indirectly while ER deactivation leads to a global decrease in active enhancer marks and to a major down-regulation in gene expression.

RESULTS

69

Figure 36. Pairwise overlaps of DEGs and genes that map to nearby DMPs or DB H3K27ac peaks. For fold-enrichment significance, a Fisher test against only genes present in both differential datasets were considered as background.

Identification of ER key cofactors in E2-dependent gene regulation

The results obtained so far show that E2 deprivation leads to DNA hypermethylation and a decrease in histone acetylation of enhancers and that these changes are partially associated with gene expression down-regulation in cis. As it is widely accepted that E2 deprivation leads to reduced ER binding and time of residence, it is conceivable that E2 deprivation could also affect the binding of ER cofactors to DNA [START_REF] He | c-Jun/AP-1 overexpression reprograms ERα signaling related to tamoxifen response in ERα-positive breast cancer[END_REF][START_REF] Swinstead | Steroid Receptors Reprogram FoxA1 Occupancy through Dynamic Chromatin Transitions[END_REF][START_REF] Zwart | Oestrogen receptor-co-factor-chromatin specificity in the transcriptional regulation of breast cancer[END_REF]. This is indeed the case for p300, SRC1/2/3 and JUN whose binding can decrease after E2 deprivation but it is not for FOXA1 whose binding is barely affected in the same context [START_REF] Glont | Comprehensive Genomic Analysis Reveals that the Pioneering Function of FOXA1 Is Independent of Hormonal Signaling[END_REF][START_REF] He | c-Jun/AP-1 overexpression reprograms ERα signaling related to tamoxifen response in ERα-positive breast cancer[END_REF][START_REF] Swinstead | Steroid Receptors Reprogram FoxA1 Occupancy through Dynamic Chromatin Transitions[END_REF][START_REF] Theodorou | GATA3 acts upstream of FOXA1 in mediating ESR1 binding by shaping enhancer accessibility[END_REF]. Based on the above, we hypothesised that the observed epigenetic changes are the result of ER's altered cistrome and we sought to explore the TFs that may coregulate these events.

For this, I first conducted an enrichment analysis of transcription factor binding sites in the proximity of DMPs, DB H3K27ac and DEGs using the ChEA database [START_REF] Kuleshov | Enrichr: a comprehensive gene set enrichment analysis web server 2016 update[END_REF]. For each type of dataset, a list of TF that significantly bind to the associated genes was obtained. The overlap of the enriched TF lists revealed that in addition to expected enrichment of ERα and ERβ binding events, binding of other TFs (including ZNF217, TFAP2C and GATA3) were found enriched among the genes that map near DMPs and DB H3K27ac regions (Figure 37). These findings support functional and physical protein interactions between these TFs and ER that have been reported in previous studies [START_REF] Mohammed | Endogenous Purification Reveals GREB1 as a Key Estrogen Receptor Regulatory Factor[END_REF][START_REF] Nguyen | A functional interplay between ZNF217 and Estrogen Receptor alpha exists in luminal breast cancers[END_REF].

Interestingly, when considered separately, the up-regulated fraction of DEGs was enriched for binding sites of PRC2 members SUZ12 and EZH2, consistent with the fact that when ER is active, it can also actively repress gene expression (Figure 38) [START_REF] Ariazi | A New Role for ERα: Silencing via DNA Methylation of Basal, Stem Cell, and EMT Genes[END_REF]. TF binding sites within the down-regulated fraction of DEGsthat are active when ER is induced by E2 -were also enriched in genes containing DMPs and DB H3K27ac. This suggests that common TFs are involved in the ERdependent maintenance of these genes' epigenetic marks and their expression. Although informative, the ChEA database is limited by the number of datasets it contains and provides a broad picture of TF enrichment based on gene names rather than specific genomic sequences. Therefore, to detect putative DNA-binding TFs in a more unsupervised approach, we searched for enriched TF motifs in proximity of E2 deprivation-mediated DMPs and DB H3K27ac by performing a de novo motif analysis. Our first observation was the lack of the majority of the motifs corresponding to the above-identified TFs, including the canonical ER-binding motif, suggesting that ER and these particular TFs do not bind directly on or in close proximity of DMPs and DB H3K27ac regions. Instead motifs matching AP-1 and FOX TFs, known to be main ER cofactors, were enriched in proximity to both DMPs and DB H3K27ac regions (Figure 39). Based on these findings and the knowledge of reported ER cofactors, I overlapped a set of publically available TF ChIP-seq datasets in MCF-7 cultured in conditions equivalent to CTR (Table 4) with DMPs, DB H3K27ac marks and DEGs identified in this study. The significantly enriched TFs were then split in three categories: 1) TFs that are both highly enriched and overlapping with top hits, 2) TFs that were highly enriched but overlapped with low to middle fraction of top hits and 3) TFs that were lowly enriched but overlapped a large fraction of top hits.

FOS, a major AP-1 component, was the only TF that matched the criteria of the first category as it was enriched by 9.1x and overlapped with nearly 1/3 of DMPs. FOS enrichment in DB HK27ac regions and DEGs was lesser (<2.0x) but overlapped a large fractions of regions. ER, FOSL2, p300, AHR and NCOA3 met the criteria of the second category as they were also importantly enriched across all differential datasets but with moderate and low percentages of overlap (Figure 40). Inversely, FOXA1, GATA3 and ZNF217 matched the third category as they were moderately enriched in DMPs and DB HK27ac regions although they highly overlapped all differential datasets. This is in line with the preferential binding of these factors in proximity to active marks or active genes. At last, a de novo motif search in the regions bound by ER and its cofactors (namely FOXA1, FOS, TFAP2C, p300, ZNF217, GATA3, and AHR) revealed an important occurrence of AP-1 and FOX motifs in all cofactors' binding regions (Figure 41).

Together, these findings emphasise the importance of ER tethering and suggest that AP-1 plays a central role in the maintenance of DNAm and H3K27ac marks, thereby orchestrating chromatin changes through enhancer regulation. It is also noteworthy that the expression of all the above-mentioned TFs remained unchanged after E2 deprivation, which hints towards a regulation in cis rather than trans. 4.

C. Partial reversibility of epigenetic and expression changes after ER re-activation

We mentioned earlier that E2 deprivation and ER antagonists are used as an in vitro proxy to application of aromatase and ER inhibitors in the clinical setting. Such treatments are applied for a defined period of time after which ER signalling can resume. Nevertheless, some patients were found to fail responding to the treatment or even relapse that is frequently associated with development of endocrine resistance to the initial treatment. This led to the hypothesisis that some irreversible molecular changes may occur during this process. Having this in mind, I asked whether the observed epigenetic and transcriptional changes could be reversed or precluded from further increasing if the E2 deprivation interrupted or ER inhibitors removed.

Partial reversibility of epigenetic and transcription changes in response to E2 re-stimulation i. DNA methylation

To test potential reversibility of E2 deprivation specific epigenetic changes, MCF-7 that were E2-deprived for 4 days were re-exposed to E2 for 10 days (ReSt) before being harvested and processed on the 850k platform (Figure 43A). Data was preprocessed and normalised concomitantly to CTR and E2D (as described on p.43). I analysed the data using three approaches of linear regressions models: 11). For the second model, we assumed there is no effect of time in culture and we compared ReSt against CTR and E2D only within d14 and for all probes on the array (n=3 for each group). Again no DMP was found between the comparison between ReSt and CTR but 679 DMPs were found between ReSt and E2D (Table 12). This suggests that ReSt methylome is closest to the methylome of CTR samples rather than E2D samples. Finally, when I compared ReSt to the other two groups using only the 952 hyperDMPs obtained from the E2 deprivation, a minor fraction of DMPs was hypomethylated in CTR compared to ReSt. This suggests that there was a gain in DNAm at d4 that persisted in several CpG sites after 10 days of E2 re-stimulation. On the other hand, the large majority of DMPs was statistically different between ReSt and E2D, indicating that the gain of DNAm was largely precluded after E2 restimulation. Globally all three models detected no or few DMPs between the CTR and the ReSt.

On the other hand, the comparison of ReSt to E2D did reveal DMPs in all models, but never as many as were found in the CTR vs E2D comparison, suggesting the ReSt group was not as different from the E2D group as was the CTR group. In addition, overall DNAm in the hyperDMPs was seemingly higher in the ReSt group as visible on Figure 43B

ii. Histone acetylation

Next I aimed to compare the ReSt H3K27ac peaks to those of the CTR by focusing on the DB peaks that were common to the DB peaks of the CTR vs E2D comparison. We found that 347 out of 3053 E2-responsive DB H3K27ac peaks were significantly lower in ReSt compared to the CTR group and their signal was equivalent to that of the E2D group (Figure 45). Despite the 10-day E2 re-stimulation, the histone acetylation levels in certain regions did not return to initial acetylation levels. We then overlapped these regions with publically available H3K37me3 peaks and found virtually no overlaps. The non-recovery of H3K27ac signal in certain regions and the complete absence of H3K27me2 marks suggest a possible shift from active to primed chromatin that cannot be inversed by solely by the addition of E2. iii.

RNA expression

At the level of gene expression, I found that after E2 re-stimulation (ReSt) the majority of DEGs that were down-regulated in E2D globally returned to their expression levels with the exception of 6 DEGs whose expression was significantly different from CTR (CTR vs ReSt, FDR < 0.05, |log 2 FC| > 1). Two of the DEGs that did not recover their initial expression level were TRAF2 and NCK Interacting Kinase (TNIK) and Matrix Metallopeptidase 13 (MMP13) that are an inducer and a target, respectively, of AP-1 (Figure 46). The expression of canonical ER targets GREB1, TFF1 and PGR fully recovered (Figure 47).

Validation of partial reversibility of DNAm and transcription following removal of ER antagonist

To further validate the reversibility of DNAm changes using an independent approach, I inhibited ER activity with antagonist ICI 182,780 for 4 days and 14 days (ICI), while the re-activated group was treated with ICI for 4 days, after which the drug was removed and DNAm analysed 10 days thereafter (ReAc) (see Figure 9). For this, 4 CpG sites were selected that were among the top DMPs from the E2deprivation analysis that met the criteria for enhancer-like genomic features, and DNAm was measured in a targeted way using pyrosequencing. Similar to E2deprivation, I observed that ICI 182,780 treatment increased DNAm, and in addition DNAm changes were of larger amplitude in comparison to the response to E2deprivation (Figure 48). Interestingly, DNAm levels between ReAc and CTRs were not significantly different, although the restoration of initial DNAm was not complete for all tested sites, indicating that reactivation of ER following the antagonist removal leads to partial and slow reversal of DNAm changes. In addition, expression of TNIK, but not that of MMP13, decreased with ICI treatment and did not fully recover with the removal of the inhibitor (Figure 49). This result further supports that AP-1 can be indirectly deregulated in an ER-dependent manner. Just like for the E2 deprivation, the expression of HDAC9, PADI4 and TET2 studied earlier (Figure 33) fully recovered after ICI removal (Figure 50). Globally, these results show that epigenetic changes are largely reversible or precluded from increasing following E2 restimulation. Nevertheless, a smaller fraction of the observed changes does not 6 CTR and E2D samples are the same as those used in Figure 32 recover, even if the expression of ERFs does, suggesting the existence of an estrogendependent epigenetic memory.

D DISCUSSION

In this study, we bring novelty to the field by adapting a cell culture protocol for studying E2 deprivation and re-stimulation in stricto sensu in and by using the latest methylation array that allows a genome-wide interrogation of methylation states, including a comprehensive panel of enhancers [START_REF] Moran | Validation of a DNA methylation microarray for 850,000 CpG sites of the human genome enriched in enhancer sequences[END_REF]. We saw that prolonged E2 deprivation and re-stimulation result in time-dependent DNAm changes and in histone acetylation modifications across diverse genomic regions, many of which occur within enhancer elements. This is the first study that comprehensively characterized DNAm and histone acetylation at enhancers in response to a sequential inactivation and re-activation of ER signalling. These observations are consistent with and extend previous studies showing that ligand-mediated ER activity has an impact on the epigenome, through mechanisms that likely involve tethering of ER by other TFs [START_REF] Fleischer | DNA methylation at enhancers identifies distinct breast cancer lineages[END_REF][START_REF] Kininis | Genomic Analyses of Transcription Factor Binding, Histone Acetylation, and Gene Expression Reveal Mechanistically Distinct Classes of Estrogen-Regulated Promoters[END_REF].

The cell line model

In order to study the sequential effect of estrogen deprivation and re-stimulation, we aimed to implement a cell culture model that can address strictly the effect of E2. So far the existing protocols of E2 deprivation consist in depriving cells from all components that are removed during the process of FBS charcoal stripping. In addition, our protocol offers an appropriate compromise in terms of duration of treatment that was sufficiently long to cause noticeable molecular changes and also short enough to prevent cells from losing their sensitivity to estrogen. This latter issue is linked to the nature of the cell line itself. MCF-7 is a cancer cell line that was established 1973 and has become ever since the most widely used cell model in the context of ER-positive breast cancer research [START_REF] Soule | A Human Cell Line From a Pleural Effusion Derived From a Breast Carcinoma2[END_REF]. Being a cancer cell line, it is not surprising that plenty of clones have been derived from the initial culture ever since its establishment, one of which, HTB-22, is commercially available by ATCC in passage no less than 147. The "age" of the cell line, as well as the large variety of MCF-7 subclones contribute at least in part to reported discrepancies between studies and lack of reproducibility. Unfortunately, this is also reflected at the level of large consortia such as ENCODE, that used 2 different sources of MCF-78 , although the HTB-22 clone was used in most cases. All these considerations underscore the necessity to use consensus models and protocols within the scientific community, especially for widely used cell lines such as MCF-7.

Because of all the limitations related to the current models, we believe that our approach to culturing MCF-7 offers a robust and reproducible alternative that addresses the effect of one single molecule and that it consolidates the knowledge of the molecular biology of ER. Indeed, the direction of the changes in expression of ER target genes in response to E2 deprivation, ER inhibition and siRNA indicates that our protocol recapitulates the expected changes that are caused from the deregulation of ER signalling. Our protocol approach could easily be applied in the context of other molecules.

Hypermethylation and histone deacetylation of ER-dependent enhancers

We showed that DNAm levels consistently increase whereas H3K27ac signal decreases in a particular set of enhancers using different molecular methods and in total three different cell culture approaches that effectively down-regulate ER activity:

• E2 deprivation which reduces dramatically ER binding events without affecting the protein expression

• siRNA targeting ER's transcript which reduces the amount of ER available at the protein level

• treatment with ICI that leads to the reduction of ER at the protein level through promoting its degradation

The common point to all of the above treatments is that ER binding decreases in a consequent number of sites. ER's reduced binding in response to these treatments is extensively shown in literature and this is also evident in this study through the deregulation of canonical and newly discovered target genes. Depending on the approach, epigenetic changes were detectable as early as 4 days and became FOX motifs were enriched by 6.5x and 3.9x in proximity of DMPs and by 3.9x and 1.7x in DB H3K27ac regions. These motif enrichments were associated to important overlaps of FOS binding -an AP-1 component -and DMPs while both FOS and FOXA1 highly overlapped with DB H3K27ac regions (Figure 40). Based on these observations, it is conceivable that AP-1 regulates DNAm of CpG sites that are proximal to its binding sites in the context of an E2-induced system. Remarkably, AP-1 and FOX motifs were consistently and highly enriched at all binding sites of all established ER cofactors investigated (namely FOXA1, FOS, TFAP2C, p300, ZNF217, GATA3, and AHR) (Figure 41). FOXA1 and GATA3, followed by FOS binding, were also significantly enriched at functional enhancers (defined by an expression-methylation quantitative trait loci analysis) that distinguish ER-positive from ER-negative BC subtypes [START_REF] Fleischer | DNA methylation at enhancers identifies distinct breast cancer lineages[END_REF]. This raises the possibility that the specific mechanisms are conserved between cell line models and clinical samples.

These results are consistent with the frequency by which ER regulates genomic regions through a tethering mechanism with other cofactors, rather than through direct binding [START_REF] Carroll | EJE PRIZE 2016: Mechanisms of oestrogen receptor (ER) gene regulation in breast cancer[END_REF][START_REF] Mohammed | Endogenous Purification Reveals GREB1 as a Key Estrogen Receptor Regulatory Factor[END_REF].

ii. Cross-talk of ER, ER cofactors and epigenetic regulatory factors

However, all these ER and cofactors dynamic binding events cannot directly modify epigenetic marks. The edition of the latter is in fact catalysed by ERFs, including different enzymes from the DNMT, TET, HAT and HDAC families. Previous studies have established that DNA hypermethylation and histone deacetylation at gene promoters are associated with gene silencing [START_REF] Weber | Distribution, silencing potential and evolutionary impact of promoter DNA methylation in the human genome[END_REF]. In contrast, there is less data on the role of epigenetic changes at enhancer regions, and the specific mechanisms that underlie specific targeting of enhancers for hypermethylation is unclear. It is known that specific chromatin interactions, involving DNAm dynamics, are abolished and are replaced by others in response to ER activation and inactivation [START_REF] Fullwood | An oestrogen-receptor-α-bound human chromatin interactome[END_REF][START_REF] Hsu | Estrogenmediated epigenetic repression of large chromosomal regions through DNA looping[END_REF][START_REF] Taberlay | Reconfiguration of nucleosome-depleted regions at distal regulatory elements accompanies DNA methylation of enhancers and insulators in cancer[END_REF]. Indeed, in some instances DNAm deposition is directed by TF-mediated DNMT recruitment, whereas in others DNAm can be the result of the lack of TF occupancy which exposes these sites to DNMT complexes that have an affinity for unmethylated CpGs [START_REF] Hervouet | Specific or not specific recruitment of DNMTs for DNA methylation, an epigenetic dilemma[END_REF][START_REF] Stadler | DNA-binding factors shape the mouse methylome at distal regulatory regions[END_REF]. The disengagement of ER and its cofactors from chromatin induced by E2 deprivation, could make unmethylated CpG sites susceptible to methylation, which is consistent with the observation that methylated DNA represents a default state [START_REF] Singh | De Novo DNA Methylation in the Male Germ Line Occurs by Default but Is Excluded at Sites of H3K4 Methylation[END_REF].

In an alternative, although not mutually exclusive scenario, E2 deprivation-associated hypermethylation may be the result of reduced demethylating activity at enhancers involving TET methylcytosine dioxygenases. This could happen through the reduction of TETs' expression and/or through a decreased recruitment of TET2 on enhancers. Interestingly, we observed both after E2 deprivation and ICI inhibition a decrease in expression of DNA demethylase TET2 (Figure 33). The possibility of a double mechanism explaining hypermethylation is supported by a recent study showing that ER not only recruits TET2 on ER-dependent enhancers but it also positively regulates TET2 expression in response to E2 induction, and that TET2 knock-out leads to hypermethylation of enhancers [START_REF] Wang | TET2 coactivates gene expression through demethylation of enhancers[END_REF]. Given that AP-1 peaks and motif were enriched in proximity to E2-dependent DMPs, it is possible that AP-1 could also interact with TETs as it was already shown for FOXA1 [START_REF] Yang | FOXA1 potentiates lineage-specific enhancer activation through modulating TET1 expression and function[END_REF][START_REF] Zhang | Nucleation of DNA repair factors by FOXA1 links DNA demethylation to transcriptional pioneering[END_REF]. Considering that ER and AP-1 often interact as part of the same complex [START_REF] He | c-Jun/AP-1 overexpression reprograms ERα signaling related to tamoxifen response in ERα-positive breast cancer[END_REF][START_REF] Jeffy | An estrogen receptoralpha/p300 complex activates the BRCA-1 promoter at an AP-1 site that binds Jun/Fos transcription factors: repressive effects of p53 on BRCA-1 transcription[END_REF][START_REF] Lupien | Growth factor stimulation induces a distinct ERα cistrome underlying breast cancer endocrine resistance[END_REF], our results support the idea that AP-1 may play a bridging role between ER and TET2 at the level of enhancers, therefore maintaining DNAm levels stable under E2stimulated conditions.

Similarly, the occurrence of AP-1 and FOX motifs and the enrichment of acetyltransferase p300, FOS and FOXA1 binding sites at DB H3K27ac regions reinforce previously reported ER-mediated maintenance of histone acetylation recruitment of p300, AP-1 and FOXA1 to enhancers [START_REF] Allison | Affinity and dose of TCR engagement yield proportional enhancer and gene activity in CD4+ T cells[END_REF][START_REF] Jeffy | An estrogen receptoralpha/p300 complex activates the BRCA-1 promoter at an AP-1 site that binds Jun/Fos transcription factors: repressive effects of p53 on BRCA-1 transcription[END_REF][START_REF] Vierbuchen | AP-1 Transcription Factors and the BAF Complex Mediate Signal-Dependent Enhancer Selection[END_REF][START_REF] Zwart | Oestrogen receptor-co-factor-chromatin specificity in the transcriptional regulation of breast cancer[END_REF]. Furthermore the decrease of H3K27ac levels after E2 deprivation is consistent with a decrease of enhancer activity. This finding, together with the absence of overlap with H3K27me3, which is characteristic of repressed and poised chromatin states, suggests that the identified enhancers may switch from an active to a primed state in response to E2 deprivation [START_REF] Caglio | Keep Them Close: PRC2 Poises Enhancer-Promoter Interactions at Anterior Neuronal Genes[END_REF][START_REF] Calo | Modification of enhancer chromatin: what, how, and why?[END_REF]. Priming of ER-dependent enhancers is further supported by the intermediate levels of DNAm which seem to play also a role in this mechanism [START_REF] Elliott | Intermediate DNA methylation is a conserved signature of genome regulation[END_REF][START_REF] Mahé | Cytosine modifications modulate the chromatin architecture of transcriptional enhancers[END_REF][START_REF] Sharifi-Zarchi | DNA methylation regulates discrimination of enhancers from promoters through a H3K4me1-H3K4me3 seesaw mechanism[END_REF].

In this regard, we propose that the reduced binding of ER could open the opportunity to HDACs to act. It is in fact plausible that the decrease of H3K27ac signal could be related to both an increased accessibility to HDACs but also the observed increase of HDAC9 expression, although there is little evidence of a direct relationship between this class IIa HDAC and H3K27ac deacetylation [START_REF] Delcuve | Roles of histone deacetylases in epigenetic regulation: emerging paradigms from studies with inhibitors[END_REF][START_REF] Petrie | The Histone Deacetylase 9 Gene Encodes Multiple Protein Isoforms[END_REF].

Based on these findings we propose that the decreased activity of ER and potentially that of its cofactors renders ER-dependent enhancers prone to DNA hypermethylation and histone deacetylation, ultimately tipping the balance from an active to a primed enhancer configuration [START_REF] Hervouet | Specific or not specific recruitment of DNMTs for DNA methylation, an epigenetic dilemma[END_REF][START_REF] Skvortsova | DNA Hypermethylation Encroachment at CpG Island Borders in Cancer Is Predisposed by H3K4 Monomethylation Patterns[END_REF]. The decrease of TET2 and the increase of HDAC9 expression that was observed in response to both E2 deprivation and ICI inhibition could contribute to the above mechanism, although further studies are needed to test this hypothesis.

A putative mark of endocrine resistance i. Partial reversibility of epigenetic changes

Our data also showed that the majority of DNA hypermethylation induced by E2 deprivation could be prevented or even reversed when the cells were re-stimulated with E2 after a period of deprivation (Figure 43, Figure 44). This trend was even more pronounced after ER inhibition by ICI treatment (Figure 48). However, we cannot rule out the possibility that the duration of deprivation followed by re-stimulation was not sufficiently long to observe a full reversal of DNAm changes, or on the opposite, to generate statistically detectable changes pointing towards a non-reversible mechanism. Indeed, we found that DNAm differences were not statistically different at d4 (Table 8). Nevertheless, DNAm levels of the ER-dependent hyperDMPs were in average higher in the re-stimulated group compared to controls (Figure 43C) which cast the doubt on whether DNAm levels are fully reversible or simply prevented from increasing after E2 re-stimulation. As for H3K27ac, there was a subset of DB peaks in the ReSt group that carried a significantly lower H3K27ac signal compared to CTRs and that was in fact as low as the signal of E2Ds (Figure 45). This strengthens the possibility that there is a fraction of enhancers that loses their activity in a long-term manner although longer E2 deprivation and re-stimulation should be applied to confirm a permanent switch-off.

ii. All lights on AP-1

Independently of the enrichment of AP-1 binding within ER's cistrome, the E2 restimulation and ICI removal assays both lead us on a second track involving AP-1.

Indeed, we observed that the expression of two AP-1-related genes, TNIK and MMP13, remained significantly lower following E2 re-stimulation. TNIK is a kinase that activates c-Jun N-terminal kinase (JNK1) in primary human B-cell which in turn activates the AP-1 complex [START_REF] Shkoda | The Germinal Center Kinase TNIK Is Required for Canonical NF-κB and JNK Signaling in B-Cells by the EBV Oncoprotein LMP1 and the CD40 Receptor[END_REF]. Interestingly, it was shown that phosphorylation and JNK1 genomic binding on both on ER and AP-1 binding sites can occur partially in response to E2-induced ER signalling [START_REF] Sun | Estrogen regulates JNK1 genomic localization to control gene expression and cell growth in breast cancer cells[END_REF]. Our results show that TNIK expression is regulated by ER, likely through E2-dependent binding on the only putative enhancer of TNIK (Figure 51). In addition, no changes in the expression of AP-1 component were observed which is in favour of a nongenomic, yet ER-dependent, regulation of AP-1 (data not shown). This opens the door to a missing link between ER and AP-1 that is worth exploring. Little is known about the upstream regulation of MMP13 in the context of breast cancer cells besides the fact it can be a target of AP-1 [START_REF] Datar | RKIP Inhibits Local Breast Cancer Invasion by Antagonizing the Transcriptional Activation of MMP13[END_REF]. High protein levels of MMP13

were detected in more aggressive and metastatic breast tumours that is neither in in direct relation to endocrine resistance nor in line with the very low expression levels measured in MCF-7 [START_REF] Kotepui | Differential expression of matrix metalloproteinase-13 in association with invasion of breast cancer[END_REF].

In agreement with a scenario of non-reversible mechanisms, some studies showed E2 deprivation and ER inhibition can lead to endocrine resistance and this has often been attributed to an overexpression or a redistribution of AP-1 binding sites [START_REF] He | c-Jun/AP-1 overexpression reprograms ERα signaling related to tamoxifen response in ERα-positive breast cancer[END_REF][START_REF] Lupien | Growth factor stimulation induces a distinct ERα cistrome underlying breast cancer endocrine resistance[END_REF][START_REF] Malorni | Blockade of AP-1 Potentiates Endocrine Therapy and Overcomes Resistance[END_REF]. In simple words, AP-1 is not recruited on ER-dependent sites following ER inactivation and that leads to the partial redistribution of the pool of AP-1 complexes on other binding sites. We think that the epigenetic changes we observed occur in these common ER/AP-1 binding sites and that when cells are E2-re-stimulated, a part of these sites is no more accessible to the complex. A study by Fujimoto et al supports this possibility by showing that DNAm on CpG sites that are adjacent to the binding motif of AP-1 can prevent AP-1 from binding in vitro [START_REF] Fujimoto | Methylation adjacent to negatively regulating AP-1 site reactivates TrkA gene expression during cancer progression[END_REF]. Taken together with our findings on AP-1 enrichment in proximity to DMPs may add a new piece to the puzzle of endocrine resistance, further studies with longer durations of treatment are needed to substantiate these findings. One would need to apply our protocol for longer periods of TET2, HDAC9 and TNIK expression levels could contribute indirectly to the enhancers' shift from an active to a primed state. et al., 2012;[START_REF] Swinstead | Steroid Receptors Reprogram FoxA1 Occupancy through Dynamic Chromatin Transitions[END_REF][START_REF] Wang | TET2 coactivates gene expression through demethylation of enhancers[END_REF].

The necessity to re-evaluate long term E2 deprivation and ER inhibition models in the context of research on endocrine resistance

In the optimisation phase of the antagonist exposures, I discovered that ER drug inhibition did not work on cells cultured in CSS-medium in presence of E2 (Figure 8).

At first, I thought that the inhibitors' concentration could not antagonise E2 at the ligand-binding domain of ER but it has been shown that both ICI and TAM can bind ER with a similar affinity to that of E2 [START_REF] Kuiper | Comparison of the Ligand Binding Specificity and Transcript Tissue Distribution of Estrogen Receptors α and β[END_REF]. Furthermore, if the concentration of the antagonists was indeed a problem, the concentration of E2 in non-stripped serum -which can also range around 10nM [START_REF] Sikora | Endocrine Response Phenotypes Are Altered by Charcoal-Stripped Serum Variability[END_REF] should also prevent the inhibitors' action. Yet this was not the case which brought me to follow a second lead. During serum stripping most lipophilic molecules are removed and as a result the polarity of the solution changes and becomes more hydrophobic [START_REF] Sikora | Endocrine Response Phenotypes Are Altered by Charcoal-Stripped Serum Variability[END_REF]. In addition to this, bio-tech providers also point that and it is known that they are highly lipophilic molecules [START_REF] Jordan | New insights into the metabolism of tamoxifen and its role in the treatment and prevention of breast cancer[END_REF][START_REF] Park | Local delivery of hormonal therapy with silastic tubing for prevention and treatment of breast cancer[END_REF]. Based on all the above, it is conceivable that these drugs are not chemically available to the cells when they are added in CSS medium. A solution to test this would be to quantify the amount of the antagonists in the medium 24h after their addition in four groups as following: CSS+E2 medium and MCF-7 vs CSS+E2 medium but no cells and equivalently regular-serum medium and MCF-7 vs regularserum medium and no cells.

In 2015, Stone et al. showed that 6 months of E2 deprivation as well as individual ICI and tamoxifen treatments of MCF-7 lead to the hypermethylation of >15000 CpG sites in common that were importantly located in enhancers [START_REF] Stone | DNA methylation of oestrogen-regulated enhancers defines endocrine sensitivity in breast cancer[END_REF].

However, the cells treated with the two inhibitors were in fact cultured in CSS medium. According to the above observation, it would not be impossible that the extensively inhibited cells underwent in fact the same treatment as E2 deprived cells.

Furthermore, the acquisition of these DMPs was done by comparing the treated cells to a parental MCF-7 that seem not to have been cultured in parallel for an equivalent number of passages. Yet it has been shown repeatedly shown that MCF-7 gradually and spontaneously develop a loss of estrogen responsiveness which is translated into endocrine resistance and for biologically consistent results, MCF-7 should not be cultured for longer than 25 passages [START_REF] Devarajan | Human breast cancer MCF-7 cell line contains inherently drug-resistant subclones with distinct genotypic and phenotypic features[END_REF][START_REF] Kleensang | Genetic variability in a frozen batch of MCF-7 cells invisible in routine authentication affecting cell function[END_REF][START_REF] Nelson | Detection of Endogenous Selective Estrogen Receptor Modulators such as 27-Hydroxycholesterol BT -Estrogen Receptors: Methods and Protocols[END_REF]. For all the reasons stated in this section and given the "popularity" of MCF-7 in the BC scientific community, it would be important to establish and follow standard guidelines of cell culture and why not turn to new normal-like ERexpressing cell lines.

C CONCLUSIONS

• This study capitalises on our optimised cell culture protocol for studying estradiol deprivation and re-stimulation in stricto sensu combined with the latest methodologies for epigenomic profiling and extensive expertise for analysis and interpretation of omics and mechanistic data available in Epigenetics group and collaborating labs.

• We show that E2 deprivation and re-stimulation result in time-dependent DNA methylation changes and in histone modifications across diverse genomic regions, many of which overlapped with enhancers.

• Our RNA-seq and ChIP-seq data reinforced the notion that ER activity may orchestrate a broad transcriptional programme through regulating a limited panel of critical enhancers and that the E2-mediated epigenetic changes are likely involved in priming of enhancers.

• ER activity regulates epigenetic marks on distant enhancers in an indirect way that involves the recruitment of a multitude of transcription factors and epigenetic regulatory factors.

• The maintenance of these epigenetic marks at ER-dependent enhancers are likely to be the result of a balance between an active recruitment of activators and the passive repulsion of inhibiting factors.

• While observed epigenetic changes induced by estrogen deprivation are largely reversible, some changes persist following a sequential E2 deprivation and re-stimulation consistent with the existence of "epigenetic memory" of estrogen-dependent gene activity states.

• Our enrichment analysis of TF binding and motif occurrence highlighted the importance of ER tethering mainly through two partner TF families, AP-1 and FOX, in the proximity of enhancers that undergo epigenetic reconfiguration.

• Based on our results, we proposed a new mechanistic model in which ligandactivated ER and its cofactors orchestrate an intricate interaction in cis and maybe trans between promoters and a panel of epigenetically dynamic enhancers.

• Collectively our study provides entirely new information that constitutes a major advance in our understanding of the events by which ER and its cofactors mediate changes in DNA methylation and chromatin states at enhancers.

• This new insight gained by our study opens new avenues for studying the role of the deregulated estrogen signalling in both breast cancer development and the mechanism underlying the "roots" of endocrine resistance that commonly develops in response to anti-estrogen therapy

Figure 1 .

 1 Figure 1. Worldwide age-standardised incidence rates of breast cancer in 2018 (ASR: Agestandardised rate) 1 .

Figure 2 .

 2 Figure 2. Schematic representation of the adult mammary gland and its main cellular subtypes. Schematic depiction of potential regulatory cross-talk between different ductal mammary cells in response to steroid hormone stimulation. Steroid hormones (red circles in lumen) activate ER+ epithelial cells (either mature cells or progenitors). BCs are classified in molecular and histopathologic subtypes. (Adapted from[START_REF] Visvader | Mammary stem cells and the differentiation hierarchy : current status and perspectives[END_REF]

Figure 3 .

 3 Figure 3. Model depicting the main molecular mechanisms through which ER can regulate gene expression and cell signalling.(From Heldring et al., 2007)

Figure 4 .

 4 Figure 4. A simplified summary of histone modifications and DNAm used as traditional predictors of chromatin states. Methylated and unmethylated CpG sites are represented as full or empty lollipop-circles.

Figure 5 .

 5 Figure 5. Mechanistic hypothesis of the maintenance of epigenetic marks and gene regulation in relation to E2-dependent ER activity.

Figure 6 .

 6 Figure 6. Experimental design for the study of E2 deprivation and re-stimulation

Figure 7 .

 7 Figure 7. Experimental design for the study of E2 deprivation and re-stimulation

Figure 8 .

 8 Figure 8. Expression of ESR1 and ER target gene GREB1 in response to 1µM of ICI 182,780 (ICI) and 1µM of tamoxifen (TAM) after 5 days of treatment. (A) MCF-7 were cultured in phenol-re-free medium, CSS and 10nM of E2 (CTR) in addition to an antagonist. (B) MCF-7 were cultured in standard medium and FBS in addition to an ER antagonist. Results are shown as the average expression of triplicates ± SD.

Figure 9 .

 9 Figure 9. Experimental design for the study of ICI 182,780 (ICI) inhibition and re-activation.

Figure 10 .

 10 Figure 10. DNA fragment size distribution after sonication of chromatin for 6 to 16 cycles (30sec ON, 30sec OFF Bioruptor Pico)

.

  DNAm and transcription correlation FDR-significant DMPs and DEGs were correlated within CTR and E2D samples (n=12) at with and without any differential cut-off (Δβ ≥10% for DNAm and |log 2 FC| >1 for gene expression), using an FDR-adjusted Spearman correlation. DNAm and transcription data were obtained from the same samples. DMPs were matched to the closest expressed gene, upstream or downstream. iv. Statistics for targeted experiments For all targeted experiments (pyrosequencing, RT-qPCR and ChIP-PCR) I performed pairwise non-parametric tests, namely Mann-Whitney test. For DNAm data and H3K27ac-ChIP, only timepoint d14 was considered (n=3 per group) and for relative mRNA expression, both d4 and d14 were considered (n=6 per group). P-values <0.05 were considered as significant. Statistical comparisons were performed with R version 3.5.1.

Figure 12 .

 12 Figure 12. Quality control of samples based on array signal intensity.

Figure 13 .

 13 Figure 13. Distribution of DNAm levels before normalisation

Figure 14 .

 14 Figure 14. Distribution of beta and M-values after Functional Normalisation

Figure 15

 15 Figure 15 Proportion of variability explained from 10 first principal components

Figure 17 .

 17 Figure 17. Overlap of DMPs associated to treatment and those associated to time.

Figure 18 .Figure 19 .

 1819 Figure 18. QQ plot showing the inflation of p-values of the comparison E2D vs CTR (A) at d4 (λ=0.64) and (B) at d14 (λ=2.15)

Figure 21 .

 21 Figure 21. (A) Genomic features distribution and (B) CpG density of hypermethylated DMPs (hyperDMPs, n= 952)), all 850k CpG probes (n=866836) and hypomethylated DMPs (hypoDMPs, n=62). (C) Fold-enrichment of DMPs overlapping with publicly available histone ChIP-seq data and (D) ChromHMM annotations originated from MCF-7 cells (Taberlay et al., 2014). Asterisks (*) mark significant enrichments (p< 10e-5, Fisher test). Distributions and enrichments are shown for DMPs with ≥10% Δβ.

Figure 22 .

 22 Figure 22. Distribution of DNAm levels of DMPs before E2-deprivation (CTR-d0 replicates 1 to 3)

Figure 24 .

 24 Figure 24. Expression of ESR1 and ER targets after 5 days after siESR1 transfection. Data shown as average ± SD of triplicates. Significance was tested using a Mann-Whitney test.

Figure 27 .B

 27 Figure 27. PCA of normalized read counts (A) for all the binding sites and samples, and (B) for DB H3K27ac regions (FDR <0.05) specifically for CTR and E2D contrast

Figure 30 .

 30 Figure 30. Genome browser snapshot of H3K27ac normalised ChIP-seq reads in CTR and E2D conditions, ER and FOXA1 ChIP-seq reads in E2-treated and untreated MCF-7 (*GSE72249, Swinstead et al., 2016) and FOS ChIP-seq (°GSE105734, ENCODE Project Consortium et al., 2012).

Figure 31 .

 31 Figure 31. (A) Distribution of -log 10 (p-values) of differentially expressed genes (DEGs) according to log 2 fold-change of expression. Coloured dots represent down-and up-regulated DEGs with an absolute log 2 (FC) >1 (blue and red) based on a FDR <0.05 (dashed horizontal). Genes that were differentially expressed with an absolute log 2 (FC) <1 are coloured in dark gray. (B) Heatmap showing top 50 DEGs (FDR <0.05) for CTR and E2D at d4 and d14.

Figure 32 .

 32 Figure 32. Technical validation of selected genes from the RNA-seq by quantitative RT-qPCR. The effect of treatment was tested using a Mann-Whitney test (* p <0.05; ** p <0.01). Results are shown as the average of triplicates ± SD.

Figure 33 .

 33 Figure 33. Gene expression of epigenetic remodeling factors HDAC9, PADI4 and TET2 in response to (A) E2 deprivation found (RNA-seq, n=3, FDR < 0.05) and (B) to ICI inhibition (RT-qPCR, Mann-Whitney CTR vs ICI, * p <0.01), for 4 and 14 days. Data of ICI treatment are shown as the average of triplicates ± SD.

Figure 34 .

 34 Figure 34. Correlation of differential methylation and expression. Significantly correlating DMPs and DEGs (Spearman |ρ| > 0.67, FDR < 0.05) are depicted in pink and red points indicating correlating an absolute |log 2 FC|>1 and |Δβ| >0.1. DMPs matching nondifferentially expressed genes are shown in grey.

Figure 35 .

 35 Figure 35. Correlation of DNAm of cg10612997 and cg26207870 and gene expression correlation to their closest gene (GREB1 and PKIB, respectively).

Figure 37 .

 37 Figure37. Enrichment of TF binding on the nearest genes to DMPs (n=829 genes), the nearest genes to differentially bound H3K27ac marks (n=577, ≥2 peaks per gene name) and DEGs (n=547). Numbers represent the amount of TFs whose binding was significantly enriched across the different sets of genes (FDR <0.05). Enrichment of TFs originating from experiments performed in MCF-7 cells are highlighted in red in the upper right panel(ChEA 2016[START_REF] Kuleshov | Enrichr: a comprehensive gene set enrichment analysis web server 2016 update[END_REF].

Figure 38 .

 38 Figure 38. Enrichment of TF binding down-and up-regulated genes (log 2 FC <-1 and log 2 FC >1, FDR <0.05, n=381 and n=166). The enrichment score is the combination of log(p-value) and Z-score of expected rank (ChEA 2016, Kuleshov et al., 2016). Bars surrounded with a black line indicate datasets that are originated from MCF-7 cells.

Figure 39 .Figure 40 .

 3940 Figure 39. De novo motif analysis performed on DMPs (top 4 hits shown) and DB H3K27ac (|FC| ≥ 1, top 2 hits shown) obtained after 14 days of E2 deprivation. The enrichment is the result of the percentage of the motif occurrence in target sequences over random background genomic sequences. The motif search was expanded over a 500bp-window for DMPs and a 1kb-window for DB H3K27ac peaks.

Figure 41 .

 41 Figure 41. De novo motifs detected in regions bound by ER and known cofactors.

Figure 42 .

 42 Figure 42. Average DNAm distribution in CpG sites bound by TFs before E2 deprivation at d0.

  a) a model taking both into account treatment and time as separate variables, comparing ReSt to CTR and E2D across all normalised probes, pdata$treatment <-relevel(pdata$treatment, ref = "ReSt") pdata$time <-relevel(pdata$time, ref = "d14") design <-model.matrix(~as.factor(treatment) +as.factor(time), data=pdata) fit <-lmFit(mval, design) b) a model where ReSt was compared against CTR and E2D only within timepoint d14 and pdata$group <-relevel(pdata$group, ref = "ReSt_d14") design <-model.matrix(~as.factor(group), data=pdata) fit <-lmFit(mval, design) c) a model where ReSt was compared against CTR and E2D on d14, and where we included only probes corresponding to the 10% Δβ hypermethylated DMPs obtained in the CTR vs E2D analysis. pdata$group <-relevel(pdata$group, ref = "ReSt_d14") design <-model.matrix(~as.factor(group), data=pdata) fit <-lmFit(mval_hyperDMPs, design) Briefly, for the first model, the ReSt group (n=3) was compared to CTR (n=9) and E2D (n=6), taking time into account in the linear regression model. No DMPs were found between CTR and ReSt but there were 51 DMPs when comparing ReSt and E2D data (Table

  and C. Validation with pyrosequencing recapitulated the trend with the majority of CpGs in ReSt group showing a return to DNAm levels equivalent to those of CTR with the exception of the CpG matching the cg04210444 probe where DNAm levels of ReSt were closer to E2D group (third panel, Figure 44).

Figure 43 .

 43 Figure 43. Hypermethylation due to E2 deprivation is in majority reversed or precluded after E2 re-stimulation (A) MCF-7 were deprived of E2 for 4 days, after which they were restimulated with E2 for 10 days (ReSt, blue dashed line). (B) Heatmap and (C) boxplot of hypermethylated DMPs (FDR < 0.05; > 10% Δβ) in response to E2-deprivation for CTR, E2D and ReSt at d14.

Figure 44 .

 44 Figure 44. Technical validation of DNAm levels. Pyrosequencing analysis of indicated DMPs in CTR, E2D and ReSt (green, orange and blue, respectively) .

Figure 45 .

 45 Figure 45. Partial recovery of H3K27ac signal after 4 days of E2 deprivation followed by 10 days of E2 re-stimulation (A) H3K27ac signal in log2 normalised reads of 347 peaks that were significantly different in both CTR vs ReSt and CTR vs E2D comparisons. (B) Genome browser snapshot of H3K27ac normalised ChIP-seq reads in CTR, E2D and ReSt conditions at d14 (chr1:24701801-24702801)

Figure 46 .Figure 47 .

 4647 Figure 46. Gene expression of AP-1 inducer TNIK and target gene MMP13 in E2 deprived and re-stimulated conditions on day 14 (RNA-seq, CTR vs ReSt, FDR < 0.05, |log 2 FC| > 1)

Figure 48 .

 48 Figure 48. DNAm dynamics of top DMPs (title: 850k probe id) following DMSO treatment for 4 and14 days (CTR, green), ICI 182,780 treatment for 14 days (ICI, orange), and 4 days of ICI 182,780 treatment followed by 10 days of DMSO (ReAc, blue dashed). Results are shown as the mean of triplicates with 95% of C.I. Significance of changes was calculated at d14 using a Mann-Whitney test (*p < 0.05).

Figure 49 .

 49 Figure 49. Gene expression of AP-1 inducer TNIK and target gene MMP13 in following treatment (ICI) and removal (ReAc) of ICI (RT-qPCR, Mann-Whitney * p < 0.05 and ** p < 0.01).

Figure 50 .

 50 Figure 50. Gene expression of epigenetic remodelling factors HDAC9, PADI4 and TET2 in response to ER inhibition and reactivation by ICI 182,780, for 4 and 14 days 7 . Data of ICI treatment are shown as the average of triplicates ± SD.

Figure 52 .

 52 Figure52. Proposed mechanistic model for shift from active to primed chromatin state in response to E2 deprivation and re-stimulation. In this model the reduction of liganddependent ER activity leads to an AP-1 and FOXA1 mediated increase of DNAm and decrease of H3K27ac and transcription. The reversibility of these changes is partial. This model include findings of this thesis work and integrates prior knowledge of ER-dependent chromatin and TF interactions(ENCODE Project Consortium et al., 2012;[START_REF] Swinstead | Steroid Receptors Reprogram FoxA1 Occupancy through Dynamic Chromatin Transitions[END_REF][START_REF] Wang | TET2 coactivates gene expression through demethylation of enhancers[END_REF].

  serum composition can affect cell membrane permeability and therefore the transfer of lipids into the cell 9 . Besides this, the chemical structures of E2, ICI 182,780 and tamoxifen are different enough to have different hydrophilic features as the three molecules have low, yet different water solubility (Figure 53, E2: 3.6 mg/L 10 ; ICI: 9.53x10 -6 mg/L 11 , TAM: NA). Finally, as ICI 182,780 and tamoxifen are commonly used drugs in the clinical setting, their chemical and physical features are well studied

Figure 53 .

 53 Figure 53. Chemical structures of E2, TAM and ICI.

Table 1 .

 1 Primer sequences to measure DNA methylation by pyrosequencing. Primers carrying a biotine modification at their 5' contain "BIO" in their name

			Sequencing		Amplicon	Annealing
	Primer	Sequence	primer		Genomic range	(bp)	Temp (°C)
		AGTTATGATGGGGT	Same	as	chr11:85,884,137-
	cg12667152_F	GGAGAAA	forward		85,884,246	110	<55.8
		AAAATCCCTACTAA
	cg12667152_BIO_R	TCCACACAAA
		TGTGTTATTGATGA	TTAGGGGAGA	chr2:11,673,902-
	cg10612997_F	AGTTTTTGAG	ATATG		11,674,119	218	61
	cg10612997_BIO_R	ACTAAAATCTACAC

Table 2 .

 2 Primer sequences used to measure gene expression

	Primer	Sequence	Amplicon (bp)	Source
	qRT_GREB1_F	ATGGGAAATTCTTACGCTGGAC	171 designed in this study
	qRT_GREB1_R	CACTCGGCTACCACCTTCT
	qRT_TFF1_F	TCGAAACAGCAGCCCTTATT	102 designed in this study
	qRT_TFF1_R	GGCCCAGACAGAGACGTGTA
	qRT_PGR_F	TTATGGTGTCCTTACCTGTGGG	112 designed in this study
	qRT_PGR_R	GCGGATTTTATCAACGATGCAG
	qRT_FOS_F	AAGCGGAGACAGACCAACTA	134 designed in this study
	qRT_FOS_R	CAGGTCATCAGGGATCTTG
	qRT_JUN_F	AACAGGTGGCACAGCTTAAAC	77 PrimerBank ID 44890066c2
	qRT_JUN_R	CAACTGCTGCGTTAGCATGAG
	qRT_HDAC9_F	AGTAGAGAGGCATCGCAGAGA	141 PrimerBank 116284378c1
	qRT_HDAC9_R	GGAGTGTCTTTCGTTGCTGAT
	qRT_PADI4_F	CAGGGGACATTGATCCGTGTG	130 PrimerBank 216548486c1
	qRT_PADI4_R	GGGAGGCGTTGATGCTGAA
	qRT_MMP13_F	ATTTCTCGGAGCCTCTCAGT	103 designed in this study
	qRT_MMP13_R	AGTTTGCAGAGCGCTACCT
	qRT_RPLP0_F	TGGCAGCATCTACAACCCTGAA	90 designed in this study
	qRT_RPLP0_R	ACACTGGCAACATTGCGGACA
	qRT_TNIK_F	TAAGGGTCGTCATGTCAAAACG	173 PrimerBank 239735578c1
	qRT_TNIK_R	CCATGCCTGGTGGGTTCTTT
	qRT_TET2_F	CTTTCCTCCCTGGAGAACAGCTC	146 designed in this study
	qRT_TET2_R	TGCTGGGACTGCTGCATGACT

Table 3 .

 3 Primer sequences used to measure histone mark enrichment by qPCR

			Amplicon
	Primer	Sequence	(bp)	Source
	Peak1_F	GCTTGCAGGCAGAAATGTAAACAC	136 designed in this study
	Peak1_R	CTGTTGCCCTAAGCCTGTATT
	Peak2_F	GTTGTGTGGACCATATTTTCTAAAC	102 designed in this study
	Peak2_R	TGTTTACATTTCTGCCTGCA
	Peak3_F	TTAAACAGTTCTCCTCCACTCC	74 designed in this study
	Peak3_R	GCAGATTACAAAACTCACTGAAGAC
	Peak4_F	TGATGTGTCAGGAACCCTCC	92 designed in this study
	Peak4_R	CCCAGACCTACTGCATCAGAAAC
	GREB1_pos_F	CACGTCCCCACCTCACTG	56 Jozwik, 2016, Cell Reports
	GREB1_pos_R	TGTTCAGCTTCGGGACACC
	SPATA16_neg_F	GACGCCAGGACTTGAAAC	102 designed in this study
	SPATA16_neg_R	TCAAGTGCAACTCAAGAGGA

Table 4 .

 4 List of ChIP-seq datasets produced in MCF-7 used in this work.

Table 5 .

 5 Summary of samples and different variables used in the subsequent analysis (referred to in the script as pdata)

	Sample name barcode	treatment	time group	sentrix ID	sentrix position
	CTR_D0_1	201533580011_R01C01	CTR	d0	CTR_d0	201533580011 R01C01
	CTR_D0_3	201533580011_R08C01	CTR	d0	CTR_d0	201533580011 R08C01

Table 6 .

 6 Significance of variables explaining 5 first principal components (p-values of Wilcoxon and Kruskal-Wallis tests for one level and multilevel variables respectively).

	A		B
	Variable	PC1	PC2	PC3	PC4	PC5
	sentrix	0.161	0.477	0.144	0.018	0.526
	sentrix_position	0.297	0.576	0.752	0.931	0.243

Table 7 .

 7 DMPs in response to E2-deprivation over time (reference: CTR and d0) (FDR=0.05) Figure 16. QQ plot observed against expected p-values of DMPs resulting from the comparison of CTR vs E2D (λ=1.84)

		E2D	d14	d4
	Hypo	745	241	0
	NotSig	819715	822293	822693
	Hyper	2233	159	0

Table 8 .

 8 DMPs in response to E2-deprivation on d4 (reference: CTR_d4 in contrast to E2D_d4) (FDR=0.05)

		CTR_d0	CTR_d14	E2D_d14	E2D_d4
	Hypo	0	0	5661	0
	NotSig	822693	822693	804857	822692
	Hyper	0	0	12175	1
	Table 9. DMPs in response to E2-deprivation on d14 (reference: CTR_d14 in contrast with
	E2D_d14) (FDR=0.05)
		CTR_d4	CTR_d0	E2D_d14	E2D_d4
	Hypo	0	48	802	0
	NotSig	822693	822493	818843	822689
	Hyper	0	152	3048	4

Table 10 .

 10 Overlap of DMPs with H3K27ac peaks

		DMPs (n=1014)	850k (n=866,836)
	H3K27ac	CTR d14	DB peaks	CTR d14	DB peaks
	window size +	(n=31,090)	(n=3053)	(n=31,090)	(n=3053)
	1 kb	10.9%	0.3%	10.7%	0.3%
	10 kb	34.9%	5.0%	27.1%	1.8%
	100 kb	84.9%	23.0%	65.1%	11.0%

+ Total genomic window evenly centred on peak summit.

Table 11 .

 11 DMPs in response to E2 re-stimulation (reference: ReSt and d14) (FDR=0.05)

	CTR	E2D	d0	d4
	Hypo	0	0	159	266
	NotSig	822693	822642	822293	821825
	Hyper	0	51	241	602

Table 12

 12

	. DMPs in response to E2 re-stimulation on d14 (reference: ReSt_d14) (FDR=0.05)
		CTR_d14	E2D_d14
	Hypo	0	23
	NotSig	822690	822014
	Hyper	0	656

Table 13 .

 13 DMPs in response to E2 re-stimulation only among E2D-responsive hypermethylated DMPs with 10% Δβ on d14 (FDR=0.05).

	CTR_d14	E2D_14
	Hypo	20	0
	NotSig	932	7
	Hyper	0	945

Data source: GLOBOCAN

2018; graph production: International Agency for Research on Cancer (http://gco.iarc.fr/today)

Percentages of mutations and copy number alterations were summed from http://www.cbioportal.org data on 817 tumours.

http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/

Method for calculation found here : https://stats.stackexchange.com/questions/110755/calculate-inflation-observed-andexpected-p-values-from-uniform-distribution-in

ChromHMM is a software based on a hidden Markov model using multiple variables to assign functional chromatin states to genomic regions[START_REF] Ernst | ChromHMM: automating chromatin-state discovery and characterization[END_REF]. For this, different epigenomic datasets are used as input, mainly ChIP-seq but also DNAse-seq and NOMe-seq.

CTR and ICI samples are the same as those used in Figure33

In majority ATCC, HTB-22 subtype was used for standard protocols. ECACC MCF-7 subtype (catalog#: 86012803) was used for E2 deprivation and re-stimulation protocols. Information found on www.encodeproject.org/documents/ under growth protocols, filtered for MCF-7.

https://www.thermofisher.com/fr/fr/home/references/gibco-cell-culture-basics/cellculture-environment/culture-media.html

10 National Center for Biotechnology Information. PubChem Database. Estradiol, CID=5757, https://pubchem.ncbi.nlm.nih.gov/compound/5757 (accessed on May 20,

2019) 11 National Center for Biotechnology Information. PubChem Database. Fulvestrant, CID=104741, https://pubchem.ncbi.nlm.nih.gov/compound/104741 (accessed on May 20, 2019)

The PhD student Athena Sklias was supported by the Fonds National de la Recherche, Luxembourg (grant number: 10100060). The experiments were partially supported by grants from the Institut National du Cancer (INCa, France) and the European Commission (EC) Seventh Framework Programme (FP7) Translational Cancer Research (TRANSCAN) Framework and the Fondation ARC pour la Recherche sur le Cancer (France) to Zdenko Herceg.

List of abbreviations

BC:

Breast Cancer CGI:

progressively more pronounced pointing towards an active adaptation of the epigenetic landscape.

i. Epigenetic changes occur in enhancers and reflect TF binding

Thanks to the genome-wide coverage of 850k BeadChip array and H3K27ac ChIPsequencing, we were not only able to measure numerous significant changes of epigenetic modifications but also to observe that these changes were significantly enriched in enhancers. The latter were identified mostly with the use of publically binding of ER at a nucleus-wide scale using a newly developed microscopy technique that allows measuring the time of residence of a TF on DNA [START_REF] Swinstead | Steroid Receptors Reprogram FoxA1 Occupancy through Dynamic Chromatin Transitions[END_REF].

Together, these observations led to postulate that these epigenetically variable enhancers could be shaped by dynamic E2-mediated ER binding.

In addition, the overlap between genes with expression and epigenetic changes in proximal cis-regulatory regions was limited (Figure 36). We therefore turned to an analysis of the TF network of ER given that their binding to ER-dependent enhancers is also affected by E2 deprivation, and consequently could impact the epigenetic landscape. In line with the limited overlap between ER binding and differential epigenetic marks, our analysis of enriched binding motifs of TFs did not identify ER motifs in proximity to DMPs and DB H3K27ac peaks. On the contrary, AP-1 and and with regular sample collections to capture the switch towards endocrine resistance setting of the methylome, chr3:170,971,567-170,988,176) iii.

Proposition of a novel molecular mechanism

Based on all our results, we propose a model in which ligand-activated ER and its cofactors orchestrate an intricate interaction in cis between promoters and a panel of epigenetically dynamic enhancers (Figure 5).