Histone H3K4me3 breadth in hypoxia reveals endometrial core functions and stress adaptation linked to endometriosis.

Trimethylation of histone H3 at lysine 4 (H3K4me3) is a marker of active promoters. Broad H3K4me3 promoter domains have been associated with cell type identity, but H3K4me3 dynamics upon cellular stress have not been well characterized. We assessed this by exposing endometrial stromal cells to hypoxia, which is a major cellular stress condition. We observed that hypoxia modifies the existing H3K4me3 marks and that promoter H3K4me3 breadth rather than height correlates with transcription. Broad H3K4me3 domains mark genes for endometrial core functions and are maintained or selectively extended upon hypoxia. Hypoxic extension of H3K4me3 breadth associates with stress adaptation genes relevant for the survival of endometrial cells including transcription factor KLF4, for which we found increased protein expression in the stroma of endometriosis lesions. These results substantiate the view on broad H3K4me3 as a marker of cell identity genes and reveal participation of H3K4me3 extension in cellular stress adaptation.

Introduction

Promoters are key components in gene regulation, and trimethylation of histone H3 at lysine 4 (H3K4me3) is a known marker of active promoters (Barski et al., 2007; Chen et al., 2015; Karlić et al., 2010; Park et al., 2020). Broad H3K4me3 promoter domains have been associated with cell identity and transcriptional consistency (Benayoun et al., 2014) and are thought to guard cell type-specific functions. However, the dynamics of broad H3K4me3 promoter domains under stress conditions have not been studied in mammals, and it is unclear if broad domains that are postulated to buffer environmental fluctuations are maintained or modified in stress. Here we study H3K4me3 dynamics in hypoxic stress, which is a major cellular stress condition that induces profound cellular adaptation mechanisms (Pugh and Ratcliffe, 2017; Wilson et al., 2020). Hypoxia-induced functionally relevant chromatin changes have been reported to be caused by modification of normoxic H3K4me3 promoter marks (Prickaerts et al., 2016). These occur at least partly via lysine demethylases 5 (KDM5) enzymes that demethylate H3K4me3 in an oxygen dependent manner (Batie et al., 2019) and are also known to participate in protective hypoxic or ischemic preconditioning (Pilz et al., 2019).

Hypoxic conditions are especially relevant in the uterus where endometrial stromal cells are exposed to hypoxia upon menstruation and during placentation (Martínez-Aguilar et al., 2021; Soares et al., 2017). Furthermore, repeated menstrual cycles are viewed as a stress preconditioning process (Brosens et al., 2009) that has parallels to hypoxic stress preconditioning observed in other tissues (Gidlöf et al., 2016; Pilz et al., 2019). Thus, endometrial stromal cells can be used to investigate epigenetic mechanisms involved in cycling stress conditions and cellular adaptation to stress, such as hypoxia. In addition, hypoxia is involved in the etiology of endometrial disorders such as endometriosis (Li et al., 2021).

To determine the dynamics of H3K4me3 promoter marks under hypoxic stress conditions, we studied hypoxic H3K4me3 marks in cultured endometrial stromal fibroblasts (ESF) and differentiated decidual stromal cells (DSC). We conducted H3K4me3 ChIP-seq and investigated the breadth and height of the H3K4me3 promoter marks and correlated these to transcription by assessing both absolute levels and fold-changes. In addition, we tested if the discovered subsets of H3K4me3 promoters were enriched within known endometrial pathogenesis genes and validated these findings using immunohistochemistry in endometrial biopsies.

Results and discussion

H3K4me3 demethylases are hypoxia regulated in endometrial stromal cells

To evaluate the contribution of the main histone methylation modifying enzymes in hypoxia responses, we inspected their mRNA expression patterns in hypoxia-treated endometrial stromal cells (Figures 1A and S1). We found that several histone lysine demethylases KDMs and histone lysine methyltransferases KMTs (that add methyl groups) targeting histone H3K4 had high transcription levels (>125 TPM in normoxia) and were downregulated in hypoxia. Notably, H3K4me3 demethylases (KDM5) were the most abundantly transcribed demethylases, and KDM5A displayed the most notable hypoxic downregulation (Figure 1A). Specifically, KDM5 enzymes demethylate H3K4me3 in an oxygen dependent manner, which suggests that chromatin H3K4me3 levels would constitute a switch that can be readily regulated in hypoxia. Our findings are consistent with those implicating a role for H3K4me3 demethylases in hypoxia regulation (Batie et al., 2019), and suggest that H3K4me3 modification is important for hypoxic stress responses in endometrial stromal cells.

Figure 1

Broad H3K4me3 domains correlate with high transcription of core endometrial genes that are maintained in hypoxia stress

(A) Hypoxia induced changes in the transcription of histone lysine demethylases (KDM) (transcripts per million, TPM).

(B) Promoter H3K4me3 signal (100 bp bins) for endometrial stromal fibroblasts (ESF) in normoxia, ESF in hypoxia (ESFhyp) (1% O2, 16 h), and decidual stromal cells in hypoxia (DSChyp). H3K4me3 signal density displays the number of ChIP-seq reads per million mapped reads.

(C) Spearman correlations of H3K4me3 peak breadth or height to the gene transcription in the studied conditions. X axis displays included gene promoters with increasing cut-offs on RNA abundance (>TPM). Arrow marks the TPM >2 cut-off. See Figure S3 for dot plots.

(D) Gene Set Enrichment Analysis (GSEA) of transcriptionally hypoxia-upregulated or hypoxia-downregulated genes (FC > 2, FDR <0.01, >2 TPM) on the H3K4me3 promoters ranked with breadth (broadest on the left). The GSEA enrichment score describes the overrepresentation of a gene set in the ranked list.

(E) Proportions of shared top 500 broadest and highest H3K4me3 promoter peak signals for the three conditions.

(F) Hierarchical clustering of the p values of the enriched of GO terms of the gene promoters with the 500 broadest and highest H3K4me3 marks. The color scale depicts the enrichment with -log10(P).

(G) GSEA using siPGR-downregulated genes (Mazur et al., 2015) (FDR < 0.01 and > 2 TPM) on the H3K4me3 promoters (ESFhyp) ranked by breadth.

(H) Hypoxia induced transcriptomic changes in a subset of siPGR-downregulated (Mazur et al., 2015) and top H3K4me3 broad genes that are shared in all the three conditions.

See also Figures S1–S4 and Tables S1, S2, and S3.

Breadth of H3K4me3 promoter domains correlates with transcription

To characterize H3K4me3 dynamics in a major cellular stress condition, we conducted H3K4me3 ChIP-seq in hypoxia-treated (16 h, 1% O2) undifferentiated endometrial stromal fibroblasts (ESF) and in differentiated decidual stromal cells (DSC) and compared the MACS2 called H3K4me3 peaks at gene promoters to those in normoxia. We observed substantial variation in H3K4me3 peak read intensity and breadth (Figure 1B and Table S1), but we did not find any significant differences in the numbers of marked promoters between hypoxia and normoxia (Figure S2A) and neither found strongly enriched functional categories among the H3K4me3 marked (TPM > 2) gene promoters specific to each condition (Figure S2B). This suggests that relevant histone changes in hypoxia are caused by modification of existing normoxic promoter marks, which also has been observed in other cell types, such as in MCF7 breast cancer cells (Prickaerts et al., 2016).

The absence of a clear functional subset of genes with gained and lost marks motivated us to examine H3K4me3 dynamics by correlating H3K4me3 signal height and breadth with transcription (Figures 1C and S3). In all studied conditions, both the H3K4me3 peak height and breadth correlated with the genes being turned on or off (TPM > 2), but breadth correlated better with the intensity of the transcription (Figure 1C). Our results substantiate previous observation from CD4+ T cells (Chen et al., 2015) that H3K4me3 peak height at the promoter partially reflects the poised or preprogrammed state rather than active transcription, whereas H3K4me3 peak breadth better correlates with transcriptional levels.

Broad H3K4me3 promoters are maintained in hypoxia and have functional enrichment for genes relevant for endometrium

We observed that the broadest H3K4me3 marked promoters were enriched for transcriptionally hypoxia-upregulated but not for hypoxia-downregulated genes (GSEA p < 0.001 and p = 0.27, respectively) (Figure 1D and Table S2). To further inspect the H3K4me3 dynamics we inspected the top 500 ranks of the H3K4me3 peak height or breadth. In hypoxic conditions, the 500 promoters with broadest H3K4me3 peaks retained a substantially higher transcriptional output than all H3K4me3 marked and transcribed (TPM > 2) promoters, whereas the level of transcription of the 500 promoters with the highest H3K4me3 peaks did not differ from the rest of the H3K4me3 marked and transcribed promoters (Figure S4). Furthermore, the same 500 broadest H3K4me3 promoters in normoxia were maintained in hypoxia, whereas the 500 highest ones were not maintained (Figure 1E). These results suggest that the broad H3K4me3 promoters constitute a functionally important and maintained subset in stress. This is further supported by gene ontology enrichment analysis where these genes were highly enriched in functional categories relevant for endometrial physiology (p < 10−15), such as blood vessels and urogenital development (Figure 1F and Table S3). In contrast, the gene promoters with the 500 tallest peaks displayed a modest enrichment (10–15 < p < 10–5) mainly in housekeeping functions, especially mRNA processing and metabolism (Figure 1F and Table S3).

Broad H3K4me3 promoters are enriched with core endometrial regulators that are transcriptionally induced by hypoxia

Progesterone receptor targets constitute core gene regulatory networks that allow the decidualization of the stromal cells from undifferentiated EFSs (Mazur et al., 2015). Both upregulated and downregulated progesterone receptor targets that were previously detected by siPGR (Mazur et al., 2015) were enriched in broadest H3K4me3 promoters (Figure 1G). We intersected the genes maintaining the broad H3K4me3 promoters in hypoxia with genes upregulated by progesterone (siPGR down). We identified genes such as homeobox A10 (HOXA10) and heart and neural crest derivatives expressed 2 (HAND2) that are transcriptionally upregulated in hypoxia (Figure 1H) and are core components in progesterone induced decidualization (Gellersen and Brosens, 2014). This suggests that hypoxic transcriptional upregulation associates with promoters present with broad H3K4me3 peaks, enforcing the maintenance of the cell type directing functions such as the expression of core decidualization genes. This also supports the earlier notion that a broad H3K4me3 promoter marks guard cell identity in stress conditions (Benayoun et al., 2014).

Decrease in H3K4me3 signal associates with general hypoxic repression and extension of H3K4me3 breadth with targeted hypoxic upregulation

We then focused on the fold-changes of the H3K4me3 intensity or breadth. When inspecting the promoters of all the expressed (TPM > 2) genes, we observed lower average H3K4me3 levels in hypoxia than normoxia (Figure 2A). We investigated the fold-change in the peak intensities using DiffBind (with MACS2 peaks) and noticed that hypoxia led to a decrease in the intensity of H3K4me3 peak signals at 1700 expressed promoters (DiffBind, q < 10−10) (Table S4). The genes of these promoters were enriched with functions related to RNA metabolism (R-HSA-8953854, p = 1.1 × 10−89) that we also detected for the genes with the highest peaks (Figure 1F). Downregulation of RNA metabolism is part of general metabolic repression, which is one of the most defining features of hypoxia responses across all cell types (Semenza, 2020). When we contrasted these hypoxic induced intensity changes to those that affect the H3K4me3 domain breadth using the Chromtime tool (detecting changes in peak breadth), we detected 112 expressed promoters (q < 0.05) with hypoxia-extended (increased breadth) H3K4me3 peaks in ESF. These peaks displayed higher H3K4me3 signals downstream of TSS in hypoxia compared to normoxia (Figure 2B). Majority of these detected hypoxia-extended H3K4me3 extend only downstream the ORF from TSS, with some examples extending to both directions or only upstream from TSS (Figure S5).

Figure 2

Hypoxic extension of H3K4me3 domains reveals endometrial hypoxia adaptation genes

(A and B) H3K4me3 promoter signal of (A) all expressed (>2 TPM) and (B) hypoxia H3K4me3-extended genes. Red, normoxic ESF; dark blue, hypoxic ESF; and light blue, hypoxic DSC. The promoter signal is diplayed as mean log2FC of H3K4me3 vs input.

(C) Relationship of promoter H3K4me3 signal fold change (DiffBind log2FC) and transcription fold change (TPM log2FC) in hypoxia treated versus normoxic ESF. The ChromTime (breadth) subset is marked with red.

(D) Hierarchical clustering of the p values of the functional enrichment terms of the genes with hypoxia increased H3K4me3 signal intensity (DiffBind) promoters and genes with hypoxia-extended promoters (ChromTime). The functional categories are from GO Biological Processes (GO), Hallmark gene sets (H), canonical pathways (CP), and Reactome (RE) accessed via metascape.org. For gene lists see Table S4.

(E) Examples of hypoxia extended H3K4me3 promoter genes for the “Hypoxia” and “Epithelial mesenchymal transition” terms: Vascular endothelial growth factor (VEGFA), Phosphoglycerate Kinase 1 (PGK1), and Basic Helix-Loop-Helix Family Member E40 (BHLH40).

See also Figures S5 and S6 and Table S4.

Hypoxic extension of H3K4me3 domain breadth reveals hypoxia adaptation genes relevant for endometrium

To investigate the functional consequences of the changes in H3K4me3 peak height and breadth, we correlated these variables with changes in transcription. The intensity changes had weak positive correlation with transcriptional changes (Pearson r = 0.25) (Figure 2C), whereas breadth changes displayed a marked positive correlation (Pearson r = 0.49). When intensity change correlation was repeated using only the peaks detected by Chromtime (breadth), a notable positive correlation was detected (Pearson r = 0.47, Figure 2C red line, see Figure S6 for extension fold change vs transcription fold-change). Thus, upon stress, especially changes in H3K4me3 breadth robustly predict functional changes. Concordantly, the detected 112 promoters with hypoxia-extended H3K4me3 peaks had robust enrichment for hypoxia related gene sets (Hallmark Hypoxia p = 1.3 × 10−11, Epithelial mesenchymal transition (EMT) p = 8.5 × 10−13) (Figure 2D). These included classical hypoxia response genes such as vascular endothelial growth factor (VEGF) that is key driver of angiogenesis, glycolytic enzyme phosphoglycerate kinase 1 (PGK1), and BHLHE40 that regulates endometrial EMT (Asanoma et al., 2019) (Figure 2E). Thus, hypoxic H3K4me3 extension is involved in angiogenesis that is crucial for placentation and for endometrial reconstruction after menstruation (Martínez-Aguilar et al., 2021; Soares et al., 2017), ensuring the energy demands of endometrial growth (Kommagani et al., 2013; Zuo et al., 2015) and the regulation of the epithelial/mesenchymal status that is altered during decidualization (Gellersen and Brosens, 2014).

Hypoxic extension of H3K4me3 domains reveals genes associated with endometriosis

To compare our results on the H3K4me3 status with clinical datasets, we intersected our H3K4me3 data and transcriptomic data of the 500 broadest genes of the hypoxic ESFs with data from endometriosis and endometrial cancers (see STAR Methods for details) using Fisher’s exact test (Figure 3A). Upregulated (DiffBind) and extended (Chromtime) H3K4me3 promoter sets were intersected with RNA-seq upregulated genes whereas downregulated (DiffBind) and compressed (Chromtime) H3K4me3 promoter sets were intersected with RNA-seq downregulated genes. Genes that showed hypoxia-induced H3K4me3 peak extension and those found to be upregulated in stromal cells isolated from endometriosis lesions (Rekker et al., 2017) displayed the highest significant overlap (p = 3.4 × 10−7, 10.4-fold), followed by broad H3K4me3 promoter mark and endometriosis-upregulated stromal genes (p = 3.0 × 10−6, 4.5-fold) (Figure 3A). Several of the genes found with both H3K4me3 extension and upregulation in endometriosis are involved in cell survival (KLF2, NUAK1, NDRG1, and SGK1) and angiogenesis or fibrosis (RARRES2, SERPINE2, S100A10, TNFSF4, and ACTA2) (Table S5). For instance, KLF2 is known to protect cells from stress in multiple conditions (Fang et al., 2017; Gao et al., 2019), RARRES2 (chemerin) is a marker of endometriosis in the peritoneal fluid (Meng et al., 2019), and ACTA2 is a marker for fibrotic state, which is a general hallmark of endometriosis lesions (Vigano et al., 2018). Overall, these results suggest that hypoxic H3K4me3 extension at gene promoters associates with gene pathways that are relevant in the progression of endometriosis.

Figure 3

Hypoxic extension of H3K4me3 domains reveals genes that are upregulated in endometriosis

(A) A heatmap of Fisher’s exact tests for overrepresentation across the overlaps of H3K4Me3 promoter and transcriptomic states (ESF hypoxia) with representative endometrial disease datasets. The color scale depicts overrepresentation with -log10(p-value). The H3K4me3 data (X axis) was intersected with corresponding hypoxia RNA-seq-upregulated or RNA-seq-downregulated genes before the tests, and the 500 broadest ESF H3K4me3 in hypoxia (Broad) were intersected with both regulatory directions of RNA-seq. The disease data (y axis) include stroma specific endometriosis upregulated genes (Rekker et al., 2017), a cell type heterogeneous endometrium data from endometriosis patients versus healthy controls (Tamaresis et al., 2014) and two Cancer Genome Atlas Uterine Corpus Endometrial Carcinoma (TCGA-UCEC) sets, mutation, and copy number altered (CNA). Blue and green boxes, most significant overlaps.

(B) Hypoxia induced transcriptomic changes in ESF hypoxia H3K4me3-extended genes (∗) or top H3K4me3 broad genes that also are transcriptionally upregulated in endometriotic stromal cells (Rekker et al., 2017).

(C) Genomic view of H3K4me3 signal at Kruppel-like Factor 2 (KLF2) locus.

(D) Immunohistochemical staining of KLF2 in three types of endometriosis lesions (deep infiltrating sacrouteral, deep infiltrating intestinal and ovarian) and healthy control endometrium. E = epithelium, S = stroma.

(E) Proportion of positive nuclear stromal staining of KLF2 in three types of endometriosis lesions (IN = intestinal (n = 3), OV = ovarian (n = 3), SA = sacrouteral (n = 2)) and healthy control endometrium (n = 3). P value from Welch Two Sample t-test.

See also Table S5.

Finally, we assessed the genes with hypoxia extended and broad H3K4me3 promoters that were transcriptionally upregulated in hypoxia and endometriosis (Figure 3B). Of these, the transcription factor Kruppel like factor 2 (KLF2) displayed the most notable hypoxic transcriptomic induction (17-fold in ESF hypoxia) (Figure 3B) with hypoxia extended H3K4me3 promoter mark (Figure 3C). For translational insight, we stained immunohistochemically KLF2 in endometriosis lesions and endometrial biopsies from healthy controls. In three types of endometriosis lesions (deep infiltrating sacrouteral, deep infiltrating intestinal and ovarian) the stromal staining of KLF2 was stronger compared to the control endometrium (combined P = 0.019) (Figures 3D and 3E). The strongest stromal staining was observed for deep infiltrating sacrouteral lesions, and we did not observe significant differences in epithelial cells of the various sample types. Together, our results indicate that hypoxic extension of H3K4me3 at the KLF2 promoter identified in vitro is associated with upregulation of KLF2 transcription and protein level in endometriosis in vivo lesions. Notably, in addition to its known stress protective roles (Fang et al., 2017; Gao et al., 2019), KLF2 is a repressor of T cells/monocyte activation and a repressor of the NFKB pathway (Jha and Das, 2017), suggesting that its upregulation may contribute to the ability of endometriotic stromal cells to escape clearance by the immune system.

Conclusions

We investigated H3K4me3 dynamics in normoxia and hypoxia stress in cultured endometrial stromal cells and observed that promoter H3K4me3 domain breadth, rather than H3K4me3 height, correlates with transcription. Our results show that broad promoter H3K4me3 domains are maintained in hypoxia and mark core cell type-specific regulators in endometrial stromal cells. Thus, our results support the earlier notion that a broad H3K4me3 promoter marks guard cell identity in stress conditions (Benayoun et al., 2014).

We discovered hypoxic H3K4me3 extension at promoters of stress response genes that may be part of hypoxia preconditioning and prepare endometrial cells to tolerate hypoxic periods in normal reproductive processes. From previous studies, it is known that broad H3K4me3 domains associate with open nucleosome regions, increased Pol II-mediated transcriptional elongation, and chromatin 3D structures involved in dynamic transcription (Benayoun et al., 2014; Chen et al., 2015; Li et al., 2018; Thibodeau et al., 2017). The potential mechanism for the changes in H3K4me3 breath may involve several regulatory layers, such as mRNA and protein levels of both the histone lysine demethylases (KDM) and histone lysine methyltransferases (KMT) targeting histone H3K4. Notably, because KDM5 enzymes demethylate H3K4me3 in oxygen dependent manner, H3K4me3 levels may constitute a rapidly inducible switch in hypoxia regulation (Batie et al., 2019). Recently it was reported that broad H3K4me3 domains are regulated by KDM5-specific inhibitor (Belhocine et al., 2021), further suggesting that KDM5 enzymes regulate H3K4me3 domain breadth and that changes in the breadth may link with cellular stress adaptation.

Based on our results, the genes with H3K4me3 extension are potentially part of the responses induced in the ectopic endometrium by hypoxic niche after retrograde menstruation induced endometriosis (Sampson, 1927; Zondervan et al., 2018). In this process, the stromal cells that lose endometrial blood supply are shed out of the uterus and constitute endometriotic lesions. Our results add to the previous finding (Horne et al., 2019) on the oxygen-dependent gene regulation as a potential nonhormonal therapeutic intervention for endometriosis. More generally, the described dynamics of H3K4me3 breadth upon hypoxia may present a relevant example of epigenetic regulatory dynamics in cellular adaptation to environmental stress.

Limitations of the study

In this study, we conducted H3K4me3 ChIP-seq in hypoxia-treated endometrial stromal fibroblasts (ESF) and differentiated decidual stromal cells (DSC) and combined the results with previous hypoxic and endometriosis transcriptome data. This study was mainly limited to the analysis of H3K4me3 histone modification and transcriptomics, and no protein levels or enzyme activities of KDM5 s were analyzed. In the future, for more accurate characterization of in vivo endometrial or endometriosis-related chromatin states, ChIP-seq experiments using extracted in vivo tissue samples are required.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Kalle T. Rytkönen (katury@utu.fi).

Materials availability

This study did not generate new unique reagents or materials.

Experimental model and subject details

Cell culture

Human immortalized endometrial stromal fibroblasts (ESF) (T HESC, Mor lab, Yale University, corresponding to ATCC CRL-4003) were grown in normoxia (21% O2) in Dulbecco’s Modified Eagle’s medium (DMEM) (Sigma-Aldrich), supplemented with 10% charcoal-treated fetal bovine serum (Hyclone), 1% antibiotic/antimycotic (Gibco), 1 nM sodium pyruvate (Gibco), 0.1% insulin-transferrin-selenium (BD Biosciences), and 0.12% sodium bicarbonate. To generate DSC, ESFs were decidualized by adding of 0.5 mM 8-bromoadenosine 3′,5′-cyclic monophosphate (8-Br-cAMP) (Sigma) and 0.5 μM of the synthetic progestin medroxyprogesterone acetate (MPA) in DMEM supplemented with 2% charcoal-treated fetal bovine serum.

For ESF hypoxia exposure was conducted for 16 h with ProOx C21 nitrogen-induced hypoxia system (BioSpherix, Red Field, NY) at 1% O2, 5% CO2 and compared to corresponding normoxic cells from the same cell batch. For DSC, ESF were first decidualized for 40 h and then similarly exposed to hypoxia for 16 h and compared to normoxic DCS from the same cell batch that were decidualized for 48 h. As hypoxia generally slows cellular processes and differentiation, we reasoned that 40 h in normoxia and 16 h in hypoxia would constitute an appropriate comparison to 48-h normoxic differentiation.

Patient samples

Patient samples were collected and processed as described previously (Heinosalo et al., 2018). Briefly, the Joint Ethics Committee of Turku University and Turku University Hospital approved collections (ETMK310/2005 ENDOMET, ETMK34/180/2012 PROENDO) and all study subjects provided written informed consent prior to recruitment. Endometriosis samples were collected during laparoscopy or laparotomy. Samples from three different types of endometriosis were studied, three from stage IV intestinal deep infiltrating endometriosis (donor ages 31, 38 and 48), two from stage IV sacrouteral deep infiltrating endometriosis (donor ages 31 and na), and three from stage IV ovarian endometriosis (women ages 27, 32 and na). As a control group, endometrial biopsies from three healthy, endometriosis-free women (donor ages 37, 43 and 46) undergoing laparoscopic tubal ligation were collected. All the samples studied were collected during the proliferative phase of the menstrual cycle from women that had not used any hormonal medication. Women with other significant disease or medication, suspicion of malignancy, pregnancy or acute infection were excluded. Histological evaluation was performed to confirm the presence of normal endometrial histology in controls and patients and to diagnose endometriosis in patients.

Method details

ChIP-seq experiment

ChIP with H3K4me3 antibody (49-1005 Invitrogen) was conducted as previously described (Lynch et al., 2015) for two biological replicates in each condition. Briefly, cells were harvested, resuspended and cross-linked in suspension by adding 1% fresh formaldehyde with 10 min incubation (RT, rotator). Cross-linking was quenched and cells were harvested by centrifugation and washed twice in cold PBS. Nuclei were extracted and sonicated with Misonix Sonicator S-4000. Fragment sizes between 150 bp and 600 bp were confirmed by gel electrophoresis and subsequently with Bionalyzer. For immunoprecipitation, 10–150 mg of isolated chromosomal DNA was incubated with 10 μg of antibody coupled to the proteinG Dynabeads (Invitrogen) that were prepared following manufactures instructions. Chromatin-antibody-bead complex was incubated overnight at 4°C in ChIP Dilution buffer (0.02% SDS, 2.2% Triton X-100, 2.4 mM EDTA, 33.4 mM Tris pH 8.1, 334 mM NaCl) with protease and phosphatase inhibitors. After incubation the complex was washed two times with IP wash buffer (NaCl) (100 mM Tris pH8.0, 500 mM NaCl, 1% NP-40, 1% deoxycholic acid) followed by two times wash with IP wash buffer (LiCl) (100 mM Tris pH8.0, 500 mM LiCl, 1% NP-40, 1% deoxycholic acid) with 2 min rotation, and once with 1 mL TE buffer. Chromatin was eluted (50 mM Tris pH 8.0, 10 mM EDTA, 1% SDS) by incubation at 65°C, 10 min with agitation. The eluted DNA was incubated at 65°C overnight to reverse the cross-linking. Immunoprecipitated DNA was treated sequentially with RNase A (Invitrogen), Proteinase K (American Bio) and QIAquick PCR purification kit (Qiagen). DNA library preparation and sequencing were performed for two biological replicates on Illumina HiSeq 2000 following the protocol recommended by Illumina with a 75 bp read length at the Yale Center for Genome Analysis, and resulting sequences were submitted to GEO under the series accession number GSE167946.

ChIP-seq peak calling, differential peak calling and peak extension/compression

The quality of the sequenced reads was checked using FastQC tool (Andrews, 2010). Bowtie2 2.2.6 (Langmead and Salzberg, 2012) was used to align the reads against the human reference genome hg19 (UCSC; downloaded from Illumina iGenomes website). The resulting SAM files were converted to BAM files, sorted and merged using samtools 1.2 (Li et al., 2009) The peaks were called with MACS 2.1.0 (Zhang et al., 2008) with the options “--broad --nomodel --extsize 147 --broadcutoff 0.1” enabled and combining the two replicates. We re-computed the height of each peak by scanning each peak for read pile-up, the highest pile-up value was used as the peak height. Peaks were then processed with bedtools 2.17.0 (Quinlan and Hall, 2010) to merge peaks closer than 3 kb and GREAT (McLean et al., 2010) to find promoter associated peaks (TSS −/+5 kb). We used DiffBind 2.14.0 (Ross-Innes et al., 2012) with separate MACS2 peaks for the two biological replicates, default parameters and with DESeq2 (Love et al., 2014) enabled to call differential peaks, and included peaks with q < 1 × 10−10 for the correlations, pathway analysis and Fisher’s exact tests. Chromtime 0.11.4 (Fiziev and Ernst, 2018) with default parameters (q < 0.05) was employed to determine extension and compression between the conditions using the merged BAM files. The Chromtime detected H3K4me3 extensions at gene TSS were inspected for upstream extension (>300 bp) directly from the MACS2 calls. For downstream analysis we filtered DiffBind differentially bound and extended/compressed gene promoters to include only the expressed genes (TPM > 2) (see below). Processing, statistics and visualization were done using R, IGV and METASCAPE (Zhou et al., 2019).

Correlation analysis with RNA-seq data

For the correlation analysis we used the previously published RNA-seq data of normoxic and 24-h hypoxia treated ESF and DSC (GSE111570 (GSM3034449, GSM3034450, GSM3034451, GSM3034452)) and GSE63733 (GSM1556296, GSM1556297, GSM1556298, GSM1556299). Because hypoxia mediated changes in H3K4me3 are faster than transcriptomic changes (Batie et al., 2019) we matched H3K4me3 data from 16-h hypoxia exposure with RNA-seq data from 24-h hypoxia exposure. The RNA-seq data were processed as previously (Rytkönen et al., 2020). For the absolute correlation analysis of K4 height and breadth we used absolute TPM values.

For the fold-change correlations, the fold change values of the differential peaks called by DiffBind were used with RNA-seq fold change values. As ChromTime does not output fold changes, we approximated the fold change of the peak width directly from corresponding MACS2 peaks called for normoxia and hypoxia, and these values were used in the correlation analysis with RNA-seq FC values.

Gene set enrichment analysis and over-representation analysis

Custom gene set enrichment analysis was done using GSEA (www.gsea-msigdb.org/gsea) with user defined RNA-seq datasets instead of databases, including hypoxia RNA-seq (Rytkönen et al., 2020) and PGR targets that were previously detected using siRNA by Demayo lab (Mazur et al., 2015). For standard gene enrichment over-representation analysis and heatmaps, gene lists of top 500 broadest and top 500 tallest H3K4me3 marked promoters, Diffbind H3K4me3 upregulated and H3K4me3 extended genes were used as input for METASCAPE (https://metascape.org/) where all available pathway databases were selected for analysis.

Disease datasets and statistical tests of the intersections

To test the significance of our cell culture results on the H3K4me3 status with clinically relevant and hypoxia associated data we overlapped specific H3K4me3 subsets with endometriosis and endometrial cancer data from literature and databases. Selected endometriosis data sets included a FACS isolated stromal cell dataset of differentially expressed genes between endometriosis lesions versus control healthy endometrium (Rekker et al., 2017) and a heterogeneous tissue dataset of differentially expressed genes in the endometrium of endometriosis patients (endometriosis and abnormal) versus healthy controls (Tamaresis et al., 2014). For cancer we selected Uterine Corpus Endometrial Carcinoma TCGA (UCEC_TCGA) PanCancer (https://www.cbioportal.org/study/summary?id=ucec_tcga_pan_can_atlas_2018) data that included two datasets, mutation and copy number alterations (CNA).

We then overlapped H3K4me3 promoter subsets for hypoxic ESF using Fisher’s exact test. Before the disease overlap the H3K4me3 data was first intersected with differentially expressed genes from RNA-seq hypoxia experiment (Rytkönen et al., 2020). For the 500 broadest H3K4me3 promoters we tested subsets that overlapped with both transcriptomic hypoxia up- and down-regulated genes. For DiffBind H3K4me3 upregulated or extended promoters (genes) we tested subset with genes that were also upregulated in the hypoxic RNA-seq and for DiffBind H3K4me3 downregulated or compressed promoters (genes) we tested subset that were also downregulated in the hypoxic RNA-seq, an visualized the p values in with pheatmap R package.

Immunohistochemistry

Tissue samples were fixed in formalin and embedded in paraffin for histological analysis. Antigen retrieval of hydrated 5 μm thick sections was performed in 10 mM sodium citrate buffer (pH 6.0), followed by immunohistochemistry with primary antibodies against KLF2 (NB100-1051, Novus Biologicals) with 1:100 dilution. Additionally, for ovarian endometriosis samples, to confirm the presence of endometriotic stroma staining with CD10 (NCL-L-CD10-270, Leica Biosystems) dilution 1:40 was conducted. Negative control samples not incubated with primary antibodies did not show any staining. Sections were scanned for analyses with the panoramic 250 Flash series digital slide scanner (3DHISTECH, Hungary). Then for each slide the proportions of stained cells were calculated as from mean from three representative locations (100× magnification), and Welch Two Sample t-test was applied to compare endometriosis samples with the healthy control endometria.

Quantification and statistical analysis

All analytical and statistical tests for ChIP-seq data, RNA-seq data and co-analysis with the publicly available data are outlined above. Any analysis step not specified was performed using R 3.6.1.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
H3K4me3	Invitrogen	Cat# 49-1005; RRID:AB_2533856
KLF2	Novus Biologicals	Cat# NB100-1051; RRID:AB_2132571
CD10	Leica Biosystems	Cat# NCL-L-CD10-270; RRID:AB_563478
Biological samples
Endometriosis biopsies from patients	Turku University Hospital	Heinosalo et al. (2018)
Endometrial biopsies from healthy women	Turku University Hospital	Heinosalo et al. (2018)
Chemicals, peptides, and recombinant proteins
Dulbecco’s Modified Eagle’s medium (DMEM)	Sigma-Aldrich	Cat# D2906
HyClone Charcoal/Dextran treated fetal bovine serum	HyClone	Cat# SH30068
Antibiotic/antimycotic (ABAM)	Gibco	Cat# 15240
Insulin-transferrin-selenium (ITS premix)	BD Biosciences	Cat# 354351
8-bromoadenosine 3′,5′-cyclic monophosphate (8-Br-cAMP)	Sigma-Aldrich	Cat# B7880
Medroxyprogesterone acetate (MPA)	Sigma-Aldrich	Cat# M1629
PureLink RNase A	Invitrogen	Cat# 12091021
Proteinase K	American Bio	Cat# AB00925-00100
Sodium pyruvate	Gibco	Cat# 11360070
Deposited data
H3K4me3 ChIP-seq	This study	GEO: GSE167946
Publicly available RNA-seq data	GEO	GEO: GSM3034449
Publicly available RNA-seq data	GEO	GEO: GSM3034450
Publicly available RNA-seq data	GEO	GEO: GSM3034451
Publicly available RNA-seq data	GEO	GEO: GSM3034452
Publicly available RNA-seq data	GEO	GEO: GSM1556296
Publicly available RNA-seq data	GEO	GEO: GSM1556297
Publicly available RNA-seq data	GEO	GEO: GSM1556298
Publicly available RNA-seq data	GEO	GEO: GSM1556299
Experimental models: Cell lines
Human immortalized endometrial stromal fibroblasts (T HESC) (ESF)	Mor lab, Yale University	T HESC corresponding to ATCC CRL-4003
Software and algorithms
FastQC tool	(Andrews, 2010)	https://www.bioinformatics.babraham.ac.uk/projects/fastqc/
Bowtie2 2.2.6	(Langmead and Salzberg, 2012)	https://github.com/BenLangmead/bowtie2
Samtools 1.2	(Li et al., 2009)	https://github.com/samtools/samtools/
R 3.6.1	CRAN	https://cran.r-project.org/
MACS 2.1.0	(Zhang et al., 2008)	https://pypi.org/project/MACS2/2.1.0.20140616/
GREAT	(McLean et al., 2010)	http://great.stanford.edu/public/html/
DiffBind 2.14.0	(Ross-Innes et al., 2012)	https://bioconductor.org/packages/release/bioc/html/DiffBind.html
DESeq2	(Love et al., 2014)	https://bioconductor.org/packages/release/bioc/html/DESeq2.html
Chromtime 0.11.4	(Fiziev and Ernst, 2018)	https://github.com/ernstlab/ChromTime
IGV	Broad Institute	https://software.broadinstitute.org/software/igv/
METASCAPE	Zhou et al., 2019	https://metascape.org/
GSEA	Broad Institute	www.gsea-msigdb.org/gsea
ProOx C21 nitrogen-induced hypoxia system	BioSpherix	https://biospherix.com/proox-c21/
Misonix Sonicator S-4000	Fisher	Cat# 13101353
ProteinG Dynabeads	Invitrogen	Cat# 10004D
QIAquick PCR purification kit	Qiagen	Cat# 28106
Pannoramic 250 Flash III digital slide scanner	3DHISTECH	https://www.3dhistech.com/research/pannoramic-digital-slide-scanners/pannoramic-250-flash-iii/