Introduction

Infertility is an emerging health issue that affects approximately 15% of couples1. One in five women aged 15 to 49 years old with no prior births suffers from infertility in the United States2. One important factor affecting fertility is failed embryo implantation or subsequent post-implantation loss due to endometrial defects. This is evident from the high number of failed pregnancies, with as many as 15% of pregnancies resulting in early pregnancy losses3. Understanding the molecular mechanism of how the maternal endometrium becomes suitable for embryo implantation and eventual decidualization, will be the key to eradicating global concerns related to infertility and early pregnancy losses.

The transforming growth factor β (TGFβ) family plays diverse roles in development, physiology, and pathophysiology4, 5, and in particular, signaling pathways downstream of the bone morphogenetic protein (BMP) subfamily are essential for decidual formation6, 7. There are more than 30 TGFβ family ligands, and these ligands signal through complexes of transmembrane type 1 Activin-Like Kinases (ALK) receptors (ALK1, ALK2, ALK3, ALK6) and transmembrane type 2 receptors (BMPR2, ACVR2A, ACVR2B) and then phosphorylate downstream SMAD1 and SMAD5 proteins. Phosphorylated SMAD1/5 form heteromeric complexes with SMAD4 and translocate to the nucleus to induce specific transcriptional programs. Our laboratory and others have used genetically engineered mouse models with deletions of ligands, receptors, and downstream effectors of BMP signaling pathways to establish that BMP signaling pathways are major regulators of early pregnancy612.

A successful pregnancy begins with reciprocal crosstalk between the maternal endometrium and the new blastocyst during the peri-implantation window. Effective implantation requires precise synchronization between the development of the blastocyst and the transformation of the maternal endometrium into a functional decidua. Endometrial stromal fibroblasts undergo the decidualization process in which they differentiate into unique secretory decidual cells that offer a supportive and immune-privileged microenvironment required for embryo implantation and placental development. Decidualizing stromal cells can react to individual embryos in a way that either supports the implantation and subsequent embryonic development or exerts early rejection13, 14. Aberrant decidualization processes are observed in patients with recurrent pregnancy loss (RPL), displaying a disordered pro-inflammatory response, decreased induction of decidual marker genes, and abnormal responses to embryonic human chorionic gonadotropin13, 15, 16. In addition to affecting early pregnancy outcomes, defective decidualization is also involved in the maternal etiology of preeclampsia causing abnormal placental phenotype17, 18. The process of decidualization is tightly regulated by hormone signaling pathways (estrogen, E2, and progesterone, P4), as well as by BMP signaling pathways. Our recent studies found that endometrial Smad1 deletion had no significant effect on fertility, Smad5 conditional deletion resulted in subfertility, while double Smad1/5 conditional deletion led to infertility due to implantation and decidualization defects9. The uteri of mice with double conditional Smad1/5 deletion also displayed decreased response to P4 during the window of implantation, suggesting synergy between the two pathways. However, the mechanistic genomic actions of SMAD1 and/or SMAD5 in the uterus have not been explored, partly because there are no specific antibodies that distinguish phospho-SMAD1 versus phospho-SMAD5.

In this study, we define how SMAD1/5 instructs the decidualization process using genomic approaches in newly generated transgenic mouse lines. We inspect the potential crosstalk between P4 and BMP signaling pathways mediated by SMAD1/5. Together, our study demonstrates that SMAD1 and SMAD5 exhibit shared and unique genomic binding features and further reveals that SMAD1/5 contributes to the P4 response through transcriptional reprogramming during decidualization.

Materials and methods

Generation of Knock-in Mouse Lines

Smad5PA/PA knock-in (KI) mice were generated using a similar approach as previously described19 . Briefly, single-guide RNA (sg-RNA) was designed to target the regions close to the start codon (Fig 1B) and the sgRNA sequence was inserted into the pX459 V2.0 plasmid (#62988, Addgene). The reference plasmids containing PA tag sequence were constructed in pBluescript II SK (+) vector (Agilent, Palo Alto, CA, USA). One μg of guide RNA inserted vector and 1.0 μg of reference plasmid were co-transfected into EGRG01 embryonic stem (ES) cells. Twelve ES clones out of 48 had the expected knock-in allele.

Mouse models with global HA tagged SMAD1 and PA tagged SMAD5 proteins.

A-B) Schematic approaches for generating Smad1HA/HA and Smad5PA/PA knock-in mouse lines. Sanger sequencing of the genotyping results are included as validation of knock-in sequence. Black and blue boxes indicate untranslated and coding regions, respectively. C-D) Immunoblot (IB) analysis of the immunoprecipitation (IP) of HA tagged SMAD1 and PA tagged SMAD5 proteins from different tissues of the tagged mouse lines. Wild type (WT) mice were used as negative controls. Antibodies used for IB and IP are as labeled. Targeted bands of SMAD1 and SMAD5 are indicated by red arrows.

ES cell clones that possessed the proper KI allele were injected into ICR embryos and chimeric blastocysts were transferred into pseudopregnant females. Chimeric male mice were mated with B6D2F1 female mice to obtain the PA-tagged SMAD5 KI heterozygous mice. Homozygous Smad5PA/PA mice were maintained in the C57BL/6 J × 129S5/SvEvBrd mixed genetic background. To generate Smad1HA/HA mice, Cas9 protein (Thermo Fisher, A36497), sg-RNA, and a repair oligo of homology-directed repair (HDR) containing HA-tag and linker sequences were electroporated into zygotes harvested from in vitro fertilization using B6D2F1 male and female mice. An ECM830 electroporation system (BTX, Holliston MA) was used for electroporation. Subsequently, embryos were cultured overnight to the 2-cell stage and then transferred to the oviducts of pseudopregnant CD-1 mice (Center for Comparative Medicine, Baylor College of Medicine). Pups were further screened for successful heterozygous or homozygous knock-in alleles by PCR using primers spanning across the HA tag.

Sequences of sgRNA, the single-stranded repair oligo for HDR and primer used for genotype are listed in Supplement table 1.

Animal ethics compliance and tissue collection

All mice were housed under standard conditions of a 12-hour light/dark cycle in a vivarium with controlled ambient temperature (70°F ± 2°F and 20-70% relative humidity). All mouse handling and experimental procedures were performed under protocols approved by the Institutional Animal Care and Use Committee of Baylor College of Medicine. All experiments were performed with female mice aged between 7 to 12 weeks a C57BL/6 J × 129S5/SvEvBrd mixed genetic background. All mice were euthanized using isoflurane induction followed by cervical dislocation, and tissues were snap-frozen in liquid nitrogen.

Cleavage Under Targets and Release Using Nuclease (CUT&RUN) Approach

Nuclei from uterine tissues were purified following a previously published protocol20. The experiments were performed using pooled biological replicates from two mice that were processed as technical replicates throughout the CUT&RUN procedure and analysis. In short, uteri were harvested from pregnant mice at 4.5 days post coitus and washed with cold swelling buffer (10 mM Tris-HCl pH 7.5, 2 mM MgCl2, 3 mM CaCl2, 1X Protease Inhibitor Cocktail (PIC, Roche, 11836170001)) immediately after collection. Then tissue was cut into small pieces (∼2-3mm) using scissors, while submerged in cold swelling buffer. Nuclear extract was prepared by dounce homogenization in cold swelling buffer (using a size 7 dounce) and filtered using the cell strainer (100 μm, BD Biosciences). Lysate was centrifuged at 400 g for 10 min, then resuspended in lysis buffer (swelling buffer with 10% glycerol and 1% CA-630, 1X PIC) using an end-cut or wide-bore tips and incubated on ice for 5 min. Nuclei were washed twice with lysis buffer and resuspended in lysis buffer. Next, CUT&RUN procedure largely follows a previous protocol21. Briefly, around 500,000 nuclei were used per reaction. 10 μl of concanavalin-coated beads (Bangs Labs, BP531) were washed twice in Bead Activation Buffer (20 mM HEPES, pH 7.9, 10 mM KCl, 1 mM CaCl2, 1 mM MnCl2) for each reaction. Then, beads were added to nuclei resuspension and incubated for 10 mins at room temperature. After incubation, bead-nuclei complexes were resuspended in 100 μl Antibody Buffer (20 mM HEPES, pH 7.5, 150 mM NaCl, 0.5 mM Spermidine, 1X PIC, 0.01% digitonin, and 2mM EDTA) per reaction. One μg of IgG antibody (Sigma, I5006), HA antibody (EpiCypher, 13-2010) and PA antibody (Fuji Film, NZ-1) were added to each group respectively. After overnight incubation at 4 °C, bead-nuclei complexes were washed twice with 200 μl cold Dig-Wash buffer (20 mM HEPES pH=7.5, 150 mM NaCl, 0.5 mM Spermidine, 1X PIC, 0.01% digitonin) and resuspended in 50 μl cold Dig-Wash buffer with 1 μl pAG-MNase (EpiCypher, 15-1016) per reaction. After incubation at room temperature for 10 min, bead-nuclei complexes were washed twice with 200 μl cold Dig-Wash buffer and resuspended in 50 μl cold Dig-Wash buffer, then 1 μl 100 mM CaCl2 was added to each reaction. The mixture was incubated at 4 °C for 2 hours and the reaction was stopped by adding 50 μl Stop Buffer (340mM NaCl, 20 mM EDTA, 4 mM EGTA, 0.05% Digitonin, 100 ug/mL RNase A, 50 mg/mL glycogen, 0.5 ng E. coli DNA Spike-in (EpiCypher, 18-1401)) and incubate at 37 °C for 10 min. 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, E7645) following manufacture’s protocol. Paired-end 150 bp sequencing was performed on a NEXTSeq550 (Illumina) platform.

Bioinformatic Analysis for CUT&RUN data and re-analysis of published single-cell RNA seq data

For CUT&Run data, raw data were de-multiplexed by bcl2fastq v2.20 with fastqc for quality control. Clean reads were mapped to reference genome mm10 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. Duplicated reads were removed, and only uniquely mapped reads were kept. 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 SECAR22 with the parameters of -norm -stringent -output. Track visualization was done by bedGraphToBigWig23, bigwig files were imported to Integrative Genomics Viewer for visualization. For peak annotation, common peaks were identified with ’mergePeaks’ function in HOMER v4.1124 and then genomic annotation was added by ChIPseeker 25. Motif analysis was conducted through HOMER v4.11 with parameter set as findMotifsGenome.pl mm10 -size 200 -mask24. For single-cell RNA seq, raw data was obtained from EMBL-EBI under accession No. E-MTAB-10287. Cells with low coverage (less than 500 genes detected) were filtered, then gene counts were normalized for each cell by converting counts to quantiles and obtaining the corresponding values from a normal distribution. Then normalized cell vectors are concatenated along the gene panel. Plot visualization was conducted through CELLXGENE platform.26

Western Blot Analysis of Immunoprecipitation (IP-WB)

Tissues were pulverized in liquid nitrogen and then lysed using NETN buffer (20 mM Tris-HCl, pH 8.0, 150 mM NaCl, 0.5 mM EDTA, 10% Glycerol, and 0.5% NP-40). Protein concentration was determined by BCA Protein Assay Kit (Thermo Fisher, 23225). 1.5 mg of total protein lysate was used for IP. IP was performed by adding HA antibody (Cell Signaling, C29F4) or PA antibody (Fuji Film, NZ-1) to the lysate and incubate for 1 hour at 4 °C. Subsequently, protein G magnetic beads (Thermo Fisher, 88847) for an additional 1 hour at 4 °C. Then, the beads were washed five times with NETN buffer, and denatured in sample buffer (Thermo Fisher, NP0007) for further analysis by Western blot. For western blot procedures, briefly, denatured protein lysates were run on the 4 to 12%, Bis-Tris protein gels (Thermo Fisher, NP0321BOX) followed by electrophoretic transfer to nitrocellulose membrane. The membrane then went through blocking by 5% milk in Tris-buffered saline with Tween20 (TBST), followed by incubation overnight at 4 °C in the primary antibodies anti-HA (Cell Signaling, C29F4), anti-PA (Fuji Film, NZ-1), anti-SMAD1 (Life technologies, 385400) and anti-SMAD5 (ProteinTech, 12167-1-AP) at 1:1,000 dilution. The next day, membranes were washed three times with TBST, then incubated with horseradish peroxidase–conjugated secondary antibody for 1 h at room temperature, then washed three times with TBST, developed and imaged on iBright Imaging System (FL1500).

Primary endometrial stromal cells isolation/RNAi/decidualization

Studies using human specimens were conducted as indicated in a protocol approved by the Institutional Review Board at Baylor College of Medicine, H-51900. Human endometrial stromal cells were collected from healthy volunteers’ menstrual effluent as previously reported2729. (N=3) In brief, samples were collected by participants in a DIVA cup during the 4-8 hours on the first night of menses and stored in DMEM/F12 with 10% FBS, antibiotic/antimycotic and 100µg/ml Primocin in a cold insulated pack until processing in the laboratory on the day of collection. The effluent was digested with 5 mg/ml collagenase and 0.2 mg/ml DNase I for 20 min at 37 °C, then cell pellet was collected by centrifuging at 2,500 rpm for 5 min at room temperature. Next, red blood cell lysis was performed by resuspending the cell pellet in 20 ml of 0.2% NaCl for 20 seconds and neutralized with 20ml of 1.6% NaCl. Then the solution was then centrifuged at 2,500 rpm for 5 min. Five ml complete medium (DMEM/F12 supplement with 10% FBS, 1X Antibiotic-Antimycotic + 100 μg/ml Primocin) was used to resuspend the pellet and the solution was passed through 100 µm and 20 µm cell strainer sequentially. The flowthrough containing the stromal cells was centrifuged for 5 min at 2,500 rpm and the pellet was resuspended in 10 ml complete medium and plated in a 10 cm dish. siRNA knockdown was performed using Lipofectamine RNAiMAX following the manufacturer’s protocol. In brief, 0.2 million stromal cells were plated in 12-well plate one day before transfection. On the day of transfection, 2 µl siRNA (20 µM, Dharmacon, D-001810-10, L-012723-00-0005, L-015791-00-0005) and 3 µl Lipofectamine RNAiMAX were diluted in 50 µl Opti-MEM respectively and then mixed to incubate at room temperature for 15 min. Then, the complex was added dropwise onto the cells. 24 h after transfection, medium was changed to DMEM/F12 supplement with 2% charcoal stripped FBS. Decidualization was induced by the addition of 35 nM estradiol (Sigma, E1024), 1 µM medroxyprogesterone (Sigma, 1378001), and 0.05 mM cyclic adenosine monophosphate (Axxora, JBS-NU-1502L) for 4 days with media changes every 48 hours.

RNA Extraction and RT-qPCR

For mRNA extraction from stromal cells, cells were lysed with TRIzol and processed using the DirectZol kit (Zymo, R2051) following manufacturer’s procedures. Approximately 100 ng of mRNA was reverse transcribed into cDNA using iScript cDNA Supermix (Bio-Rad, 1708890) and amplified using specific primers listed in Supplementary Table 1. Primers were amplified using 2X SYBR Green Reagent (Life Technologies, 4364346) using a BioRad CFX384 Touch Real Time PCR Detection System. Data analysis was performed by calculating ΔΔCT value towards GAPDH and then normalized to siCTL. P-value was determined by One-Way ANOVA test using PRISM 9. * P ≤ 0.05, ** P ≤ 0.01, *** P ≤ 0.001, **** P ≤ 0.0001.

Results

Generation of Mouse Models with Global HA-Tagged SMAD1 and PA-Tagged SMAD5 Proteins

Activation of BMP signaling pathways has been established as one of the hallmarks of the decidualization process30, 31. Canonically, SMAD1/5 are regarded as downstream effectors of BMP2 signaling pathways to regulate decidual-specific gene expressions6, 32. However, our recent findings demonstrated that SMAD1/5 can also affect the sensitivity of the endometrium towards E2 and P4 stimulation9. Since we observed phenotypical differences between uterine-specific single SMAD1 and single SMAD5 deletion mice, it is beneficial to delineate the role of SMAD1 and SMAD5 in mediating P4 responses during early pregnancy. We used CRISPR technology to generate genetically engineered knock-in mice with and HA-tagged Smad1 allele (herein called Smad1HA/HA) and PA-tagged Smad5 allele (herein called Smad5PA/PA) as shown in Figure 1A, B. The HA tag and the PA tag were inserted into the N-terminus of the SMAD1 and SMAD5 proteins, respectively. Sanger sequencing was used to confirm genomic insertion (Figure 1A, B). To validate the global detection of tagged proteins, we performed immunoprecipitation followed by western blot analysis on different tissues from Smad1HA/HA and Smad5PA/PA mice. We confirmed the HA and PA antibodies can readily detect HA-tagged SMAD1 and PA-tagged SMAD5 proteins at the predicted size (Figure 1C, D). We also demonstrated the molecular size and expression pattern of HA antibody detected SMAD1 protein was comparable to the SMAD1 antibody detected SMAD1 protein across different tissue types. Similarly, PA antibody showed comparable signal intensity to the SMAD5 antibody in detecting SMAD5 protein across different tissue types (Figure 1C, D). Thus, we successfully generated viable mouse models with global HA-tagged SMAD1 and PA-tagged SMAD5 proteins.

SMAD1 and SMAD5 Exhibit Shared and Unique Genomic Binding Sites During Decidualization

The BMP signaling pathway regulates multiple key events during early pregnancy5, mediated through receptor-regulated SMAD proteins, including SMAD1 and SMAD5. As transducers of the BMP signaling pathway, phosphorylated SMAD1 and SMAD5 form homomeric complexes and then couple with SMAD4 to assemble hetero-oligomeric complexes in the nucleus to execute transcription programs. Our previous studies revealed that conditional ablation of SMAD1 and SMAD5 in the uterus decreased P4 response in during the peri-implantation period, suggesting that the transcriptional of PR depends on BMP/SMAD1/5 signaling9. Furthermore, previous genome-wide PR binding studies show that SMAD1 and SMAD4 binding motifs are enriched in PR binding sites in the uterus33.

To determine the shared and unique transcriptional regulomes of SMAD1 and SMAD5 contributing to the diverse effects of BMP and P4 signaling pathways during decidualization, we first utilized Cleavage Under Targets & Release Using Nuclease (CUT&RUN)21 coupled with next generation sequencing to profile genomic loci bound by SMAD1, SMAD5 and PR from mouse uterine tissues. We performed CUT&RUN on the uterine tissues collected at 4.5 days post coitus (dpc), the time when the fertilized embryo reaches the uterus physically and initiates the decidualization program34 (Figure 2A). After aligning CUT&RUN reads to the mm10 mouse genome, we called peaks using Sparse Enrichment Analysis for CUT&RUN (SEACR)22. To identify high confidence peaks, background noise was normalized to IgG and the stringent criteria for peak calling in SEACR was used. After merging common peaks from two biological replicates, we identified 118,778 peaks for SMAD1 and 166,025 total peaks for SMAD5. We visualized the enrichment of SMAD1 and SMAD5 peaks to the overall aligned chromatin regions as shown in Figure 2B. We found that 7.55% of SMAD1 peaks and 9.53% of SMAD5 peaks were located within the ± 3 kb of the promoter regions (Figure 2C, D). This corresponded to 10,368 genes that were directly bound by SMAD1 at the promoter regions (± 3 kb), whereas 18,270 genes were directly bound by SMAD5 at the promoter regions (± 3 kb). Among these, 4,933 genes were found in common between SMAD1 (47.5%) and SMAD5 (27.0%), while 2,744 and 7,427 genes were found to be uniquely bound by SMAD1 and SMAD5, respectively, providing evidence for the shared and unique functions of SMAD1 and SMAD5 at the transcriptional level (Supplement Figure 1B). To date, only a limited amount of transcription factors have been investigated using the CUT&RUN-seq technique from the tissue samples due to antibody compatibility issue. We recognize that the binding sites and gene number identified here are quite high; however, the high density of binding events was also observed in the ENCODE35 chromatin immunoprecipitation followed by sequencing (ChIP-seq) data for SMAD1 and SMAD5 in the human K562 cells, detecting an average of 63,563 peaks for SMAD1 and 109,682 peaks for SMAD5. (Data accessed through GSE95876 and GSE127365 from Gene Expression Omnibus) Such observations suggest that the SMAD1/5 transcription factors may be dwelling on the chromatin and are poised to drive transcription upon stimulus or following co-factor recruitment as previously shown36. Hence, interpreting how the binding events correlate to biological activity requires comparisons with gene expression profiling in a tissue specific manner.

Genomic profiling of SMAD1 and SMAD5 binding sites during decidualization in vivo.

A) Diagram outlining experimental approaches for tissue collection, processing, and CUT&RUN. B) Heatmaps and summary plots showing the enrichment of SMAD1 and SMAD5 binding peaks from one exemplary replicate. C-D) Feature distribution of the annotated peaks for the SMAD1 (C) binding sites and SMAD5 (D) binding sites.

Identification of Direct Target Genes of SMAD1 and SMAD5 During Early Pregnancy

To pinpoint the direct target genes of SMAD1 and SMAD5, we integrated transcriptomic data from previously published9 SMAD1/5 double conditional knockout mice using progesterone receptor cre (SMAD1/5 cKO) (GSE152675) with SMAD1 and SMAD5 genomic data from this paper. We cross-compared the differentially expressed genes in the transcriptomic data to the SMAD1 and SMAD5 bound genes, respectively. Among the 805 significantly up-regulated genes, we identified 449 genes that were both significantly up-regulated upon SMAD1/5 depletion and were directly bound by SMAD1 and SMAD5, whereas 187 of the up-regulated genes were bound by SMAD5 only and 30 were bound only by SMAD1. (Figure 3A) Among the 683 significantly down-regulated genes, we identified 523 genes that were both significantly down-regulated upon SMAD1/5 depletion and were directly bound by SMAD1 and SMAD5, whereas 83 of the down-regulated genes were bound by SMAD5 only and 13 were bound by SMAD1 only. (Figure 3B, Supplement Table 2)

SMAD1 and SMAD5 show unique direct target genes during early pregnancy.

A-B) Venn diagrams showing the shared and unique direct up-target genes (A) and down-target genes (B) of SMAD1, SMAD5 Numbers indicate genes numbers. C-D) Motif enrichment analysis from the up-targets and down-targets for SMAD1 (C) and SMAD5 (D). E-F) Dot plot showing Gene Ontology enrichment analysis of shared direct target genes of SMAD1/5 from the up-targets (E) and the down-targets (F), respectively. Dot size represents the gene ratio in the enriched categories compared to background genes, dot colors reflect P-value.

Next, we utilized Binding and Expression Target Analysis (BETA) algorithm37 to perform motif enrichment analysis of the direct target genes to identify putative co-factors working together with SMAD1 and SMAD5 in controlling gene expression (Figure 3 C,D). “Up-targets” represent genes that were up-regulated in the SMAD1/5 cKO mouse uteri and showed either a SMAD1 or a SMAD5 binding site in the genomic profiling data. Similarly, “down-targets” represent genes that were down-regulated in the SMAD1/5 cKO mouse uteri and displayed either a SMAD1 or a SMAD5 binding site. Thus, motifs enriched in the “up-targets” indicate potential repressive SMAD1/5 co-factors while motifs enriched in the “down-targets” indicate potential SMAD1/5 co-activators. Among the “up-targets” of SMAD1, MYB Proto-Oncogene (Myb)/MYB Proto-Oncogene Like 1(Mybl1) motif was the most highly enriched with a P-value of 1.85E-02. Myb and Mybl1 transcription factors belong to MYB gene family, which has been well-defined in controlling cell survival, proliferation and differentiation in cancer38. In addition, they have also been reported to be E2 induced in human uterine leiomyoma samples39. Homeobox containing 1 (Hmbox1) and Krüppel-like factor (Klf) family members (Klf4/Klf1/Klf12) were also identified as potential repressive co-factors of SMAD1 with P-values of 2.85E-02 and 3.75E-02 respectively. (Figure 2C) Of note, KLF4 has been reported to inhibit the binding activity of estrogen receptor α (Erα) to estrogen response elements in promoter regions40. Among the “up-targets” of SMAD5, EBF Transcription Factor 1 (Ebf1) motif was the most enriched with a P-Value of 1.57E-02. Interestingly, Ebf1 can directly repress the transcription of Forkhead box protein O1 (Foxo1)41. It is also recognized as downstream effectors of steroid hormone receptors in the mouse uterus42. Additionally, motifs from transcription factors Zfp128 and Otx1 were also significantly enriched in the up-regulated genes bound by SMAD5 (Figure 2D). Taken together, our enrichment analysis provided robust evidence for identifying novel co-factors of SMAD1/5, and such co-regulating mechanisms are in line with the unopposed E2 response observed in the SMAD1/5 cKO mice9. Furthermore, odd-skipped– related genes (Osr1 and Osr2) were identified as potential co-activators for SMAD1. Osr2 has been reported to be highly expressed in the human endometrium43, and it was also abundantly detected at the protein level in the human decidual tissues44. Decreased OSR2 level was observed in the patients with recurrent spontaneous abortion and knockdown of OSR2 impairs the decidualization process in the human endometrial stromal cells44. Moreover, OSR1 has been reported to suppress BMP4 expression, which in turn reduced the Wnt/β-catenin signaling pathways during lung development in xenopus45. Apart from Osr family, motifs in the Homeobox genes (HOX) were found to be enriched in the “down-targets” from both SMAD1 and SMAD5 datasets. Specifically, Hoxa11/Hoxd12/Hoxc10 were predicted to be co-activators for SMAD1 while Hoxd10 was indicated to be closely interacting with SMAD5. Indeed, HOX genes are critical for endometrial development in normal and disease conditions and are essential during the establishment of pregnancy4649.

With direct targets genes of SMAD1 and SMAD5 identified, we then analyzed the Gene Ontology enrichment for the SMAD1/5 shared up-targets and down-targets, respectively. We found that “up-targets” genes exhibit enrichment for regulation of cell-cell adhesion, cell junction organization and desmosomes organization (Figure 3E). Moreover, among the “down-targets” genes, we found the enrichment for blood vessel / vasculature development and extracellular matrix organization categorizes (Figure 3F). Indeed, during early pregnancy, the stimulation from corpus-luteum derived P4 enabled the endometrium to be transformed to a receptive state, which allows subsequent embryo attachment and develop through the epithelium into the stromal sections30. During this process, apportioned direct cell-cell contacts are ensured by tight and adherent junctions and such interactions are key in facilitating implantation and embryo invasion. In accordance with our findings, desmosomes and adherens junctions were extensively described to decline in the early pregnancy period, which facilitates the invasion of trophoblast through the epithelial layer5053. In addition, the stromal compartment of the endometrium also undergoes profound vascular remodeling. Precise regulations of the angiogenesis are required to establish extensive vascular network, which is essential to ensure blood supply and successful embryonic development54, 55. Collectively, our findings present evidence that emphasizes the shared roles of SMAD1 and SMAD5 in facilitating the endometrial transitions during early pregnancy.

Direct Target Genes of SMAD1 and SMAD5 Maintain the Homeostasis of Uterine Function

To discover novel direct target genes of SMAD1/5, we visualized keys genes of interest from the up-targets and down-targets. As shown in Figure 4A, data from RNA-seq represents the decrease of several "down-targets” in the SMAD1/5 cKO mouse uteri, including Retinoic Acid-Related Orphan Receptor B (Rorb), Follistatin (Fst), Lymphoid Enhancer Binding Factor 1 (Lef1), and Insulin Like Growth Factor 1 (Igf1). Integrative Genomics Viewer (IGV) track view shows the exemplary SMAD1/5 binding activities near the promoter regions of Rorb and Fst Figure 4B, demonstrating that these genes are bona fide direct target genes of SMAD1/5. Rorb belongs to the nuclear receptor families in the retinoic acid (RA) signaling pathways56 and is considered as a marker for mesenchymal progenitor cells in the stroma compartment of the endometrium57. In murine models, deficient RA signaling through the perturbation of RA receptor in the uterus leads to implantation and decidualization failure58. Fst binds several TGFβ family ligands and thereby inhibits TGFβ family signaling extracellularly59. Under physiological conditions, Fst is up-regulated in the decidua during early pregnancy. Conditional deletion of Fst in the mouse uterus results in severe subfertility with a phenotype of non-receptive epithelium and poor-differentiated stroma60. Notably, RA signaling deficiency also decreases Fst levels in the uterus and systematically administration of FST can fully rescue the deficient-decidualization phenotype but not the non-receptive phenotype observed in the RA receptor mutant mice58. Our results suggest a direct relationship between BMP and RA signaling pathway, accomplished by SMAD1/5 at the transcriptional level, likely establishing a positive signaling feedback loop. Apart from being a crucial transcriptional activator, SMAD1/5 also plays a role in repressing key gene expression pathways. Shown in Figure 4C, upon the deletion of SMAD1/5 in the mouse uteri, several E2-responsive genes were significantly up-regulated, including Fibroblast Growth Factor Receptor 2 (Fgfr2), Matrix Metallopeptidase 7 (Mmp7) and Wnt Family Member 7B (Wnt7b). In addition, Inhbb, a downstream target of Fst60, is also a target gene of SMAD1/5 that resulted in transcriptional repression. SMAD1/5 binding on the Fgfr2 and Mmp7 genes are exemplified in an IGV track view in Figure 4D. Fgfr2 and its ligands regulate epithelial cell proliferation and differentiation. Components of the Fibroblast Growth Factor (Fgf) signaling pathway are cyclically expressed in the uterus and act as paracrine and/or autocrine mediators of epithelial-stromal interactions61, 62. During early pregnancy in mice, P4 inhibits expression of Fgf2 in the stromal cells, which is critical to counteract the E2-driven epithelial proliferation61. Similar observations are reported in gilts, where the expression of Fgfr2 decreased alongside with increased parity of the sows63. It is also noteworthy that loss of function of Fgfr2 in the mouse uterus leads to luminal epithelial stratification and peri-implantation pregnancy loss62. Moreover, Mmp7 and Wnt7b are up-regulated upon E2 stimulation and participate in the re-epithelialization of the endometrium and implantation process, respectively6466. In accordance with the phenotype of hyperproliferative endometrial epithelium during early pregnancy observed SMAD1/5 cKO mice, we demonstrated that the suppression of key E2-responsive genes, such as Fgfr2 and Mmp7, by SMAD1/5 maintains the precise balance between E2 and P4.

Direct target genes of SMAD1/5 mediate uterine homeostasis.

A) Histogram of normalized Fragments Per Kilobase of transcript per Million mapped reads (FPKM) of downregulated transcripts in the Control and SMAD1/5 cKO groups as indicated by the label. Histograms represent average +/- SEM of experiments uteri from Control mice (N=3) and SMAD1/5 cKO mice (N=4). Analyzed by a unpaired t-test. B) Integrative Genomics Viewer (IGV) track view of SMAD1, SMAD5 binding activities. Gene loci are as indicated in the figure, genomic coordinates are annotated in mm10. C) Histogram of FPKM of up regulated transcripts in the Control and SMAD1/5 cKO groups as indicated by the label. D) IGV track view of SMAD1, SMAD5 binding activities. Gene loci are as indicated in the figure, genomic coordinates are annotated in mm10. E) Dot plot showing the gene expression pattern of the key SMAD1/5 direct target genes in different cell types from published human endometrium single-cell RNA-seq dataset.

To explore the major cell types regulated by SMAD1/5 direct targets in human, we profiled the expression levels of the key “up-targets” and “down-targets” in the different cell types of the human endometrium. Using previously published single-cell RNA seq data of human endometrium67, we visualized the expression patterns of suppressive targets and activating targets of SMAD1/5. Apart from the major epithelial and stromal compartments, SMAD1/5 target genes are also widely expressed in the immune cell populations. Such observation reinforced the importance of the BMP signaling pathways in establishing an immune privileged environment at the maternal-fetal interface68.

SMAD1 and SMAD5 Co-regulate PR Target Genes

SMAD1/5 cKO mice were infertile due to endometrial defects and displayed decreased P4 response during the peri-implantation period9. Hence, we hypothesized that SMAD1 and SMAD5 act as co-regulators of P4-responsive genes during the window of implantation and are required for endometrial receptivity and decidualization. By determining the genomic co-occupancy of SMAD1, SMAD5 and PR, we aimed to clarify the transcriptional interplay between the BMP and P4 signaling pathways. To this end, we performed additional PR CUT&RUN experiments on the uteri of mice collected at 4.5 dpc and identified 134,737 peaks showing PR binding activities (Figure 5A). We identified 7,393 genes that were directly bound by PR at the promoter regions (± 3 kb), among which, 2596 genes were also concurrently bound by both SMAD1 and SMAD5 at the promoter regions (± 3 kb) (Supplement Figure 1B).

SMAD1 and SMAD5 co-regulate PR target genes.

A) Heatmaps and summary plots showing the enrichment of PR binding peaks from one exemplary replicate. B) Dot plot showing KEGG pathway enrichment analysis for shared genes bound by SMAD1, SMAD5, and PR. C) IGV track view of SMAD1, SMAD5 and PR binding activities. Gene loci are as indicated in the figure, genomic coordinates are annotated in mm10. D) Table of motif analysis results for shared peaks between SMAD1, SMAD5 and PR, with P-value and motif annotation specified for each motif.

Next, we performed KEGG pathway enrichment for the genes co-bound by SMAD1, SMAD5 and PR. As expected, pathways critical for decidualization such as relaxin signaling pathways and WNT signaling cascade were identified in the enrichment results (Figure 5B). We visualized exemplary genes co-regulated by SMAD1, SMAD5 and PR and presented in the normalized IGV track view. (Figure 5C) We demonstrated SMAD1, SMAD5 and PR showed co-occupancy at the loci of the SRY-Box Transcription Factor 17 (Sox17), Inhibitor of DNA binding 2 (Id2), Forkhead box protein O1 (Foxo1), Insulin-like growth factor 1 (Igf1), Transforming growth factor beta receptor 2 (Tgfbr2) and RUNX family transcription factor 1 (Runx1) (Figure 5C). Sox17 has been reported as one of the direct target genes of PR33 and is essential for uterine functions during implantation and early pregnancy69, 70. More recent studies also showed the importance of Sox17 in regulating uterine epithelial–stromal crosstalk and its indispensable role in female fertility71. We provided evidence that Sox17 is also directly regulated by SMAD1/5 complexes. Our results indicated that Id2, considered as canonical direct transcriptional targets of BMP-SMAD signaling72, 73 is also regulated by PR. We also confirmed that known P4-responsive genes such as Tgfbr274 and Runx175, as well as decidual markers such as Foxo176 and Igf177, were co-regulated by SMAD1, SMAD5 and PR (Figure 5C).

To identify additional transcription factors that are associated with the regulatory interplay between SMAD1/5 and PR during decidualization, we performed unbiased motif analysis on the shared CUT&RUN peaks between SMAD1/5 and PR. We reported the top 10 transcription factors harboring the enriched motifs, including NANOG, Homeobox A protein family (HOXA11 and HOXA9), NK6 homeobox 1(NKX6.1), TGFB induced factor homeobox 2 (TGIF2), FOS, RUNX family transcription factor 2 (RUNX2), Androgen receptor (AR), Sox17 and Lymphoid enhancer-binding factor 1 (LEF1) (Figure 5D). Many of these putative interactors have been reported to interact with the SMAD proteins in other biological process. For example, NANOG interacts with SMAD1 during mesoderm differentiation78.

HOXA9 forms heterodimers with SMAD4, leading to BMP-driven initiation of transcription from the mouse Opn promoter in vitro79, 80. Transcription factor AP-1 family (FOS) and RUNX2, as well as β-catenin/Lef1 complex, increase the effectiveness and specificity of DNA binding activities of SMAD1/5 in response to BMP ligand stimuli8183. Overall, our analyses demonstrate that the transcriptional activity of SMAD1, SMAD5 and PR coordinate the expression of key genes required for endometrial receptivity and decidualization.

Decidualization of Human Endometrial Stromal Cells Requires SMAD1/SMAD5

We next sought to functionally characterize the role of SMAD1/5 during decidualization in human endometrial stromal cells. To do so, we examined the effect of SMAD1/5 perturbations on the decidualization of primary human endometrial stromal cells (EnSCs). EnSCs were transfected with short interfering RNAs (siRNAs) targeting each gene (SMAD1 and SMAD5) and subjected to in vitro decidualization by treatment with E2-cAMP-and MPA (EPC) for 4.5 days (Figure 6A). We hypothesized that the combined SMAD1/5 knockdown would impair the decidualization process significantly compared to cells treated with non-targeting siRNAs. Our results demonstrated that SMAD1/5 knockdown affected decidualization and led to significantly decreased expression of the canonical decidual markers, PRL and IGFBP1 in EnSCs (Figure 6B). The PR co-regulator, FOXO184, exhibited a decreasing trend in the siSMAD1/5 group although with a P-value of 0.07 due to variance derived from different individual samples. We also examined the expression level of the RA pathway regulator gene, RORB, and of the SMAD4-PR target gene, KLF1512, following SMAD1/5 perturbation. We observed a significant decrease in both RORB and KLF15 expression upon SMAD1/5 knockdown during in vitro decidualization treatment (Figure 6C). Taken together, our findings indicate SMAD1/5 can modulate PR activity during decidualization and that this transcriptional cooperation is required for the in vitro decidualization of primary human endometrial stromal cells.

SMAD1 and SMAD5 are required for PR responses during decidualization of human endometrial stromal cells.

A) Schematic approach and timeline outlining in vitro decidualization for endometrial stromal cells (EnSCs). B-C) RT-qPCR results showing mRNA levels of PRL, IGFBP1, FOXO1, RORB and KLF15 after SMAD1/5 perturbation using siRNAs. Data are normalized to siCTL-Veh for visualization. Histograms represent average +/- SEM of experiments on cells from three different individuals with technical triplicates. Analyzed by a One-Way ANOVA test.

Discussion

SMAD proteins are canonical transcription factors that are activated in response to TGFβ family signaling and mediate the biological effects of these pathophysiologically critical ligands82. While SMAD2 and SMAD3 are downstream of TGFβs, activins, and multiple other family ligands, SMAD1 and SMAD5 preferentially transduce BMP signaling pathways and are regarded as pivotal activators for many physiological processes, including bone development, cardiac conduction system development, and embryonic pattern specification8587. Importantly, SMAD1 and SMAD5 are implicated in diverse female reproductive physiology and pathophysiology processes5, 9, 8890

Due to high structural similarity, SMAD1/5 have been suggested to be redundant from the studies in ovarian biology and chondrogenesis89, 91. However, other studies clearly demonstrated that SMAD1/5 have different roles in governing hematopoiesis and uterine functions9, 92. The DNA binding activities of SMAD1 and SMAD5 have not been readily distinguished between each other due to anti-phospho antibody limitations. To robustly define the roles of SMAD1/5 in regulating transcriptional programs in vivo, we produced two genetically engineered mouse models with global knock-in of an HA tag and a PA tag in the Smad1 and Smad5 loci, respectively. We showed that SMAD1 and SMAD5 not only have shared transcriptional activities but also have unique roles in uterine physiology. In agreement with previous studies showing that SMAD1/5 function is partially redundant89, 91, we confirmed that SMAD1/5 share a total of 972 direct target genes in the uterus. Furthermore, we demonstrated that 43 genes were uniquely regulated by SMAD1 whereas 270 genes are specifically regulated by SMAD5 only. Our motif analysis also revealed distinct potential co-factors between SMAD1 and SMAD5, providing evidence at the molecular level to mechanistically delineating the distinct roles of SMAD1 and SMAD5 in directing cellular processes in the uterus.

Apart from directly regulating target gene expression, our data demonstrate that SMAD1/5 present as dense genomic occupancies. Multiple aspects can contribute to this observation. First, transcription factors (TFs) tend to dwell or “search and bind” throughout the genome36. Such events may not yield actual biological effects but rather are due to differences in motif binding affinities93. Second, apart from robust binding activities, TFs may not initiate transcription programs owing to the lack of co-factors or favorable conditions to exert their functions94. Additionally, TF binding sites and target genes are unlikely a one-to-one relationship. TFs could be positioned from the proximal promoter regions to hundreds of kilobases afar to modulate gene expression. In the meantime, the same binding site could regulate multiple genes by interacting with different promoters in different subpopulations of cells. Lastly, TFs usually direct target genes expression in a cell-type specific manner95. Our genomic profiling samples were collected from whole uterus at the time of 4.5 dpc, containing a great range of cell populations, including but not limited to the epithelium (luminal and glandular), stroma (progenitors and differentiated cells), myometrium, endothelium, and immune cell populations. The data is therefore expected to depict the dynamic and complex activities of SMAD1/5 in the entire uterus. Together, the stringent filtering and normalization criteria, comparable peak number to the published dataset and IGV track view visualization collectively validate our CUT&RUN experiments and uncover the enriched regions as robust SMAD1/5 binding events.

Although our studies herein confirm that SMAD1 and SMAD5 proteins have distinct transcriptional regulatory activities, our previous studies demonstrated that while SMAD5 can functionally replace SMAD1, SMAD1 cannot replace SMAD5 in the uterus9. How this epistatic relationship is established in the tissue-specific manner still needs to be determined by further biochemical investigations. In addition, further studies are needed to uncover whether SMAD1 and SMAD5 response differently upon ligand stimulation in the uterus, and if so, how the preference is achieved. Our study provides versatile in vivo genetic tools for these questions and can advance the toolbox for the field studying BMP signaling pathways. Because our mouse models are global knock-in mice, they will not only serve as a powerful tool for studying BMP signaling pathways in the reproductive system but will also promote the study of BMP signaling in other organs and tissues.

BMP signaling pathways are involved in a plethora of cellular processes and appropriate functioning of the BMP pathway depends on the precise crosstalk with other signaling pathways. Coordinated communication with other pathways can yield synergistic effects and leads to a complex regulatory network of biological processes. To be specific, SMAD1/5 mediates the crosstalk with the WNT/β-catenin pathway. WNT signaling inhibits glycogen synthase kinase 3β (GSK3β) activity and prevents SMAD1 from degradation which governs the embryonic pattern formation96. Also, SMAD1/5 can physically interact with T-cell factor (TCF) or lymphoid enhancer factor (LEF) transcription factors to form transcriptional complexes to activate the transcription of many WNT-and BMP-responsive genes97. In addition, SMAD1 and SMAD5 can directly associate with Notch intracellular domain and enhance known Notch target gene expression by binding to their regulatory DNA sequences98. Intriguingly, in prostate cells, SMAD1 physically interacts with the androgen receptor (AR) and halts the androgen-stimulated prostate cell growth99. Moreover, we provide first-hand evidence showing that BMP signaling pathways converge with RA signaling pathways through the regulation of RORB by SMAD1/5. Further studies will grant a more detailed mechanism of the positive feedback loop between BMP and RA signaling.

Our previous studies suggest that the mouse endometrium presents decreased P4 responsiveness following conditional deletion of SMAD1/5 in the uterus9. In accordance with the phenotypical observation, we offer compelling support in our current study that SMAD1/5 work collectively with PR to regulate their target genes and that SMAD1/5 mediate the crosstalk between BMP and P4 signaling pathways during decidualization, a key process to ensure a successful pregnancy, and ultimately direct the biological transformations of the uterus during early pregnancy. We provide genomic evidence that SMAD1/5 are co-bound at around 35% of PR target genes in the mouse uterus during decidualization. We also identified nuclear receptor motifs (i.e., PR sequence motifs) enriched in the SMAD1/5 binding sites (Supplement Figure 1C,D). Correspondingly, in a previously published study where they performed PR ChIP-seq in the mouse uterus after P4 stimulation, the SMAD1 motif was the 5th most significantly enriched sequence motifs identified33.

SMADs are known to recruit co-repressors (i.e., Ski100) or co-activators (i.e., p300101) to inhibit or activate target gene transcription, less is known about their cell-specific co-factors that confer the precise spatial-temporal control over binding activities to target genes. Our study highlights the potential co-factors by integrating both genomic and transcriptomic data to delineate signaling crosstalk that are responsible for maintaining tissue homeostasis, especially in the female reproductive tract.

In summary, our findings and those of others indicate that SMAD1 and SMAD5 not only are signal transducers for BMP signaling pathways, but also engage extensively in the crosstalk with PR signaling pathways. While P4 responses are critical for early pregnancy establishment, abnormal P4 responses are implicated in diseases such as endometriosis and endometrial cancers102105. Hence, our results which show that BMP and P4 signaling pathways synergize within the endometrium; these key pathways can shed light on the endometrial contribution to conditions that impact reproductive health in women, including early pregnancy loss, endometriosis, and endometrial cancer. Furthermore, we anticipate that the SMAD1/5 knock-in tagged transgenic mouse models developed herein will be useful for studying BMP/SMAD1/5 signaling pathways in other reproductive and non-reproductive tract tissues in the body.

Data availability

Sequencing data and analyses are deposited in the Gene Expression Omnibus under accession number GSE237975. (Reviewer token: gludeooirhczxcj)

Supplementary data

Supplement Table 1: DNA sequences used in the study.

Supplement Table 2: Direct target genes of SMAD1/5.

Supplementary Figures: Supplementary Figures 1-3

Author contributions

Zian Liao: Conceptualization, Formal analysis, Methodology, Validation, Writing—original draft. Suni Tang: Data curation, Investigation, Methodology, Visualization, Writing—review & editing. Kaori Nozawa: Data curation, Investigation, Methodology, Writing—review & editing. Keisuke Shimada: Data curation, Methodology, Writing—review & editing. Masahito Ikawa: Project administration, Resources. Diana Monsivais: Conceptualization, Investigation, Supervision, Funding acquisition, Writing – review & editing. Martin M. Matzuk: Conceptualization, Supervision, Funding acquisition, Writing – review & editing.

Acknowledgements

This research is supported by Eunice Kennedy Shriver National Institute of Child Health and Human Development grants HD105800 (D.M.), HD096057 (D.M.), HD032067 (M.M.M.), and HD110038 (M.M.M.). DM is supported by a Next Gen Pregnancy Award from the Burroughs Wellcome Fund (NGP10125).

Conflict of interest

No conflict of interest is declared by the authors.

REFERENCE

Table and figures legends

Gene numbers with SMAD1/5 promoter binding activities and motif analysis of SMAD1/5 peaks.

A) Venn diagrams showing the shared and unique genes bound by SMAD1 or SMAD5 in the +/- 3kb region of the promoter regions. B) Venn diagrams showing the shared and unique genes bound by SMAD1, SMAD5 or PR in the +/- 3kb of the promoter regions. C-D) Table of motif analysis results for unique peaks for SMAD1(C) and SMAD5 (D), with P-value and motif annotation specified for each motif.

Knockdown effect validation of SMAD1/5 perturbation.

A-B) RT-qPCR results showing mRNA levels of SMAD1(A) and SMAD5 (B) after siRNA treatments in the both Veh and EPC conditions. Data are normalized to siCTL for visualization. Histograms represent average +/- SEM of experiments on cells from three different individuals with technical triplicates. Analyzed by a One-Way ANOVA test.

Genotype of the knock-in mouse lines.

A) Schematic design of the genotype primers for Smad1HA/HA and Smad5PA/PA mouse lines. B) Exemplary gel electrophoresis of PCR products derived from homozygous knock-in mice, heterozygous mice, and WT mice using genotyping primers.

Uncropped western Blots for figure 1