Introduction

The myometrium is the muscular component of the uterus. Throughout gestation, the myometrium maintains contractile quiescence until labor, where uterine contractility facilitates parturition. Pre-term birth, often characterized by aberrant myometrial activity before term, results in the termination of pregnancy prior to 37 weeks gestation. Complications associated with pre-term birth account for 35% of neonatal death and contribute to lifelong physical and socioeconomic multi-morbidities to surviving neonates and mothers (Chawanpaiboon, Vogel et al. 2019). It is therefore critical to understand the molecular mechanisms underlying the maintenance and shift of myometrial quiescence throughout gestation.

The myometrium undergoes extensive structural and functional remodeling preparing for parturition through genomic regulatory mechanisms which influence gene expression throughout pregnancy. Major transcription factor families have been identified to contribute throughout the remodeling process. Studies on the role of steroid hormone receptors in myometrial remodeling have provided evidence that the withdrawal of functional progesterone signaling at term is due to a stoichiometric increase of progesterone receptor (PGR) A to B isoform related estrogen receptor (ESR) alpha expression activation at term. (Mesiano, Chan et al. 2002) (Merlino, Welsh et al. 2007) (Nadeem, Shynlova et al. 2016). Activator protein 1 (AP-1) complex subunits have been observed to act as a PGR coregulators (Dai, Litman et al. 2003) and also have dynamic expression patterns throughout gestation in both humans and rodent models, for example, FOS:JUN heterodimers are implicated to be critical for the initiation of labor through transcriptional regulation of gap junction proteins such as Cx43 (Nadeem, Farine et al. 2018) (Balducci, Risek et al. 1993).

Contributing to the dynamic nature of the myometrial transcriptome during term is the epigenome and its reprogramming. Studies investigating epigenetic markers related to gene activation in the mouse myometrium have revealed that the promoters for contractility-driving genes, such as Cx43, are epigenetically activated well before the onset of labor (Shchuka, Abatti et al. 2020). Considering this, the myometrial epigenome’s role in parturition disorders, such as preterm birth, has thus been under investigation. For example, it is known that altered DNA methylation is linked with preterm birth (Barcelona de Mendoza, Wright et al. 2017), including at the promoters of genes associated with myometrial contraction (Erickson, Myatt et al. 2023) and fetal membrane rupture (Wang, Ogawa et al. 2008). However, the role of differential myometrial DNA methylation at contractility-driving gene promoter CpG islands in preterm birth is not thought to be major (Mitsuya, Singh et al. 2014), but given that DNA methylation mediated gene regulation often occurs outside of CpG islands (Irizarry, Ladd-Acosta et al. 2009), there is still work to be done at this interface. Regardless, whether through epigenetic or transcriptomic gene regulatory programs, the molecular mechanisms underlying the regulation of genes critical for myometrial reprogramming and quiescent maintenance are still poorly understood.

In this study, we aim to improve our understanding of the cellular processes involved in shaping the myometrium’s dynamic state of quiescence throughout pregnancy through the epigenetic and transcriptomic profiling of the human myometrium. Because of the relevance of genomic regulatory mechanisms in coordinating myometrial activity, we identified candidate cis-acting regulatory regions using chromatin immunoprecipitation sequencing (ChIP-Seq) assays of surrogate histone enhancer markers H3K27ac (Creyghton, Cheng et al. 2010) and H3K4me1 (Heintzman, Stuart et al. 2007) alongside chromatin conformation capture assays in term pregnant myometrial tissues. Putative enhancers upstream of the PLCL2, a gene encoding for the protein PLCL2 which has been implicated in the modulation of calcium signaling (Uji, Matsuda et al. 2002) and maintenance of myometrial quiescence (Peavey, Wu et al. 2021), transcriptional start site were subject to functional assessment using CRISPR activation based assays. Here, we identified a genomic region 35-kilobases upstream of the PLCL2 transcriptional start site involved in PLCL2 transcriptional regulation. Furthermore, we aimed to characterize this genomic region through the identification of candidate PLCL2 transcriptional regulators co-localizing with this cis-acting element using integrative cistrome and transcriptome analysis and have demonstrated PGR to be a direct regulator for PLCL2 expression. These findings build upon our understanding of myometrial remodeling throughout gestation and will be pertinent for the development of medical interventions aiming to address pre-term birth.

Results

Epigenomic Landscape in Term Pregnant Myometrial Specimens

The chromatin immunoprecipitation (ChIP) assay identified an average of 44238 H3K27ac and 74325 H3K4me1 positive genomic regions in three term-pregnant, nonlabor myometrial biopsies (Table 1). An RNA-Seq assay revealed an average of 12157 active genes in these specimens that manifested expression levels of fragments per kilobase of transcript per million mapped reads (FPKM) greater or equal to 1 (Table 1). A high throughput chromatin conformation capture (Hi-C) assay further found a total of 27162 chromatin loops from two of the three myometrial specimens (Table 1, subject 1 and 3). We failed to identify chromatin loops in the second subject’s biopsy due to limited sample availability. Together, these results reveal an epigenomic landscape in the human term pregnant myometrial tissue before the onset of labor which we use as a resource to investigate the molecular mechanisms that prepare the myometrium for subsequent parturition.

Table 1.

When comparing the present study to previous findings (Dotts, Reiman et al. 2023), 8563 genomic regions carry common H3K27ac-positive histone marks. However, significant variations on the number and location of H3K27ac-positive intervals are present across the six samples among these two studies (Figure S1). Depending on the number of mapped intervals, the 8563 common regions constitute between 20.7% to 71.5% of total H3K27ac-positive intervals across the six myometrial specimens (Figure S1). Notably, the fraction of intervals commonly identified in the present study’s three specimens ranges between 56.4% and 70.0% of each of their total H3K27ac-positive regions, while the Dotts dataset has a wider range between 28.0% and 74.0%. With a less stringent criterium, most the H3K27ac-positive intervals are found in at least two samples, ranging between 75.4% and 98.3% (Figure S1). These data together highlight the degree of variation on mapping the epigenome among specimens and datasets.

Putative Enhancers for Gene Regulation in the Human Myometrium

Since H3K27ac and H3K4me1 marked genome regions are associated with enhancers for gene regulation (Creyghton, Cheng et al. 2010) (Heintzman, Stuart et al. 2007), we define the regions that have overlapping H3K27ac and H3K4me1 marks as putative myometrial enhancers at the term pregnant nonlabor stage (Dataset S1). Among the three specimens, 13114 putative enhancers are commonly present (Figure 1A). A significant variance on the location of myometrial putative enhancers is also seen among biopsies as evident by the observation that these common putative enhancers make up less than half of the total putative enhancers in each individual specimen (Figure 1A). More than one-third of these 13114 common putative enhancers are located within a 100-kilobase vicinity of actively high-expressing genes (FPKM >= 15), while a much smaller fraction of them (13.3%) are associated with actively low-expressing genes (Figure 1B). Motif enrichment analysis on these common putative enhancers further reveals an overrepresentation of binding motifs for myometrial transcription factors AP-1, STAT5, and NFκB, steroid hormone receptors PGR and NR3C1, and smooth muscle transcription regulators SRF and ELK1 (Figure 1C and Dataset S2). These findings collectively suggest that the putative myometrial enhancers bring together smooth muscle and hormonal control programs for the regulation of myometrial gene expression at term pregnancy.

Figure 1.

Super enhancers house hormone-dependent gene regulatory programs for female reproductive tract homeostasis (Shin, Willi et al. 2016) (Hewitt, Grimm et al. 2020). Across the three human subjects, 540 putative super enhancers are commonly identified in the term pregnant nonlabor myometrial specimens (Figure 2A and Dataset S3). More than 40% of the 540 putative super enhancers are located within a 100-kilobase distance to high-expressing genes, while only 7.3% of putative myometrial super enhancers are found near low-expressing genes (Figure 2B). Compared with the regular putative enhancers, the putative myometrial super enhancers are found more frequently near active genes that are expressed at relatively higher levels (Figure 1B and Figure 2B). Notably, 76% of the putative super enhancers are co-localized with known PGR-occupied regions in the human myometrial tissue (Figure S2). This is significantly higher than the 20% co-localization in the regular enhancer group (Figure S2). Further examining the myometrial active genes that are associated with putative super enhancers revealed an enrichment of gene functions in intercellular and cell-matrix interactions, hormonal responses, and mRNA stability (Figure 2C). Taken together, these results establish the association among the putative myometrial enhancers, the potential regulatory program within these super enhancers, and the cellular functions of the genes they may control.

Figure 2.

Cis-acting elements for the control of the contractile gene PLCL2

Next, we take advantage of this myometrial epigenome map to identify cis-acting elements that regulate the PLCL2 gene expression. The Hi-C assay identified a chromatin loop with one end at the PLCL2 transcription start site (TSS) and the other end at approximately 90-kilobases upstream of the TSS (Figure 3A). Within this chromatin loop, ChIP-Seq results revealed seven genomic regions that are marked with H3K27ac in all three myometrial specimens (Figure 3A). Six of the seven regions are also co-localized with previously published genome occupancy of transcription regulators curated by the ReMap Atlas (Hammal, de Langen et al. 2022) (Figure 3A). For the purpose of functional screening, we focus on H3K27ac signals instead of using H3K27ac/H3K4me1 double positive criterium to cast a wider net. Eight guide RNAs (gRNAs) were designed to target these seven candidate regions for CRISPR activation (CRISPRa) based screening of PLCL2 gene regulation (Figure 3A). The two previously reported gRNAs (PLCL2-7 and PLCL2-8) were able to elevate the endogenous PLCL2 mRNA levels significantly higher than the non-targeting gRNA in immortalized human myometrial cells (hTERT-HM) (Figure 3B) (Peavey, Wu et al. 2021). gRNAs PCLC2-1, PLCL2-3, PLCL2-5, and PLCL2-6 also significantly increased PLCL2 expression above the levels of the non-targeting group (Figure 3B). Moreover, results from the Perturb-Seq assay found that the PLCL2 gene ranks high in the differentially expressed gene list of the PLCL2-5 gRNA expressing cells compared with the non-targeting gRNA expressing cells (Dataset S4), in line with the qRT-PCR assay finding in which the gRNA PLCL2-5 induced PLCL2 expression to the greatest extent among the upstream non-coding genome targeted gRNA (Figure 3B). Moreover, gRNA PLCL2-5 was capable of mediating the increase of the PLCL2 transcript abundance in another uterine mesenchymal cell type T-HESC (Figure 3C). Collectively, these findings support the PLCL2-5 targeted genomic region as a cis-acting element for regulation of the PLCL2 gene.

Figure 3.

PGR as an upstream regulator for the PLCL2 gene in the human myometrium

After determining its cis-acting element role, the PLCL2-5 targeted genomic region was further examined to identify upstream regulators that control myometrial PLCL2 expression. This was achieved by first identifying transcription factor occupancy in this genome region in any tissues/cells that are documented in public databases (ReMAP Atlas) (Hammal, de Langen et al. 2022), followed by filtering for factors that have known mRNA and/or protein expression in hTERT-HM cells and tissues using publicly available transcriptomic and histological data including the Human Protein Atlas (Human Protein Atlas proteinatlas.org) (Uhlen, Fagerberg et al. 2015) and published RNA-Seq datasets (Stanfield, Amini et al. 2019) (Wu, Anderson et al. 2020). The resulting candidate upstream regulators for the myometrial PLCL2 gene are documented in Table 2. Since the PLCL2 mRNA abundance is much lower in the hTERT-HM cells than in the myometrial tissues, these candidate upstream regulators were further grouped into candidate activators and repressors based on an assumption that a candidate activator would be expressed higher in human myometrial tissue, and a candidate repressor being expressed higher instead in the hTERT-HM cells, holding the PLCL2 mRNA levels down. Candidate activators for the PLCL2 gene include known myometrial regulators such as PGR, ESR1, and AP-1 (Table 2). Interestingly, the mediator complex protein subunit MED1, cohesion complex members RAD21 and SMC3, and many epigenomic modifiers are among the list of the candidate repressors (Table 2). These findings not only provide candidates for subsequent functional assessment, but also highlight potential pathways for future investigations on the regulation of the contractile gene PCLC2 expression in the myometrium.

Table 2.

We previously demonstrated the regulation of mouse Plcl2 gene by the myometrial PGR (Wu, Wang et al. 2022). To further test whether such a regulatory mechanism is also present in the human, the CRISPRa technology was employed to increase the expression of PGR mRNA and proteins in hTERT-HM cells (Figure 4A and Figure 4B). Given that the PGR640 gRNA led to higher PGR protein production (Figure 4B), this gRNA was utilized to stimulate myometrial PGR expression for subsequent experiments. Indeed, RT-qPCR results showed that PGR overexpression increased PLCL2 mRNA abundance in hTERT-HM cells (Figure 4C), demonstrating a consistent finding with the mouse model.

Figure 4.

To further explore this regulatory relationship in the human myometrial tissue, we examined the degree of correlation between the inferred PGR activity and the PLCL2 transcript abundance in myometrial biopsies. Gene expression profiles of myometrial specimens from 13 proliferative-phase, 6 secretory-phase, 14 post-menopausal, and 10 Provera-treated human subjects were determined by the RNA-Seq assay (Dataset S5) to take advantage of the wide spectrum of progesterone signaling activities in these samples. Inferred myometrial PGR activities, represented as T-scores (Dataset S5), were derived from this normalized gene expression matrix of 43 human myometrial specimens with the previously published mouse myometrial Pgr gene signature (Wu, Wang et al. 2022) (Li, Bushel et al. 2021). A Pearson correlation analysis found a positive correlation (r = 0.47) between the T-scores and normalized PLCL2 mRNA levels in these 43 specimens (Figure 4D). Together, these in vitro and in vivo findings collectively support that PGR is a PLCL2 activator, likely acting through the 35-kb upstream cis-acting element.

Discussion

Temporally and spatially regulated gene expression is the foundation of tissue phenotype manifestation. Shchuka and colleagues utilize the mouse model to show that the myometrial epigenome is readily programed for parturition-associated gene induction days ahead of the term pregnancy (Shchuka, Abatti et al. 2020), laying the ground for subsequent laboring event. We chose to study the human myometrium at the term pregnant nonlabor stage in order to examine the association between the epigenome landscape and the gene expression pattern not only for preparing the myometrium for parturition, but also for maintaining myometrial quiescence. We also performed a mullti-omic study on a myometrial specimen from each individual human subject to better align all the datasets together, reflecting specimen and subject individuality. While the variation among specimens is substantial, we were able to identify myometrial cis-acting elements commonly shared among three individuals. As a proof of principle, we subsequently used this information to functionally screen for genomic regions that may mediate the regulation of the contractile-suppressing gene PLCL2. The H3K27ac marked 35-kilobase upstream region as a result of the screening then serves as a bridge to identify PGR as one of the upstream regulators for modulating PLCL2 gene expression in human myometrium. These findings showcase the value of these datasets as a resource for future investigations on the mechanisms of myometrial gene regulation.

Variations across cistromic datasets for the human myometrium are significant (Figure 1 and Figure S1) (Dotts, Reiman et al. 2023). Contributing factors likely include but are not limited to subject-to-subject differences and batch effects from sample preparations. Taking the H3K27ac marked genomic regions as an example, regions commonly identified in the specimens of this study’s three subjects account for 55.9%, 67.7%, and 57.9% of total H3K27ac-positive intervals in each individual sample, despite our best attempt to reduce the batch effect. Heterogeneity on cell type compositions, likely due to the sampling location, could be a major contributing factor to the variation. It is also speculated that the genomic regions with subject-to-subject variances could be transiently present and thus the signal could not be captured across samples. In addition, the underlying health conditions, medications, and environmental challenges among others may also affect the epigenomic profile. Variations could also arise from differences in the detection threshold. We identified an average of 44238 H3K27ac marked regions per specimen, while Dotts and colleagues found nearly 19000 per sample (Dotts, Reiman et al. 2023). Reagent and sample handling could confound with biological variances of specimens leading to this difference. Future investigations using standardized sample preparation protocols at the single cell resolution may help to reduce this variation.

Results from the CRISPRa-based Perturb-Seq assay demonstrate the capability of this technology to screen for cis-acting elements that are in topological vicinity of the gene of interest. This assay is particularly useful on genes that are expressed at low levels in the screening platform, i.e., the cells. The dCas9VPR activator can generate satisfactory signal-to-noise ratios to be detected by single cell RNAseq. The advantage of assaying with this system is that each individual cell serves as a container to generate data points instead of relying on multiple wells on cell culture plates. Multiplexing with multiple gRNAs can be achieved by using the gRNA sequence as a unique barcode of the cell, enabling the simultaneous collection of numerous data points for each individual gRNA. In the present study, 56.7% of tested cells carried one species of detectable gRNAs, permitting the study of the resulting effect under a single gRNA. Moreover, another 27.0% of assayed cells have more than one species of detectable gRNAs, which open the possibility to study interactions between different cis-acting elements. However, we failed to use this approach for the purpose of identifying target genes of a putative enhancer of interest (Dataset S4). Challenges such as difficulties on defining the criteria of the associated genes and unsatisfactory signal-to-noise ratios on genes already at high baseline levels should be considered in future experimental designs. Taking together, results of the present study support the use of the CRISPRa-based Pertub-Seq assay under a defined scope.

PLCL2 and its paralog PLCL1 are catalytically inactive members of the phospholipase C family functioning as negative regulators of calcium signaling (Kanematsu, Takeuchi et al. 2005) (Takenaka, Fukami et al. 2003) (Uji, Matsuda et al. 2002). Both of them are expressed in human myometrial tissues and are able to attenuate oxytocin-induced muscle cell contraction (Peavey, Wu et al. 2021). PLCL1 mRNA abundance is higher in the myometrium at term pregnancy than the non-gravid stage, while PLCL2 transcript levels remain comparable between the two stages (Wu, Anderson et al. 2020). Data from cultured human endometrial stromal cells and the genetically engineered mouse model reveal progesterone signaling as the upstream regulator of PLCL1 and Plcl2, respectively (Muter, Brighton et al. 2016) (Peavey, Wu et al. 2021). The present study further demonstrates the PGR-dependent regulation of the PLCL2 gene in the human myometrial cells, suggesting that this pathway is conserved between human and mouse. Importantly, the identification of the cis-acting element 35-kilobase upstream of the PLCL2 gene opens a venue to investigate the impact of interactions between PGR and other myometrial transcription regulators on the mechanism of action that controls the myometrial contractile machinery. The future findings in this space could provide insight on progesterone signaling modification and identify progesterone signaling modifiers for the development of novel therapy that targets parturition disorders.

Materials and Methods

Collection of myometrial specimens

Permission to collect human tissue specimens was prospectively obtained from subjects undergoing hysterectomy or cesarean section for benign clinical indications (H-33461). Matched sSpecimens of gravid myometrium were collected at the margin of hysterotomy from women undergoing clinically indicated cesarean section at term (>38 weeks estimated gestation age) without evidence of labor. Specimens of healthy, non-gravid myometrium were also pecimens were collected from uteri removed from pre-menopausal women undergoing hysterectomy for benign clinical indications. Specimens of the core of leiomyomas as well as healthy myometrium were obtained at least 2 cm from the closest adjacent tumorany adjacent leiomyomas. Any specimens from gravid women receiving active treatment for pre-eclampsia, eclampsia or pregnancy-related hypertension or pre-term labor were excluded.

Human Cell Lines

Human immortalized myometrial cells (hTERT-HM) (Condon, Yin et al. 2002) cells and human immortalized endometrial cells (T-hESC) (ATCC, CRL-4003) were cultured in DMEM/F-12 (Invitrogen, Grand Island, NY, USA) supplemented with 10% Fetal Bovine Serum (FBS, Gibco) and antibiotics (10 000 IU/mL penicillin, 10 000 IU/ mL streptomycin; Life Technologies, Grand Island, NY). Cell culture media was filtered using the 0.22um Rapid-Flow™ Sterile Disposable Filter Units (Nalgene). Cells were incubated at 5% CO2 and 37°C.

RNA Isolation and RT-qPCR

RNA was isolated from cells using TRIzol Reagent (Invitrogen) and the RNeasy mini RNA isolation kit (Qiagen, Hilden, Germany). cDNA was prepared using Moloney Murine Leukemia Virus reverse transcriptase (Thermo Fisher Scientific) with Random Hexamers (Invitrogen, Waltham, MA, USA) according to manufacturer protocol. For quantitative analysis of mRNA, SsoAdvanced Universal SYBR Green Supermix (BioRad, Hercules, CA, USA) was used according to manufacturer instructions. Each reaction was performed in technical duplicates using the standard curve-based method (Larionov, Krause et al. 2005). Briefly, reaction samples were prepared to a volume of 20ul with 5uM of both forward and reverse primers, cDNA, and a final 1x concentration of the SYBR Green Supermix. The reaction was heated to 95°C for 30 seconds, followed by 39 cycles of denaturation at 95°C for 15 seconds and annealing and elongation at 60°C for 30 seconds. Temperature cycles were performed on the CFX Connect^TM Real-Time PCR Detection System (Bio-Rad).

Western Blot Assay

Protein was isolated from cells using the RIPA Lysis and Extraction Buffer (Thermo Scientific) with the following specifications: approximately 300,000 cells were pelleted and lysed with 100ul of complete Pierce RIPA buffer. Protein was quantified using the BCA Protein Assay Kit (Thermo Scientific, Waltham, MA, USA) as instructed by the manufacturer, and 40ug of protein were loaded per lane on a Mini-PROTEAN TGX Pre-cast Protein gel (Bio-Rad, #4568094) with the Precision Plus Protein Dual Color Standards Ladder (Bio-Rad #1610374). Protein lysates from MCF7 cells were used as positive controls, while protein lysates from HEK293T cells were used as negative controls. Gels were transferred to nitrocellulose membrane using the Turbo-Blot transfer system (BioRad) according to the manufacturer’s instructions. The membrane was blocked with 5% milk (Santa Cruz Biotech, Santa Cruz, CA, USA) in TBST (20 mM Tris, pH 7.4 (Lonza, Morrisville, NC, USA), 140 mM NaCl (Lonza), 1% TWEEN-20 (Sigma). PGR protein was detected using Monoclonal Mouse Anti-Human Progesterone receptor clone PgR 1294 diluted 1:1000 in milk overnight at 4C. Bands were detected using Donkey Anti-Mouse 926-32212 (Lot # C60524-15) diluted 1:20,000 in milk for 45 minutes.

Blots were imaged using Odyssey Fc Imager (LI-COR Biosciences) using 800nm channel for 10 minutes, and 700nm channel for 30 seconds (to image ladder). The control protein (B-actin) was detected similarly as above using the following antibodies: Actin (I-19)-R sc-11616-R (Lot#DO406) Rabbit polyclonal IgG Santa Cruz (1:1000 dilution) and Goat anti-rabbit 926-32211 (Lot#DOO304-15) (1:20,000 dilution).

RNA-Seq Library Preparation

RNA-Seq libraries were prepared according to the Illumina TruSeq Stranded mRNA protocol Document # 1000000040498 v00. Libraries were sequenced using the Illumina NextSeq 500 and NovaSeq 6000 systems.

ChIP-Seq Library Preparation

The ChIP-Seq assay was performed by the Active Motif service laboratory using snap-freeze human myometrial specimens. The H3K27Ac ChIP reactions were conducted with 10 μg of tissue chromatin and 4 μg of H3K27Ac antibody (Active Motif, cat # 39133). The H3K4me1 ChIP reactions were carried out by using 10 ug of tissue chromatin and 4 ul of H3K4me1 antibody (Active Motif, cat # 39297). Libraries were prepared by a custom Illumina library with the standard Illumina PE adaptors (Bentley, Balasubramanian et al. 2008). Libraries were sequenced using the Illumina NextSeq 500 and NovaSeq 6000 systems.

Hi-C Library Preparation

Snap-freeze human myometrial specimens were shipped to the Arima Genomics and the Active Motif service laboratories for preparation of the Hi-C library using the Arima-HiC Kit (Arima Genomics A510008), protocol version A160132 v00. Libraries were sequenced in the NIEHS on the Illumina NovaSeq 6000 platform.

CRISPR Activation (CRISPRa) assay

Acquisition of guide RNA expression vectors for CRISPRa assay

All gRNA expression vectors were synthesized by and acquired from VectorBuilder. gRNA expressing vectors for PGR targeting (PGR640A and P6SM), non-targeting control, and the empty backbone were purchased from VectorBuilder under the catalog numbers VB191117-1498rhp (PGR640), VB210824-1268wmh (P6SM), VB191117-1500trj, and VB190918-1522gwq, respectively.

Acquisition of dCas9-VPR-mCherry for CRISPRa assay

IGI-P0492 pHR-dCas9-NLS-VPR-mCherry was a gift from Jacob Corn (Addgene plasmid # 102245 ; http://n2t.net/addgene:102245 ; RRID:Addgene_102245).

Production of Lentiviruses

All lentiviruses were packaged in HEK293T/17 cells (ATCC # CRL-11268) according to published method in Current Protocols in Neuroscience (Chen, Haam et al. 2019). Briefly, 293T cells were transiently transfected with pMD2G, psPAX2 and a transfer vector containing the desired gene using Lipofectamine 2000. Supernatant was collected 48 hours post transfection and concentrated by centrifugation at 50,000 g for 2 hours. Pellets were resuspended in PBS and used for infection. To determine titer, HEK293T/17 cells were infected with lentiviral samples. Five days post infection, Qiagen DNeasy Blood & Tissue Kits was used to isolate chromosomal DNA from infected cells. All titers were determined by performing droplet digital PCR (ddPCR) to measure the number of lentiviral particles that integrated into the host genome.

CRISPRa viral transduction

hTERT-HM cells were infected with lentiviral gRNA expression vectors with a multiplicity of infection of 4 (MOI = 4) and with dCas9-VPR-mCherry at an MOI = 4. gRNA expression vector and dCas9-VPR-mCherry expression vectors, following lentiviral transduction, confer GFP and mCherry fluorescent markers, respectively.

Fluorescence Activated Cell Sorting for CRISPRa assay

Cells were examined using a BD FACSAriaII cell sorter (Becton Dickinson Biosciences, San Jose, CA) equipped with FACSDiVa software. Initially, a “scatter” gate was set on a forward scatter (FSC-A) versus side scatter (SSC-A) dot plot to isolate the principal population of cells free of debris. Subsequently, cells were consecutively gated on a side scatter height (SSC-H) versus width (SSC-W), then a forward scatter height (FSC-H) versus width (FSC-W) dot plot to isolate single cells. GFP expressing cells were excited off a 488 nm laser and detected using a 525/50 nm filter. mCherry expressing cells were excited off a 561 nm laser and detected using a 610/20 nm filter. Double positive cells were collected based off non-transfected controls and used for further cultural/biochemical analysis.

Perturb-Seq library preparation

The cells in suspension were counted and examined for viability with trypan blue staining using a TC-20 cell counter (Bio-Rad). Approximately 16,500 live cells at 1×10⁶ cells/ml concentration with 65% or above viability were loaded into the Single Cell Chip to generate single cell emulsion in Chromium Controller with Chromium Single Cell 3’ Library & Gel Bead Kit v3.1 (10x Genomics, Cat. 1000268). Reverse transcription of mRNA and cDNA amplification were carried out following the manufacture’s instruction (10x Genomics, Cat. 1000268, Cat. 1000262 with 10x Genomics protocol CG000316). The amplified cDNA was separated into CRISPR sgRNA derived cDNA and transcriptome derived cDNA. The CRISPR sgRNA derived cDNA was used to make NGS sequencing libraries. The transcriptome derived cDNA was further fragmented to construct NGS libraries. Both libraries were then sequenced together with the molar ratio of 1 to 4 by the NIEHS Epigenomics and DNA Sequencing Core Laboratory with the parameters recommended in the manufacture’s instruction.

Human enhancer perturb-seq data analysis

The raw sequencing FASTQ files generated from both the transcriptome and CRISPR screening libraries were processed together by Cell Ranger software (version 4.0.0, 10× Genomics). The “cellranger count” pipeline used STAR for aligning the reads to the human reference, GRCh38 “refdata-gex-GRCh38-2020-A” (10X Genomics), and associated gene expression profile with guide RNA (gRNA) identity by unique barcode in each cell. Seurat software (version 3.6.3) was utilized to perform clustering analysis on the combined dataset (Satija 2015, PMID: 25867923). We applied the SCTransform package to normalize gene expression counts across cells (Hafemeister 2019, PMID: 31870423). The cells were clustered based on number of unique gRNA type detected. The cell populations containing more than one type of gRNA were excluded from further analyses.

We applied two approaches to determine the association between gRNAs and target genes:

1. For each cell cluster having gRNA with known target genes, such as PLCL1 and PLCL2, we compared the expression of target genes in these clusters to the cluster with scramble control gRNA. We defined individual gRNA activation as a fold change greater than 1.5.

2. For cell clusters having gRNAs with unknown target genes, we identified differentially expressed genes by performing pair-wise comparisons between each cluster and the cluster with scramble control gRNA. This analysis was conducted using the “FindAllMarkers” function in the Seurat package. Target genes were defined as those with a fold change greater than 1.5 and an adjusted p-value less than 0.05.

ChIP-Seq

The raw ChIP-seq reads (51-bp, single-end) were filtered with average quality scores greater than 20. Adaptor sequences were trimmed from reads using cutadapt (v1.12). Then the reads were aligned to the human reference genome (hg38) using Bowtie (v1.1.2) (Langmead, Trapnell et al. 2009), requiring unique mapping and allowing up to 2 mismatches per read (-m 1 -v 2). Duplicated reads with identical sequences were removed using Picard tools. To visualize the read coverage, BigWig files were generated from the bedgraph files of each sample using bedGraphToBigWig. These bigWig files were displayed as custom tracks on the UCSC genome browser.

The normalization of sequencing depth across all ChIP-seq datasets were achieved by down sampling to 20 million uniquely mapped reads per sample. The initial peaks of each sample were identified using MACS2 with a cutoff of adjusted p-value < 0.0001 (Zhang, Liu et al. 2008). The gene associated with each peak was predicted by searching the transcription start site (TSS) of nearby genes within a 100Kb range using HOMER (Heinz, Benner et al. 2010).

Identification of union peaks between H3K27ac and H3K4me1 peaks

The initial union peaks between H3K27ac and H3K4me1 in each sample were identified using the “merge” function of BedTools. The common union peaks from multiple samples were used for downstream analysis.

Identification of union peaks for term pregnant myometrial specimens

Additional H3K27ac ChIP-seq and related input data were obtained from GEO (GSE202027; patient 1, 4, and 7). The same ChIP-seq analysis pipeline described above was applied to those samples. To obtain a comprehensive set of peaks across all H3K27ac ChIP-seq datasets, the peak intervals from each dataset were merged using BedTools.

Identification of super enhancers

H3K27ac ChIP-seq peaks were defined as H3K27ac-positive enhancers in each sample. The enhancers within 12.5Kb were merged. The combined enhancer regions were called super enhancers if they were larger than 15Kb. The common super enhancers from multiple samples were used for downstream analysis.

RNA-Seq analysis of Term Myometrial Specimens

The raw paired-end reads (51-bp, paired-end) were initially processed by applying a filter with an average quality score of 20. Adaptor sequences were subsequently removed from reads using cutadapt version 1.12. The processed reads were then aligned to the human reference genome (hg38) using the STAR aligner version 2.5.2b. The gene expression levels in each sample were determined by counting the total number of paired-end reads mapped to each gene using DESeq2 R package version 1.12.4 (Love, Huber et al. 2014). The gene expression count matrix was normalized based on ratio of total mapped read pairs in each sample to 56.5 million. Only the genes with average normalized count across all samples larger than one were used for further analysis. The Bioconductor package edgeR was applied to the gene expression count matrix to detect differentially expressed genes between groups of interest (Robinson, McCarthy et al. 2010). The false discovery rate (FDR) of differentially expressed genes was estimated using the Benjamini and Hochberg method.

RNA-Seq analysis of Provera Treated Human Myometrial Specimens

The raw paired-end reads (76-bp, paired-end) were processed by applying a filter with an average quality score of 20. Adaptor sequences were removed from reads using cutadapt (v1.12). The processed reads were then aligned to the human reference genome (hg38) using the STAR aligner (v2.5.2b). Fragments Per Kilobase Million (FPKM) was calculated for each gene in each sample by using Cufflinks version 2.0.2.

Hi-C

The raw reads (51-bp, paired-end) generated from HiC library were mapped to human reference genome (hg38; Genome Reference Consortium Human GRCh38 from December 2013) using HiCUP version 0.7.1 (Wingett, Ewels et al. 2015). The uniquely mapped di-tag passing quality filtering with distance larger than 10Kb were used for downstream analysis. Chromatin loops were identified using “hiccups” function of Juicer version 1.8.9 with default parameters (Durand, Shamim et al. 2016). The A/B compartments were predicted at 100Kb resolution using “eigenvector” function of Juicer.

DNA-binding factor motif enrichment analysis

Enriched motifs were identified by HOMER (Hypergeometric Optimization of Motif EnRichment) v4.11 (Heinz, Benner et al. 2010).

Inferred myometrial PGR activities and the correlation analysis

The differentially expressed gene list between control and myometrial Pgr knockout groups at the mid-pregnancy stage was used as the mouse myometrial Pgr gene signature (Wu, Wang et al. 2022). T-scores of the 43 human myometrial specimens were computed by the SEMIPs application with the mouse myometrial Pgr gene signature and the normalized gene expression counts of the human biopsies (Li, Bushel et al. 2021). The Pearson correlation analysis was performed on Excel.

Data Availability

All the Nest Generation Sequencing Data were deposited to the National Center for Biotechnology Information Gene Expression Omnibus under the accession number GSE244735 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE244735).

Author Contribution

SPW, LL, and FJD conceived the study. MLA provided crucial reagent. SPW, EQ, SMR, and XX performed experiments. SPW, TW, EQ, and SMR conducted data analyses. SPW, TW, EQ, XX, MLA, and SMR wrote the initial draft.

Acknowledgements

We thank the following National Institute of Environmental Health Sciences (NIEHS) Core facilities for exceptional technical supports: Epigenomic and DNA Sequencing Core, Viral Vector Core, Flowcytometry Core, and Integrative Bioinformatics Supportive Group. We thank Ms. Olivia Emery for her technical assistance. This study is partly supported by the Intramural Research Program of the National Institute of Environmental Health Sciences Z1AES103311 (FJD) and by the Gaine Research Foundation GRF-2018-01 (LL). EQ is the awardee of the NIH Intramural Training Program Fellowship. SMR is the recipient of the NIEHS Scholar Connect Program Fellowship.

Figure S1. Subject variations on H3K27ac-positive histone marks.

Figure S2. PGR occupancy in myometrial enhancers. (A) Percentages of enhancers and super enhancers that manifest PGR occupancy. PGR genome occupancy data was previously published in NCBI GEO accession numbers GSM4081683 and GSM4081684.

Dataset S1. H3K27ac and H3k4me1 double positive enhancers in term pregnant not in labor human myometrial specimens.

Dataset S2. Enrichment of known transcription factor binding motifs in putative myometrial enhancers.

Dataset S3. Super enhancers in term pregnant not in labor human myometrial specimens.

Dataset S4. CRISPRa dependent gene expression patterns in hTERT-HM cells and gRNA information. DEG, differentially expressed genes between denoted gRNAs.

Dataset S5. Normalized gene expression counts across the 43 human myometrial specimens and the PGR T-Scores of individual specimens.