Gene ontology classification and upstream factor analyses reveal enrichment of TGFβ signaling and oxidative stress response in endometriosis

To understand the pathways that were overrepresented among the differentially expressed genes, we performed a gene ontology analysis of all the differentially regulated genes (>1.4, <0.4-fold change, FDR <0.05) in the normal or endometriosis datasets (Supplementary Data File 2). Among the top categories in the normal donors, we found that “VEGFA/VEGFR2 signaling” (Day 8) (Supplementary Fig. 1B) corresponded to increased levels of SH2D2A, encoding an adaptor protein that binds to VEGFR and amplifies its signal⁴³. Other enriched categories were “mRNA Protein and Metabolite Inducation Pathway By Cyclosporine A” and “Amino Acid Metabolism in Triple Negative Cancer”, which represented genes such as SLC1A5 and SLC7A5, whose gene expression decreased during decidualization (Supplementary Fig. 1B and Supplementary Data File 2). These transporters have been previously shown to be dynamically controlled in the endometrium during the peri-implantation period of pregnancy in other animal species^44,45. Enriched pathways in the endometriosis dataset included “Neuroinflammation and Glutamatergic Signaling”, with genes such as TNFRSF1B showing increased expression during decidualization. TNFRSF1B encodes the TNFα receptor, TNFR2; both TNFR1 and TNFR2 have been investigated as potential biomarkers for endometriosis, given the pro-inflammatory roles of TNFα signaling^46,47. The “Nuclear Receptors Meta Pathway,” listed that vitamin D receptor, VDR, as being represented in our dataset of decidualized cells from donors with endometriosis. Our time course decidualization results (Supplementary Data File 1) showed that VDR was consistently decreased during decidualization (EPC Day 8, −2.00 Log FC, FDR=0.0004). VDR signaling is crucial for stromal cell decidualization and early pregnancy progression, suggesting that defects in this pathway may contribute to fertility defects in women with endometriosis⁴⁸. Another overrepresented pathway in our analysis included the “BMP2 WNT4 FOXO1 Pathway In Primary Endometrial Stromal Cell Differentiation,” where we observed that the gene encoding the canonical transcription factor for BMP signaling, SMAD1, was decreased during decidualization. We conclude that this time course decidualization profiling identified several gene categories and signaling pathways that are differentially regulated in decidualizing stromal cells from donors with and without endometriosis.

We performed an upstream transcription factor analysis of the differentially regulated genes in the endometriosis or normal stromal cells to identify master regulatory networks driving the differential transcriptional response (Supplementary Data File 3). By mining the consensus gene targets in the ENCODE and ChEA Transcription Factor Targets dataset^49,50, we identified that genes regulated by the CCAAT Enhancer Binding Protein Beta (CEBP/β) and Transcription Factor 3 (TCF3) were enriched in the normal stromal cells (Fig. 1B). CEBP/β has been shown to be a key factor in endometrial stromal cell decidualization that controls the transcription of the PR^51,52 TCF3 is also shown to control endometrial stromal cell proliferation and decidualization⁵³.

On the other hand, regulation of genes by the NFE2 Like BZIP Transcription Factor 2 (NFE2L2) and SMAD4 transcription factors was highly enriched in the endometriosis dataset (Fig. 1C). NFE2L2, encodes NRF2, and is an important regulator of oxidative stress response that controls the expression of genes that contribute to ferroptosis resistance^54–56. SMAD4 is the downstream activated transcription factor controlling expression of genes downstream of bone morphogenetic proteins (BMPs, through SMAD1 and SMAD5) or TGFβ/activin ligands (through SMAD2 and SMAD3)⁵⁷. Some of these differential responses could also be observed at the gene level (Fig. 1D), as indicated by the expression of IGFBP1⁵⁸, ZBTB16⁵⁹ (decidualization markers), SLC40A1, GPX3, PTGS2^60,61 (markers of oxidative stress and ferroptosis), or LEFTY2, SMAD1⁶² (BMP/SMAD-signaling pathways). To further expand on the expression patterns of the TGFβ signaling family, we visualized the time-dependent expression patterns of all members of this family (Supplementary Fig. 1E).

In summary, we observed that conventional transcriptional programs driving decidualization were found to be overrepresented in the endometrial stromal cells from individuals without endometriosis, indicating a normal decidualization response. However, endometrial stromal cells from individuals with endometriosis displayed an impaired response to oxidative stress and defective BMP/TGFβ signaling pathways.

Identification of perturbed BMP/TGFβ signaling pathways in the decidualizing stromal cells from individuals with endometriosis

We examined the dynamic profiles in the early and late decidual cell transcriptomes and visualized the differentially expressed genes between individuals with and without endometriosis in volcano plots using a >2 or <½ fold-change and FDR <0.05 (Fig. 2A–D). At baseline, we observed that 612 transcripts were differentially expressed between the normal and endometriosis groups, with 261 being upregulated and 351 downregulated in the endometriosis group compared to the normal group (Supplementary Data File 4). Two days after EPC administration, 425 genes were downregulated in the endometriosis samples compared to the normal counterparts. We then observed 364, 207, and 179 genes to be downregulated on Day 4, Day 6, and Day 8, in endometriosis samples compared to the normal samples, respectively. Meanwhile, 335, 263, 147 and 97 genes were upregulated in endometriosis compared to normal at Day 2, Day 4, Day 6, and Day 8, respectively (Supplementary Data File 4).

Open in a separate window

Fig. 2

BMP/TGFβ signaling pathways are defective in the decidualizing stromal cells from individuals with endometriosis.

Differentially expressed genes between the endometrial stromal cells of individuals without (n=4) or with endometriosis (n=3) at each time point during the decidualization treatment were identified and visualized as volcano plots. Differentially expressed genes were determined using a cut off (|log2 fold-change|>1 and FDR<0.05, red denotes increased genes, blue denotes decreased genes, gray indicates no significant change) and displayed following Day 2 (A), Day 4 (B), Day 6 (C), or Day 8 (D) of treatment with the EPC decidualization cocktail (35nM estradiol, 1µM medroxyprogesterone acetate, 50µM cAMP). E Gene ontology analysis of the differentially expressed genes was performed at each time point and visualized as a dot plot. Genes in the TGFβ signaling pathway were identified to be enriched at each of the time points after EPC treatment.

To understand the pathway dependent differences in the endometriosis groups, we implemented Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis in the differentially expressed genes spanning from Day 0 to Day 8 timepoints between normal and endometriosis donors (Fig. 2E, and Supplementary Data File 4). Our enrichment analysis showed distinct enrichment patterns at each of the different time points we assessed. Of notice, the TGFβ signaling pathway was the only category shared in the timepoints during decidualization. Alterations in this category included mostly the ligands of the TGFβ family, such as BMP4, BMP2, LEFTY2, GDF6, INHBA, and LTBP1, in which BMPs typically signal through the SMAD1/5/4 transcription factors while TGFβ and activins signal through the SMAD2/3/4 transcription factors⁶². We also observed persistent downregulation of the classical BMP/SMAD1/5/4 target genes, ID2 and ID3, suggesting that impaired BMP signaling was resulting in decreased transcriptional activation by SMAD1/5/4 in stromal cells derived from individuals with endometriosis.

Interestingly, we observed that while the TGFβ signaling pathway was persistently overrepresented during the time course analysis, the difference was less pronounced at day 6 (Fig. 2E). We observed that fewer genes categorized as members of the TGFβ signaling pathway were present at this timepoint (4 genes at day 6, versus 12 genes at day 2, 9 genes at day 4, and 6 genes at day 8) (Supplementary Data File 4). Three of the genes differentially regulated at day 6 were conserved during all the other timepoints, LRRC32 (0.142-fold, FDR<0.0001), FST (0.109-fold change, FDR<0.001), and BMP6 (0.027-fold change, FDR<0.01). Notably, fibrillin (FBN1), was the only gene to show differential expression between the normal and endometriosis groups at day 6 of decidualization. Compared to decidualized cells from normal donors, FBN1 expression was elevated (2.09-fold, FDR<0.0001) in the decidualized cells of donors with endometriosis (Supplementary Data File 4). Fibrillin is an extracellular glycoprotein that associates with elastic fibers of the extracellular matrix and controls the bioavailability of TGFβ, leading to increased signaling and excessive extracellular matrix deposition⁶³. Fibrillin autoantibodies are commonly found in patients with autoimmune conditions characterized by excessive extracellular matrix deposition, such as systemic sclerosis, connective tissue disease and scleroderma^64–66. Increased anti-fibrillin antibodies have also been found in women with recurrent pregnancy loss⁶⁷, suggesting elevated fibrillin has a pathogenic role in pregnancy progression.

After overlapping the differentially expressed genes from different time points, we identified that 48 genes were consistently down-regulated and that 20 genes were consistently up-regulated regardless of the EPC treatment length (Supplementary Fig. 2A, B). Further enrichment for genes associated with human diseases by DisGeNET⁶⁸ revealed that the 48 consistently down-regulated genes in the individuals with endometriosis were significantly associated with fertility complications such as early pregnancy loss, miscarriage, and spontaneous abortion (Supplementary Fig. 2C, Supplementary Data File 4).

We also observed that genes related to retinoic acid synthesis and metabolism were persistently decreased in the endometriosis group compared to the normal group during decidualization (Supplementary Data File 4). For example, the retinoic acid receptor responder 1 (RARRES1) was significantly decreased in endometriosis across all the time points. Aldehyde dehydrogenase 1 family member B1 (ALDH1B1) was also significantly decreased in endometriosis relative to normal decreased across all time points. Retinoic acid receptor beta (RARB) was decreased on days 2 and 4 of EPC treatment in endometriosis relative to normal. Reprogramming of the endometrium by retinoic acid signaling is critical for endometrial decidualization and early pregnancy^69,70. Furthermore, altered retinoic acid metabolism also affects endometriotic stromal cell decidualization⁷¹. Therefore, our results are in line with previous observations, and support the hypothesis that alterations in retinoic acid metabolism drive fertility defects in individuals with endometriosis.

We also observed that the GATA protein binding 6 (GATA6) was significantly decreased in the endometriosis group relative to normal cells on days 0, 2, 4, 6, and 8 of EPC treatment. GATA6 is a PR direct target gene⁷² and is upregulated in ectopic endometriotic stromal cells from lesions due to methylation defects^73,74. Overall, the time course decidualization analysis between normal and endometriosis donors highlighted various pathways that are differentially expressed between the two groups, providing additional therapeutic or diagnostic opportunities for endometriosis-associated infertility. Furthermore, results from our transcriptomic profiling emphasized the critical roles of the BMP signaling pathways in driving the decidualization processes of the normal endometrium.

Time course analysis identifies expression of key branching points during decidualization

Using time course analysis of decidualization at the single cell level, previous studies identified the trajectory of decidualizing cells over time and showed the presence of branching points, with cells diverging into decidual or senescent decidual cells⁷⁵. In their study, it was shown that non-senescent decidual cells clustered into groups expressing FTL (ferritin light chain) and SCARA5 (scavenger receptor class A member 5), while senescent decidual cells expressed markers of extracellular matrix (ECM) remodeling factors CLU (clusterin), ADAMTS5 (ADAM metallopeptidase with thrombospondin type 1 motif 5), CEMIP (cell migration inducing hyaluronidase 1, also known as KIAA1199), and ABI3BP (ABI family member 3 binding protein). To determine how our time course analysis of decidualization aligned with their results, we correlated the top 50 genes corresponding to the decidual/senescence branchpoint with our set of DEGs from the time course analysis at EPC Day 8 (from Supplementary Data File 1). We identified that genes involved in immune surveillance, such as TIMP3 and CXCL14 were increased in the decidualizing stromal cells from donors with endometriosis relative to cells from normal donors (Supplementary Fig. 1D). These genes are implicated in immune recognition of stressed or senescent cells, suggesting that the expression of these genes rises in decidualizing cells of endometriosis possibly due to decidual senescence. Indicative of a higher state of cellular stress, we also found that the decidualizing cells of donors with endometriosis displayed higher expression of stress response genes, such as CRYAB, HSD11B1 and GLRX (Supplementary Fig. 1D). In the study by Lucas ES et al.⁷⁶, high expression of DIO2 was identified as a key marker in the divergence into decidual cells along with a cluster of cells expressing high CLU, ALDH1A1, ADAMTS5, ABI3BP, and CEMPI. Our analysis showed that the decidualizing cells from patients with endometriosis displayed decreasing expression of DIO2, CLU and ALDH1A1, but elevated expression of ADAMTS5, ABI3BP, CEMPI and COL14A1 (Supplementary Fig. 1D). Thus, our results suggest that decidualizing cells from patients with endometriosis display expression of genes associated with elevated immune surveillance, stress response, and that their transcriptomic signature partially correlates with the signature of decidual senescence. The discordance in some of the markers of senescence could be attributed to the fact our analysis is based on the results of bulk sequencing versus the more sensitive single cell analysis performed by Lucas ES et al.⁷⁶.

Altered BMP signaling impairs decidualization in the endometrium of individuals with endometriosis

Our time course transcriptomic analyses revealed that BMP ligands and SMAD1/5/4 transcriptional targets, ID2 and ID3, were consistently downregulated in stromal cells from individuals with endometriosis during in vitro decidualization (Fig. 2). Previous studies in mouse models indicate that BMP signaling via the SMAD1/5 transcription factors is essential for endometrial receptivity, embryo implantation and decidualization^77–80. To explore the specific roles of BMP signaling pathways in the decidualization defects observed in individuals with endometriosis, we first assessed the activation of BMP signaling pathways at different time points during time course decidualization. As shown in Fig. 3A, in the samples from individuals without endometriosis (n=5), SMAD1/5, the signal transducers of the BMPs, were activated in a time-dependent manner. Upon EPC stimulation, phosphorylated SMAD1/5 (pSMAD1/5) gradually increased along with the EPC treatment length. However, in the samples derived from individuals with endometriosis (n=4), activation of the BMP signaling pathway was impaired, as manifested by the decreased levels of pSMAD1/5 during the time course EPC stimulation (Fig. 3B, C). Gene expression analysis using RT-qPCR demonstrated decidual markers such as BMP2 and IGFBP1 showed increasing levels throughout Day 2 to Day 8 EPC treatment in the stromal cells with blunted induction in the endometriosis cells (Fig. 3D, E). These results indicate that while normal decidualizing stromal cells successfully engage BMP/SMAD1/5 signaling, stromal cells from endometriosis did not induce BMP/SMAD1/5 signaling and failed to decidualize.

Open in a separate window

Fig. 3

Impaired BMP signaling perturbs decidualization in the endometrium of individuals with endometriosis.

A, B Lysates from endometrial stromal cells of individuals without (A, “normal”) with endometriosis (B) after 2, 4, 6, or 8 days of EPC treatment were probed with antibodies to detect phosphorylated SMAD1/5 (pSMAD1/5), total SMAD1, total SMAD5, or GAPDH expression. C Densitometric analysis of pSMAD1/5 in the EPC-treated stromal cells from individuals without (n=5) or with endometriosis (n=4). The different symbols represent individual patient trajectories per sample. One way ANOVA with a Tukey’s posttest. D, E Quantitative reverse transcriptase PCR (qRT-PCR) was used to determine the expression of BMP2 and IGFBP1 in the endometrial stromal cells from individuals without (C, n=4) or with endometriosis (D, n=4). The different symbols represent individual patient trajectories per sample. Data in (D, E) were log transformed and analyzed using a 2-Way ANOVA with a Tukey’s multiple comparison posttest. F, G Time course analysis from the RNAseq analysis comparing the increasing gene expression patterns of ID1 and ID3 in normal and decreasing gene expression pattern in endometriosis stromal cells. EPC, 35nM estradiol, 1µM medroxyprogesterone acetate, 50µM cAMP.

The discrepancy in response to the EPC treatment observed between individuals with and without endometriosis was also identified from the RNA-seq data. Figure 3F, G shows exemplary genes that present contrasting trends during our time course decidualization analysis. Similar to BMP2, the ID1 and ID3 genes were progressively increased over the time course treatment in the samples derived from donors without endometriosis while they were conversely downregulated in the endometriosis cohort. Inhibitor of DNA-binding (ID) genes are not only known downstream BMP responsive genes^81,82 but also are important for endometrial remodeling and decidua formation^83–85. Such inverted trends substantiated the dysfunctional BMP signaling pathways in the endometriosis groups. Our data indicated that in individuals with endometriosis, an impaired BMP signaling pathway is accompanied by dysfunctional endometrial decidualization.

Genome-wide binding of SMAD4 reveals differential binding patterns in the endometrium of individuals with and without endometriosis

Upon ligand binding, canonical BMP signaling pathways use SMAD1/5/4 proteins to initiate transcriptional regulation. SMAD1/5 forms heterodimers and translocate into the nucleus together with common SMAD4⁶². Given that BMP/SMAD1/5/4 signaling is essential for implantation and decidualization^{77–80,85,86}, our goal was to investigate the mechanisms that underpin defective BMP signaling in individuals with endometriosis during decidualization at the transcriptional level. To do so, we utilized the Cleavage Under Targets & Release Using Nuclease (CUT&RUN) method to profile the genome-wide SMAD4 binding sites in the EPC-treated (4 days) endometrial stromal cells derived from both individuals with and without endometriosis. We observed a distinct pattern of SMAD4 binding activities between the two groups (Supplementary Fig. 3A, B). We exemplified the binding activities by showing the Integrative Genomics Viewer (IGV) track view of the ID1 and ID3 loci. We observed that SMAD4 binding was diminished in the endometriosis groups at the ID1 and ID3 loci (Fig. 4A). In total, we identified 2060 peaks showing differences in signal intensity between normal and endometriosis groups (Supplementary Data File 5). Among these, 1190 peaks showed decreased enrichment in the endometriosis group, while 870 peaks had increased enrichment in the endometriosis group.

Open in a separate window

Fig. 4

SMAD4 and H3K27ac CUT&RUN reveals differential binding events in the endometrial stromal cells of individuals with endometriosis.

CUT&RUN was performed for SMAD4 and H3K27ac in endometrial stromal cells from individuals with or without endometriosis induced to decidualize for 4 days to identify differences or similarities in their genome-wide distribution. A Genome track views for the ID1 and ID3 genes displaying the enriched SMAD4 peaks obtained from the normal cells (blue) that are decreased in the endometriosis cells (pink). B Peak annotation of the SMAD4 peaks in 2060 peaks showing differences in signal intensity between normal and endometriosis samples, showing many of the differential peaks (72.29%) were located proximal to the promoter region (within ± 3kb promoter region). C Reactome analysis showing classification of genes that were differentially bound by SMAD4 in the endometriosis samples. Chromatin modifications and signaling by TGFβ family members were in the top three categories. D H3K27ac CUT&RUN was performed in endometrial stromal cells from individuals without endometriosis (“Normal”) or with endometriosis (“E-Osis”) after 4 days of EPC treatment. The heatmap shows the peak signal obtained for H3K27Ac in normal versus endometriosis stromal cells. E Gene ontology classification of the 1122 peaks that were more enriched in the endometriosis samples.

Peak annotation revealed that the majority of these peaks were located within the ±3kb promoter region (72.29%) (Fig. 4B, Supplementary Data File 5). Additionally, we performed the Reactome pathway enrichment⁸⁷ for the genes that were differentially bound by SMAD4 in the endometriosis group (Fig. 4C, Supplementary Data File 5). We found that categories related to ‘signaling by TGFβ family members’, ‘signaling by TGFβ receptor complexes’, and ‘TGFβ activated SMADs’ were enriched. These findings agreed with our transcriptomic results and indicated that a defective TGFβ/SMAD4 signaling program is abnormal in the endometrium of individuals with endometriosis. Because SMAD4 is the “common” SMAD that forms a complex with both BMP-SMAD1/5 and TGFβ/activin-SMAD2/3, our CUT&RUN experiments cannot specifically differentiate between the two pathways or differentiate between SMAD1/5 and SMAD2/3-mediated signaling. However, given that SMAD4 enrichment was decreased at the ID1 and ID3 promoter regions, and this corresponds with decreased ID gene expression, it was likely that BMP signaling pathways also impair SMAD1/5/4 binding activities in the stromal cells from individuals with endometriosis. To identify the potential for a direct transcriptional regulation of SMAD4 on canonical decidualization genes (IGFBP1, PRL, SPP1, FOXO1), we searched for SMAD4 peaks mapping to these genes. While no significant peaks were identified for IGFBP1, PRL or SPP1, we did observe SMAD4 enrichment mapping to the FOXO1 TSS and gene body (Supplementary Fig. 3C). The data suggest that BMP2/SMAD4 signaling directly control the transcriptional activation of FOXO1 during decidualization, but not of other decidual markers. However, this does not rule out that the expression of key decidual markers was affected by a secondary cascade.

Interestingly, apart from TGFβ related categories, pathways involved in ‘chromatin modifying enzymes’, ‘signaling by NTRKs,’ and ‘NGF-stimulated transcription,’ were also among the top enriched categories (Fig. 4C). Neurotrophic tyrosine receptor kinases (NTRKs) are well-documented for their roles in pain and inflammation in endometriosis and are elevated in the endometriotic lesions of affected patients^88–90. Nerve growth factors⁹¹ signal through the NTRKs and were recently shown to be associated with endometriosis through genome-wide association studies⁹². Hence, our studies suggest that abnormal NTRK signaling may also impact the eutopic endometrium and affect receptivity in patients with endometriosis.

To further delineate the chromatin level differences between the normal and endometriosis groups, we profiled the depositions of histone mark H3K27 acetylation (H3K27ac) in the EPC-treated stromal cells (Supplementary Data File 6). H3K27ac modification on the chromatin has been well-defined in the enhancer and promoter regions and is usually accompanied by active transcription activities^93–95. Similar to SMAD4 binding patterns, H3K27ac marks also showed distinct patterns between the normal and endometriosis groups following a 4-day EPC treatment (Fig. 4D, Supplementary Fig. 4A). We identified 1439 peaks that had increased enrichment in the normal group and 1122 peaks that had more enrichment in the endometriosis group (Supplementary Data File 6).

For the genes that preferentially have more H3K27ac peaks in the endometriosis group, Gene Ontology (GO) enrichment analysis indicated positive regulation of cell adhesion being the most enriched category (Supplementary Fig. 4B, Supplementary Data File 6). This includes genes such as vascular cell adhesion molecule 1 (VCAM1), CD44 molecule (CD44), and Wnt family member 5A (WNT5A) (Supplementary Data File 6). Higher levels of VCAM1 are found in ectopic endometriotic lesions and in the eutopic endometrium of individuals with endometriosis and may contribute to disease establishment and progression⁹⁶. CD44 was previously shown to be elevated in the eutopic endometrium of patients with endometriosis and is involved in the attachment and invasion of endometrial cells into the peritoneum^97,98. WNT5A controls endometrial mesenchymal stem cell renewal by activating WNT/β-catenin signaling⁹⁹. Overall, these results corroborate previous studies indicating that the increased cell adhesion abilities in the eutopic endometrium of those with endometriosis facilitate lesion establishment at ectopic sites¹⁰⁰ (Supplementary Fig. 4B, and Supplementary Data File 6).

GO enrichment on the genes that have less H3K27ac peaks in the endometriosis group indicated that categories involving transcription factor binding, extracellular matrix structural constituent and transcription corepressor activity were deficient during decidualization in the endometrium of individuals with endometriosis (Fig. 4E, Supplementary Data File 6). Additionally, genes in the SMAD binding category, such as SMAD3, SMAD6, the SMAD specific E3 ubiquitin protein ligase 2 (SMURF2), and the transforming growth factor beta receptor 1 (TGFBR1/ALK5), had fewer H3K27ac peaks in the endometrial stromal cells from individuals with endometriosis, corroborating our previous transcriptomic results in which a dysfunctional TGFβ/BMP signaling pathway was identified in the decidualizing stromal cells from individuals with endometriosis (Fig. 4E, Supplementary Data File 6). Additional defects in H3K27ac deposition in individuals with endometriosis were observed in the well-known progesterone-responsive genes RARB and CEBPA loci^101,102 (Supplementary Fig. 4C). These results show that defective endometrial transcriptional responses driven by TGFβ and BMP signaling in individuals with endometriosis are detected at the chromatin level, as evidenced by different genome wide SMAD4 and H3K27ac binding patterns in normal versus endometriosis groups.

We also calculated the overlapping numbers of endometriosis-associated DEGs (determined from Day 4 EPC, FC>1.4, <0.4, FDR<0.05), SMAD4, and H3K27Ac bound genes (Supplementary Fig. 4D, E). The genes are visualized as Venn diagrams and grouped by Up- and Down-regulated (Supplementary Fig. 4D) or as a Venn diagram without distinguishing between Up or Down-regulated genes (Supplementary Fig. 4E). Overall, we found that out of the 512 DEGs in the endometriosis dataset at Day 4 EPC, 90 (17.5%) also shared both a SMAD4 and H3K27Ac peak, while 15 (2.9%) shared only a SMAD4 peak, and 282 (55%) shared only an H3K27Ac peak (Supplementary Data File 6).

Silencing of SMAD1 and SMAD5 perturbs endometrial stromal cell decidualization

Previously published studies in mouse models indicate that the BMP/SMAD1/5 signaling pathways are critical for decidualization and endometrial receptivity^{77–80,85,86}. In our present study, we found that the decreased decidualization potential of stromal cells from individuals with endometriosis correlated with defective BMP/SMAD1/5 activation. To functionally examine the role of BMP signaling pathway in mediating decidualization, we perturbed the SMAD1/5 complex using small interfering RNA (siRNA) in endometrial stromal cells from individuals without endometriosis and treated them with EPC to induce in vitro decidualization (Supplementary Fig. 5A). The knockdown effect was validated at the transcript and protein level (Fig. 5A, Supplementary Fig. 5B, Supplementary Data File 7). Upon the knockdown of SMAD1/5, we observed that canonical decidualization markers such as IGFBP1 and WNT4 were significantly downregulated (Fig. 5A). KEGG pathway enrichment on the differentially expressed genes revealed that the TGFβ and FOXO signaling pathways were also enriched within the downregulated group of genes (Fig. 5B). FOXO family plays a critical role in regulating progesterone-dependent differentiation and decidualization¹⁰³ and is indispensable for implantation and decidua formation¹⁰⁴. We highlighted several key genes changes in the heatmap format to visualize the effect of SMAD1/5 knockdown (Fig. 5C).

Open in a separate window

Fig. 5

Knockdown of SMAD1 and SMAD5 perturbs decidualization in endometrial stromal cells.

A Volcano plot showing the expression of differentially expressed genes in siCTL + EPC vs. siSMAD1/5+EPC treated endometrial stromal cells (using a cutoff of Log 2 FC>0.30,<−0.30, FDR<0.05). Blue indicates genes that are down-regulated in siSMAD1/5+EPC vs. siCTL + EPC, red indicates genes that are increased. (n=1 individual without endometriosis). B Dot plot showing the enrichment of genes in key signaling pathways after SMAD1/5 knockdown. C Heatmap showing the expression and functional classification of key genes following SMAD1/5 knockdown + EPC versus siCTL + EPC treatment (n=3 individuals without endometriosis). D The Binding and expression target analysis program was used to integrate SMAD4 binding peaks with the transcriptional changes after SMAD1/5 knockdown in EPC-treated endometrial stromal cells. Dotplot displays the gene ontology classification of genes that were activated by SMAD1/5 (i.e., were downregulated by SMAD1/5 and have a SMAD4 binding site). E Motif analysis was performed on the group genes identified to be SMAD1/5/4 direct targets and displayed as “uptargets” (genes that were increased after SMAD1/5 knockdown and had a SMAD4 peak) or as “downtargets” (genes that were downregulated after SMAD1/5 knockdown and had a SMAD4 peak).

To further map the direct target genes and potential co-factors of SMAD1/5 during decidualization, we used the Binding and expression target analysis program¹⁰⁵ to consolidate our genomic profiling of SMAD4 and the transcriptomic profiling of SMAD1/5 perturbation. Among the direct targets that were activated by SMAD1/5 (which were down regulated upon SMAD1/5 perturbation and were bound by SMAD4, labeled as Down-targets), were the TGFβ signaling pathway and pathways regulating the pluripotency of stem cells (Fig. 5D). We also performed motif analysis on the direct target genes to provide mechanistic insight to the SMAD1/5 mediated gene expression during decidualization. We uncovered potential SMAD1/5 co-repressors such as NFATC2 and T-box family (TBX1/TBX15). NFATC2 is involved in cGMP-PKG signaling pathways and has a role in regulating immune, inflammatory responses¹⁰⁶, it is also reported to be elevated in the thin-endometrium patients who usually have deficient implantation and lower pregnancy rate¹⁰⁷. Acting mainly as repressors, TBX family genes are crucial for embryonic development and tissue differentiation and formation¹⁰⁸. A recent study has shown that TBX15 was elevated in patients with adenomyosis¹⁰⁹. As for transcriptional co-activators, apart from the canonical pan-tissue co-activator NFIC¹¹⁰, the PR motif was enriched in the SMAD1/5 direct target genes, confirming our previous finding that SMAD1/5 may regulate progesterone-responsive genes at the transcriptional level (Fig. 5E).

We also identified two genes, MALAT1 and HDAC4, that contained a SMAD4 binding site and were decreased by SMAD1/5 siRNA knockdown in EPC-treated stromal cells, suggesting that they are direct target genes that are activated by BMP/SMAD1/5/4. The direct SMAD4 binding activities in MALAT1 and HDAC4 loci were visualized as genome track views in Supplementary Fig. 5C, D. MALAT1 and HDAC4 are both involved in facilitating decidualization and in the pathogenesis of endometriosis^111–116. In summary, our studies combine datasets from SMAD1/5 siRNA-mediated knockdowns with SMAD4 binding studies in endometrial stromal cells during decidualization. The results from these not only validated the indispensable roles of SMAD1/5 during human decidualization, but also provided additional layers of regulation in the downstream networks of SMAD1/5 mediated BMP signaling pathways.

Establishment and decidualization of primary endometrial stromal cells

Primary endometrial stromal cells were isolated from surgically resected endometrial biopsies, which were immediately placed in stromal culturing media, DMEM/F12 (Gibco #11330032) supplemented with 10% FBS, 1% Antibiotic-Antimycotic (Gibco #15240062), and 100µg/mL of Primocin (InvivoGen, Cat # MSPP-ANTPM2). Endometrial biopsies were cut into small pieces, digested in Hanks’ Balanced Salt Solution (HBSS) containing 5mg/mL of collagenase (Sigma, #C0130-1G) and 0.2mg/mL of DNase I (Sigma, Cat #DN25-100MG), and then incubated at 37°C for 20min on an orbital shaker at 120rpm. After incubation, the digested tissues were spun down, and the pellets were resuspended in stromal culturing media. The cell suspension was passed through 100µm cell strainers and then 20µm cell strainers. The cell fraction from the flowthrough after the 20µm cell strainers contains the stromal cells and was cultured in stromal culturing media. Stromal cells were passaged once they reached 90% confluency. For decidualization, stromal cells were seeded on 12-well plates at 2×10⁵ cells/well and 10cm dishes at 1×10⁶cells/dish and treated with phenol red-free DMEM/F12 (Gibco #11039021) supplemented with 2% charcoal-stripped FBS and EPC cocktail (1µM MPA, Sigma Cat#1378001-200MG, 0.05mM cAMP, Axxora Cat #JBS-NU-1502-50, and 35nM E2, Sigma Cat #E1024-1G).

Establishment and decidualization of endometrial assembloids

Establishment of endometrial assembloids was performed according to a published protocol^117,134 with minor modifications. In brief, primary endometrial stromal cells and glandular epithelial organoids were established from human endometrial samples. After culturing epithelial organoids and stromal cells separately for two passages, the two were mixed gently at a ratio of 1:2 (v/v) and resuspended in 20 times ice-cold Collagen (Sigma, Cat # C0130-1G). Cells were then aliquoted in 20µl volumes into a 48-well plate and allowed to polymerize at 37°C for 45min, after which the collagen assembloid droplets were overlayed with 500µl of Expansion Medium and maintained in a 37°C cell culture incubator. The collagen droplets were maintained under constant shaking at 90rpm for 48h. (Expansion Medium consists of: Advanced DMEM/F12 (Invitrogen, Cat #12634010), supplemented with 1X N2 supplement (Invitrogen, Cat # 17502048), 1X B27 supplement (Invitrogen, Cat # 12587010), 100μg/ml Primocin (Invivogen, Cat # MSPP-ANTPM2), 2mM L-glutamine (Invitrogen, Cat # 25030-024), 500nM A83-01 (Sigma, Cat #2939), 10% WNT3a, 10% R-Spondin conditioned media, 10mM Nicotinamide (Sigma, Cat # N0636-100G), 1.25mM N-acetyl-L-cysteine (Sigma, Cat # A9165-5G), 10% Noggin conditioned media, 100ng/ml FGF10 (Peprotech, Cat # 100-26), 50ng/ml HGF (Peprotech, Cat #100-39), 50ng/ml EGF (Peprotech, Cat # AF-100-15) and 10nM E2 (Sigma, Cat. #E1024-1G). Conditioned media was produced in HEK293 cells and obtained from the Center for Digestive Diseases and Organoid Core Facility at Baylor College of Medicine. To induce decidualization of the assembloids, the culturing media was changed to decidualization media in minimal differentiation media (MDM) (Advanced DMEM/F12 supplemented with 1X N2 supplement, 1X B27 supplement, 100μg/ml Primocin, 2mM L-glutamine, 1.25mM N-acetyl-L-cysteine, 1µM MPA, 0.5mM cAMP, and 1µM E2) supplemented with or without rhBMP2 (R&D, Cat #355-BM-010/CF) at 25ng/ml. Media was refreshed every 48h.

Histological assessment of endometrial assembloids

Assembloids were fixed in 4% PFA at room temperature for 15min and then immobilized in the Histogel (Thermo Fisher, Cat #22-110-678). After processing the assembloids in Histogel, assembloid blocks were dehydrated through a series of ethanol washes and processed for paraffin embedding at the Human Tissue Acquisition and Pathology Core at Baylor College of Medicine. Paraffin blocks were sectioned using 5µm thick sections. Assembloid sections were deparaffinized in Histoclear and rehydrated in a series of 100%, 95%, 80%, and 70% ethanol washes, followed by washing in dH2O. For identifying the morphological structures, assembloid sections were sectioned and stained with hematoxylin and eosin. For immunofluorescence staining, assembloid sections were heated in boiling 10mM sodium citrate, pH 6.0 for 20min for antigen retrieval and quenched in 3% hydrogen peroxide for 10min. After blocking with 3% BSA for 1h, assembloid sections were incubated with primary and fluorescent secondary antibodies according to the manufacturer’s instructions and nuclei were stained with 1mg/ml DAPI (1:100- dilution, ThermoFisher, Cat # D1306). The stained slides were mounted in VECTASHIELD antifade mounting medium (Vector Laboratories #H-1000-10). Fluorescence images were taken on a Zeiss LSM780 confocal microscope at the Optical and Vital Microscopy Core at Baylor College of Medicine.

RNA extraction and quantitative PCR from endometrial stromal cells and endometrial assembloids

The RNAs of endometrial stromal cells were extracted by using QIAGEN RNeasy micro kit (QIAGEN, Cat #74004) according to the manufacturer’s instruction. The assembloids were lysed in Trizol reagent (Life Tech, Cat #10296010) and the RNAs were extracted by using Direct-zol RNA microprep kit (Zymo Research, Cat #R2062). A total of 50–200ng of RNA from each sample was transcribed into cDNA by using qScript cDNA supermix (Quantabio, Cat #95048-100). Real time quantitative PCR was performed on Bio-Rad CFX384 Touch Real-Time PCR Detection System. Fold changes of target genes were calculated using delta Ct method and normalized GAPDH¹³⁵. Primer sequences are listed in Supplementary Table 1.

Protein extraction and western blotting

Cells were washed with ice-cold 1×DPBS and lysed in M-PER mammalian protein extraction reagent (ThermoFisher, Prod#78505) supplement with protease inhibitor cocktail (ThermoFisher, Cat #78437) and phosphatase inhibitor cocktail (ThermoFisher, Cat #78426). Protein concentrations were quantified by using Pierce BCA protein assay kit (ThermoFisher, Cat #23225). A total of 20µg of protein lysate was loaded onto 4–12% Bis-Tris Plus Mini protein gels (ThermoFisher, Cat # NW04122BOX) and transferred onto nitrocellulose membranes (Bio-Rad, Cat #1704270). The membranes were blocked with 5% non-fat milk in TBST buffer for an hour at room temperature and then incubated with primary antibodies at 4°C overnight. Antibody information is listed in Supplementary Table 2. The next day, the membranes were probed with HRP-conjugated secondary antibodies (Jackson ImmunoResearch) for two hours at room temperature and protein bands were visualized by using SuperSignal West chemiluminescent substrate (Pierce) on Bio-Rad Chemidoc Touch Imaging system. Protein bands were quantified by using Image Lab software (Bio-Rad). All the uncropped blots are included in Supplementary Figs. 7–9.

Statistics and reproducibility

All experiments were performed using >3 biological replicates using technical triplicates in each experiment. Data in graphs are presented as mean +/− standard error of the mean, or as described in each of the figure legends. Statistical tests were performed using t-tests or one-way ANOVA tests with multiple comparison post-tests, as indicated in the figure legend for each experiment. Multiple testing for the large datasets was performed using cut-offs and FDR<0.05. GraphPad Prism was used for statistical analyses, *, 0.033; **, 0.002; ***, <0.001

Gene expression profiling using RNA sequencing

All sequencing data are available in the NCBI Gene Expression Ominibus under SuperSeries GSE243158.

Time course EPC studies

Endometrial stromal cells from 4 normal and 4 endometriosis samples were treated with 35nM estradiol, 1µM medroxyprogesterone acetate and 50µM 8-Br-cyclic AMP for 0, 2, 4, 6, or 8 days. RNA expression profiles were obtained at each timepoint using RNA sequencing analyses (20–30 million paired-end reads using NovaSeq System from Novogene Corporation Inc. Reads were trimmed with fastp v0.23.2 and aligned using STAR 2.7.10a to human genome assembly GRCh38.p13. Differentially expressed genes between normal patients and endometriosis EPC-treated cells were obtained by comparing to the baseline samples (Day 0). Significantly changed genes during the time course treatment were obtained using an ANOVA F-test using an FDR<0.05 from patients without (n=3) and with endometriosis (n=4). All DEGs are presented in Supplementary Data File 1. EnrichR^136–138 was used to identify the gene ontology classifications (Supplementary Data File 2), as well as ENCODE and ChEA consensus transcription factors known to regulate differentially expressed genes in the normal and endometriosis EPC-treated stromal cells (Supplementary Data File 3). Differentially expressed genes between normal and endometriosis stromal cells were obtained by comparing transcripts at each time point of treatment. Differentially expressed genes between normal (n=4) and endometriosis (n=3) were identified using a Wald test with cutoff values of fold-change >2 or <1/2 and FDR<0.05 (Supplementary Data File 4). Gene and pathway enrichment analysis was conducted using R package Cluster Profiler¹³⁹.

SMAD1/5 siRNA studies

Endometrial stromal cells from 3 individuals without endometriosis were treated with either 80nM negative control (siCTL, Horizon Cat #D-001810-10-20) or 40nM of SMAD1 plus 40nM of SMAD5 (siSMAD1/5, Horizon Cat # L-012723-00-0005 & L-015791-00-0005) siRNA followed by 4 days’ EPC treatment. Cells were transfected using Lipofectamine RNAiMAX (LifeTechnologies, Cat #13778500). RNAs were isolated by using QIAGEN RNeasy micro kit and subjected to RNA sequencing analysis to identify differentially regulated transcripts. Samples were normalized through effective library sizes and DESeq2 was used identify differentially expressed genes between the siCTL and siSMAD1/5 EPC-treated cells using an FDR<0.05 (Supplementary Data File 7). Gene and pathway enrichment analysis was conducted using R package Cluster Profiler¹³⁹.

SMAD4 and H3K27 genome-wide binding studies using CUT & RUN

CUT&RUN experiments were performed according to a protocol by Skene and Henikoff¹⁴⁰. After endometrial stromal cells were treated with EPC for 4 days, they were collected by digesting with 0.25% Trypsin (ThermoFisher, Cat #25200056) for 3min. After the digestion, cells were pelleted down at 300×g for 3min and viably frozen down in the freezing medium (90% FBS with 10% DMSO) until experiment day. On the day of the experiment, cell vials were quickly thawed and washed 3 times with washing buffer (20mM HEPES pH 7.5, 150mM NaCl, 0.5mM Spermidine, 1 X Protease Inhibitor). For each reaction, 1.3×10⁶ cells were used for the subsequent Concanavalin A bead binding step. After 10min incubation with Concanavalin A beads, bead-cell complexes were resuspended in 100μl antibody buffer (washing buffer supplemented with 0.01% digitonin, and 2mM EDTA) per reaction. 1μl of IgG antibody (Sigma, Cat #I5006), H3K27ac (Cell Signaling, Cat #8173) and SMAD4 antibody (Abcam, Cat #ab40759) were added to each reaction respectively. After overnight incubation at 4°C, bead-cell complexes were washed twice with 200μl cold dig-washing buffer (washing buffer supplemented with 0.01% digitonin) and resuspended in 50μl cold dig-washing buffer with 1μl pAG-MNase (EpiCypher, Cat #15-1016). After incubation at room temperature for 10min, bead-cell complexes were washed twice with 200μl cold dig-washing buffer and resuspended in 50μl cold dig-washing buffer, then 1μl 100mM CaCl₂ was added to each reaction. The mixture was incubated at 4°C for 2h and the reaction was stopped by adding 50μl stop buffer (340mM NaCl, 20mM EDTA, 4mM EGTA, 0.05% Digitonin, 100ug/mL RNase A, 50mg/mL glycogen, 0.5ng E. coli DNA Spike-in (EpiCypher, Cat #18–1401) and incubated at 37°C for 10min. The supernatant was collected and subjected to DNA purification with phenol-chloroform and ethanol precipitation. Sequencing libraries were prepared using NEBNext Ultra II DNA Library Prep Kit (New England BioLabs, Cat #E7645) following manufacture’s protocol. Paired-end 150bp sequencing was performed on a NEXTSeq550 (Illumina) platform and each sample was targeted for 10 million reads.

Sequencing raw data were de-multiplexed by bcl2fastq v2.20 with fastqc for quality control and then mapped to reference genome hg19 by Bowtie2, with parameters of -end-to-end -very-sensitive -no-mixed -no-discordant -phred33 -I 10 -X 700. For Spike-in mapping, reads were mapped to E. coli genome U00096.3. Spike-in normalization was achieved through multiply primary genome coverage by scale factor (100000 / fragments mapped to E. coli genome). CUT&RUN peaks were called by Model-based Analysis of ChIP-Seq (MACS/2.0.10)¹⁴¹ with the parameters of -f BAMPE -g 1.87e9 -q 0.05 (H3K27ac) or -q 0.1 (SMAD4). Track visualization was done by bedGraphToBigWig, bigwig files were imported to Integrative Genomics Viewer for visualization. For peak annotation, genomic coordinates were annotated by ChIPseeker¹⁴². Differential binding analysis and clustering were conducted using DiffBind¹⁴³. Direct targets motif analysis was conducted through Binding and Expression Target Analysis (BETA)¹⁰⁵ with parameter BETA plus –p –e –k LIM –g hg19 -gs hg19.fa -bl. Gene and pathway enrichment analysis was conducted using R package Cluster Profiler¹³⁹. Annotated peak files were included in Supplementary Data File 5 (SMAD4) and Supplementary Data File 6 (H3K27ac).

We are grateful to Dr. Martin M. Matzuk (M.M.M) for support and guidance on this project. Studies were supported by Eunice Kennedy Shriver National Institute of Child Health and Human Development grants R00-HD096057 (D.M.), R01-HD105800 (D.M.), R01-HD032067 (M.M.M.) and R01-HD110038 (M.M.M). Diana Monsivais, Ph.D. holds a Next Gen Pregnancy Award (NGP10125) from the Burroughs Wellcome Fund.