The transition of endometrial stromal cells from a mesenchymal to an epithelial phenotype (MET) is an essential prerequisite to blastocyst implantation (16, 17). hESCs were isolated from endometrial tissue biopsies and the response to a decidualization stimulus determined by monitoring the levels of prolactin. Decidualization was induced by incubating hESCs in the presence of cAMP (0.5mM) or cAMP (0.5 mM) plus MPA (1μM) for 48 h. WT1 mRNA levels were significantly (p < 0.05) increased in hESCs isolated from fertile donors in response to cAMP +/- MPA, but not in endometrial PCOS cells ( Figure 1A ). Decidual prolactin (dPRL) levels were determined from samples cultured ex vivo and found to be increased in cells obtained from fertile individuals, whereas there was no response from PCOS derived stromal cells ( Figure 1B ).

Having established that stromal cells derived from fertile donors were functionally active, the genome-wide targets of WT1 were evaluated to determine how the ‘gatekeeper’ role that we have previously proposed may manifest itself (7).

We undertook the mapping of genome-wide WT1 DNA-binding (deposited as GSE240055 in the NCBI GEO repository) in hESCs obtained from fertile patients treated with cAMP (0.5 mM). This de novo data was analysed and identified 19,417 called peaks (regions bound by WT1) (FDR < 0.05) revealing that WT1 DNA binding occurs across the genome ( Figure 2A ). DNA binding site motif analysis on genomic sequences contained within the WT1 binding peaks was used to find any consensus binding sites. As expected, MEME-ChIP identified a motif (E-value = 9.1e-47) that mapped to the WT1 motif (E-value = 5.35e-6) (18) ( Figure 2B ). Interestingly a second highly significant motif (E-value = 1.5e-81) was identified that mapped to the ANDR (AR) motif (E-value = 3.76e-2) ( Figure 2D ). The observation that AR binding sites are found in the immediate proximity of a high proportion of WT1 binding sites suggested that the two transcription factors could colocalise or that competitive exclusion of WT by AR could occur.

A more in-depth analysis identified 34% of peaks located within gene promoter regions close to the transcription start site (TSS, Figure 2C ), with the majority of WT1 binding (29% of peaks) occurring within 1kb of the TSS suggesting that WT1 plays a role in transcription regulation through interactions with the core transcription machinery that assembles at these proximal promoter locations. WT1 was also seen to bind extensively within gene bodies, primarily within intron 1 (13% of peaks) and other introns (17% of peaks) suggesting a potential role in intron retention (19). There was also extensive WT1 binding at distal intergenic regions (23% of peaks, as defined by ChIPseeker) suggesting WT1 may function at enhancer sites.

WT1 binding was found to occur in the promoter regions of a large number of FOX genes (20 from 21 FOX genes, P-value < 0.001, Figure 2E and Supplementary Table S1 ) and homeobox (HOX) genes (32 from 34 HOX genes, P-value < 0.001, Figure 2F , Supplementary Table S1 ). HOX genes are essential for endometrial development and endometrial receptivity (20), and loss of FOX gene regulation has been linked PCOS (21), which suggests that the ‘gate keeper’ function of WT1 could be in the regulation of forkhead and homeobox genes during decidualization. Unsurprisingly therefore, over-representation analysis (ORA) (22) revealed ‘activation of HOX genes during differentiation’ (FDR < 0.001) as a major process, as well as ‘Estrogen dependent gene expression’ (FDR < 0.001), and Gene Ontology (GO) (23) biological process pathways revealed ‘Wnt signalling pathway’ (FDR < 0.001), ‘angiogenesis’ (FDR< 0.01) and ‘myeloid leukocyte cytokine production’ (FDR < 0.1) as significantly enriched ( Figure 2G ). During decidualization, endometrial stromal cells undergo morphological and functional changes, becoming mesenchymal, to allow implantation of the embryo (17, 24). Progression into the secretory phase of the menstrual cycle is governed by increases in the expression of estrogen (25). Wnt signalling is implicated in a number of implantation and decidualization events during mammalian pregnancy, and aberrant Wnt signalling negatively effects these processes (26, 27). Angiogenesis occurs during the secretory phase of the menstrual cycle, forming new blood vessels to supply nutrients in the case of a potential implantation event (28, 29). Many studies have shown the importance of interactions between endometrial cells including hESCs and leukocytes to regulate the processes associated with decidualization (30–32) including vascular remodelling of the decidua and angiogenesis (33).

The promoter regions of a large number (5,341) of DNA coding regions for mRNA, miRNAs and lncRNAs were bound by WT1 in our hESC decidualization model. We sought to determine if this corresponded to differential gene regulation using comparative RNA sequence data analysis. We accessed publicly available transcriptomic data RNA-seq (NCBI GEO repository accession GEO86491) from endometrium samples obtained from fertile patients and identified genes that were differentially expressed during decidualization (34). Comparing differentially expressed genes (DEGs) with WT1 gene targets determined that 598 were significantly upregulated and 836 significantly downregulated during decidualization (FDR < 0.05) ( Figure 2H , Supplementary Table S2 ). These genes included several WT1 HOX genes targets, which displayed significant differential expression between proliferative and secretory phases. HOXB2, HOXB4, HOXB5 and HOXB7 were significantly up-regulated and HOXC8, HOXD8, HOXD9, HOXD10 and HOXD11 were significantly down-regulated in the secretory phase compared with the proliferative phase, and FOXO3 was downregulated in the secretory phase ( Figure 2I , Supplementary Table S1 ). Consistent with our analysis, HOXD10 and HOXD11 have previously been found to be downregulated in the proliferative phase of the menstrual cycle (20). Additionally artificial induction of decidualization in mice elicits an increase in FOXO3 expression that is more pronounced at the implantation site in mouse uteri (35).

We determined the genome-wide localisation of AR (deposited as GSE240055 in the NCBI GEO repository) in hESC isolated from PCOS patients treated with DHT (15) to recapitulate events driven by hyperandrogenemia in vivo. Here we identified 12,017 significant AR peaks (FDR < 0.05%), of which only 5% were localised in gene promoter regions ( Figures 3A, C ), whereas 36% of peaks were located in distal intergenic (potential enhancer) regions and 51% were located in the first (16%) or other (35%) introns ( Figure 3A ). This suggests that the regulation of gene expression by AR in PCOS hESCs in response to elevated androgen levels is primarily via binding at putative enhancer sites, and that there could also be a role in intron retention (19) due to the other major proportion of AR binding being in introns. MOTIF analysis of DNA sequences within AR peaks identified a motif (E=2.7e-268) which mapped to both androgen receptor motif ANDR (E = 1.26e-2) ( Figure 3B ) and progesterone receptor motif PRGR (1.28e-6), and we have previously shown that the progesterone pathway is active in ovulatory PCOS patients (7).

The motif with second highest enrichment (E = 5.1e-89) mapped to the consensus binding site for JUN proteins (E = 2.89e-08) ( Figure 3D ), components of the AP-1 transcription activator complex that are known enhance AR activity and stimulate cellular proliferation in prostate cancer cells (36). ORA analysis of genes with AR peaks in promoter regions using the GO biological processes gene set, identified enrichment of ‘regulation of endothelial cell proliferation’, ‘extracellular matrix organization’ and ‘angiogenesis’ (FDR < 0.3) ( Figure 3E ), pathways known to be associated with the decidualization process (25). Comparing genes with AR peaks in their promoter regions with DEGs between proliferative and secretory phases revealed 83 genes which are up-regulated and 66 genes which are down-regulated, suggesting a specific AR mediated gene set profile that is dysregulated in PCOS ( Supplementary Figure S1 ).

The discovery of an enriched AR motif in the WT1 peak set suggested that WT1 and AR binding sites could be co-located across the genome ( Figure 4A , Supplementary Figure S2 ). We therefore measured the distance between each AR peak to the nearest WT1 peak ( Figure 4B ). We found that 7.0% of AR ChIP peaks directly overlapped with WT1 peaks (P < 2.2e-16) suggesting that the two transcription factors could mutually exclude one another from binding. Merging the WT1 and AR peak sets and keeping only overlapping WT1/AR peaks identified in 826 merged peaks (including multiple overlaps i.e. two WT1 peaks overlapping a single AR peak were considered as a single merged peak). Only 13% of merged peaks were localised in gene promoter regions, corresponding to 106 peaks and 100 genes (considering promoters may contain multiple peaks) ( Figure 4C , Supplementary Table S3 ) including 15 lncRNAs and 10 miRNAs ( Supplementary Table S4 ). The remaining 741 (87%) of WT1/AR merged peaks were located outside gene promoter regions ( Figure 4D ), with 39% located in distal intergenic (potential enhancer) regions ( Figure 4D ).

ORA analysis of shared WT1/AR targets revealed enrichment of GO biological function, including ‘blood vessel endothelial cell proliferation involved in sprouting angiogenesis’, ‘regulation of cytokine production involved in immune response’, ‘regulation of Notch signalling pathway’ and ‘positive regulation of angiogenesis’ ( Figure 4E ). Examination of RNA levels for DEGs with WT1/AR merged peaks in promoter regions using the Sigurgeirsson data set revealed 20 significantly up-regulated and 15 significantly down-regulated DEGs ( Figure 4F ) (34). FKBP5, which had the greatest increase in differential expression, has been shown to regulate decidualization through Ser473 dephosphorylation of AKT (37), and its dysregulation in rats is associated with aberrant PGR-targeted gene expression (38).

STRING database analysis (39) was then used to identify protein/protein interactions of significantly differentially expressed candidate WT1/AR target genes, and using networks visualization we revealed that several of these DEGs were present in the GO biological function pathways identified above ( Figures 4E, G ). This detailed network analysis identified interactions between proteins encoded by STAT3, BCL6, ILR1, HMOX1, dysregulation of which has been implicated in diseases of the female reproductive system (see discussion). Together this reinforces the likely importance of these proteins in the decidualization process, that their dysregulation, due to aberrant AR/WT1 binding in PCOS patients, may affect implantation success.

Finally, we considered the involvement of candidate miRNA and lncRNA WT1/AR targets ( Supplementary Table S4 ) in biological processes and identified pathways upregulated by miRNAs to again include ‘Estrogen signalling pathway’ (P = 0.0015), ‘endometrial cancer’ (P = 0.0012), ‘Wnt signalling pathway’ (P = 0.0195), ‘Focal adhesion’ (P = 0.0002), ‘Adherens junction’ (P = 1.11e-16) and ‘FOXO signalling pathway’ (P = 0.008) ( Supplementary Figure S3 ). Gene products targeted by lncRNAs with merged WT1/AR peaks include HOXA9, IGF1R, WT1, BCL6, FOXO1 and HOXA16 ( Supplementary Figure S4 ).

The proposed functional interplay between WT1 and AR occurs predominantly at distal intergenic regions, and therefore suggests a link to certain epigenetic marks that occur at gene enhancer sites. Enhancer regions are often located many 1000s of bp from gene TSS and are involved in the regulation of these genes through direct interactions with the core promoter via enhancer regions looping over to allow proteins recruited to that site to interact with the core transcription machinery (40). The histone marks H3K4me3 and H3K27ac are purported indicators of promoters and active enhancers respectively (41), with both marks enriched in regions of open chromatin, resulting in active gene transcription. To enable comparative meta-analysis between WT1 and AR binding events and the active histone marks, we exploited publicly available H3K4me3 and H3K27ac ChIP-seq data sets from endometrial stromal cells treated with cAMP and MPA to induce decidualization (GSE61793) (42). Of the 6,676 WT1 peaks located in gene promoter regions, 5,180 were co-located with H3K4me3 peaks (77.6%, P-value < 0.001 one-sided Fishers exact test) (see Figure 5A , Supplementary Table S5 ) indicating that WT1 binding is associated with active gene transcription in decidualization. Performing the same analysis on merged WT1/AR peaks, we found 59 were co-located with H3K4me3 peaks (63.2%, P-value < 0.001 one-sided Fishers exact test), ( Figure 5A , Supplementary Table S5 ), further substantiating that any competitive binding between WT1 and AR in PCOS patients is likely to result in aberrant regulation of genes that are actively transcribed during decidualization. Performing ORA analysis on genes with both H3K4me3 and WT1 peaks within promoter regions identified a similar set of pathways as the previous analysis for genes with WT1 peaks in promoter regions. In contrast only 16.1% of WT1 peaks were located in H3K27ac enriched distal intergenic regions (P-value < 0.001 one-sided Fishers exact test) ( Figure 5B , Supplementary Table S5 ), indicating that WT1 plays a limited function in regulating decidualization via these putative enhancers. From the WT1/AR merged peaks, 165 occupied H3K27ac enriched regions (22.9%, P-value < 0.001 one-sided Fishers exact test) ( Figure 5B , Supplementary Table S4 ), implying that competitive binding of AR in WT1 binding regions affects putative enhancer activity in PCOS patients, resulting in aberrant gene regulation. Finally, we explored the GeneHancer database (43) to understand the biological functions regulated by WT1 and WT1/AR binding in distal intergenic regions and determined that 1,157 WT1 occupied locations (of the 4,437 sites identified in distal intergenic regions) were at defined enhancer locations known to regulate 4,117 unique genes (GeneHancer interaction score > 10), and potentially revealing 3,280 novel enhancer sites. Known enhancer sites corresponded to regulator regions for ‘Notch signalling’ (FDR < 0.01) and ‘FoxO signalling’ (FDR < 0.1) identified using pathway analysis (44) ( Figure 5C ), supporting direct ChiP-seq and comparative RNA-seq meta-analysis. Applying the same analysis to WT1/AR occupied distal intergenic regions, we identified a unique set of 153 annotated enhancers locations (from the 322 overlapping sites identified in distal intergenic regions), regulating 565 genes (GeneHancer interaction score > 10) that identified enriched pathways include ‘response to cAMP’ (P < 0.05) and ‘interferon signalling’ (P <0.001) ( Figure 5D ).

The loss of WT1 in PCOS patients coincides with an increase in the levels of activated AR (7, 15). Here we began to decipher how AR could dysregulate processes in PCOS patients that are regulated by WT1 in fertile women. We identified 106 promoter and 322 enhancers regions that contained overlapping WT1/AR peaks, and infer this is indicative of a competitive transcription factor binding between AR and WT1 within cis-regulatory regions in PCOS patients, governed by the relative abundance of AR within the cells. We identified a small number of putative WT1/AR regulated genes which exhibit significant differential expression between secretory and proliferative phases including FKBP5, which has the highest log₂ fold-change. FKBP5 has been shown to regulate decidualization through Ser473 dephosphorylation of AKT (37), and its dysregulation in rats is associated with aberrant PR-targeted gene expression (38), which coincides with the recruitment to AR to PR sites in PCOS.

Network analysis of WT/AR targets identified interactions between STAT3, BCL6, ILR1, HMOX1, NFKBIA and THBS1; dysregulation of STAT3, BCL6, ILR1 and HMOX1 has been implicated in diseases of the female reproductive system, again demonstrating the validity of our approach. Ablation of STAT3 in murine models results in dysregulation of uterine epithelial remodelling, resulting in aberrant embryonic attachment (51), infertility (52) and implantation failure (53). BCL6 is involved in endometrial cell differentiation, migration, and invasion of trophoblastic cells (54). IL1R1 imbalance in ectopic endometrial tissue of women suffering from / with endometriosis results in heightened sensitivity to IL1 stimulation, which affects their ability to develop into host tissues for implantation (55). HMOX1 dysregulation has been implicated in PCOS (56), endometriosis (57), and in-vitro is down-regulated in endometrial stromal fibroblasts obtained from healthy patients in late secretory phase (58). Finally, PCOS patients exhibit significant differences in their immune cell population distribution and associated cytokine production, and display defective vascular remodelling of the endometrium (59). Our results indicate that WT1/AR crosstalk may affect genes that regulate immune-related and vascular remodelling-related pathways, resulting in reduced fertility in PCOS patients.

Primary hESCs were isolated from endometrial biopsies as previously described from patients of proven fertility and infertile patients diagnosed with anovulatory PCOS (47). Confluent hESC cells were washed twice in PBS and maintained in DMEM/F12 medium supplemented with 10% charcoal stripped media for 24 h prior to the start of the experiment. To induce decidualization, cells cultured in charcoal stripped media were treated with cAMP (0.5 mM), or cAMP (0.5 mM) and medroxyprogesterone acetate (MPA (1 μM) for 48 h. Decidualization was evaluated by analysing changes in cell morphology, and by measuring the secreted levels decidual PRL (dPRL) in the cell culture media, or with DHT (10⁻⁸ M), or DHT (10⁻⁸ M) + cAMP  (0.5 mM) for 48 h to simulate PCOS conditions.

Concentrations of secreted dPRL in cell culture media were measured using commercial ELISA kits according to the manufacturer’s instructions (DY682; R&D systems). Measurements were performed in triplicate.

Chromatin immunoprecipitation (ChIP) was carried out as described previously following treatment with cAMP (0.5 mM) for fertile controls (15, 60). Following treatment, hESC cells were fixed using 1% formaldehyde solution (Sigma^®), quenched with 2.5M Glycine (Sigma^®) and centrifuged following Active Motif’s Epigenetic Services ChIP Cell Fixation Protocol instructions. The pellet was then sent to Active Motif for sequencing. An anti-WT1 antibody (ab89901 Abcam, UK) and anti-androgen Receptor antibody (ab9474, Abcam, UK) was used to probe for WT1 and AR-target region enrichment respectively.

Fastq files were received from Active Motif^®. Adaptor sequences were removed by Active Motif^®. Prior to genome mapping, reads with more than 5 bases with Phred score less than 30 or reads containing undefined bases were removed (61). The sequencing quality of the remaining reads was determined using FASTQC (version 0.11.05, https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) ( Supplementary Figure S5 ). Trimmed sequencing reads were then mapped against the reference genome (hg38, GRCh38) using Bowtie 2 (version 2.2.9, http://bowtie-bio.sourceforge.net/bowtie2/index.shtml) with default parameters (62). Normalized strand coefficient (NSC) and relative strand correlation (RSC), indicators of ChIP-seq experiment signal to noise ratio, were assessed using phantompeakqualtools (63) ( Supplementary Figure S6 ). SAM files were converted to BAM files, sorted and indexed using SAMTools (64). Reads mapping to DAC consensus excluded regions were removed (65). Peak calling for WT1 and AR enrichment above the input control using MACS software (1.4.2) (66) with p-value cut-off 1e-5 (66).

Fastq files of H3K4me3 and H3K27ac ChIP-seq samples on human endometrial tissue samples treated with cAMP were downloaded from the NCBI Geo repository from dataset GSE61793 using the NCBI SRA-toolkit. Quality control and mapping to the genome were performed in the same way as with transcription factor ChIP-seq samples. Peak calling of H3K4me3 samples was performed using MACS2 software, with broad peak option, and qvalue cut-offs of 5e-2 and 1e-1 for narrow peak and broad peak regions respectively. Peak calling for H3K27ac samples was performed using epic2 (version 0.041, https://github.com/biocore-ntnu/epic2) with FDR cut-off of 0.05 (67). Bigwig files were generated from BAM files using Deeptools (version 3.3.1) (68).

P-values of overlapping peaks from WT1, AR and histone mark datasets were calculated using one-sided Fisher’s exact test. For each pairwise comparison between peak sets, the genome was divided into equal length bins of 654 bp long (corresponding to the mean length of WT1 and AR peaks). A two-by-two contingency table was constructed with upper left containing the number of bins without any peaks, lower left containing the number of bins exclusively containing peaks from the first peak set, the upper right bin containing the number of bins exclusively containing peaks from the second peak set and the lower right cell containing the number of bins with peaks from both peak set 1 and peak set 1. P-values were calculated from the contingency table using the fisher.test R function (69).

RNA-seq raw count data, containing sequencing results from seven paired healthy endometrial tissue samples taken in proliferative and secretory phases of the menstrual cycle were downloaded from the NCBI GEO repository GSE86491 (34). Differential analysis was performed comparing mRNA expression in secretory to proliferative phases using DESeq2 (70).

Data analysis was performed using R (version 3.6.3) (69). Overlapping peak between samples were evaluated and merged peak sets were generated using the findOverlapsOfPeaks function in the ChIPPeakAnno (71) Bioconductor R package. Pairs of peaks defined as overlapping contained at least one common genomic coordinate. Peaks were assigned to genomic regions using the ChIPPeakAnno function assignChromosomeRegion. Genes were associated with peaks using the findOverlaps function in the GenomicRanges R Bioconductor package (72), to determine genes within 5000bp of peaks. Gene overlaps were assessed using the intersect function in R. Gene coordinates were obtained from EnsDb.Hsapiens.v86 R Bioconductor package (73). Enrichment heatmaps were created using the EnrichedHeatmap Bioconductor package (74). Genome browser plots were created using GViz R package (75). Venn diagrams were created using the ChIPPeakAnno function makeVennDiagram for overlapping peaks and the R VennDiagram package for overlapping genes. All other graphics were created using ggplot2 (76) in R. Motif analysis was performed using the MEME Suite MEME-ChIP web application (http://meme-suite.org/tools/meme-chip) (77). Multiple Expectation maximization for Motif Elicitation (MEME) (version 5.1.1) and Discriminative Regular Expression Motif Elicitation (DREME) software’s (version 5.1.1) were used for motif discovery, with threshold of E ≤ 0.05. HOCOMOCO Human (v11 CORE) input motif set was used. Over-representation analysis was performed using WebGestaltR R package (78). The ‘BH’ FDR method was used, with an FDR threshold of 0.05. ORA analysis was performed for GO Biological process, GO cellular component (23), and Reactome (79) gene sets. Network analysis to identify protein/protein interactions was performed using the STRING DB plugin (39) for Cytoscape (80). Enhancer locations and enhancer/gene regulatory relationships were obtained from the GeneHancer database (43). Genes and pathways associated with miRNAs were identified using DIANA-mirPath v.3 web application (81). Genes and pathways associated with lncRNAs were identified using LncSEA web application (82).

Gene expression analysis was performed according to the MIQE guidelines (83). Total RNA was isolated from cells following lysis in RLT buffer using the RNeasy Mini kit (Qiagen, Manchester, UK), according to the manufacturer’s instructions. Reverse transcription of 1 μg mRNA was performed in a 20 μl reaction volume using the High Capacity cDNA reverse transcription kit (Thermo Fisher), according to the manufacturer’s instructions. Quantitative PCR primers were designed using the Primer-BLAST primer design software (http://www.ncbi.nlm.nih.gov/tools/primer-blast/) and validated by BLAST analysis against the Homo sapiens (taxid:9606) Refseq mRNA database. Primers and GAPDH reference gene primers were obtained from Sigma–Aldrich: WT1 forward, 5’-CTATTCGCAATCAGGGTTACAG-3’, reverse, 5’-CATGCTTGAATGAGTGGTTGG-3’; GAPDH forward, 5’-GTCCACTGGCGTCTTCAC-3’, reverse, 5’-CTTGAGGCTGTTGTCATACTT-3’. Quantitative PCR was performed in a 10 μl reaction volume comprising 1 × iTaq Universal SYBR Green Supermix (Bio-Rad) with primers added in nuclease-free water to a final concentration of 0.4 mM and 2 μl of cDNA. Thermal cycling parameters were as follows: one cycle of 95°C for 5 min, followed by 40 cycles of 95°C for 10 s and 60°C for 60 s. The relative quantification method was employed to quantify target gene mRNA within samples (84). To generate standard curves, total RNA extracted from cells was reverse transcribed to cDNA, as described. Ten-fold serial dilutions of this reference cDNA were prepared (neat to 1 × 10^-3) in nuclease-free water (Qiagen). For each sample, target and reference gene mRNA abundance was determined from the appropriate standard curve (quantification cycle, Cq). Changes in mRNA abundance between samples were then determined from the ratio of the target gene Cq to the reference gene Cq.

Statistical analyses were performed using IBM SPSS Statistics 22 with biological replicate as the experimental unit. Initially the data were tested for homogeneity, and log or square root transformed if appropriate. Parametric data were analysed by analysis of variance (ANOVA) using Dunnett’s pairwise multiple comparison t-test for individual group comparisons. Data are presented as mean with standard deviation (SD), p < 0.05 was considered statistically significant, and the number of independent experiments is stated in the figure legends.

Statistical analysis of enrichment overlapping gene or peak sets was performed in R using the phyper function (which performs one-sided Fisher’s exact test).