ERβ is expressed in granulosa cells throughout folliculogenesis

ERβ is known to be expressed in granulosa cells, but its expression in other cells including theca cells has been debated [19]. We performed immunohistochemistry (IHC) with an ERβ antibody (PPZ0506) that has been thoroughly validated across different tissues in humans and rodents [17, 19,20,21]. We confirmed that ERβ was strongly expressed in the nucleus of granulosa cells. Some nuclear staining was further observed in cells surrounding the theca layer and in a few cells in the stroma, but not within the theca cell layer itself (Fig. 1A, upper panel). In ERβKO ovaries, ERβ expression (IHC) was absent from the granulosa cells (Fig. 1A, lower panel). However, some cytoplasmic staining, mostly of stroma cells, was still present and was deemed to be non-specific. We corroborated the ERβ expression using ERβ RNA in situ hybridization. This confirmed its high expression in granulosa cells and showed clear evidence of ERβ in the layer surrounding the theca cells of the follicle, which may correspond to the follicular microvasculature, and in a minority of other stromal cells, while the theca cells themself appeared blank (Fig. 1B). ERβ expression was further corroborated by western blot (WB) using the same antibody and by qPCR directed towards the exon deleted in the knockout (Fig. 1C-D). The WB of the WT ovary reveals two bands of similar intensity near the expected size. These are likely to correspond to the two murine splice variants of ERβ, the 567 amino acid (aa) ERβ_ins (isoform 1 / NM_207707.1) and the 549 aa ERβ (isoform 2 / NM_010157.3) which differ by an 18-aa long sequence (approx. 2 kDa) of the ligand binding domain. The epitope for the antibody is at the N-terminal part of the receptor (corresponding to the first coding exon), and the antibody thus recognizes both splice variants. The expression of the two splice variants was confirmed at the gene expression level (Additional file 1: Fig. S1). To be noted, the knockout mice have exon 3 deleted which creates a frameshift and premature termination during translation. Had any truncated proteins or peptides been expressed, these should have been recognized by the antibody. Since the granulosa cells of the knockout were primarily negative by IHC, we determined that ovarian ERβ is fully absent following knockout. Finally, we investigated if ERβ protein expression differed depending on the follicular stage (primary, small preantral, large preantral, antral), but did not find a significant change (Fig. 1E). We conclude that ERβ protein is highly and consistently expressed in granulosa cells during folliculogenesis and that ERβ is also expressed in cells surrounding the theca cells.

Fig. 1

ERβ is expressed in granulosa cells and some stromal cells of the ovary. A IHC with the validated antibody PPZ0506 in WT (upper panel, bar indicates 20 μm in the left panel and 50 μm in all other, black arrows indicate granulosa cells and white theca cells, n = 4) and ERβKO (lower panel, n = 4) mouse ovary. B In situ hybridization with a probe against Esr2 in WT female ovary (n = 2 mice, one representative ovary shown in image) confirms ERβ expression in granulosa cells, in cells surrounding the outer layer of theca cells, and in some stromal cells in mouse ovary. C Western blot (one of 2 replicates, each with 2 pooled ovaries), and D qPCR analysis (n = 5, paired t-test) corroborates the loss of ERβ expression in ERβKO ovaries. Both mouse ERβ isoforms (549 aa, calculated molecular weight 61.2 kDa and 567 aa ERβ_ins with 18 inserted aa,, calculated molecular weight 63.2 kDa) are visible near the expected sizes, with both bands being absent in the knockout ovary. E Assessment of ERβ protein expression by IHC in WT mice (n = 11), separated by follicular stage and each follicle scored according to staining intensity and area

The endogenous ovarian ERβ cistrome

ERβ is located within the nucleus (as can be noted in Fig. 1A) and functions as a transcription factor. To explore its endogenous genome-wide chromatin binding in the ovary, we performed ChIP-seq of WT mouse ovaries with the same ERβ antibody (PPZ0506) that was previously optimized for ChIP of human ERβ in cell lines [33, 34]. The experiment was performed in biological triplicates and compared to inputs. More than 80% of the produced sequencing reads were of high quality, and between 25 and 31 M reads per ChIP sample were aligned to the genome (GRCm38, Additional file 2: Table S1). A heatmap and Venn diagram illustrate that the majority of binding sites were detected in all three replicates (3175 sites, Fig. 2A-B). As many as 4875 ERβ-binding sites were detected in at least two ChIP replicates and were used for further analysis. The ERβ ChIP-seq was further validated by using ChIP-seq data from ERβKO ovaries instead of input (> 70% of sites confirmed, Additional file 1: Fig. S2A-B and Additional file 3: Table S2). In accordance with ERβ cistromes from other cell types [34,35,36,37], the majority of binding sites were located in introns and intergenic regions (87%), and a minority (5%) in promoter regions (-1 kilobase (kb) to +100 base pairs (bp) from the gene transcription start sites, TSS, Fig. 2C). Notably, the intronic binding sites were most frequently located in intron 1 (36%) or intron 2 (18%).

Fig. 2

Genome-wide landscape of ERβ chromatin binding in mouse ovary. A ERβ chromatin-binding sites in mouse ovary per ChIP-seq replicate (n = 3 replicates, each with 14 pooled ovaries from 7 mice) in relation to corresponding input samples, visualized in a heatmap. B Venn diagram of detected ERβ-binding sites in ChIP-seq triplicates normalized against input. Sites present in at least two replicates, colored in blue, were used in further analysis. C The genomic distribution of ERβ-binding sites. All regions, except intergenic regions, were significantly enriched (p < 0.005). D Top-enriched DNA motifs among ERβ-bound genomic sequences, identified using HOMER de novo motif analysis (sorted by p-value). % sequences represent the percentage of ChIP:ed sequences (bound sites) that have a particular motif. The percentages will not add up to 100% as each sequence can have more than one motif. E Enriched biological functions among genes nearest ERβ chromatin-bound sites (within −50 kb to +2 kb)

We next determined which transcription factor binding sequence motifs were significantly enriched in the ERβ-bound DNA. Using HOMER, we identified that the estrogen response element (ERE) was the most enriched (p = 1e−1814) and abundant motif (present in 45.9% of the bound sequences whereas the background frequency is 3.6%), as would be expected (Fig. 2D). Further, motifs of the well-known ER-tethering and pioneering factors AP-1, GATA, and FOXO were highly enriched (Fig. 2D). This confirmed the specificity of the antibody and the accuracy of the experiment. Interestingly, the NR5A motif was the second most enriched motif and as abundant (44.8% of bound sequences, significant enrichment over the 16.6% background frequency, p = 1e−459) as the ERE. The NR5A motif can be bound by two nuclear receptors: the steroidogenic factor 1 (SF-1/Nr5a1) and the liver homolog 1 (LRH-1/Nr5a2). Both are expressed in the mouse ovary and essential for ovulation [38, 39]. While SF-1 has been observed in relation to ERβ previously, an interaction with LRH-1 has not been reported. We found no difference in the types of motifs bound in promoter regions versus more distant (enhancer) regions (Additional file 1: Fig. S2C).

Finally, we mapped which genes’ TSS were located closest to each ERβ-binding site (Additional file 3: Table S2). Biological pathway analysis on those genes’ functions revealed enrichments related to lipid metabolism, cell proliferation, cell differentiation, multicellular organism development, cell adhesion, transcription regulation, and apoptosis (Fig. 2E). These are functions that ERβ is known (transcription regulation, cell proliferation, apoptosis, cell adhesion) or proposed (cell differentiation, multicellular organism development, lipid metabolism) to modulate in various cell types. This significant enrichment suggests that ovarian ERβ has the potential to impact these functions, although not all bound genes may be de facto regulated. Thus, we here describe the complete endogenous ovarian ERβ cistrome for the first time.

ERβ influences the ovarian transcriptomic landscape

To map the consequences of ERβ deletion on the ovarian transcriptional landscape, we performed RNA-seq on ovaries from WT (n = 5) and ERβKO (n = 4) mice. We identified 803 differentially expressed genes (DEGs) following the loss of ERβ, with a relatively uniform distribution between up- and downregulated genes (375 upregulated and 427 downregulated, Fig. 3A, Additional file 4: Table S3). Thrombospondin 4 (Thbs4, adhesion) was the most significantly downregulated gene following loss of ERβ (Fig. 3A). Thbs4 is abundant in normal ovary and related to the polycystic ovary syndrome-associated gene Thbs1 [40]. Further, as would be expected, genes involved in response to estrogen were significantly enriched and specifically downregulated (incl., Star, Igfbp2, Tgfbr1, Pdgfrb) in the absence of ERβ (Fig. 3B, Additional file 4: Table S3). Other enriched functions among the downregulated genes included cell adhesion (e.g., Thbs4, Itga1, Itga2, Itga4, Itga5, Cdh11), oxidation-reduction (e.g., Cyp11a1, Hsd17b7), and ion transport (incl. potassium channels Kcnab3, Kcnd2, Kcnj16, Kcnk3, Kcnma1) (Fig. 3B). The ovarian expression of these genes is thus indicated to be upregulated as a consequence of ERβ expression. Among the genes upregulated in absence of ERβ (i.e., downregulated as a consequence of ERβ expression), we find FSH and LH targets (incl. Fshr, Cyp11a1, Gata4, Runx2) (Fig. 3A), along with genes related to neural crest cell migration (Sema3b, Sema3c, Sema3e, and Sema3g), fatty acid metabolism (e.g., Fasn), TGFβ receptor signaling (e.g., Smad6), and male gonad development (e.g., Inha, Kitl, Gata1, Gata4) (Fig. 3B). Moreover, the male sex determination gene desert hedgehog (Dhh) was upregulated. We further noted upregulation of plexin C1 (Plxnc1) in ERβKO ovaries (Fig. 3A). Plxnc1 is related to Plxnb1 which is involved in mouse follicular development [41]. Interestingly, Greb1, which is upregulated both by ERα [42] and ERβ (engineered expression, [43]) in human breast cancer cells, was found among these genes. These functions and genes are thus potentially repressed by ERβ in the ovary. Genes with functions in lipid metabolism (e.g., Fads6 down, Lep up) and angiogenesis (e.g., Angpt1 down, Angpt2 and Angpt4 up) were enriched among both up- and downregulated genes. When all DEGs (regardless of direction) were analyzed for biological process enrichment, lipid metabolism, glucose homeostasis, response to hypoxia, response to stimulus, and multicellular organism development were among the enriched functions (Fig. 3C, Additional file 4: Table S3). We conclude that deletion of ERβ impacts several pathways essential for normal ovarian function, including the repression of male gonad development.

Fig. 3

Impact of ERβ on the ovarian transcriptional landscape. A Volcano plot of ERβ-regulated genes. Genes were considered differentially expressed when FDR < 0.05 and log2FC >|0.4| (n = 5 WT, n = 4 ERβKO). B Top-10 enriched biological pathways of the down- (blue) and upregulated (red) DEGs. The size of the bubbles corresponds to the enrichment score (-log10(p-value)). C All significantly enriched biological pathways visualized in semantic space. The size and color of the bubbles correspond to the enrichment score (-log10(p-value)) and gene count, respectively. D Bar chart presenting the projected cell type abundance in the ovary of WT and ERβKO mice (Wilcoxon rank-sum exact test)

ERβ impacts the ovarian cell composition

Since the analysis above identified the gene expression of the complete ovary, we next used the RNA-seq data to investigate whether the loss of ERβ impacted the ovarian cell composition. By applying digital cytometry, using gene signatures from published mouse ovary single-cell RNA-seq data [44] along with CIBERSORTx [45], on our bulk RNA-seq data, we could estimate the abundance of different cell types. This identified in total 10 different cell types encompassing two types of granulosa cells (Inhbahigh and Kctd14high), two types of luteal cells (regular and large), theca cells, two types of cumulus cells (Nupr1high and Ube2chigh), ovarian surface epithelial cells, macrophages (Lyz2high), and stromal cells. Following the loss of ERβ, a significant decrease in luteal cells, along with an increase of theca and Nupr1high cumulus cell populations were indicated (Fig. 3D). The decrease in luteal cells confirms previous histological observations and can be directly linked to the impaired ovulatory phenotype that results in a reduced formation of corpus luteum in the ERβKO mice. Overall, our transcriptome analysis demonstrates a clear impact of ERβ on the cellular composition of the ovary.

Characterizing ovarian ERβ target genes

In an effort to identify the direct ERβ transcriptional target genes in the ovary, we integrated our RNA-seq and ChIP-seq data. This revealed that as much as a third (30%, or 245 out of 803 genes) of the ERβKO ovarian DEGs were located nearest (their TSS) to one or multiple ERβ-binding chromatin sites (Fig. 4A). Assuming a connection between the ERβ chromatin binding and the subsequent transcript regulation upon its loss, we here denote these 245 genes as direct transcriptional targets of endogenous ERβ in the ovary. A slightly larger proportion of these binding sites (52%) contained an ERE motif compared to all ERβ-bound sites (46%). The direct targets were regulated in both directions following the deletion of ERβ (49% up and 51% down). Again, the direct target genes were enriched for functions related to response to estrogen (e.g., Dhh, Pdgfrb, Gata4), lipid metabolism (e.g., Fasn, Cyp11a1), positive regulation of angiogenesis (e.g., Angpt2, Angpt4), cell differentiation (e.g., Etv6), and multicellular organism development (e.g., Cebpa, Pak3, Fzd1, Greb1), (Fig. 4B). Also, response to hypoxia (e.g., Angpt2, Endra) was among the most enriched functions of the target genes (Fig. 4B). Related ERβ chromatin binding is exemplified in Fig. 4C. It can be noted that overall, few of the directly regulated genes had a binding site in the promoter region (13 out of the 245 genes, or 5%), encompassing Inha, Hsd17b1, Gata4, Neat1 (long non-coding RNA), Cpm, Epb41l2, Mylk3, Mdfic, Zfp219, Epb41l1, Pik3cd, Skil, and Fosl2. Notably, the majority of the direct targets (158 out of 245 or 64%) had an ERβ binding site in an intron, most often (86 genes, or 35%) in intron 1. This included Greb1 (4 binding sites in intron 1), Angpt4, Pak3 and Bcl2 (each 2 binding sites in intron 1), and Dhh and Pdgfrb (each 1 binding site in intron 1). Thus, a primary mechanism whereby ERβ regulates genes appears to involve its binding to intron 1.

Fig. 4

Identification of direct ERβ-regulated genes and biological pathways in mouse ovary. A Venn diagram representing the genes overlapping between ERβ ChIP-seq in WT ovaries (n = 3 replicates, each with 14 pooled ovaries from 7 mice) and RNA-seq performed in WT (n = 5) and ERβKO ovaries (n = 4). B Circle plot representing the top-10 enriched biological functions of the 245 overlapping genes between ChIP-seq and RNA-seq (from A). For each biological function, up- (red) and downregulated (blue) genes are represented. C ChIP peaks show enrichment of ERβ at chromatin by genes involved in specific functions. D ChIP enrichment at binding regions nearest to Esr2 and Nr5a2. All tracks are set to the same Y-axis height for the ChIP-seq and input. E The top-10 enriched biological functions of the regulated genes identified by both ChIP-seq and microarray (before and after ovulatory signal, from [29]). The size of the bubbles corresponds to the enrichment score (-log10(p-value)). F ERβ ChIP signal at Fshr and Tgfbr1

Since nuclear receptors may crosstalk at several levels (regulation, interaction, shared transcriptional targets), we next explored our datasets to identify potential nuclear receptor regulation by ERβ. We found that the TSS of multiple nuclear receptor genes were located closest to several ovarian ERβ-binding sites and could thus be potential direct targets (Additional file 5: Table S4). This included ERβ itself (Esr2, three binding sites visualized in Fig. 4D) and LRH-1 (Nr5a2, five binding sites, Fig. 4D), indicating a potential regulatory loop. Similarly, SF-1 (Nr5a1) had two intronic ERβ-binding sites. However, neither LRH-1 nor SF-1 transcripts were significantly regulated in the ERβKO ovary. The remaining nuclear receptor genes that were located by ERβ-binding sites (ERα, PR/Pgr, SHP/Nr0b2, PXR/Nr1i2, GR/Nr3c1, NUR77/Nr4a1, COUP-TF-I/Nr2f1, and COUP-TF-II/Nr2f2), as well as several nuclear receptor coregulators (incl. Nrip1, SRC1/Ncoa1, and Ncor), were also not detected as significantly regulated in the knockout ovaries. This lack of transcriptional regulation is in line with the result generated by Binder and colleagues using isolated granulosa cells from WT and ERβKO ovaries, where only Pgr was differentially expressed following ERβ knockout after ovulatory stimuli [29, 46]. To conclude, although ERβ can bind to chromatin regions close to several nuclear receptors and potentially regulate them, we did not find evidence that it does so in the ovary.

Finally, to explore which direct targets may be most consequential for ovarian function, we searched for genes where ERβ would be essential for their expression in the ovary. That is, genes whose transcripts were near absent (< 2 Fragments Per Kilobase of transcript per Million mapped reads, FPKM) in the ERβKO ovary but at least 2-fold upregulated in the WT ovary, and with high confidence (FDR < 0.0001) across the individual samples. The expression of 14 genes fulfilled these strict requirements (Table 1). Eight of these had an ERβ chromatin-binding site, nearly all (7/8) located in an intron. Interestingly, a majority (9 out of 14, or 64%) encoded for proteins located in membranes such as cell surface receptors for low density lipoprotein (Lrp8, Lrp11) and catecholamine (Adrb2). Thus, the proteins that appear to be most dependent on ERβ for their ovarian expression are located in cellular membranes with roles in cell signaling, and they are regulated by ERβ through intronic chromatin binding. Taken together our study identifies direct targets of ERβ which furthers our understanding of its role in the ovary.

Table 1 Ovarian genes that appear dependent on ERβ for their in vivo expressionERβ gene regulation in granulosa cells during the ovulatory process

Previously, microarray comparisons of granulosa cells collected from large antral follicles of WT and ERβKO mice, before and after ovulatory signal, have identified 1361 genes related to ERβ expression [29, 46]. Only 5% (70 genes) of those were detected both before and after ovulatory signaling (Additional file 1: Fig. S2D). Since ERβ is primarily expressed in granulosa cells, we compared this data with our ERβ results from the whole ovary. While only 12% (164 of 1361) of the genes detected as regulated in the microarray were differentially expressed also in the WT and ERβKO ovarian RNA-seq (Additional file 1: Fig. S2D), the exact same proportion of genes (30%, 413 of 1361 vs 245 of 803 genes detected by RNA-seq) had an ERβ chromatin-binding site per our cistrome (Fig. 4A and Additional file 1: Fig. S2D). Exploring the pathway enrichment of all ERβ-bound (per our ChIP-seq) microarray-detected genes (413) showed that, similarly to our ERβ-bound genes in RNA-seq, developmental and differentiation genes were regulated both before and after ovulatory signal whereas other pathways were restricted to either condition (Fig. 4E). For example, genes involved in cell migration and cell adhesion were primarily regulated before the induction of ovulation, while target genes with functions in response to insulin, and cholesterol or lipoprotein metabolic processes were regulated after ovulatory signal. Genes with an ERβ chromatin-binding site, whose expression were detected as regulated in both studies (77 genes), included the key ovarian genes Dhh (Fig. 4C), Fshr, and Tgfbr1 (Fig. 4F). Thus, several of the same genes and functions were detected in both studies, but some (including response to estrogen and angiogenesis) were only detected in our in vivo study.

Conserved ERβ regulation between species

To explore a potential relation to human health, we compared the genes located by ERβ-bound sites in the mouse ovary (endogenous ERβ) with those previously assessed in human cell lines (only ERβ ChIP-seq data from exogenous ERβ in other cell types than ovary were available), and performed with the same validated antibody [34, 47]. This analysis showed that over one-third (1332 of 3472, or 38%) of the genes bound by ovarian ERβ (i.e., whose TSS was the closest located to the bound sites) in the mouse genome, were also bound by ERβ in human cells (Additional file 6: Table S5). Moreover, most of these genes were bound in their intronic regions in both species. In contrast, the genes bound in their promoter, transcription termination site (TTS), exon, 3’ or 5’ UTR regions, rarely had such conserved location. We also note that most of the binding sites in genes bound by ERβ in both species (1332 genes, corresponding to 2038 linked binding sites in mice) had an ERE motif (1138 of 2038, or 56%), and a nearly as large proportion (923 of 2038, or 45%) contained the NR5A motif, most commonly in combination with an ERE. Thus, also in the conserved genes, both ERE and NR5A were more frequent than the occurrence of the pioneering GATA motif (561 sites). When we specifically looked at the ovarian direct target genes (the 245 genes that were also regulated at the transcript level in the mouse), more than half of them (54% or 132 genes) also exhibited an ERβ binding site in the human cells, again most commonly in the introns (72 out of 132 genes). Here, however, the conservation was strong for the few promoter-bound target genes (mouse ovary) where nearly all (11/13 genes or 85%) were bound by ERβ in humans too, albeit only three (INHA, EPB41L2, SKIL) were bound at the promoter also in humans. These three direct targets with conserved promoter binding were either upregulated (Skil and Epb41l2) or repressed (Inha) by ERβ in the mouse ovary (Additional file 6: Table S5). Finally, we compared our cistrome and transcriptome to human granulosa-enriched genes from the Human Protein Atlas [48]. Out of the genes defined as granulosa-enriched in humans, over a quarter (27%, or 133 of 496 genes) were either bound by ERβ (mouse or human) and/or regulated following ERβ knockout (mouse ovary) in our data sets (Additional file 6: Table S5). Twenty-four of these genes were bound and regulated (here denoted direct ERβ targets) in the mouse ovary (incl. DHH, FSHR, GATA4, GREB1, HSD17B1, INHA, INHBB, LRP5, and LRP8). Thus, although there are major species differences between human and murine fertility, the chromatin binding of ERβ on key targets is relatively conserved. This indicates that key ovarian genes (e.g., INHA, EPB41L2, and SKIL) and transcription factors (e.g., GATA6, LRH-1, and PPARG) are of particular importance for ERβ-mediated female fertility also in humans.

ERβ may co-regulate targets with LRH-1 in granulosa cells

Two species-conserved ERβ-bound genes of particular importance, LRH-1 and SF-1, are nuclear receptors that bind the DNA motif (NR5A) that we found enriched at a considerable fraction (44%, p = 1e−459, Fig. 2D). This suggests that ERβ may regulate these nuclear receptors and also function in the same chromatin-bound complex as either of them. In other tissues, these transcription factors are known to recruit, and thereby enable the function of, other nuclear receptors. For example, LRH-1 enables the function of FXR in the liver [49], and SF-1 can recruit DAX-1 (Nr0b1) to promoters in the adrenal gland [50]. Exploring all ovarian ERβ-bound sequences for presence of ERE, NR5A, or GATA (pioneering factor) motifs, we found that most (4105 out of 4875, 84%) harbored at least one of these motifs (Fig. 5A). Among the sequences harboring only one of these three motifs, ERE was predominant (1199), followed by NR5A (747) which was about twice as common as GATA (313). The combination ERE and NR5A at one single binding site (1092) was more prevalent than ERE and GATA (672). That the NR5A motif was more common than the pioneering factor GATA motif could indicate that an NR5A nuclear receptor may co-regulate ovarian genes with ERβ and/or bind first and bring ERβ in. To investigate if the identified NR5A motifs were de facto bound by an NR5A nuclear receptor (LRH-1 or SF-1) in the ovary, we used a publicly available LRH-1 ChIP-seq data set from isolated mouse granulosa cells (GSE119508, [51, 52]). We investigated whether LRH-1 chromatin-binding sites overlapped with ERβ-binding sites within the ChIP-seq peak sizes (200 nucleotides, nt). We identified that over a third of all ERβ-bound sequences (1740 out of 4875, 36%) were indeed de facto also bound by LRH-1 in granulosa cells (Fig. 5B, Additional file 7: Table S6). This demonstrates that LRH-1 and ERβ, to a large extent, bind the same (200 nt) regulatory chromatin sites in the ovarian genome. This supports the interpretation that LRH-1 is important for ERβ function in the ovary, or vice versa.

Fig. 5

ERβ shares ovarian cistrome with LRH-1. A Venn diagram comparing distribution and co-occurrence of ERE, GATA, and NR5A motifs in all identified ERβ-bound sequences. B Venn diagram representing the overlapping chromatin-binding sites between ERβ and LRH-1. C Enriched biological pathways, summarized by REVIGO, comparing genes located nearest to chromatin bound by both ERβ and LRH-1 or specifically bound by either ERβ or LRH-1. The size of the bubbles corresponds to the enrichment score (-log10(p-value)). D Venn diagram comparing distribution and co-occurrence of ERE and NR5A motifs in the the sequences bound by both ERβ and LRH-1. E The distance between ERE and NR5A motifs in the sequences containing dual motifs and bound by both ERβ and LRH-1 (from D) plotted within +/- 200 bp. F, G Density plots representing motif occurrence within +/- 1500 bp distance of F all common ERβ and LRH-1 binding sites, and G all ERβ-binding sites. H Venn diagram comparing ERβ ChIP-seq with accessible chromatin (genes) using publicly available FAIRE-seq data [51] from granulosa cells during ovulatory stimulation. I, J ChIP enrichment and corresponding gene regulation (RNA-seq) of critical genes that require ERβ. All tracks are set to the same Y-axis height for the ChIP-seq and input

Next, we explored potential differences in function (enriched biological pathways) between genes co-bound by ERβ and LRH-1 and those bound by only one receptor. To condense the list, all biological pathways were run through REVIGO (top-ten functions illustrated in Fig. 5C). We determined that genes bound by either nuclear receptor (whether alone or together) were strongly enriched for genes involved in apoptosis. The ERβ-only bound genes were specifically enriched for PI3K signaling, extracellular matrix organization, fertilization, ERK1/ERK2 cascade, male gonad development, and positive regulation of NFκB transcription factor activity. The LRH-1-only pathways, on the other hand, were strongly enriched for chromatin organization and cell differentiation. Finally, the co-bound genes had multicellular organism development as the most enriched function (ERβ-only and LRH-1-only were also enriched for various developmental genes) and were uniquely enriched for lipid metabolic process and the TGFβ receptor signaling pathway. We thus note clear differences in functions of the genes that ERβ and LRH1 both bind (at close distance), compared to those bound only by one of the receptors.

ERβ-LRH-1 chromatin interactions

The co-binding of two ovulation regulators, ERβ and LRH-1, to the same ovarian DNA locations, points in the direction that this is a key ovarian molecular mechanism. To better understand this mechanism, we questioned whether LRH-1 might recruit ERβ (bound DNA harbors only NR5A motif), if ERβ might recruit LRH-1 (only ERE motifs), if both may bind simultaneously, or if they compete. In an effort to explore this, we first investigated to which extent the de facto dual-bound sequences (approx. 200 nt) harbored both NR5A and ERE motifs. We found that the largest group (702 out of 1740 de facto dual bound sites, or 40%) harbored an NR5A but not an ERE motif. The second largest group (532 sites, 31%) held both an ERE and an NR5A motif (Fig. 5D), whereas the smallest group (274 sites, 16%) had an ERE without an NR5A motif (the remaining 13% lacked both two motifs). Altogether, this may indicate that LRH-1 first binds to the NR5A motif and brings ERβ into the complex.

However, a second mechanism is possible for a sizeable fraction (31%) of chromatin events where both motifs are present: LRH-1 and ERβ either bind adjacent to each other or compete for binding. To explore the latter question, we first investigated the distance between the dual ERE and NR5A motifs (532 sites de facto bound by both receptors, although not necessarily at the same time). We found that in nearly all cases (92%), the two motifs were within 100 bp of each other (Fig. 5E). Further, the density plot clearly shows that both motifs had an equally high probability occurrence within a 500 bp distance of the binding site for these dual-bound sequences (Fig. 5F). This is to be compared with the density plot for all ERβ-binding sites where the ERE has the highest probability occurrence, and NR5A comes a distant second (Fig. 5G). A technical caveat is that the NR5A motif sequence overlaps with the complementary sequence of the ERE half-site motif with all but 3 additional nucleotides, making this analysis sensitive to artifacts. For this reason, we separated the sites where the NR5A motif directly or within 10 bases overlapped with an ERE motif and found that this only occurred for 88 out of the 532 sequences and did thus not impact the overall results. Evidently, the NR5A and ERE binding sites were commonly clearly separated, supporting that binding of the two nuclear receptors to their respective motifs near each other may occur.

Finally, we investigated to what extent the identified ERβ-binding sites were present in active chromatin regions. We compared our ERβ ChIP-seq data with Formaldehyde-Assisted Isolation of Regulatory Elements (FAIRE)-seq data performed in the granulosa cells under the same settings as the LRH-1 ChIP-seq [51]. We identified that 69% of the ERβ-bound genes were located nearby active chromatin regulatory regions (Fig. 5H). According to the FAIRE-seq data, the top motifs in the active chromatin regions were LRH-1 and CTCF, while ERE motifs represented 10% of all binding sites. Thus, as this shows that ERβ primarily binds