Steroid hormones regulate genome-wide epigenetic programming and gene transcription in human endometrial cells with marked aberrancies in endometriosis

Endometrium is a dynamic tissue whose cellular components undergo cyclic proliferation and differentiation, preparing for embryo implantation by highly coordinated spatiotemporal actions of ovarian-derived estradiol (E 2 ) and progesterone (P 4 ) [ 1 , 2 ]. These hormones bind cognate receptors [estrogen receptor (ER) and progesterone receptor (PR)], whose activities are tightly regulated by post-translational modifications and interactions with cell- and tissue-specific co-regulators [ 3 – 5 ]. Binding ER and PR leads to their nuclear translocation, complexing with nuclear response elements, remodeling chromatin by co-modulator recruitment [ 3 ], and alteration of the transcriptional machinery. In endometrium, dynamic circulating E 2 and P 4 levels drive the normal functionality of the tissue. Moreover, environmental and inflammatory signals can alter steroid hormone-driven endometrial gene transcription and cellular function resulting in tissue dyshomeostasis [ 6 ], including endometrial hyperplasia and cancer, endometrial-based infertility, endometriosis, and poor pregnancy outcomes [ 7 , 8 ]. While changes in chromatin accessibility, PR targets and changes in histones and gene expression in eSF decidualization by E 2 , cAMP and MPA have been shown [ 9 – 11 ], how E 2 and P 4 individually interact with the endometrial epigenome normally or in inflammatory disorders that compromise endometrial function, e.g., as in the disorder endometriosis, are incompletely understood. We hypothesized that these steroid hormones induce unique genome-wide signatures in normal human endometrial stromal fibroblasts with aberrant signatures in endometrial cells from women with endometriosis, and their effect on the epigenome is directed by specific genomic sequences and locations. Endometriosis is a common, chronic disorder wherein endometrial tissue, shed into the pelvis at menses, elicits an inflammatory response, neuroangiogenesis and fibrosis, resulting in infertility and chronic pelvic pain [ 12 ]. Hallmarks of the disorder are its dependence on E 2 for growth, disrupted P 4 signaling caused by chronic inflammation in endometriosis lesions and in the eutopic endometrium (uterine lining) [ 13 ], and epigenetic chromatin changes that determine endometrial cellular responses to mitogenic and differentiative signals [ 6 , 12 , 14 ]. Normally, the eutopic endometrial DNA methylome varies according to the hormonal milieu, with greatest differences in the E 2 -dominant (proliferative) versus P 4 -dominant (secretory) phase of the cycle and associated with gene expression changes [ 10 , 13 – 16 ]. Chronic inflammation affects the chromatin landscape in endometrium of women and animal models of endometriosis [ 6 ]. These observations on bulk tissue raise fundamental questions about steroid hormone-epigenome interactions in cellular components of the endometrium normally and in women with disease, how steroid hormones affect the epigenome, how the epigenome affects steroid hormone response and action, and if there are epigenetic differences in endometrium of women with endometriosis and different stages of disease, and if so, what role they play in these processes. Herein, we studied responses of endometrial stromal fibroblasts (eSF) isolated from normal women (controls) and those with endometriosis. eSF comprise a major endometrial cell type whose programmed responses to E 2 and P 4 are essential for pregnancy success and whose responses are compromised in inflammatory disorders [ 17 ], including endometriosis [ 18 , 19 ]. Given the centrality of a normal eSF P 4 response for pregnancy and that women with endometriosis have infertility and poor pregnancy outcomes reported by some to be due, in part, to altered eSF basal gene expression and abnormal response to P 4 [ 18 – 21 ] , understanding steroid hormone-signaling and regulation of transcription in this cell type is paramount. Moreover, controls for this study were women with no known gynecologic or systemic disorders, thereby enabling establishment of a normative platform for steroid hormone effects on the epigenome and gene transcription in this cell type in endometrium, the tissue that is the anatomic pre-requisite for continuation of the species.

We assessed the genome-wide effect of individual steroid hormones, E 2 and P 4 , and their combination (E 2 + P 4 ) on endometrial stromal fibroblasts (eSF) after 14 days of exposure mimicking the timeframe in the menstrual cycle for maximal hormone responsiveness. Since the effects of E 2 and P 4 individually and together on the hormone-epigenome interplay in normal endometrial cells was unknown and was a main aim of this study, we applied stringent criteria and utilized only endometrial samples from extensively screened volunteers without any gynecologic disorders and no uterine pathology (NUP), with confirmed in vitro eSF progesterone responsiveness (see methods and S1 Fig ). Interrogation of 485,577 methylation targets across the genome revealed that E 2 and P 4 and their combination affected the DNA methylome to different extents and with distinct patterns in eSF normal . (Note: throughout the Results section steroid responses of eSF DNA methylomes are compared to untreated (vehicle) cells for each group.) E 2 induced the most extensive changes in the eSF normal DNA methylome (2047 CpG sites), followed by combined E 2 +P 4 (569 CpG sites), and P 4 alone had the least effect (505 CpG sites) ( Fig 1A and 1B , S1 Table ). Importantly, combined E 2 +P 4 resulted in dramatically reduced numbers of differentially methylated loci compared to E 2 alone ( Fig 1A ). While individual hormone treatments elicited hormone-specific DNA methylome changes, the simultaneous presence of both hormones altered their individual effects. The pattern of loss and gain of methylation is also distinct for each hormone. E 2 , and E 2 +P 4 induced more loss than gain of methylation ( Fig 1A , yellow: gain of methylation, blue: loss of methylation vs vehicle), while P 4 induced similar numbers of loss and gain of methylation. Concordant with differential patterns and extents of methylation changes, we found minimal overlap in the differentially methylated CpG sites affected by each hormone and the majority was unique ( Fig 1C , S1 Table ). In particular, loci differentially methylated in response to E 2 +P 4 were mostly unique compared to those in response to P 4 or E 2 alone and were not a combination of the response to each hormone individually ( Fig 1C ). 1A. Differentially methylated CpG sites induced by E 2 , P 4 and E 2 +P 4 versus vehicle. Each heatmap reflects differential methylation of each sample in each hormone treatment versus its corresponding non-treated vehicle control (Δβ: Hormone treated minus vehicle control). Yellow heatmaps above the X-axis reflect gain of methylation vs. vehicle; blue heatmaps below the X-axis reflect Δβ loss of methylation. In each heatmap, rows show Δβ of differentially methylated loci, columns indicate samples. Y-axis shows the number of differentially methylated loci for either gain or loss of methylation for each hormone treatment. 1B. Number of differentially methylated CpG sites and in gain/loss of methylation for each hormone treatment. 1C. Venn diagram of unique and common differentially methylated CpG sites for each hormone shows little overlap between differentially methylated loci in each hormone treatment 1D. Enrichment of intergenic regions in % in each hormone treatment for all differentially methylated loci (All Loci), those with gain or loss of methylation (Gain, Loss) and by individual hormones (E 2 , P 4 , E 2 +P 4 ). Enrichment is assessed by Z-test and p<0.05 are shown in parentheses. Black bar represents percentage of intergenic loci of total interrogated CpG sites. 1E. Statistically significant involvement of enhancers by hormones and gain or loss of methylation. Enrichment is assessed by Z-test and p<0.05 are shown in parentheses. Black bar represents percentage of enhancers of total interrogated CpG sites. 1F. Genomic distribution of all differentially methylated CpG sites in each hormone and by gain or loss of methylation, assessed at TSS1500, TSS200, 5’UTR, 1 st exon, gene body, 3’UTR, and intergenic regions. Black line represents the percentage of interrogated CpG site at each location, green line (top panel) shows all differentially methylated loci in NUP for all hormones, yellow line (middle panel) shows all loci with gain of methylation in NUP, and blue line (bottom panel) shows all loci with loss of methylation for all hormones. Enrichment is assessed by Z-test and p<0.05 are shown in parentheses for each genomic location. 1G. Distribution of differentially methylated CpGs by CpG islands (CGI), CGI north/south shores and shelves for all loci with gain of methylation in all hormone treatments (orange lines) or with loss of methylation (blue line) in comparison to the distribution of the interrogated CpG sites in each of these genomic locations (black line). N Shelf: North Shelf; S Shelf: South Shelf; N Shore: North Shore; S Shore: South Shore. NUP: normal (no uterine pathology). Pathways and biofunctions ( Table 1 ) as well as functional enrichment clustering ( Table 2 ) were also unique to each hormone with E 2 pathways enriching for gap junctions, melanogenesis, and glutamatergic and dopaminergic synapses pathways, and zinc and ion binding, cell membrane, glycoprotein and signal peptide functional clusters, with fewer and unique statistically significant pathways and functional clusters affected by P 4 and E 2 +P 4 . Together these data indicate that each hormone affects different regions and E 2 +P 4 targets are not a combination of E 2 and P 4 . Differentially methylated loci in all hormonally treated eSF normal involved several pathways, many important in normal endometrial function and dysfunction. Important pathways affected by each hormone and the differentially methylated genes in each pathway are shown in S1 Data . The data were further mined for differences in DNA methylation patterns, profiles, and genomic locations, regulatory elements, transcribed genes and biofunctions induced by each hormonal treatment in cells from normal and endometriosis women (see below). Interestingly, while the patterns, profiles, differentially methylated CpG sites, pathways and biofunctions were unique to each hormone, the genome-wide distribution of their affected CpG sites shows specific enrichments and depletions. All hormones (E 2 , P 4 and E 2 +P 4 ) were statistically significantly enriched in intergenic regions ( Fig 1D ) and in enhancers ( Fig 1E ) , albeit with different extents and in gain vs loss of methylation and some variations based on hormones. These may reflect hormone binding sites in these regions, as had been reported in breast cancer cell lines [ 22 , 23 ]. There was a marked depletion of differential methylation in close proximity to transcription start sites (TSS) up to -200 nt upstream (TSS200) for all hormonal treatments in both of gain or loss of methylation ( Fig 1F ; S2 Fig for gain/loss for each hormone ; S2 Table ). But, CpG sites with gain of methylation in all hormonal treatments exhibited low enrichment at 5’UTRs, 1 st exons and gene bodies, while in loss of methylation TSS1500 and 1 st exons were less involved but gene bodies and 3’UTRs were more enriched. Greatest differences in gain versus loss of methylation ( Fig 1F ) involved gene bodies, and 3’UTRs, and less at TSS1500, TSS200, 5’UTRs, and 1 st exons. CpG islands (CGI) , CGI shores and shelves . For all hormones, there was low involvement of CGIs and CGI shelves and shores, compared to interrogated loci on the HM450 platform ( Fig 1G and S3 Table ). Most DNA methylation changes in eSF normal did not involve CGIs. Indeed, while 31% of all interrogated loci were at CGIs and 33% at CGI shores and shelves (total CGI-related 64%), only 5–7% of the differentially methylated loci for any hormonal treatment were located at CGIs (33–40% overall in CGIs, shores and shelves, compared to 64% arrayed on the platform). The majority of differentially methylated CpG sites involved non-CGIs (59–66%) versus 36% non-CGI CpG sites on the platform. However, there were more differentially methylated CpG sites at CGI shores and shelves than in CGIs, in both gain and loss of methylation. Notably, more CGIs, less CGI shelves (north and south shelves), more CGI shores (north and south shores), and less non-CGI CpG sites were involved in gain versus loss of methylation. Steroid hormones affect their target genes through various mechanisms and as such, changes in DNA methylation may not fully reflect their effect on changes in gene expression, particularly in the case of those loci whose transcriptional regulation does not involve chromatin modifiers. To elucidate a more complete effect of hormones, transcriptomic profiles were determined in the same steroid hormone-treated eSF used for DNA methylation analysis and was compared to its corresponding gene expression profiles in untreated eSF. E 2 induced more up- than down-regulated genes ( Table 3 , top up- and down-regulated loci; S4 Table , full gene list), and P 4 elicited similar numbers of up- and down-regulated genes, ( S4 Table ). However, more genes were differentially expressed when E 2 and P 4 were combined, with more genes up- than down-regulated ( S4 Table ). E 2 +P 4 induced the largest and P 4 the smallest changes in gene expression. Almost all the genes differentially induced by P 4 were shared with E 2 +P 4 and some were shared with the E 2 treatment. Notably, half of the E 2 -induced and the majority of E 2 +P 4 induced differentially expressed genes were unique. However, in commonly up-regulated genes between E 2 and E 2 +P 4 , the variable fold changes indicate inhibitory or stimulatory effects when E 2 is combined with P 4 . For example, PGR is upregulated by both E 2 and E 2 +P 4 (FC = 4.5 vs 1.7, respectively), indicating that the addition of P 4 limited up-regulation of PGR compared to E 2 alone. Among other up-regulated gene in common in E 2 , P 4 , E 2 +P 4 , are IGF1 and SPARCL1 with known roles in endometrial biology. But, the FCs were different ( IGF1 : E 2 = 21.5, P 4 = 4.6, E 2 +P 4 = 13.6; SPARCL1 : E 2 = 4.6, P 4 = 13.2; E 2 +P 4 = 53.1) indicating potentially different mechanisms for up- or down-regulation and for different genes, potentially affected by genomic location and other regulatory factors, modifiers and gene/region-specific mechanisms involved in hormone-induced gene expression regulation. E 2 increased tissue and cellular development, growth and maintenance, and downregulated cell-to-cell signaling, immune cell trafficking, inflammatory response, apoptosis and cellular migration ( S5 Table ). P 4 elicited down-regulation of cellular regeneration and proliferation and cell-cell signaling and adhesion. E 2 +P 4 upregulated cell death and molecular transport and downregulated cell growth and proliferation, carbohydrate metabolism and molecular transport. The genes commonly upregulated by E 2 , P 4 , E 2 +P 4 involved catalytic activity, receptor and signal transduction, binding, transporter and structural molecule activity. The main biofunctions of differentially expressed genes that were shared with differentially methylated loci involved cell membrane and signaling in response to E 2 . Hormonally-induced differentially methylated CpG sites were assessed for association with differential gene expression for each corresponding locus, noting that not all transcribed loci are included in both platforms and many intergenic regions in DNA methylation platform were not represented on the gene expression array used in this study. Only loci with a strong positive or negative association (by Spearman rho, and corrected p<0.05, see Methods ) were considered. There was a large number of functional gene clusters with strong association of DNA methylation and gene expression for E 2 in eSF normal ( Table 5 ), which was not observed for P 4 or E 2 +P 4 treatments. We next aimed to determine the effect of hormones on the endometrium of endometriosis patient, known to have abnormal P 4 response. We applied strict criteria using eSF from patients with only endometriosis and no other uterine, pelvic or gynecologic disorders and those that show P 4 -resistance confirmed by microscopy and IGFBP1 assay ( S1 Fig ). Furthermore, to understand the effect of disease stage on the hormone-epigenome interplay, we used early (stage I) and late stage (stage IV) disease. Similar to normal, E 2 induced the most and P 4 the least DNA methylation changes in eSF stage-I and eSF stage-IV ( Fig 2A and 2B , S6 and S7 Tables , respectively). But, in eSF endo , the extent of E 2 -induced changes was less than in eSF normal specifically in eSF stage-IV (418 CpG sites) exhibiting significantly less methylation alterations compared to eSF normal and eSF stage-I (2047 and 1633 CpG sites, respectively) ( Fig 2B ) . Opposite to that of eSF normal the majority of changes in both stages were gain of methylation ( Fig 2A heatmap). The extent of E 2 -induced differentially methylated loci differed considerably between the two stages ( Fig 2A and 2B ): stage I showed extensive changes induced by E 2 , much reduced in stage IV. The considerable difference in the extent of E 2 - induced methylation in eSF stage-IV and in the gain/loss pattern ( Fig 2A and 2B ), indicate an aberrant response to E 2 in both stages of disease and more extensively in stage IV, not previously reported. Progesterone, similar to eSF normal induced the least DNA methylome alterations in both eSF stage-I and eSF stage-IV despite the difference in the robust decidualization response to P 4 in eSF normal and the refractory decidualization response to P 4 in eSF endo ( S1 Fig ) . While P 4 induced similar numbers of loss and gain of methylation in eSF normal in eSF stage-I there was more loss than gain of methylation and in eSF stage-IV more gain than loss of methylation. Interestingly, despite complete lack of decidualization, eSF stage-I exhibited more E 2 +P 4 -induced differentially methylated loci versus eSF normal , and in both disease stages there was more loss than gain of methylation ( Fig 2A and 2B ). eSF stage-IV showed the fewest methylome changes in all three hormone treatments, suggesting extensive aberrancies in hormone-methylome interactions in late stage disease. Particularly important is the novel observation of an aberrant response to E 2 and not just to P 4 and E 2 +P 4, as previously believed [ 24 ]. 2A. Differentially methylated CpG sites induced by E 2 , P 4 and E 2 +P 4 versus vehicle. Heatmaps reflect the differential methylation of each sample in each hormone treatment versus its corresponding non-treated vehicle control (Δβ: Hormone treated minus vehicle control) (see Fig 1A legend for details). 2B. Number of differentially methylated CpG sites and in gain/loss of methylation for each hormone treatment in Endo I and Endo IV. C. Unique and common differentially methylated CpG sites for each hormone in Endo I (left) and Endo IV (right) indicating mostly unique loci for each hormone. 2 D . Unique and common differentially methylated CpG sites across normal (NUP), Endo I and Endo IV, for each hormone: E 2 : left, P 4 : middle and E 2 +P 4 : right. 2E. Enrichment of intergenic regions; Endo I, top and Endo IV bottom charts (see Fig 1D legend for details). 2F. Enrichment of enhancers for each hormone and based on loss or gain of methylation in Endo I (left panel) and Endo IV (right panel) (see Fig 1E legend for details). 2G. Genomic distribution of all differentially methylated CpG sites in each group (Endo I, left, Endo IV right panel) and by gain or loss of methylation (see Fig 1F legend for details). 2H. Distribution of differentially methylated CpGs by CpG islands (CGI), CGI north/south shores and shelves for Endo I (left) and Endo IV (right) based on gain or loss of methylation for all hormones. For details see Fig 1G legend. N Shelf: North Shelf; S Shelf: South Shelf; N Shore: North Shore; S Shore: South Shore. Endo I: stage I; Endo IV: stage IV. Differentially methylated loci were unique in response to different hormones in each disease stage and between the two stages. ( Fig 2C ). As in normal, E 2 +P 4 induced methylation were mostly unique and not a combination of the response to E 2 or P 4 individually ( Fig 2C ), reaffirming that E 2 and P 4 interact differently with the epigenome when combined than when individually administered ( Fig 2C , S6 and S7 Tables ). Moreover, the majority of loci differentially methylated in response to each specific hormone was also unique in eSF normal vs eSF stage-I vs eSF stage-IV ( Fig 2D ). These data suggest that the hormone-DNA methylome dynamics differ under normal and disease conditions, and furthermore that the stage of disease affect the hormone-methylome response. Despite distinct profile differences with normal, hormone-induced differentially methylated CpGs for E 2 , P 4 and E 2 +P 4 in both stages of eSF endo were also statistically significantly enriched in intergenic regions ( Fig 2E ) although the extent of this enrichment differed among hormone treatments, and between disease stages ( Fig 2E , S2 Fig for gain/loss for each hormone ; S2 Table ). Similar to normal and in both gain and loss of methylation, there was marked enrichment of enhancers, although the extent differed with specific hormone treatments and disease stage ( Fig 2F ). In both E 2 and P 4 treatment, eSF normal involved more enhancers than eSF stage-I and eSF stage-IV , but E 2 +P 4 treatment induced involvement of more enhancers in eSF stage-IV , particularly in loss of methylation with nearly 50% of CpGs associated with enhancers ( Fig 2F ). In eSF stage-I , the genomic distribution of loci with gain or loss of methylation differed from eSF normal at 1 st exons, gene bodies, and 3’UTRs. In eSF stage-IV the genomic distribution was mostly similar in gain and loss of methylation. These differed at 5’UTRs, gene bodies and 3’UTRs compared to eSF normal and at 5’UTRs and 1 st exons compared to eSF stage-I ( Fig 2G ). Overall, these data demonstrate that: hormone treatments regardless of disease and its stage affected CpG sites more at the 3’UTR and intergenic regions and much less at proximal promoters/TSS; genomic locations of CpG sites differentially methylated in response to hormones differed based on loss/gain of methylation; while decreased proximal promoter (TSS200) involvement and increased intergenic region involvement were common in eSF normal, eSF stage-I and eSF stage-IV . Low involvement of promoters/TSS and increased involvement of 3’UTR and intergenic regions were remarkable, considering vast differences in patterns, profiles and loci differentially methylated in eSF under the different hormonal treatments and in normal versus disease. These observations underscore key roles for genomic locations and potentially chromatin configurations further directing hormonal effects. CpG islands (CGI) , CGI shores and shelves . There was low involvement of CGIs and CGI shelves and shores in both disease stages Fig 2H (and S3 Table ), similar to normal. But, in both disease stages loss of methylation involved more CGIs than gain of methylation. Interestingly, hormone treatments significantly enriched more pathways in eSF stage-I and eSF stage-IV versus eSF normal ( Table 1 ) ( S1 Data shows important pathways affected by each hormone in each eSF group and marking differentially methylated genes in those pathways ) . Thus, while there were fewer loci in eSF stage-I and eSF stage-IV compared to eSF normal , more specific pathways were affected; whereas, with eSF normal hormone effects did not particularly affect specific canonical pathways and likely involved broader targets across the genome. E 2 affected pathways in eSF stage-I involved endometrial function/dysfunction and endometriosis (e.g. MAPK, PI3K-Akt, ErbB signaling, focal adhesion, gap junctions, among others ( Table 1 )). E 2 elicited pathways in eSF stage-IV associated with proteoglycans and estrogen, ErbB, Ras, GnRH, and FoxO signaling ( Table 1 )–all relevant to endometriosis pathophysiology [ 8 , 14 ]. While P 4 did not significantly enrich specific pathways in eSF normal, indicating a more genome-wide effect instead of limited effect at specific canonical pathways, several statistically significant pathways were enriched in eSF stage-I and even more in eSF stage-IV . These data suggest an aberrant response to P 4 in eSF from women with disease, which is enhanced in stage IV ( Table 1 ) involving specific pathways including estrogen, MAPK and ErbB signaling, confirming pathways associated with transcriptomic data [ 19 , 20 ]. Similarly, in response to E 2 +P 4 , there were more enriched pathways in eSF endo than in eSF normal relevant to endometrial function (adhesion, and disease/cancer ( Table 1 )). Gene functional enrichments in eSF normal , eSF stage-I and eSF stage-IV for each hormone treatment are shown in Table 2 . While there was little overlap in genes or functional clusters in eSF endo compared to eSF normal , the greatest number of genes in the same functional cluster in all eSF groups induced by E 2 involved those with signal peptide, membrane, and glycoprotein functions. eSF stage-I had more gene functional clusters with specific functions in all hormone treatments compared to eSF stage-IV or eSF normal . eSF stage-IV had the fewest functional clusters in all treatments compared to eSF stage-I and eSF normal , and the genes affected specific pathways involved in endometriosis and cancer as observed in the pathway analysis (above). Importantly, P 4 treatment affected specific gene functions in disease (adhesion, synapse, cell junction, cadherins), different from eSF normal . E 2 +P 4 elicited, only in eSF stage-I , several distinct functional clusters with specific functions in endometrial biology and endometriosis, including EGF/EGF-like genes, ECM receptor interaction, focal adhesion, PI3K-Akt pathway, synapse, cell junctions, spectrins and others. The most enriched cluster elicited by E 2 +P 4 in eSF stage-IV included fibronectins (large glycoproteins in ECM that bind integrins and other matrix components with major roles in cell adhesion, growth, migration differentiation, fibrosis and cancer). Functional enrichment differences did not show a gradual change from eSF normal to eSF stage-I and then to eSF stage-IV , rather showed distinct enrichments, suggesting inherent differences between disease stages. These data support that stage I and stage IV belong to distinct disease subtypes. As patterns, profiles, pathways and gene functions differed in responses of eSF stage-I, and eSF stage-IV to E 2 , P 4 and E 2 +P 4 in comparison to those of eSF normal, the question arose whether these could be due, in part, to aberrant DNA methylation signatures present prior to hormone treatment (referred to “pre-existing differences” herein). The DNA methylation status of untreated (vehicle) eSF endo from women with endometriosis were assessed for loci with aberrant methylation changes in response to each hormone, compared to eSF normal and whether they differed from untreated (vehicle) eSF normal (see Methods ). Numerous aberrantly differentially methylated loci in disease were found to be due to pre-existing DNA methylation differences across the genome ( S3 Fig ), including up to 53% of aberrant E 2 +P 4 loss of methylation in eSF stage-IV , showing an aberrant methylation pattern from that of the untreated normal eSF, further resulting in aberrant response to hormone treatments. These data are supported by previous gene expression analysis of eSF normal and eSF endo at t = 0 in culture [ 19 ], demonstrating intrinsic and pre-existing abnormalities in the eSF endo cells, although unclear whether these aberrancies are due to disease or are contributing to its progression/pathogenesis. As in eSF normal , transcriptomic profiles were determined for both stages of disease for each hormone (versus control). Whether these responses were abnormal was investigated compared to normal. In eSF endo , all treatments resulted in fewer differentially expressed genes versus eSF normal ( Table 3 , S8 and S9 Tables) . P 4 alone induced the fewest differentially expressed genes in eSF stage-I, and eSF stage-IV , consistent with the DNA methylation changes, an aberrant response to P 4 and abnormal decidualization in eSF from women with endometriosis ( Fig 2A ) [ 18 , 19 ]. In general, genes differentially expressed in response to each hormonal treatment were different in each disease stage versus normal, although some were in common ( Table 4 ). Of interest to endometrial function and in disease, PGR was also upregulated in response to E 2 in eSF stage-I, and eSF stage-IV similar to eSF normal ( Table 4 ) and was also among the top up-regulated genes in all groups despite the aberrant and very limited P 4 response observed in the DNA methylation, gene expression and IGFBP1 production in eSF from women with disease ( Table 3 , S1 Fig ). E 2 , P 4 and E 2 +P 4 up-regulated IGF1 and IL1R1 in all eSF and SPARCL1 was up-regulated in all E 2 +P 4 treated eSF ( Table 4 , S10 Table (full list and Venn diagrams of unique and common genes within each group and each hormone treatment across groups)). Gene expression profiles in response to hormones, similar to the DNA methylome, demonstrated distinct and aberrant molecular signatures in eSF stage-I versus eSF stage-IV and compared to eSF normal . Moreover, these were not limited to P 4 and E 2 +P 4 treatments and, importantly, included an abnormal response of eSF stage-I and eSF stage-IV to E 2. * in both common E2 and common P4 up-regulated genes ^ in both E2 and E2+P4 up-regulated; BOLD, in E2, P4 and E2+P4 common loci ** Up regulated in E2 and down regulated in E2+P4 *** Down-regulated in E2 and E+P. NOTE : There are no P4-induced downregulated loci in common in NUP , Endo I , Endo IV . In eSF stage I , similar to eSF normal , all hormone treatments resulted in more gene up-regulation than down-regulation, but unlike eSF normal , all hormones, including E 2 affected fewer differentially expressed genes, with a marked minimal effect with P 4 ( S8 Table , full gene list; S10 Table , common genes and Venn diagrams). Also, in E 2 and in E 2 +P 4 half and the majority of the genes, respectively, were unique, but in P 4 , the majority of differentially expressed genes were in common with E 2 +P 4 . Note that the number of P 4 -induced differentially expressed genes were very limited in stage I, while the combination of E 2 +P 4 in stage I disease resulted in more differentially expressed genes than with E 2 or P 4 treatments alone ( S8 Table ). Similar to eSF normal , where addition of E 2 minimally affected P 4 target genes, E 2 combined with P 4 affected the target genes of E 2 alone. While there were commonly upregulated P 4 target genes eSF stage I and eSF normal including IGF1 , GREB1 , and PGR , many key genes were missing in eSF stage I , such as, SPARCL1 which was upregulated in normal but not in stage I disease, further indicating an aberrant E 2 response in disease. Similar to eSF stage I the number of differentially expressed genes in response to E 2 as well as to P 4 treatments were far fewer than what was observed in eSF normal ( S9 Table , full gene list). Similar to eSF normal and eSF stage I there were more differentially expressed genes by E 2 +P 4 . Among commonly up-regulated E 2 +P 4 induced eSF stageIV and eSF normal were SPARCL1 , IGF1 , and LAMA3 . Among the 103 genes commonly down-regulated were CCL2 , RGS4 , RGS5 , IL-6 , MEST , KRT19 , KRT18 and H19 . The overlap in differential expression of specific up- and down-regulated genes with E 2 +P 4 treatment of eSF from women with and without endometriosis is remarkable, since stage IV disease eSF cells did not decidualize and are not considered to be P 4 -responsive [ 24 , 25 ]. Pathways and biofunctions, derived from the gene expression data, underscored distinct differences between eSF endo and eSF normal and between stages of disease, similar to the DNA methylation data. In eSF stage-I, with more limited E 2 effects, pathways included activation of cellular proliferation and viability ( S5 Table ), and eSF stage-IV involved increased tissue and cellular development (as with eSF normal ), proliferation, cell-cell signaling and adhesion (unlike eSF normal ). Note that E 2 -induced biofunctions and pathways were different in eSF stage-I and eSF stage-IV and both differed from eSF normal ( S5 Table ), consistent with the DNA methylation data. There were no enriched pathways in eSF stage-I and moderately enriched (Z score = 1.9) up-regulation of cell growth and proliferation in eSF stage-IV in reponse to P4 (far fewer loci). Similar to eSF normal , carbohydrate metabolism and molecular transport was also seen in eSF stage-I in response to E 2 +P 4 , which also showed up-regulation of cell invasion and viability. Importantly, E 2 +P 4 increased cell survival, cell movement and invasion, cell-to-cell signaling and adhesion, and downregulated cellular proliferation and growth in eSF stage-IV . Genes involved in these pathways and their upstream regulators are shown in S5 Table . Similar to normal, the main biofunctions of differentially expressed genes that were shared with differentially methylated loci involved cell membrane and signaling in response to E 2 in eSF stage-IV . In loci with a strong positive or negative association of DNA methylation and gene expression (by Spearman rho, and corrected p<0.05) distinct differences were found in eSF stage-I and eSF stage-IV and versus eSF normal (full lists, S11 – S13 Tables; S14 Table for unique and common loci between each group). Functional enrichment analyses revealed distinct differences in numbers and types of gene functional clusters in each stage of disease versus normal. While there was a large number of functional gene clusters with strong association for E 2 in eSF normal ( Table 5 ), eSF stage-I and eSF stage-IV showed different and more limited functional clusters. This result further indicates that the E 2 response is aberrant in eSF endo compared to normal eSF. There were also multiple differences in response to P 4 , and E 2 +P 4 among the eSF groups, further highlighting distinct molecular signatures in each disease stage. Importantly, eSF stage-I showed distinct clusters in response to P 4 and to E 2 +P 4, including functions characteristic of endometrial biology and endometriosis pathophysiology ( Table 5 ). While there were no statistically significantly enriched gene functional clusters in response to P4 in stage IV disease, a moderate enrichment of fibronectins, cell adhesion and secreted proteins were noted. These are consistent with the important role for cell adhesion in stage I and stage IV disease. In response to E 2 +P 4 , eSF stage-IV showed enrichment of calcium channels and integrins. Together these data suggest that the responses to all hormone treatments are aberrant in eSF derived from women with stage I and stage IV disease and are specific to each stage, supporting distinct disease subtypes. Herein, eSF hormone treatments in vitro were chosen to approximate the hormonal milieu in vivo (E 2 -dominant proliferative phase endometrium (PE) and E 2 +P 4 -dominant mid-secretory phase endometrium (MSE)). Comparing in vitro hormone eSF transcriptome data to corresponding phases in bulk endometrial tissue in normal versus disease [ 26 ] and FACS-isolated eSF endo and eSF normal [ 20 ] revealed great overlap of differentially expressed genes ( Table 6 , S15 and S16 Tables ). GO functional analysis of genes differentially expressed in E 2 treated eSF stage-IV and PE tissue revealed many genes in common involved in regulation of cell migration and motility, proteolysis, negative regulation of cell death, regulation of fibroblast proliferation and others. Regulation of inflammatory response, cell migration/motility, transport, protein import to nucleus, signal transduction, wound healing, and epithelial development and others were noted in stage I ( S15 Table ). Comparing MSE tissue and E 2 +P 4 -treated eSF stage-I and eSF stage-IV vs normal, pathway analysis revealed common signaling pathways involving PI3K-Akt, Rap1 and Ras, and cancer ( S15 Table ). Comparing transcriptomes of cultured eSF endo and eSF normal and freshly isolated (uncultured) FACS-sorted eSF endo and eSF normal from human eutopic endometrium [ 20 ] revealed many shared genes ( Table 6 ; S16 Table ). Note that FACS-sorted eSF normal and eSF endo included samples from various cycle phases, different disease stages [ 20 ], and a more limited sample number compared to the bulk tissue study. Thus, the number of overlapping loci in cultured and freshly isolated eSF is expectedly smaller than those shared with whole tissue. Together, the extent of overlap with whole tissue samples and FACS isolated eSF indicates the in vitro hormonal treatment of eSF, the predominant cell type in endometrium, is a good model and reflects a persistent eSF signature in the whole tissue. Since E 2 induced the largest changes in the DNA methylome of eSF, we sought to assess its effect on the histone marks to better understand how E 2 affected the regulatory function of the epigenetic machinery. We assessed silencing and activating histone modifications, H3K27me3, and H3K27ac, using chromatin immunoprecipitation followed by deep sequencing (ChIP-Seq). Modifications of H3K27me3, and H3K27ac have been found in loci involved in eSF decidualization [ 10 , 11 , 27 ]. In response to E 2 we observed more differential peaks in H3K27ac than H3K27me3 in line with our observation of more loss of DNA methylation corresponding to a more open chromatin state induced by E 2 ( S17 Table for peaks associated with each histone mark ). GO gene functional analysis for each histone modification renriched pathways related to regulation of signaling, cell morphogenesis and differentiation, G-protein coupled receptor signaling, regulation of mitotic cell cycle and intracellular protein transport among others, many of which are shared with DNA methylation data ( S17 Table for pathways for each histone mark ). Increased binding of PGR to open chromatin was shown previously in decidualizing cells by ChIP-Seq experiments [ 10 , 27 ] and that the presence of PGR binding site and its putative co-regulator FOSL2 in a genomic location is associated with open chromatin during decidualization [ 10 , 11 ]. Moreover, direct PGR targets in eSF treated for 72hrs with E 2 +MPA+cAMP were identified by Mazur et . al ., using ChIP-Seq and RNA-Seq [ 10 ]. We assessed the overlap of the E 2 +P 4 induced differentially methylated CpG sites associated with genes in normal, stage I and stage IV disease to genes with PGR binding sites present in the Mazur et . al . study within the extended promoter region (as defined to be -7500bp and +2500bp from TSS) and intervals within ±10KB, as well as ±25KB from transcriptional start/stop site in normal eSF. We found a small subset of genes overlapping in normal eSF ( Table 7 ); however, these common genes were enriched for biofunctions that are involved in cell morphogenesis, differentiation and cell projections, endosome organization and cytoskeletal organization ( Table 7 ). These are important during decidualization as eSF decidualization is characterized by morphological changes, expansion/restructuring of extracellular matrix, surface projections and expansion of endoplasmic reticulum. Interestingly, a larger number of genes in stage I and IV overlapped with PGR binding sites than normal eSF ( Table 7 ). Stage IV and normal shared more common genes with PGR binding sites than they did with stage I ( Table 7 ). Biofunction analysis showed more biofunctions involved in stage I than normal, such as tissue morphogenesis, response to TGF-beta signaling, response to growth factor and extracellular matrix among others ( Table 7 ). In stage IV, the biofunctions involved negative regulation of Wnt signaling, and intracellular signal transduction ( Table 7 ) known to be affected in endometriosis. These data further support the notion of aberrant P 4 response, rather than P 4 -resistance in endometriosis. * Common in NUP and Endo I ** Common in Endo I and Endo IV *** Common in NUP, Endo I and Endo IV; Bold Italics , common in NUP and Endo IV

The eSF is the most abundant steroid hormone-responsive cell in endometrium and is a master regulator of tissue function and pregnancy success, and thus how the steroid hormones E 2 and P 4 regulate the epigenome and transcriptional machinery in this cell type in a timeline similar to in vivo exposure is of high priority in understanding normal and abnormal eSF function in women. Effects of E 2 and P 4 alone on the hormone-epigenome interplay has largely been studied in breast cancer cell lines, providing key insights into hormone receptor topology, epigenetic genomic alterations, transcriptional regulation, and chromatin dynamics [ 28 , 29 ]. As these complex interactions are cell- and tissue-specific, extrapolating their properties to normal endometrium is limited, although a few studies have investigated the effects of E 2 plus progestins, such as medroxyprogesterone acetate, in the presence of cAMP on eSF for 72 hrs [ 10 ] or longer with or without glucose in the culture medium on chromatin accessibility or histone marks [ 9 , 11 ]. These studies show altered chromatin accessibility in eSF decidualized with E 2 +MPA+cAMP [ 11 ] and provide significant insights into PGR binding across the genome and association with open chromatin [ 10 ]. The current study investigated the effects of estradiol, progesterone (individually) and their combination (without cAMP or other progestins) on the DNA methylome and transcriptome and their interplay in normal eSF at 14 days, mimicking in vivo exposure times. Moreover, we compared eSF from the inflammatory disorder, endometriosis, in the setting of lesser and more advanced stage disease to the normal eSF. The data herein revealed, for the first time, that E 2 and P 4 individually and together promote unique patterns and profiles in the normal DNA methylome of this cell type. E 2 alone elicited broad changes, blunted by P 4 , and mostly result in open chromatin by inducing more loss of methylation and increased H3K27ac histone mark. Progesterone alone had a limited effect on the DNA methylome, and unlike E 2, elicited loss and gain of methylation equally. E 2 +P 4 affected the epigenome less robustly than E 2 alone, but showed more loss than gain of methylation. In support of our observation Vrljicak et . al . using transposase accessible chromatin followed by sequencing (ATAC-Seq) found altered chromatin accessibility with more open than closed chromatin loci after 4days treatment with MPA and cAMP [ 11 ]. These data suggest that E 2 and P 4 interact differently with the epigenome when combined than when individually administered suggestive of different mechanisms involved in the response of eSF to E 2 and to P 4 individually and in combination. Hormone-specific patterns and profiles were abnormal in both disease stages with more severe abnormalities associated with stage IV disease. The range of differences in individual loci with differential methylation and the number of enriched clusters and gene functions induced by each hormonal treatment in disease versus normal suggest inherent differences in disease and disease stages. However, in disease, as in normal, E 2 induced more extensive alteration than E 2 +P 4 followed by P 4 . Despite these differences, hormone-induced changes overall mainly involved CpG sites at the 3’end, intergenic regions, and enhancers, limited involvement of 5’end and 1 st exons and rarely involved CpG sites in close proximity to transcription start sites (TSS200) or CpG islands. Notably, E 2 treatment of MCF-7 breast cancer cells also demonstrated minimal binding of ER to proximal promoter regions (up to 5kb) [ 30 ], despite their containing the majority of known EREs. Whereas CpG sites in CpG islands (CGI) were minimally affected, CGI shores and shelves were more involved, regardless of methylation loss or gain or the type of hormone treatment, indicating a specific genome landscape interaction of hormones in this endometrial cell type. Whether the lack of involvement of CGIs reflects regulation of genes whose functions are not regulated by direct or indirect hormone-targeted mechanisms, or whether hormone response elements are not affected in CGIs is yet to be determined. In breast cancer cell lines gene expression [ 31 ] and DNA methylation profiles [ 32 ] as well as DNA methylation at several candidate genes at their CGIs [ 33 ] depend on their ER and PR status. This observation further highlights findings herein that the majority of differentially methylated loci in eSF normal are located in the intergenic regions, 3’UTRs and enhancers, and do not involve regions in close proximity to TSS, 5’UTR and 1 st exons, where most CpG islands are located. Women with endometriosis have high prevalence of infertility with otherwise unknown etiologies and lower implantation, clinical pregnancy and live birth rates compared to those without disease [ 34 ]. Studies in humans [ 35 ] and animal models [ 36 ] suggest compromised implantation attributed, in part, to an abnormal response to P 4 and the inflammatory milieu of the endometrium. The current study confirmed abnormal P 4 -regulated decidualization marker expression in eSF endo , largely attributed to “P 4 -resistance”, although P 4 did have effects across the eSF endo genome and PGR targets. However, eSF endo additionally had different responses to E 2 compared to eSF normal , which likely also contributes to abnormal endometrial function in women with disease (as described below). Of note, aberrant lack of ERα down-regulation at the time of implantation in endometrium of endometriosis women is considered key in implantation failure in women with disease [ 37 ]. However, endometrial-based infertility and effects on pregnancy outcomes are controversial, as large studies on IVF/ICSI outcomes in women with endometriosis and ovarian endometriomas revealed no differences in pregnancy rates [ 38 , 39 ] or a significant difference in endometrial receptivity array test in women with endometriosis versus controls [ 40 ]. How the aberrancies observed herein in P 4 -, or E 2 -induced epigenetic signatures are linked with implantation outcomes in women with endometriosis warrants further investigation. Whether stage I and stage IV endometriosis are distinct disease sub-types has been the subject of debate. That eSF stage-I differ greatly from eSF stage-IV in hormone response supports distinct disease subtypes. Also, disparities in the DNA methylomes of eSF stage-I and eSF stage-IV before hormonal treatments further support distinct disease subtypes. The latter observation underscores pre-existing abnormalities in the eSF epigenome in the setting of endometriosis, and the data showed eSF stage-IV with more extensive pre-existing differences affecting its responses, compared to eSF stage-I . This is consistent with previous findings that the endometrial bulk tissue transcriptome differs between the two stages [ 18 , 26 ] and several endometriosis genome-wide analysis studies suggestive of a stronger genetically driven component for stage IV than stage I disease [ 41 ]. Of note, clinically, women with stage IV versus stage I endometriosis have significantly lower implantation rates (13.7% vs. 28.3%, respectively), pregnancy rates (22.6% vs. 40.0%, respectively) [ 42 ], and lower IVF pregnancy rates (13.84% vs. 21.12% respectively) [ 43 ]- believed due to endometrial abnormalities that reflect distinct subtypes of disease. Mapping hormone-genome interactions of these subtypes holds promise for innovative, targeted therapies to modify pre-existing and stage-specific abnormalities in endometrium of women with endometriosis and optimize endometrial receptivity for implantation and pregnancy success of women with endometriosis. Dyson et . al . [ 44 ] have observed aberrant DNA methylation in the stromal fibroblast isolated from the endometriotic ectopic lesion. It remains to be determined whether and to what extent the ectopic lesion aberrances stem from the eutopic endometrial stromal fibroblasts. Recently, Maekawa et . al . assessed the genome-wide methylome changes during decidualization and in contrast to our data reported no changes in the DNA methylome [ 45 ]. DNA methylation distribution follows a bimodal distribution with the majority of CpG sites either hypomethylated or hypermethylated [ 46 ], as also reported in their study as well as in the current study and as we have previously observed in normal and endometriosis endometrium [ 14 , 15 ] Furthermore, we have also observed that the majority of CpG sites remain unchanged in decidualized versus non-decidualized eSF. The differences in our observation could be due to different analyses methods, where Maekawa et . al limited the definition of differential methylation to >Δβ of 0.3, which would not detect smaller changes. In our analyses we considered smaller changes in the DNA methylome but with the stringency that they were observed in at least 75% of each sample group. Another reason could be due to differences in the samples, where we used normal controls while patients with myoma or cervical cancer were used in their study, or it is likely that E 2 +MPA used in that study affects the epigenome differently than E 2 +P 4 in our study. Pursuing bulk tissue transcriptomic analysis, we first described “P 4 resistance” in endometrium from women with endometriosis [ 8 , 35 ] a phenomenon also observed by others [ 47 – 49 ]. In samples obtained in the implantation window and timed to the LH surge, there was evidence for impaired expression of key epithelial and stromal fibroblast markers of embryo receptivity and decidualization, respectively [ 35 ]. Analysis of endometrium across the menstrual cycle from women with severe disease strikingly revealed persistent E 2 -regulated genes in the early secretory phase, consistent with impaired P 4 action [ 8 ]. Moreover, these data were substantiated in a larger cohort [ 26 ], that also revealed a marked pro-inflammatory phenotype within the endometrium of women with disease. Inflammation can cause epigenetic changes in endometrium as demonstrated in an animal model of the disease [ 50 ]. We and others found P 4 -resistance in eSF [ 19 , 51 ]. Notably, inflammatory cytokines (e.g., IL-1β and TNFα) epigenetically silence the eSF PR, promoting P 4 -resistance with diminished expression of decidualization markers IGFBP1 and prolactin [ 50 , 52 ] and enhanced secretion of matrix metalloproteinases, which are normally suppressed in eSF by P 4 [ 53 ]. Epigenetic mechanisms underlying P 4 -resistance in endometriosis have mostly focused on the disease itself (as opposed to the eutopic endometrium studied herein) which exhibits P 4 and progestin resistance for pain relief initially or acquired over time [ 54 ]. We suggest that the nomenclature of “P 4 resistance” be re-evaluated, since the data herein show eSF endo respond to P 4 with regard to epigenetic marks, PGR target sites, and gene expression, albeit differently from eSF normal , although they do not fully decidualize. Notably, endometrium of women with disease does not retain a proliferative phenotype throughout the cycle [ 6 , 8 ], although there is compromised implantation in women with disease [ 35 , 52 ]. Importantly, “P 4 -resistance” was observed in endometrium of non-pregnant women who previously had severe pre-eclampsia, and this was also found in the decidua at delivery of women with this disorder [ 55 ], underscoring the need to understand this process for normal pregnancy. Thus, P 4 signaling in the endometrium and aberrancies in decidualization therein in vivo are likely influenced by other cell types in the tissue, including the inflammatory status of the individual and warrant further investigation. Herein, for the first time the observations have been made that in addition to aberrant eSF endo P 4 response, eSF from women with endometriosis show vastly aberrant responses to E 2. Specifically, E 2 -induced eSF DNA methylation changes blunted in stage IV disease and were more extensive in stage I. Loci with strong associations of DNA methylation and gene expression had distinct enrichment in gene functions in stage I and stage IV, including ion channels, ATP and nucleotide binding in stage I and plasma membrane and signaling in stage IV, suggesting functional impairment of eSF from women with versus without endometriosis. Moreover, they underscore that not only is the eSF response to P 4 abnormal in women with disease, but also their response to E 2 is abnormal . The latter has received little attention in the endometriosis literature, which is surprising, as the disorder is estrogen-dependent [ 12 , 21 ]. As eSF normally require E 2 priming prior to the full decidualization P 4 response, it is not unanticipated that with abnormal E 2 signaling in eSF, P 4 signaling would also be disrupted. Aromatase, essential for E 2 production, as well as E 2 levels are highly expressed in endometriotic tissue [ 56 – 59 ] with an increased COX2 expression in turn resulting in increased E 2 production in a positive feedback loop in ectopic and eutopic tissue of endometriosis patients [ 57 ]. Whether the aberrant response to E 2 observed herein could be affected by these aberrancies in E 2 production in endometriosis remains to be determined. The binding of some hormone NRs commonly occurs at accessible regions of the chromatin before hormone induction [ 60 ] or their recruitment occurs almost equally at the nucleosome-occupied and nucleosome-free states before hormone induction [ 61 ]. E 2 (biological active estrogen) enters cells, binds to subtypes of ERα and β that have high affinity for E 2 and are encoded by different genes [ 62 ]. While both ER subtypes are expressed in human endometrium, ERα is the primary mediator of E 2 action in this tissue [ 63 ]. ERα recruitment is complex involving multiple mechanisms depending on cell type and culture conditions [ 64 ]. ERα can bind to compact chromatin while there are abundant accessible regions before E 2 induction that will further recruit ERα [ 65 ]. The DNA methylation and histone modification findings herein suggest that E 2 can increase open chromatin. Chromatin accessibility can be induced by ERα binding, as these accessible clusters are found near estrogen-target genes [ 66 ]. While increased open chromatin was found in eSF normal in response to E 2 , the opposite was found in eSF endo . This could be due to pre-existing abnormalities in disease affecting chromatin structure, a combination of transcriptional machinery preloaded across the genome, or, as found in disease, up to 50% of loci displaying pre-existing differences in epigenetic signatures influencing this response. Furthermore, the state of chromatin compaction may play an important role. About half of EREs are in regions of DNA with open chromatin prior to estrogen induction [ 67 ], but many ERα binding sites in open chromatin are associated with differentially expressed genes after estrogen induction. These data indicate that chromatin compaction can directly affect ERα recruitment and subsequently target gene transcription. Thus, pre-existing and distinct differential methylation observed in stage I and stage IV can potentially affect chromatin compaction in patients with endometriosis. This is currently under investigation in our laboratory. ERα can be activated by phosphorylation by growth factors binding to tyrosine kinase receptors such as EGFR [ 68 ], which were dysregulated in the current study. Genes targeted by phosphorylated ERα are distinct from those targeted by estrogen-induced ERα activation [ 69 ]. Signaling through EGFR is a key pathway in eSF response to E 2 , and constitutive activation of EGFR in eSF from women with endometriosis has been reported [ 70 ]. Inhibition of EGFR in eSF from women with disease restores decidualization markers [ 71 ], underscoring the complexity of the interplay between E 2 and P 4 signaling in eSF in endometrium of women with endometriosis. The data overall support phosphorylation of ERα in eSF treated with E 2 may contribute, in part, to differential DNA methylation signatures and gene expression profiles observed in E 2 versus E 2 +P 4 in women without and with disease, which remains to be proven experimentally. Interestingly, E 2 and EGF can induce ERα recruitment at three classes of enhancers [ 72 ], bound only with EGF stimulation, only with E 2 stimulation, or either. Herein, enrichment of enhancer involvement upon E 2 stimulation and with E 2 +P 4 was observed in normal eSF as well as eSF from women with both stages of disease, but with different effects on downstream target genes. We propose that even small differences in EGFR signaling pathways could greatly alter the eSF responses to hormones, as observed herein. In this study, effects of E 2 , P 4 and their combination were elucidated on genome-wide DNA methylation marks of the endometrial stromal fibroblast, the predominant cell type in human endometrium essential for establishing and continuing pregnancy. The clinical phenotyping of truly normal controls and specific, well-phenotyped disease stages is a great strength of this study. Also, using the same cells for DNA methylation and gene expression analyses also added to the robustness of the data. Moreover, comparisons of the data herein with published gene expression and DNA methylation data in bulk tissue underscore signatures in the latter due to this predominant cell type. Single cell analysis of eSF from bulk tissue by FACS further underscores the signature of this cell type in overall bulk tissue analyses and opens the door for single cell RNA-Seq and DNA-me analyses in the future. While use of an in vitro system can address whether/how steroid hormones directly affect the eSF epigenome, an in vivo model using freshly isolated, sorted endometrial cells would offer an opportunity to assess functionality of the ER and PR landscape in human eSF and other endometrial cell types. Primary epithelial cells were not viable in culture when treated with hormones and as such we did not expand the current study involving epithelial cells, but organoid systems may offer a tool for this investigation. Further analysis using chromatin immunoprecipitation followed by deep sequencing of ER and expanding on the PR binding sites by Mazur et . al [ 10 ], identification of ERE and PRE-specific to endometrium and the status important to pioneer and co-activators and in a larger sample size are required for detailed mapping of the steroid hormone landscapes and hormone-epigenome interplay in normal human endometrial cells and in disease. Transcriptome data from the same cells demonstrate extensive overlap with previously identified differentially expressed genes in whole tissue in the corresponding hormone milieu. We note that utilizing microarray instead of a more comprehensive transcriptome analysis such as RNA-Seq limited the number and the type of transcripts investigated herein. Furthermore, protein data will enable full assessment of epigenomic and transcriptomic effects of hormones in endometrial function normally and in women with endometriosis. An important limitation of this study is the small sample size in each group. We had used very strict criteria for sample selection, both in identifying normal samples without any gynecological and pelvic disease/disorder and for endometriosis to not have any other disease no matter how benign, such as uterine fibroids. A follow up study with a much larger sample size is required to confirm our observations. Overall this study has elucidated the array of responses of eSF in health and disease in hormone milieu encountered in cycling women that can also serve for comparisons with actions of pharmaceutical steroids used clinically and potentially environmental estrogens that can compromise reproductive function. Moreover, the data reveal unique responses and pre-existing epigenetic abnormalities in women with endometriosis that can benefit endometrial-based diagnostic development and novel targeted therapies for endometrial dysfunction in women with this disorder.

This study was approved by the Committee on Human Research of the University of California, San Francisco (UCSF) (IRB# 10–02786). All samples were collected after written informed consent was obtained from all subjects. Eutopic (within the uterus) endometrial tissue samples were collected through the UCSF/NIH Human Endometrial Tissue/DNA Bank. Stringent inclusion criteria were applied as follows: I) for normal controls, samples were collected from oocyte donor volunteers with no uterine or pelvic pathology (NUP, normal controls); endometriosis samples were collected from stage I and stage IV endometriosis patients ( S18 Table ). Oocyte donor volunteers (controls) were extensively screened, had no gynecologic disorders, and donated endometrial samples six months post oocyte retrieval. Endometriosis patients (stage I and IV) were surgically confirmed and had no other gynecologic abnormalities. II) All samples were collected in the proliferative phase and matched for age, BMI, no smoking history (one exception), no contraceptive steroid use three months prior to sample collection, and endometrial stromal fibroblast (eSF) passage number. Menstrual cycle phase was determined by histological evaluation [ 73 ] as well as serum levels of E 2 and P 4 . Disease stage was determined by ASRM criteria [ 74 ]. Primary eSF were isolated from endometrial biopsies by digestion with collagenase and size fractionation as described [ 75 ] and cultured as monolayers in stromal cell medium (SCM, [ 18 , 19 , 76 ]). To ensure the purity of stromal cells in culture, after digestion of endometrial sample biopsies, the digested tissue was size fractionated using a 40μM filter to separate epithelial glands, followed by stromal cells selective attachment and growth in stromal cell medium. The purity of primary eSF was monitored morphologically during the culture and the homogeneity was verified by immunocytochemical localization of vimentin for eSF, keratin for epithelial cells and actin for vascular cells [ 76 ] before further hormone experiments. Only pure primary eSF (with <0.1% other cells) were used for this study. After 24 hours serum starvation, eSF from normal women (n = 7, controls) and endometriosis women (n = 6 stage I, n = 9 stage IV) were treated with four different hormonal treatments of 10 nM E 2 , 1μM P 4 , 10 nM E 2 +1μM P 4 , or vehicle (0.1% ethanol) control for 14 days [ 76 ]) after which conditioned media and cells were collected for further analysis. Decidualization was assessed in E 2 +P 4 treated eSF from normal, stage I and stage IV disease (see below). As eSF from women without endometriosis mostly have a robust decidualization response to E 2 +P 4 , our controls were eSF that fully decidualized (n = 4) by the decidualization marker IGFBP1 by ELISA and morphologically. As rarely do eSF from endometriosis women decidualize in vitro in response to E 2 +P 4 , eSF from stage I (n = 4) and stage IV (n = 4) with non-detectable decidualization (the most common phenotype) by morphology and IGFBP1 marker by ELISA were used for further analysis. Insulin-like growth factor binding protein-1 (IGFBP1), a P 4 -induced decidualization marker [ 77 ], was measured in media conditioned by 14 day E 2 +P 4 treated cultures, by ELISA (Alpha Diagnostic International Inc., San Antonio, TX) as a marker for decidualization. Concentrations were measured in duplicate, averaged and normalized to cell number. eSF were assessed by microscopy for morphological changes corresponding to decidualization. After treatments cells were harvested, pelleted and frozen at -80C for DNA and RNA extraction as previously described [ 14 , 15 ]. Genomic DNA was extracted using QIAGEN (QIAamp DNA Tissue Kit, QIAGEN, Germantown, MD) and RNA was extracted using the Macherey-Nagel NuceloSpin Tissue Kit with DNase treatment (Macherey-Nagel Inc., Bethlehem, PA) according to manufacturers’ recommendations and stored at -80C. Genomic DNA was bisulfite converted at the University of Southern California (USC) Epigenome Center using the EZ-96 DNA Methylation Kit (Zymo Research, Irvine, CA), according to the manufacturer’s protocol, and as previously described [ 14 , 15 ]. Quality, completeness of bisulfite conversion and amount of bisulfite-converted DNA were assessed by a panel of MethyLight reactions [ 78 ]. All samples passed all quality controls (QCs) and were further assayed by the Illumina Infinium HumanMethylation450K DNA methylation platform (HM450) based on Illumina’s specifications. (All data files are submitted to GEO, under SuperSeries accession number GSE145702 ). RNA quality was assessed by Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA). RNA samples were prepared for microarray analysis according to Affymetrix (Affymetrix, Inc., Santa Clara, CA) specifications [ 15 ]. cDNA sample quality was assessed by Bioanalyzer, and samples passing QCs were hybridized to Affymetrix HU133 Plus 2.0 gene array, interrogating >38,500 genes at the UCSF Genomics Core. HM450 interrogates 485,577 methylation targets across the genome. The ratio of methylated signal over total fluorescent signal was used to calculate β values, ranging from 0 (no) to 1 (complete) methylation. 850 control bead types were used as positive and negative controls and to calculate a detection P value to assess DNA methylation measurement quality for each probe of each sample [ 79 ]. The passing threshold of P value was set at P0.05 were indicated as “missing” (no statistically significant differences from background). Probe dropout-rates (percent probes with missing values versus total number of platform probes) were calculated to exclude samples with dropout-rates >1%. All 12 samples passed these criteria. Probes with a “missing” value in >1 sample were removed. Differential DNA methylation in each hormone treatment in each group (control, stage I, stage IV) was identified compared to non-treated cells (vehicle control) in the same group. For each probe, β values of hormone treatment in an individual sample were compared to its corresponding vehicle treated β value (Δβ). Changes in β values (Δβ) < 10% were not considered as differentially methylated for each sample. Median or average changes between hormone treatment vs vehicle control were not used, to obviate limited numbers of strong signals affecting selecting differentially methylated loci. Instead, probes were considered differentially methylated if they exhibited >10% change in β value (Δβ >10%) in hormone treatment versus vehicle and in at least 3 of 4 samples within each group, and with the same direction of methylation change (gain or loss). To assess if the aberrant signatures observed in hormone induced changes in disease compared to normal were due to pre-existing aberrations in the non-treated cells, we compared the methylation signatures of the non-treated (vehicle) eSF in each disease stage to non-treated (vehicle) eSF in normal for each aberrantly methylated locus in response to each hormone in disease. Pre-existing differential methylation, loci differentially methylated in normal vs each stage of disease in untreated cells (vehicle) were determined. The percentage of the differences observed in hormonal treatment of disease were assessed. Association of CpG islands with enhancers, distribution across the genome, association with CGIs and CGI shores and shelves were extracted from Illumina Infinium HumanMethylation450K manifest. The raw.CEL gene expression data files were RMA normalized using GeneSpring (GX13.1 version, Agilent Technologies). Loci were considered differentially expressed with Benjamini-Hochberg corrected ANOVA p<0.05 and fold change (FC)≥1.5 in each comparison. Differentially methylated loci and normalized gene expression were imported into R, and corresponding probes from each platform were matched using the transcript identifier. Every DNA methylation probe for a given locus identifier was compared to all corresponding transcripts of that locus using the non-parametric Spearman’s rank-order correlation method, as bivariate normality could not be assumed (DNA methylation data are not normally distributed). Spearman's rank correlation coefficients (ρ) on gene expression and DNA methylation were computed for each probe, along with a p value testing against the null hypothesis that ρ equals zero. Ingenuity Pathway Analysis (IPA, QIAGEN) software was used to determine pathways, upstream regulators and biofunctions of differentially expressed genes, as described [ 14 , 15 , 26 ]. Pathways with Z-scores ≥|2| were considered significantly enriched. For differentially methylated loci, DAVID and KEGG databases were used to identify functional classification, functional enrichment and pathways. For pathway selection an enrichment score ≥|2| with a Benjamini-Hochberg corrected p<0.05 was considered. For genomic distribution, element enrichment analyses, the null hypothesis that the observed proportions in two groups are the same, a test of proportions was performed in R using the prop.test() function. eSF treatment with E 2 and E 2 +P 4 used herein, mimicked in vivo proliferative endometrium (PE, E 2 -dominant phase) and mid-secretory endometrium (MSE (E 2 and P 4 - dominant phase). To assess commonalities, the current data were compared to previously published [ 26 ] whole endometrial gene expression in PE and MSE in normal versus disease. E 2 -treated differential expression in eSF from stage I versus control and stage IV versus control were compared to endometriosis (all stages) versus control in bulk tissue PE, and E 2 +P 4 - treated compared to MSE. FACS-sorted eSF gene expression data [ 20 ] in disease versus control were also compared. As FACS-sorted eSF endo and eSF normal included a mixture of phases and endometriosis stages, gene expression signatures of FACS-sorted eSF were compared to stage I, stage IV eSF treated with E 2 and E 2 +P 4 . We found that E 2 affected the methylome more robustly than P 4 or E 2 +P 4 and that, unexpectedly, along with aberrant P 4 response in disease, E 2 response was also aberrant in both stages of disease. Therefore, we sought to investigate further the effect of E 2 on two repressive and open chromatin histone marks, H3K27me3 and H3K27ac. eSF cells from two independent control participant were isolated from endometrial biopsies by digestion followed by size fractionation, and primary eSF were cultured in SCM and purity of cultured eSF was assessed as described above. eSF was passaged with trypsin and 1x10 5 cell/well were seeded. Confluent cultures were serum starved for 24hrs and treated with E 2 or vehicle for 14 days. Cells were cross-linked by a final concentration of 1% formaldehyde and terminated after 10 minutes by 0.125 M final concentration glycine. Chromatin was extracted using Chromatin Extraction Kit according to manufacturer’s recommendation (ab117152, Abcam, Cambridge, UK) sonicated by Diagenode Bioruptor and the size of sheared chromatin was visualized on agarose gel (100-600bp). Chromatin Immunoprecipitation was done using Abcam ChIP Kit (ab117138, Abcam, Cambridge, UK) with antibodies for H3K27me3 (ab6002, Abcam, Cambridge, UK) or H3K27ac (ab4729, Abcam, Cambridge, UK). Input control and immunoprecipitated DNA were paired-end sequenced using Illumina NextSeq 500 after library preparation according to the manufacturer’s instructions. Data were analyzed by removing adapter sequences, then aligned to reference human genome. Peaks called using Macs2 callpeaks and were selected with q-value 3.

(PDF) Click here for additional data file. (PDF) Click here for additional data file. (PDF) Click here for additional data file. (XLSX) Click here for additional data file. (XLSX) Click here for additional data file. (XLSX) Click here for additional data file. (XLSX) Click here for additional data file. (XLSX) Click here for additional data file. (XLSX) Click here for additional data file. (XLSX) Click here for additional data file. (XLSX) Click here for additional data file. (XLSX) Click here for additional data file. (XLSX) Click here for additional data file. (XLSX) Click here for additional data file. (XLSX) Click here for additional data file. (XLSX) Click here for additional data file. Common/unique genes are shown in columns and in Venn diagrams, in normal (NUP), stage I (Endo I), and stage IV (Endo IV). (XLSX) Click here for additional data file. (XLSX) Click here for additional data file. (XLSX) Click here for additional data file. (XLSX) Click here for additional data file. (XLSX) Click here for additional data file. (PDF) Click here for additional data file.