Up-regulation of labor-associated genes involves a transcriptional mechanism

To establish a comprehensive profile of pregnant and laboring myometrial transcriptomes, we conducted total strand-specific RNA-seq on RNA isolated from the myometrium of pregnant C57BL/6 mice at gestational day 15 or day 19 while in active labor (n = 5 each, Fig 1A). Based on the RNA-seq data, we observed clustering of the same samples within each time point of collection, as expected (S1 Fig and S1 Data). Differential gene expression analysis based on exon read counts (S1 Table) revealed that a total of 956 genes showed gestational time point–varying expression levels (Fig 1B, fold change cutoff of 4, p < 0.01). Hierarchical clustering analysis of these genes confirmed similar expression trends from mice of the same gestational age (Fig 1C and S2 Data) while gene ontology (GO) term analyses highlighted the involvement of down-regulated and up-regulated genes at term in myometrial relaxation and contraction pathways, respectively (S2 and S3 Tables). In all, 578 genes exhibited a significant increase in expression during active labor compared with day 15. Apart from up-regulation of Fos (Fig 1D), these genes included (but were not limited to) prominent labor-associated players Gja1, Prostaglandin G/H synthase 2 (Ptgs2), and Oxytocin receptor (Oxtr), as well as matrix metalloproteinases (Matrix metalloproteinases 7 [Mmp7], 11 [Mmp11], and 12 [Mmp12]), signaling proteins (Chemokine C-X-C motif ligands 1 [Cxcl1], and 5 [Cxcl5]), and adhesion molecules and proteins (Vascular adhesion molecule 1 [Vcam1], Thrombospondin 1 [Thbs1], Carcinoembryonic antigen-related cell adhesion molecule 1 [Ceacam1]) known to exhibit elevated levels at term. Conversely, 378 genes were found to be significantly down-regulated during active labor compared with day 15, including (but not limited to) genes encoding proteins responsible for cell–extracellular matrix interactions (Collagens type IV alpha 6 [Col4a6], type XI alpha 1 [Col11a1], type XIII alpha 1 [Col13a1], type XV alpha 1 [Col15a1], type XXVI alpha 1 [Col26a1], and Sparc/osteonectin, cwcv and kazal-like domains proteoglycan 2 [Spock2]), proteins involved in calcium signaling (Melanin-concentrating hormone receptor 1 [Mchr1], Calmodulin-like 3 [Calml3], Calbindin 2 [Calb2]), proteins regulating myometrium response to low oxygen tension and resistance to oxidative stress (Hypoxia inducible factor 3 alpha subunit [Hif3a], Aldo-keto reductase family 1 member b7 [Akr1b7]), and voltage-dependent calcium, potassium, and water channels (Calcium channel voltage-dependent R type alpha 1E subunit [Cacna1e], Potassium voltage-gated channel subfamily G member 1 [Kcng1], and Aquaporin 8 [Aqp8], respectively).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Quiescent and term laboring myometrial transcriptomes exhibit differential expression profiles.

(A) Gestational schematic outlining days and time points at which myometrial tissues were collected for transcriptome- or genome-wide sequencing analyses. Collection days include gestational d15, d19 TNIL, d19 LAB, and pp. (B) RNA-seq volcano plot highlighting transcriptional status of genes exhibiting differential expression levels between d15 and d19 LAB myometrial tissues. Data associated with this figure can be found in S1 Table. (C) Hierarchical clustering of gene groups based on RNA expression changes between d15 and d19 LAB samples (n = 5 per gestational day). Data associated with this figure can be found in S2 Data. (D) Total RNA-seq reads (RPM) at the labor-associated Fos gene locus for d15 and d19 LAB samples mapped to the mm10 mouse genome assembly. chr12, chromosome 12; d, day; FC, fold-change; LAB, active labor; NS, nonsignificant; pp, postpartum; RNA-seq, RNA-sequencing; RPM, reads per million; TNIL, term-not-in-labor.

https://doi.org/10.1371/journal.pbio.3000710.g001

Differential exonic RNA profiles, however, do not in and of themselves reflect a regulatory mechanism change at the level of transcription for those genes. Mediation of gene expression regulation can take place at multiple stages within a gene’s expression pathway prior to formation of a final translated and functional protein product. Apart from varying of levels of primary transcript generation in the nucleus, these mechanisms can include nuclear retention of processed mRNA [24,25], alternative splicing [26,27], and increased mRNA stability [28,29]. Furthermore, multiple gene activation and repression mechanisms can compete to determine the expression outcome of a particular gene, as demonstrated by select biological contexts in which there is an imperfect correlation between the levels of a gene’s transcribed and its translated products [30].

Because sequencing of total RNA provides a transcriptome-wide profile of reads spanning both exons and introns, both spliced and unspliced RNA species can be detected [31,32]. If critical labor-driving genes were significantly up-regulated at the level of active transcription, we would expect to identify a substantial increase in reads corresponding to gene introns as well as those corresponding to their exons. Conversely, if activated genes were regulated exclusively or predominantly by posttranscriptional mechanisms such as RNA stability, we would expect that those genes’ corresponding intron reads would remain relatively constant throughout gestation. Given that a prior study has posited that mRNA stability may act as a critical player in regulation of Gja1 in particular [33], we sought to definitively ascertain whether increased contractility-associated gene expression during labor involved any regulatory input from transcriptional mechanisms. Upon inspection of the Fos and Gja1 loci, we noted striking labor-specific intronic RNA enrichment profiles at both genes (Fig 2A). Subsequently, we conducted a genome-wide intron reads–concentrated gene expression analysis (intron RNA-seq [iRNA-seq]) [31], which uncovered multiple genes that displayed increased primary transcript generation at term during active labor relative to day 15 (Fig 2B and S4 Table). In fact, the majority of genes (55%) up-regulated at labor on the basis of increased exon read accumulation also displayed a significant increase in intron reads (Fig 2C). Using exon–intron junction-spanning primers, we confirmed a significant increase in primary transcript levels of well-known labor-associated genes—Gja1, Fos, Fosl2, Oxtr, and Ptgs2—in active labor relative to day 15 and, with the exception of Fos, to day 19 term-not-in-labor (Fig 2D and S3 Data). These results demonstrate that genes promoting the contractile state in the myometrium at term act because of a rapid gestational time point–specific increase in primary transcript levels, which suggests that substantial transcriptional activity occurs in myometrial nuclei during labor.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Up-regulation of labor-associated genes involves an increase in their primary transcript levels.

(A) Exon- and intron-specific RNA-seq reads (RPM) at the labor-associated Fos gene locus from d15 and d19 LAB samples, mapped to the mm10 mouse genome assembly. (B) iRNA-seq volcano plot highlighting expression changes based on intron reads between d15 and d19 LAB samples. (C) Venn diagram displaying the number of genes significantly up-regulated during labor relative to d15 that show enrichment in exon and intron reads (region of overlap) or enrichment in either exon or intron reads (regions of nonoverlap). (D) Confirmation of laboring time point–specific up-regulation of primary transcript expression of contractility-promoting genes by RT-qPCR. Groups labeled with different letters show significant difference, with p < 0.05. Data associated with this figure can be found in S3 Data. chr12, chromosome 12; d, day; FC, fold-change; iRNA-seq, intron RNA-seq; LAB, active labor; NS, nonsignificant; P(adj), adjusted p-value; RNA-seq, RNA-sequencing; RPM, reads per million; RT-qPCR, reverse transcriptase-quantitative polymerase chain reaction; TNIL, term-not-in-labor.

https://doi.org/10.1371/journal.pbio.3000710.g002

Histone marks associated with gene activation are enriched at labor-associated gene promoters well in advance of labor onset

Having determined that the majority of genes exhibiting increased expression levels during labor were associated with increased primary transcript generation, we next investigated the active chromatin landscape surrounding these genes. After optimizing the protocol for myometrial tissue and confirming target enrichment at control regions (S2 Fig and S4 Data), we conducted chromatin immunoprecipitation with massively parallel sequencing (ChIP-seq). We targeted H3K4me3 and H3K27ac enrichment events on a genome-wide scale at day 15, day 19 term-not-in-labor, day 19 active labor, and postpartum, and found gestational time point sample replicates were highly correlated, as expected (R² > 0.79, S5 Table). When we compared the histone profiles at all promoters to the corresponding genes’ transcriptional status within a particular time point, we observed that both H3K4me3 and H3K27ac are enriched at the promoters of genes highly expressed on the designated gestational day (S3 Fig and S5 Data). However, an initial examination of the Fos locus unexpectedly revealed an enrichment of both active chromatin markers at the Fos promoter across all four time points (Fig 3A); furthermore, the same trend was observed for labor-associated genes Fosl2, Gja1, Oxtr, and Ptgs2 (S4 Fig, histone panels). We next sought to establish the active histone profile across all gene promoters at the designated gestational stages (S6/S7 Tables). Even more surprisingly, we found a similar genome-wide active histone marker enrichment pattern across gene promoters (+/− 2 kb of TSSs) in all four time points (Fig 3B). A K-means clustering analysis further revealed groups of promoters with differing H3K27ac and H3K4me3 enrichment levels but failed to identify groups of genes displaying gestational day–associated changes in these histone modifications’ profiles (S5 Fig). When we investigated the levels of these modifications at promoters of exclusively those genes that displayed a significant increase in transcription (based on an increase in intron reads) during active labor, we observed only a moderate increase in accumulation of both markers from day 15 to term (Fig 3C and 3D and S6, S7 and S8 Datas). Our analyses therefore indicate that H3K27ac and H3K4me3 enrichment premarks the promoters of labor–up-regulated genes as early as day 15 of gestation.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Gene expression activation–associated histone marks are enriched at promoters of contractility-driving genes in advance of labor onset.

(A) Anti-H3K4me3 and anti-H3K27ac ChIP-seq reads (RPM) mapped at promoter of labor-associated gene Fos in d15, d19 TNIL, d19 LAB, and pp samples. (B) Genome-wide enrichment of H3K4me3 and H3K27ac at gene promoters throughout the genome. Signal +/− 2 kb of TSSs is displayed, with genes ordered at each indicated time point according to decreasing enrichment (red → white) profile in d15 samples. (C) Plots exhibiting log₂-fold H3K4me3 or H3K27ac signal, normalized to input, at promoters of genes whose expression is enriched in laboring samples relative to d15 samples based on intron reads. Data associated with this figure can be found in S6 Data. (D) Violin plots displaying log₂-fold H3K4me3 or H3K27ac signal (RPM), normalized to input, at promoters of genes whose intron read–based expression is enriched in laboring samples relative to d15 samples. Groups labeled with different letters show significant difference, with p < 0.05. Data associated with this figure can be found in S7 and S8 Datas. ChIP-seq, chromatin immunoprecipitation with massively parallel sequencing; chr12, chromosome 12; d, day; Fos, FBJ osteosarcoma oncogene; H3K4me3, H3 trimethylation of lysine residue 4; H3K27ac, H3 acetylation on lysine residue 27; LAB, active labor; pp, postpartum; RPM, reads per million; TNIL, term-not-in-labor; TSSs, transcription start site.

https://doi.org/10.1371/journal.pbio.3000710.g003

Loci of labor-associated genes acquire RNAPII gene body occupancy and eRNA enrichment events closer to term

Given that labor-associated gene loci exhibited strikingly similar active histone mark enrichment profiles across all four gestational time points, we next sought to establish whether concomitant binding of RNAPII, an event required for transcription initiation, occurred at contractility-driving gene promoters and bodies as early in the gestational time course. We conducted ChIP-seq to identify RNAPII enrichment events (targeting serine 5 phosphorylated RNA Polymerase II Subunit A) at gestational days 15 and 19 (during active labor) and found RNAPII-enriched broad peaks in both stages (S8 Table). Again, our correlation analysis confirmed the clustering of gestational time point replicates (S6 Fig and S9 Data) while an examination of RNAPII enrichment values alongside our RNA-seq data revealed that RNAPII has a more prominent binding profile at highly expressed genes (S7 Fig and S5 Data and S9 Table). Contrary to the active histone marker trend we observed, labor-driving gene promoters exhibited differential RNAPII binding profiles at either gestational stage. We found that promoter and gene body polymerase occupancy at genes with substantially higher expression levels during labor was significantly higher at term relative to day 15 (Fig 4A and 4B and S6 and S10 Datas); such genes include, but are not limited to, Fos, Gja1, Ptgs2, and Fosl2 (Figs 4C and S4). This result correlated with our observations of these genes’ gestational time point–specific primary transcript levels and further supported the notion that transcriptional mechanisms underlie the rapid and prompt gene expression up-regulation events that underlie the myometrial state transition toward contractility.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Active labor is associated with recruitment of RNAPII to labor-driving genes and eRNA expression at gene-adjacent intergenic regions with H3K27ac peaks.

(A) Metagene plot exhibiting log₂-fold RNAPII signal, normalized to input, at bodies of genes whose intron read–based expression levels are enriched in d19 LAB samples relative to d15 samples. Data associated with this figure can be found in S6 Data. (B) Violin plots displaying log₂-fold RNAPII signal (RPM), normalized to input, at promoters of genes whose intron read–based expression levels are enriched in d19 LAB samples relative to d15 samples. Groups labeled with different letters show significant difference, with p < 0.05. Data associated with this figure can be found in S10 Data. (C) RNA-seq, RNAPII, and H3K27ac ChIP-seq reads (RPM) mapped at the Fos locus in d15 or d19 LAB samples. Regions containing RNAPII (red bars) or H3K27ac peaks (blue bars) indicated in d15 and d19 LAB samples. RNA-seq profile at labor–up-regulated eRNA-exhibiting region (situated between the dashed lines) downstream of Fos is presented at a different scale compared with the rest of the locus. (D) eRNA volcano plot highlighting intergenic H3K27ac peaks in genomic regions that exhibit significant differences in eRNA levels between d15 and d19 LAB samples. (E) All enriched transcription factor motifs identified at genomic regions that feature intergenic H3K27ac peaks and exhibit a significant increase in eRNA levels in d19 LAB samples compared with d15 samples. (F) The top 10 enriched motifs identified at promoters for genes that exhibit a significant increase in intronic RNA levels for d19 laboring samples compared with d15 samples. ChIP-seq, chromatin immunoprecipitation with massively parallel sequencing; chr12, chromosome 12; d, day; eRNA, enhancer RNA; FC, fold-change; Fos, FBJ osteosarcoma oncogene; H3K27ac, H3 acetylation on lysine residue 27; LAB, active labor; NS, nonsignificant; RNAPII, RNA polymerase II; RNA-seq, RNA-sequencing; RPM, reads per million; P(adj), adjusted p-value; TES, transcription end site; TSS, transcription start site.

https://doi.org/10.1371/journal.pbio.3000710.g004

Among the genes significantly up-regulated during labor, Fos exhibited increased RNAPII occupancy across its gene body relative to day 15, as we expected. When we further expanded our view outside the gene body, we observed a genomic region that does not encode a gene but includes both a prominent RNAPII peak (marked by red bar within dashed lines in Fig 4C) as well as an intergenic H3K27ac enrichment event 12 kb downstream of the Fos gene. As was the case with the H3K27ac signal at the Fos promoter, an intergenic H3K27ac peak was identified in both day 15 as well as day 19 labor samples; however, RNAPII association at this region was more pronounced and only identified as a peak in the labor context (Fig 4C). Because intergenic regions containing RNAPII and H3K27ac peaks have been noted to occur at active enhancer regions [6,12,34], we examined the loci of other labor-associated genes to see whether they contained intergenic regions that similarly exhibit this epigenetic profile. We found that, although we focused our earlier analyses on H3K27ac signal enrichment at promoter regions, 43% of the identified H3K27ac peaks in the laboring samples are located at distances greater than 2 kb from a gene TSS (Fig 4D and S10 Table). As was the case with H3K27ac peaks at promoters throughout the genome, intergenic H3K27ac peaks were mostly invariant across our tested gestational time points, with only 11/5,041 (0.2%) of the intergenic H3K27ac peaks displaying a significant increase in H3K27ac signal enrichment in labor samples compared with day 15 samples (S11 Table). However, despite the presence of H3K27ac enrichment at intergenic regions across all four tested time points, many gene loci contain H3K27ac-modified regions with associated differentially transcribed eRNAs (S12 Table). Furthermore, we observed that several of these regions contain H3K27ac peaks and display a significant increase in eRNA expression levels during labor (Fig 4D). Therefore, although most H3K27ac peaks across the myometrial genome are present both at day 15 and at term, several regions containing those peaks transcribe significantly higher amounts of eRNA at term.

To identify the transcription factor motifs that could underlie the changes in gene expression that occur during labor, we conducted motif enrichment analyses using Hypergeometric Optimization of Motif Enrichment (HOMER) [35]. We uncovered a significant enrichment of AP-1 motifs (TGACTCA) in labor-associated eRNA–up-regulated intergenic regions (Fig 4E), implicating this family of transcription factors in a modulatory role with regard to enhancer activity at labor onset. We made a similar observation at promoters of genes displaying increased primary transcript abundance during labor (Figs 4F, S8 and S6 Data). Additionally, apart from AP-1 motifs, the promoters of these genes contain an enrichment of several other motifs: Transcription factor 3 (TCF3), an E-protein transcription factor that has been shown to assist coactivator proteins in the induction of gene transcription events in other cell contexts [36]; CCCTC-binding factor (CTCF), a zinc finger protein prominently known to bind promoter and enhancer regulatory elements [37]; and RELA or NFkB-p65 (Nuclear factor-kappa-beta p65 subunit), a member of the labor-associated NFkB-p65/IL-6 inflammatory pathway [38,39] that has also been affiliated with inducing transcription at the Oxtr promoter [40]. Taken together, these results lead us to propose that the controlled expression of labor-associated genes is driven by transcriptional regulation mechanisms, despite the apparent epigenetic activation of labor-associated loci well in advance of labor onset.

Ethics statement

This study was carried out in accordance with standards set out by the Canadian Council on Animal Care (CCAC). Animal Use Protocol AUP# 21-0164H was reviewed and approved by the Animal Care Committee of The Center for Phenogenomics (TCP), Toronto, Canada. All animals in this AUP were maintained and used in accordance with the current recommendations of CCAC, the requirements under the Animals for Research Act, RSO 1990, and The Centre for Phenogenomics Committee Policies and Guidelines. All research using animal tissues was performed in a class II certified laboratory by qualified staff trained in biological and chemical safety protocols and in accordance with Health Canada guidelines and regulations.

Animal model

Bl6 or C57/Bl6 mice used in these experiments were purchased from Harlan Laboratories (http://www.harlan.com/). All mice were housed under specific pathogen–free conditions on a 12L:12D cycle and were administered food and water ad libitum. Female mice were mated overnight with males, and the day on which vaginal plugs were detected was designated as day 1 of gestation. Pregnant mice were maintained until the appropriate gestational time point. The day of delivery was day 19 of gestation. Our criteria for labor were based on delivery of at least one pup.

Tissue collection

Animals were euthanized by carbon dioxide inhalation, and myometrial samples were collected on gestational day 15, day 19 term not in labor term-not-in-labor, day 19 during active labor, and 2–6 hours postpartum. Tissue was collected at 10 AM on all days with the exception of the labor sample (during active labor), which was collected once the animals had delivered at least one pup. The part of the uterine horn close to the cervix from which the fetus was already expelled was removed and discarded; the remainder was collected for analysis. For each day of gestation, tissue was collected from 4–6 different animals. Uteri were placed into ice-cold PBS. Uterine horns were bisected longitudinally and dissected away from both pups and placentas. The decidua basalis was cut away from the myometrial tissue. The decidua parietalis was carefully removed from the myometrial tissue by mechanical scraping on ice, which removed the entire luminal and glandular epithelium and the majority of the uterine stroma. Myometrial tissues were flash-frozen in liquid nitrogen and stored at −80°C. When necessary, myometrial tissues were crushed into fine powder on dry ice prior to subjection to the experimental methods listed below.

Chromatin immunoprecipitation

Histone marker–targeting chromatin immunoprecipitation (ChIP) was conducted using the protocol described by Young laboratory (younglab.wi.mit.edu/hESRegulation/Young_Protocol.doc), with some necessary modifications for myometrial tissue. Crushed myometrial tissue was fixed in a 1% paraformaldehyde solution at room temperature, and the reaction quenched in a 0.125 M glycine solution. Cells were rinsed twice with 1× cold PBS, and pellets were flash-frozen and stored at −80°C until needed. Pellets were washed and samples were lysed in successive lysis buffers, followed by subjection to sonication via the Covaris sonicator with a pulse ON time of 10 seconds at 30 amps for a total of 30 cycles. Aliquots of sonicated sample were run on a 2% agarose gel to confirm chromatin was sonicated to a range of 300–500 bp in size. Sonicated samples were treated with 10% Triton X-100 and spun down at 4°C to pellet debris. Aliquots of cell lysate supernatant to be used as input were stored at −20°C.

To bind antibody to magnetic beads, Dynal Protein A and Protein G beads (added in a 1:1 ratio) were washed and resuspended in block solution. Anti-H3K27ac antibody (Abcam, ab4729) or anti-H3K4me3 (Abcam, ab8580) was added, as appropriate, to beads and incubated at 4°C with rotation. Beads were again washed and resuspended in block solution. Antibody and magnetic bead mix was added to remaining cell lysate, and samples were incubated at 4°C overnight with rotation. IP samples were washed at 4°C and eluted at 65°C. Supernatant was removed from spun-down beads, and cross-links were reversed at 65°C. Samples were treated with RNaseA at 37°C, treated with Proteinase K at 55°C, cleaned via phenol-chloroform treatment, and stored in ethanol at −20°C overnight. DNA pellets were washed with 80% EtOH and resuspended in Tris-HCl. Validation of the ChIP method was performed using ChIP-qPCR primers (sequences in S11 Table) that targeted regions expected to be enriched in our marker of interest.

RNAPII-targeting ChIP was performed as outlined in the supplemental methods in Mitchell and Fraser [59], with modifications for collection of myometrial tissue as performed in the case of the histone ChIP, and using anti-RNAPII (RPB1 serine 5 phosphorylated form) antibody (Abcam, ab5131).

ChIP sequencing and mapping

ChIP (n = 2) and input (n = 1) samples from each gestational time point—day 15, day 19 term-not-in-labor, day 19 term labor, and postpartum—were submitted for single-end 50-bp read sequencing using standard Illumina HiSeq 2500 protocols. Reads were quality checked using FastQC, trimmed with bbduk, and mapped to the GRCm38/mm10 mouse reference genome using STAR [60].

ChIP normalization and peak calling

Peaks were called for each individual replicate using MACS2 broad peak calling [61]. Significantly conserved peaks in both biological replicates were combined using irreproducible discovery rate (IDR). Differential peak analysis was performed using the diffBind package [62]. Peaks with a fold change ≥ 4 and adjusted p-value < 0.01 were considered significantly different between day 15 and term labor samples. Peaks with a significant increase in signal intensity in labor samples were linked to the closest gene TSSs using bedtools [63]. Normalized ChIP-seq reads (RPM) at promoters (+/− 2 kb of TSSs) of labor-associated up-regulated genes in H3K4me3-, H3K27ac-, and RNAPII-targeted samples were quantified using Seqmonk (https://www.bioinformatics.babraham.ac.uk/projects/seqmonk/). Kruskal–Wallis test was used to measure significant (p < 0.05) changes in enrichment values (RPM) among different time points. NGS.plot [64] was used to plot the K-means clustering heatmaps using default settings. Results were plotted using ggplot2 [65]. Sequencing data files were submitted to the Gene Expression Omnibus (GEO) repository (GSE124295).

Gene expression quantification by RNA extraction and RT-qPCR

Total RNA was extracted from crushed myometrial tissue using Trizol and further treated with DNaseI to remove genomic DNA. RNA was reverse transcribed using the high-capacity cDNA synthesis kit (Thermo Fisher Scientific). Target gene expression was monitored by qPCR using exon–intron boundary-spanning primers (S12 Table) for primary transcript detection and normalized to levels of total Hist1 mRNA because the Hist1 reference gene was consistently expressed at similar levels across gestational time points. Expression levels were calculated against Bl6 or F1 genomic DNA-based standard curve references. All samples were confirmed not to have DNA contamination because no amplification was observed in reverse transcriptase negative samples. Relative expression values were plotted using GraphPad Prism 8. Significant changes in expression were determined by one-way ANOVA with Tukey correction.

RNA-seq quantification and differential expression analysis

DNAseI-treated total RNA samples isolated from day 15 and term labor mice (n = 5 each) were subjected to paired-end sequencing using standard Illumina HiSeq 2500 protocols. Reads were quality checked using FastQC, trimmed with bbduk, and mapped to the GRCm38/mm10 mouse reference genome using STAR [60]. Exon-mapped reads were quantified using featureCounts [66]. Intron reads were quantified using SeqMonk’s active transcription quantitation pipeline (http://www.bioinformatics.babraham.ac.uk/projects/seqmonk/). Alternative transcript counts were summed together for every gene. Intron reads were then imported into DESeq2 [67] for differential expression analysis. Genes with a fold change ≥ 4 and adjusted p-value < 0.01 were considered significantly changing. Differential RNA expression data were plotted using the EnhancedVolcano package (https://github.com/kevinblighe/EnhancedVolcano). Reads were normalized for gene expression levels across replicates. Heatmaps were plotted using the pheatmap package (https://cran.r-project.org/web/packages/pheatmap/index.html). Genomic regions of interest for eRNA expression analysis were selected based on intergenic regions featuring H3K27ac peaks. RNA-seq data at these regions were subjected to differential RNA expression analysis by DESeq2. H3K27ac peaks with an eRNA fold change ≥ 4 and adjusted p-value < 0.05 were considered as peaks with significantly changing signal intensity from day 15 to day 19 (labor), and differential eRNA expression was plotted using EnhancedVolcano package. Sequencing data files were submitted to the GEO repository (GSE124295).

Motif enrichment analyses

Enrichment of transcription factor motifs at promoters and intergenic regions of labor-associated genes was performed using HOMER motif analysis tool [35]. Promoter sequences (−1 kb of TSSs) of up-regulated genes at labor were compared against an input of randomly selected 1-kb-length promoter sequences. Intergenic regions containing labor-associated H3K27ac peaks and exhibiting significant up-regulation of eRNA expression were compared against random input sequences of varying size. Significantly enriched motifs in the HOMER database were calculated with p-value < 0.05.

Metagene analyses

Genome-wide exon-normalized counts were divided into four quartiles according to the average expression of genes (transcripts per million [TPM]) across replicates in either day 15 or term labor time points. Expression quartiles were used to plot the average H3K4me3 and H3K27ac signal and RNAPII coverage at each individual time point using ngs.plot [64]. Significantly up-regulated and down-regulated intron-corresponding reads in term labor samples were used to plot H3K4me3 and H3K27ac signal and RNAPII coverage at all time points using ngs.plot.

GO term analyses

A list of genes meeting our differential expression cutoff in the exon-based RNA-seq analyses was used as input in the EnrichR database. Outputs denoting genes’ associated molecular functions, cellular components, biological processes, and pathways were sorted according to p-value.