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

Summary

Programmed cellular responses to cycling ovarian-derived steroid hormones are central to normal endometrial function. Abnormalities therein, as in the estrogen-dependent, progesterone-“resistant” disorder, endometriosis, predispose to infertility and poor pregnancy outcomes. The endometrial stromal fibroblast (eSF) is a master regulator of pregnancy success. However, the complex hormone-epigenome-transcriptome interplay in eSF by each individual steroid hormone, estradiol (E₂) and/or progesterone (P₄), under physiologic and pathophysiologic conditions, is poorly understood and was investigated herein. Genome-wide analysis in normal, early and late stage eutopic eSF revealed: i) In contrast to P₄, E₂ extensively affected the eSF DNA methylome and transcriptome. Importantly, E₂ resulted in a more open versus closed chromatin, confirmed by histone modification analysis. Combined E₂ with P₄ affected a totally different landscape than E₂ or P₄ alone. ii) P₄ responses were aberrant in early and late stage endometriosis, and mapping differentially methylated CpG sites with progesterone receptor targets from the literature revealed different but not decreased P₄-targets, leading to question the P₄-“resistant” phenotype in endometriosis. Interestingly, an aberrant E₂-response was noted in eSF from endometriosis women; iii) Steroid hormones affected specific genomic contexts and locations, significantly enriching enhancers and intergenic regions and minimally involving proximal promoters and CpG islands, regardless of hormone type and eSF disease state. iv) In eSF from women with endometriosis, aberrant hormone-induced methylation signatures were mainly due to existing DNA methylation marks prior to hormone treatments and involved known endometriosis genes and pathways. v) Distinct DNA methylation and transcriptomic signatures revealed early and late stage endometriosis comprise unique disease subtypes. Taken together, the data herein, for the first time, provide significant insight into the hormone-epigenome-transcriptome interplay of each steroid hormone in normal eSF, and aberrant E₂ response, distinct disease subtypes, and pre-existing epigenetic aberrancies in the setting of endometriosis, provide mechanistic insights into how endometriosis affects endometrial function/dysfunction.

Keywords:

Chromatin – DNA methylation – Gene expression – Methylation – Genetic loci – Gene regulation – Steroid hormones – endometrium

Introduction

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₂) and progesterone (P₄) [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₂ and P₄ 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₂, cAMP and MPA have been shown [9–11], how E₂ and P₄ 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₂ for growth, disrupted P₄ 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₂-dominant (proliferative) versus P₄-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₂ and P₄ 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₄ 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₄ [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.

Results

Distinct Effects of Ovarian Steroid Hormones in Normal Human Endometrium

We assessed the genome-wide effect of individual steroid hormones, E₂ and P₄, and their combination (E₂+ P₄) 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₂ and P₄ 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₂ and P₄ 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₂ induced the most extensive changes in the eSF_normal DNA methylome (2047 CpG sites), followed by combined E₂+P₄ (569 CpG sites), and P₄ alone had the least effect (505 CpG sites) (Fig 1A and 1B, S1 Table). Importantly, combined E₂+P₄ resulted in dramatically reduced numbers of differentially methylated loci compared to E₂ 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₂, and E₂+P₄ induced more loss than gain of methylation (Fig 1A, yellow: gain of methylation, blue: loss of methylation vs vehicle), while P₄ 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₂+P₄ were mostly unique compared to those in response to P₄ or E₂ alone and were not a combination of the response to each hormone individually (Fig 1C).

Pathways and biofunctions (Table 1) as well as functional enrichment clustering (Table 2) were also unique to each hormone with E₂ 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₄ and E₂+P₄. Together these data indicate that each hormone affects different regions and E₂+P₄ targets are not a combination of E₂ and P₄. 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).

Genomic locations, regulatory elements, CpG islands and neighborhood context

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₂, P₄ and E₂+P₄) 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.

Changes in gene expression in response to hormones

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₂ induced more up- than down-regulated genes (Table 3, top up- and down-regulated loci; S4 Table, full gene list), and P₄ elicited similar numbers of up- and down-regulated genes, (S4 Table). However, more genes were differentially expressed when E₂ and P₄ were combined, with more genes up- than down-regulated (S4 Table).

E₂+P₄ induced the largest and P₄ the smallest changes in gene expression. Almost all the genes differentially induced by P₄ were shared with E₂+P₄ and some were shared with the E₂ treatment. Notably, half of the E₂-induced and the majority of E₂+P₄ induced differentially expressed genes were unique. However, in commonly up-regulated genes between E₂ and E₂+P₄, the variable fold changes indicate inhibitory or stimulatory effects when E₂ is combined with P₄. For example, PGR is upregulated by both E₂ and E₂+P₄ (FC = 4.5 vs 1.7, respectively), indicating that the addition of P₄ limited up-regulation of PGR compared to E₂ alone. Among other up-regulated gene in common in E₂, P₄, E₂+P₄, are IGF1 and SPARCL1 with known roles in endometrial biology. But, the FCs were different (IGF1: E₂ = 21.5, P₄ = 4.6, E₂+P₄ = 13.6; SPARCL1: E₂ = 4.6, P₄ = 13.2; E₂+P₄ = 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.

Pathways and biofunctions

E₂ 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₄ elicited down-regulation of cellular regeneration and proliferation and cell-cell signaling and adhesion. E₂+P₄ upregulated cell death and molecular transport and downregulated cell growth and proliferation, carbohydrate metabolism and molecular transport. The genes commonly upregulated by E₂, P₄, E₂+P₄ 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₂.

Association of gene expression with DNA methylation

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₂ in eSF_normal (Table 5), which was not observed for P₄ or E₂+P₄ treatments.

Effects of Ovarian Steroid Hormones in Endometrium of Women with Endometriosis

We next aimed to determine the effect of hormones on the endometrium of endometriosis patient, known to have abnormal P₄ 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₄-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₂ induced the most and P₄ 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₂-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₂-induced differentially methylated loci differed considerably between the two stages (Fig 2A and 2B): stage I showed extensive changes induced by E₂, much reduced in stage IV. The considerable difference in the extent of E₂- induced methylation in eSF_stage-IV and in the gain/loss pattern (Fig 2A and 2B), indicate an aberrant response to E₂ 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₄ in eSF_normal and the refractory decidualization response to P₄ in eSF_endo (S1 Fig). While P₄ 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₂+P₄-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₂ and not just to P₄ and E₂+P_4, as previously believed [24].

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₂+P₄ induced methylation were mostly unique and not a combination of the response to E₂ or P₄ individually (Fig 2C), reaffirming that E₂ and P₄ 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₂, P₄ and E₂+P₄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₂ and P₄ treatment, eSF_normal involved more enhancers than eSF_stage-I and eSF_stage-IV, but E₂+P₄ 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.

Pathways and biofunctions associated with differentially methylated loci

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₂ 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₂ 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₄ 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₄ 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₂+P₄, 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₂ 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₄ treatment affected specific gene functions in disease (adhesion, synapse, cell junction, cadherins), different from eSF_normal. E₂+P₄ 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₂+P₄ 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.

Aberrant hormone-induced methylation in eSF_endo are due to pre-existing methylation abnormalities

As patterns, profiles, pathways and gene functions differed in responses of eSF_stage-I, and eSF_stage-IV to E₂, P₄ and E₂+P₄ 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₂+P₄ 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.

Changes in gene expression in response to hormones

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₄ 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₄ 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₂ 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₄ response observed in the DNA methylation, gene expression and IGFBP1 production in eSF from women with disease (Table 3, S1 Fig). E₂, P₄ and E₂+P₄ up-regulated IGF1 and IL1R1 in all eSF and SPARCL1 was up-regulated in all E₂+P₄ 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₄ and E₂+P₄ treatments and, importantly, included an abnormal response of eSF_stage-I and eSF_stage-IV to E_2.

Changes in gene expression in stage I disease

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₂ affected fewer differentially expressed genes, with a marked minimal effect with P₄ (S8 Table, full gene list; S10 Table, common genes and Venn diagrams). Also, in E₂ and in E₂+P₄ half and the majority of the genes, respectively, were unique, but in P₄, the majority of differentially expressed genes were in common with E₂+P₄. Note that the number of P₄-induced differentially expressed genes were very limited in stage I, while the combination of E₂+P₄ in stage I disease resulted in more differentially expressed genes than with E₂ or P₄ treatments alone (S8 Table). Similar to eSF_normal, where addition of E₂ minimally affected P₄ target genes, E₂ combined with P₄ affected the target genes of E₂ alone. While there were commonly upregulated P₄ 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₂ response in disease.

Changes in gene expression in stage IV

Similar to eSF_{stage I} the number of differentially expressed genes in response to E₂ as well as to P₄ 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₂+P₄. Among commonly up-regulated E₂+P₄ 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₂+P₄ 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₄-responsive [24,25].

Pathways and biofunctions

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₂ 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₂-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₂+P₄, which also showed up-regulation of cell invasion and viability. Importantly, E₂+P₄ 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₂ in eSF_stage-IV.

Association of gene expression with DNA methylation

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₂ 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₂ response is aberrant in eSF_endo compared to normal eSF.

There were also multiple differences in response to P₄, and E₂+P₄ among the eSF groups, further highlighting distinct molecular signatures in each disease stage. Importantly, eSF_stage-I showed distinct clusters in response to P₄ and to E₂+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₂+P₄, 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.

Comparison of in vitro versus in vivo data

Herein, eSF hormone treatments in vitro were chosen to approximate the hormonal milieu in vivo (E₂-dominant proliferative phase endometrium (PE) and E₂+P₄-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₂ 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₂+P₄-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.

Histone H3K27me3 and H3K27ac modifications in response to E₂

Since E₂ 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₂ 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₂ 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₂ (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).

Association of PGR target loci identified in E₂+cAMP+MPA decidualization and E₂+P₄ induced DNA methylation

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₂+MPA+cAMP were identified by Mazur et. al., using ChIP-Seq and RNA-Seq [10]. We assessed the overlap of the E₂+P₄ 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₄ response, rather than P₄-resistance in endometriosis.

Discussion

Unique steroid hormone effects on normal endometrial stromal fibroblasts

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₂ and P₄ 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₂ and P₄ 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₂ 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₂+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₂ and P₄ individually and together promote unique patterns and profiles in the normal DNA methylome of this cell type. E₂ alone elicited broad changes, blunted by P₄, 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₂+P₄ affected the epigenome less robustly than E₂ 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₂ and P₄ interact differently with the epigenome when combined than when individually administered suggestive of different mechanisms involved in the response of eSF to E₂ and to P₄ 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₂ induced more extensive alteration than E₂+P₄ followed by P₄. 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₂ 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.

Epigenetic signatures in endometriosis: Hormone response, disease subtypes, pre-existing abnormalities

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₄ and the inflammatory milieu of the endometrium. The current study confirmed abnormal P₄-regulated decidualization marker expression in eSF_endo, largely attributed to “P₄-resistance”, although P₄ did have effects across the eSF_endo genome and PGR targets. However, eSF_endo additionally had different responses to E₂ 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₄-, or E₂-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₂+MPA used in that study affects the epigenome differently than E₂+P₄ in our study.

Progesterone “resistance”

Pursuing bulk tissue transcriptomic analysis, we first described “P₄ 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₂-regulated genes in the early secretory phase, consistent with impaired P₄ 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₄-resistance in eSF [19,51]. Notably, inflammatory cytokines (e.g., IL-1β and TNFα) epigenetically silence the eSF PR, promoting P₄-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₄ [53]. Epigenetic mechanisms underlying P₄-resistance in endometriosis have mostly focused on the disease itself (as opposed to the eutopic endometrium studied herein) which exhibits P₄ and progestin resistance for pain relief initially or acquired over time [54].

We suggest that the nomenclature of “P₄ resistance” be re-evaluated, since the data herein show eSF_endo respond to P₄ 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₄ -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₄ 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.

Abnormal response to E₂

Herein, for the first time the observations have been made that in addition to aberrant eSF_endo P₄ response, eSF from women with endometriosis show vastly aberrant responses to E_2. Specifically, E₂-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₄ abnormal in women with disease, but also their response to E₂ 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₂ priming prior to the full decidualization P₄ response, it is not unanticipated that with abnormal E₂ signaling in eSF, P₄ signaling would also be disrupted. Aromatase, essential for E₂ production, as well as E₂ levels are highly expressed in endometriotic tissue [56–59] with an increased COX2 expression in turn resulting in increased E₂ production in a positive feedback loop in ectopic and eutopic tissue of endometriosis patients [57]. Whether the aberrant response to E₂ observed herein could be affected by these aberrancies in E₂ production in endometriosis remains to be determined.

Potential mechanisms of E₂-epigenome interactions

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₂ (biological active estrogen) enters cells, binds to subtypes of ERα and β that have high affinity for E₂ and are encoded by different genes [62]. While both ER subtypes are expressed in human endometrium, ERα is the primary mediator of E₂ 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₂ induction that will further recruit ERα [65]. The DNA methylation and histone modification findings herein suggest that E₂ 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₂, 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₂, 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₂ and P₄ signaling in eSF in endometrium of women with endometriosis. The data overall support phosphorylation of ERα in eSF treated with E₂ may contribute, in part, to differential DNA methylation signatures and gene expression profiles observed in E₂ versus E₂+P₄ in women without and with disease, which remains to be proven experimentally. Interestingly, E₂ and EGF can induce ERα recruitment at three classes of enhancers [72], bound only with EGF stimulation, only with E₂ stimulation, or either. Herein, enrichment of enhancer involvement upon E₂ stimulation and with E₂+P₄ 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.

Study strengths and limitations and future directions

In this study, effects of E₂, P₄ 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.

Materials and methods

Ethics statement

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.

Samples

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₂ and P₄. Disease stage was determined by ASRM criteria [74].

Stromal cell isolation and hormone treatment

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₂, 1μM P₄, 10 nM E₂ +1μM P₄, 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₂+P₄ 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₂+P₄, 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₂+P₄, 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.

Decidualization assessment

Insulin-like growth factor binding protein-1 (IGFBP1), a P₄-induced decidualization marker [77], was measured in media conditioned by 14 day E₂+P₄ 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.

DNA and RNA extraction

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.

DNA methylation

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).

Gene expression microarray analysis

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.

DNA methylation data analysis

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 P<0.05, and probes with P>0.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.

Gene expression analysis

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.

DNA methylation association with changes in gene expression

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.

Pathways, biofunctions and genomic distribution analyses

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.

Association of in vitro and in vivo differential gene expression

eSF treatment with E₂ and E₂+P₄ used herein, mimicked in vivo proliferative endometrium (PE, E₂-dominant phase) and mid-secretory endometrium (MSE (E₂ and P₄- 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₂-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₂+P₄- 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₂ and E₂+P₄.

ChIP-Seq for Histone H3K27me3 and H3K27ac in response to E₂

We found that E₂ affected the methylome more robustly than P₄ or E₂+P₄ and that, unexpectedly, along with aberrant P₄ response in disease, E₂ response was also aberrant in both stages of disease. Therefore, we sought to investigate further the effect of E₂ 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⁵ cell/well were seeded. Confluent cultures were serum starved for 24hrs and treated with E₂ 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 <0.05. Differential peaks were identified using Macs2 bdgdiff and log likelihood ratio >3.

Supporting information

S1 Fig [pdf]
IGFBP1 ELISA assay of normal, stage I and IV eSFs used in the study.

S2 Fig [pdf]
Genomic distribution of differentially methylated CpG sites in each hormonal treatment, by gain or loss of methylation, in normal (NUP), stage I (Endo I) and stage IV (Endo IV) eSF.

S3 Fig [pdf]
Heat-map of pre-existing aberrancies prior to hormone treatments and percentage of contribution to each hormone induced methylation in stage I and stage IV eSF.

S1 Table [xlsx]
Differentially methylated CpG sites in response to E, P, E+P (vs. vehicle) in Normal eSF, based on loss and gain of methylation.

S2 Table [xlsx]
Distribution of differentially methylated loci (in %) based on their location across the genome and in comparison to those interrogated within the platform in normal (NUP), stage I (Endo I) and stage IV (Endo IV).

S3 Table [xlsx]
Distribution of differentially methylated loci (in %) based on their association with CpG islands (CGI), and CGI north/south shelf, CGI north/south shore and not associated with CGIs (open sea), and in comparison to the interrogated CGI in the platform in normal (NUP), stage I (Endo I) and stage IV (Endo IV).

S4 Table [xlsx]
Differentially expressed genes in hormone treatments (E, P, E+P) vs. vehicle in normal eSF (NUP).

S5 Table [xlsx]
Enriched pathways and the associated genes and upstream regulators of differentially expressed gene in response to hormones (E, P, E+P) vs. vehicle, in normal (NUP), stage I (Endo I), and stage IV (Endo IV) disease.

S6 Table [xlsx]
Differentially methylated CpG sites in response to E, P, E+P (vs. vehicle) in endometriosis stage I eSF, based on loss and gain of methylation.

S7 Table [xlsx]
Differentially methylated CpG sites in response to E, P, E+P (vs. vehicle) in endometriosis stage IV eSF, based on loss and gain of methylation.

S8 Table [xlsx]
Differentially expressed genes in hormone treatments (E, P, E+P) vs. vehicle in stage I disease (Endo I).

S9 Table [xlsx]
Differentially expressed genes in hormone treatments (E, P, E+P) vs. vehicle in eSF from stage IV disease (Endo IV).

S10 Table [xlsx]
Unique and common up- and down-regulated genes in each hormone treatment and in normal (NUP), stage I (Endo I) and stage IV Endo (IV).

S11 Table [xlsx]
Association of differentially methylated (DM) loci with differentially expressed (DE) genes in normal eSF (NUP) for all hormone treatments (E, P, E+P) vs. vehicle, based on positive or negative association (positive/negative rho) and gain and loss of methylation (gain, loss).

S12 Table [xlsx]
Association of differentially methylated (DM) loci with differentially expressed (DE) genes in stage I eSF (Endo I) for all hormone treatments (E, P, E+P) vs. vehicle, based on positive or negative association (positive/negative rho) and gain and loss of methylation (gain/loss).

S13 Table [xlsx]
Association of differentially methylated (DM) loci with differentially expressed (DE) genes in stage IV eSF (Endo IV) for all hormone treatments (E, P, E+P) vs. vehicle, based on positive or negative association (positive/negative rho) and gain and loss of methylation (gain/loss).

S14 Table [nup]
Differentially methylated (DM) loci associated with differentially expressed (DE) assessed by Spearman rho, with either positive association (pos) or negative association (neg).

S15 Table [xlsx]
Overlap of with genes and biofunctions from hormonal treatments in culture and endometrial tissue cycle phases in endometriosis stage I or stage IV versus normal.

S16 Table [xlsx]
Genes commonly differentially expressed in eSF from disease vs. normal, treated or untreated in vitro with hormones and in eSF FACS-isolated in disease vs. normal.

S17 Table [xlsx]
Peaks and GO biofunctions enriched in histone marks H3K27me3 and H3K27ac in normal eSF induced by E versus vehicle.

S18 Table [xlsx]
Sample information.

S1 Data [pdf]
Differentially methylated loci and the associated pathways/biofunctions that are affected in hormonal treatments of eSF eSF and eSF with known roles/importance in normal endometrial function and dysfunction in endometriosis.