Epigenetic changes occur at decidualisation genes as a function of reproductive ageing in mice

Adaptation of the uterus to an implanting embryo is essential for a successful pregnancy. This process, termed decidualisation, is chiefly initiated by the ovarian steroid hormones oestrogen (E2) and progesterone (P4) (Wang and Dey, 2006; Adams and DeMayo, 2015). In response to this hormonal action, the uterus acquires a receptive state, ready to accept a blastocyst for implantation. Embryo attachment elicits a transformation of the primed uterine stromal cells to proliferate and then differentiate, which in the mouse typically results in the acquisition of a polyploid DNA content (Das et al., 1999). These hallmarks of the decidualisation reaction are accompanied by a swollen appearance of the uterus and an increase in uterine weight (Sharma et al., 2006; Lee et al., 2007).

We have recently shown that, in mice, this decidualisation reaction is markedly affected by maternal age (Woods et al., 2017). Thus, advanced maternal age is associated with a dramatically attenuated decidualisation response, which is the main cause of the subsequent developmental problems observed in pregnancies of older females. This defect persists in vitro under optimal decidualisation conditions; uterine stromal cells isolated from older females exhibit a significantly slower or blunted response to E2 and P4 and fail to upregulate key marker genes of decidualisation such as Bmp2, Hoxa10 and Hand2.

Ageing has been widely reported to correlate with changes in the DNA methylation landscape as a result of imperfect maintenance of DNA methylation marks (Issa, 2014). Consequently, sites that are normally heavily methylated lose methylation, and others that are unmethylated may gain some methylation (Pal and Tyler, 2016; Ciccarone et al., 2018). Generally, this loss of methylation fidelity results in the relative hypomethylation of retrovirally derived sequences, repeat elements and centromeres as well as loss of gene body methylation, whereas normally unmethylated CpG islands (CGIs) tend to acquire DNA methylation (Pal and Tyler, 2016; Ciccarone et al., 2018). Measuring DNA methylation at a few specific genomic sites can even predict chronological age in multiple tissues, including in the endometrium (Horvath, 2013; Olesen et al., 2017; Stubbs et al., 2017).

In the mouse, expression of the DNA methyltransferases Dnmt1 and Dnmt3a increases at implantation sites between embryonic day (E) 4-8 (Gao et al., 2012). Both genes are also dramatically upregulated in artificially induced deciduomas, in which the hormonally primed uterus undergoes a decidualisation reaction in the absence of an implanting embryo. This upregulation of de novo and maintenance DNA methyltransferases suggests that DNA methylation may be an important epigenetic regulator of gene expression during decidualisation. In fact, inhibition of DNA methylation by 5-azacytidine treatment does not affect implantation but dramatically dampens the decidualisation response, with reduced uterine stromal cell proliferation and a failure to adequately upregulate decidualisation markers (Gao et al., 2012). This phenotype is similar to the abnormal decidualisation response we observed in aged pregnant females (Woods et al., 2017). This correlation in phenotypes raises the question of whether maternal ageing is accompanied by DNA methylation changes in the endometrium that have a negative impact on its capacity to decidualise.

Consequently, we set out to characterise the DNA methylation changes that take place in the endometrium as it differentiates to form the decidua during normal pregnancy, and to investigate the impact of advanced maternal age on the uterine and decidual epigenome.

To gain insight into the overall expression dynamics of epigenetic modifiers during decidualisation in vivo, we first analysed RNA-seq data from uterine samples during the implantation window at E3.5, and of dissected deciduas of E9.5, E10.5, E11.5 and E12.5 conceptuses developed in young control females (Woods et al., 2017). T-distributed stochastic neighbour embedding (t-SNE) analysis clustered the samples into four main groups, namely E3.5 uterus, E9.5 and E10.5 decidua, E11.5 decidua, and E12.5 decidua (Fig. 1A), showing the progression of transcriptional changes that accompany decidual development. As anticipated, the greatest transcriptomic differences were seen between the E3.5 and E12.5 time points (Fig. 1A). Among the differentially expressed genes were a number of epigenetic modifiers, notably those involved in DNA methylation dynamics (Fig. 1B). Most of these, including the ten-eleven translocation (Tet) enzymes, the DNA methylation maintenance surveillance factor Uhrf1, and multiple members of the methyl-CpG-binding domain (Mbd) family, were most highly expressed at the early stages (E3.5 and E9.5) and declined towards the later stages of decidual maturation (E11.5-E12.5) (Fig. 1B). The DNA methyltransferases Dnmt1, Dnmt3a and Dnmt3b were transiently upregulated at E9.5, in line with previous observations (Gao et al., 2012). The most dramatic change was observed for Tet1, a 5-methylcytosine (5mC) dioxygenase enzyme involved in DNA demethylation, which was highly expressed at E3.5 but strongly downregulated thereafter.

The pronounced shift in expression of multiple factors involved in 5mC dynamics between the hormonally primed E3.5 uterus and later stages of decidualisation prompted us to compare the global DNA methylation landscapes at E3.5 and E11.5 in more detail, using DNA immunoprecipitation with a 5mC antibody followed by high-throughput sequencing (MeDIP-seq). When quantifying DNA methylation coverage in 2 kb tiling probes across the genome, the E3.5 and E11.5 samples exhibited very similar DNA methylation landscapes, as expected. However, a subset of regions were markedly differentially methylated, with the majority (90.6%) being hypermethylated at E11.5 (Fig. 1C). To better characterise the methylation changes between these stages, we used the MEDIPS program to identify genomic regions that exhibited significant differential methylation levels. This analysis revealed 8991 primary sequence stretches that exhibited a differential enrichment of 5mC, of which 8945 (99.5%) gained and 46 (0.5%) lost methylation between E3.5 and E11.5 (Fig. 1D). Merging of neighbouring tiles resulted in the identification of 3802 differentially methylated regions (DMRs) between E3.5 and E11.5 (Table S1) that were used for further analysis.

The identified DMRs were predominantly located in distal intergenic regions (∼50%), defined as >3 kb upstream or downstream of the nearest gene. Approximately 40% were in introns and 5% at promoter regions (Fig. 2A). Interestingly, analysis with the Genomic Regions Enrichment of Annotations Tool (GREAT) indicated an association of the identified DMRs with genes related to proliferation, growth and endometrial maturation (Fig. 2B). This included terms related to the regulation of stem cell maintenance, TGFβ-, NODAL-, VEGF- and WNT-signalling, and mesenchymal-to-epithelial transition, a well-characterised process in the later stages of decidualisation (Patterson et al., 2013; Zhang et al., 2013; Yu et al., 2016). Genes involved in the regulation of cholesterol transport, a substrate required for the synthesis of E2 and P4, and E2-signalling pathway components, were also enriched for DMRs and may therefore be subject to regulation by DNA methylation during decidualisation.

Integrating the methylome and transcriptome datasets revealed that a significant proportion (36%) of DMRs overlapping or upstream of a gene (within ≤5 kb upstream of a transcription start site or ≥1 bp overlap with a gene) were associated with differentially expressed transcripts between E3.5 and E11.5 (Fig. 2C). This considerable overlap may suggest that DNA methylation changes laid down during pregnancy progression contribute to the transcriptional regulation of these genes. This regulation may be directed through the action of the progesterone receptor (PGR), as a significant proportion of genes associated with a DMR were also affected by Pgr ablation (Fig. 2D; Jeong et al., 2005). Notably, hypermethylated DMRs were found at genes with known roles in placental development and decidualisation. For example, DMRs that gain methylation over gestation were detected in the promoter region of the P4-regulated gene Hif1a and within introns of the oestrogen receptor Esr1, the cyclase-associated actin cytoskeleton regulatory protein Cap2 and the immunophilin Fkbp5, which acts as a chaperone for PGR (Fig. 2E; Fig. S1). In all of these cases, the increased methylation at the identified DMRs correlated with decreased gene expression at later stages of decidualisation, suggesting that the intron methylation affects gene regulatory enhancer elements at these loci.

Overall, our approaches identified a set of uterine DMRs that overwhelmingly gain DNA methylation in the decidua during gestation. These DMRs are significantly associated with decidualisation genes, in particular those regulated by P4. This is consistent with previous reports that showed an increase in DNA methyltransferase expression after implantation (Gao et al., 2012). These data imply that a cohort of P4-responsive decidualisation genes are epigenetically regulated; after their initial transcriptional burst during early stages of decidualisation, they are methylated and repressed in mature decidual cells during later stages of gestation. This finding is in line with the reported failure in the progression of decidualisation upon exposure to the DNA methylation inhibitor 5-azacytidine (Gao et al., 2012).

The enrichment of DMRs at P4-regulated genes was of particular interest as these genes are highly sensitive to age-related changes (Woods et al., 2017). To assess whether alterations to the DNA methylation profile may affect the regulation of these genes and thereby account for the decidualisation defect of older females, we performed the MeDIP-seq analysis also on E3.5 uterine and E11.5 decidual samples from females aged 41-49 weeks. Overall, the methylation dynamics were very similar in the aged samples compared with young, with most changes entailing a gain of methylation between E3.5 and E11.5 at the vast majority of differentially enriched loci (5198/5610; 93%) (Fig. S2A). The 5mC enrichment distribution across various genomic features was also globally very similar (Fig. S2B). Contrary to the well-reported loss of gene body methylation as a function of age, no statistically significant difference was detected between young and aged uterine and decidual samples across gene bodies. However, an exception to this overall similarity was evident at CGIs that exhibited globally higher DNA methylation levels in the aged samples (Fig. S2B).

We thus turned our attention specifically to CGI methylation changes. When comparing the E3.5 and E11.5 developmental stages within each age group, we identified 302 and 1212 differentially methylated CGIs in the young and aged samples, respectively. As observed with the whole genome-covering probes, the majority of these developmentally dynamic CGIs gained methylation with the progression of decidualisation in both young and aged samples (Fig. S2C). However, the overall gain of 5mC was more pronounced in the aged cohort (Fig. S2C), resulting in the greater number of differentially enriched CGIs in aged than in young samples. CGIs are therefore the genomic feature most susceptible to methylation changes during reproductive ageing.

Intriguingly, the most conspicuous difference was observed when assessing the smaller fraction of developmentally dynamic CGIs that lose methylation during the progression of decidualisation. This cohort of CGIs was far more highly 5mC-enriched in aged than in young E3.5 uteri (Fig. S2C). Directly comparing CGI methylation of E3.5 young and aged samples revealed 41 differentially methylated CGIs, 40 of which were hypermethylated in aged samples (Fig. 3A). We independently verified the differential methylation pattern at the Fzd2 CGI by bisulphite sequencing (Fig. S2D). When assessing whether promoter-associated differentially methylated CGIs affected expression of the corresponding genes (Woods et al., 2017), we found that 75% of these genes were downregulated in aged samples (Fig. 3B). This included pivotal regulators of decidualisation, notably the WNT receptor Fzd2 and the Hoxa10/Hoxa11 gene cluster (Fig. 3C,D; Fig. S3A, Table S3) (Hsieh-Li et al., 1995; Benson et al., 1996; Lim et al., 1999; Hayashi et al., 2009). Thus, age-associated CGI hypermethylation likely prevents the induction of these genes to achieve normal expression levels during early decidualisation.

CGI hypermethylation was correlated with a significant downregulation of Tet1 expression in E3.5 uteri of aged females compared with young (Fig. 3E, Fig. S3B, Tables S3 and S4). These lower Tet1 levels corresponded with a strikingly reduced staining intensity for 5-hydroxymethylcytosine (5hmC) in aged uteri, specifically in glandular epithelial cells (Fig. 3F). Although 5mC staining intensities were not appreciably altered or only mildly increased in aged samples by visual inspection, we observed significant changes in the ground state expression of the DNA methylation machinery components Dnmt1, Dnmt3b, Uhrf1 and Trim28 (also known as Kap1), which were all expressed at significantly higher levels in uterine stromal cells from aged compared with young females (Fig. S3C). Collectively, this deregulation of epigenetic modifiers may well account for the gain in CGI methylation that is evident in the uteri of aged females, and may be a major contributor to the frequent developmental demise of embryos in these pregnancies.

Taken together, in the temporal context of reproductive ageing, early epigenetic changes seem to affect CGI methylation specifically whereas other epigenetic hallmarks of organismal ageing, such as loss of gene body methylation, are not yet evident. We find that a number of CGIs at crucial decidualisation genes are hypermethylated at E3.5 in aged female uteri, which may interfere with appropriate gene induction thus impeding the progression of early postimplantation development. The most prominent examples of genes with differentially methylated promoter CGIs were Hoxa10, Hoxa11 and Fzd2. Hoxa11 deficiency causes a hypoplastic uterine phenotype with fewer glands whereas Hoxa10 ablation causes reduced stromal cell proliferation and overall decidualisation failure, with both knockouts resulting in female infertility (Hsieh-Li et al., 1995; Satokata et al., 1995; Benson et al., 1996; Lim et al., 1999). WNT signalling through FZD receptors is also highly dynamically regulated during decidualisation and has important roles in supporting this cellular transformation process (Hayashi et al., 2009). Overall, here we show that the uterine epigenome of older females undergoes age-related alterations that likely interfere with the robust response to steroid hormones required for normal decidualisation. These epigenetic changes may explain the significant decidualisation defects and consequent developmental failures associated with declining reproduction in older female mice. Although it may be tempting to speculate that treatment with methylation inhibitors such as 5-azacytidine could prevent or reset some of these changes, the difficulty will be to direct them to the correct target loci. However, since DNA methylation changes also correlate with gene expression dynamics in the human endometrium during the endometrial cycle (Kukushkina et al., 2017), such age-induced epigenomic rearrangements may be of significant importance for the decline in uterine function and fertility linked to advanced maternal age in humans.

C57BL/6Babr mice were used throughout this study. All animal experiments were conducted in full compliance with UK Home Office regulations and with the approval of the local animal welfare committee (AWERB) at The Babraham Institute, and with the relevant project and personal licences in place. Mice were housed in standard individually ventilated cages with environmental enrichment as per best animal husbandry practice guidelines, under 12 h light-dark cycles. Mice were fed CRM (P) diet (Special Diet Services) ad libitum and received seeds (e.g. sunflower or millet) as part of their environmental enrichment. ‘Young’ females were 8-20 weeks old and ‘Aged’ females were 41-49 weeks old. Timed matings were set up between virgin females and stud males of 8-16 weeks of age, counting the morning of the vaginal plug as E0.5. Pregnant females were dissected at the gestational age indicated.

RNA and DNA were isolated simultaneously from whole frozen tissue using the AllPrep DNA/RNA Mini Kit (Qiagen, 80204). RNA from uterine stromal cells was isolated using TRI Reagent (Sigma-Aldrich, T9424), following the manufacturer's instructions. For RT-qPCR, RNA was reverse transcribed using RevertAid H Minus Reverse Transcriptase (ThermoFisher, EP0451) to generate cDNA and PCRs performed using QuantiNova SYBR Green PCR Kit (Qiagen, 208056) on a Bio-Rad CFX96 cycler. Bisulphite conversion was performed with the EpiTect Bisulfite Kit (Qiagen, 59104) according to the manufacturer's instructions. PCR products were cloned into the pGEM-T Easy Vector System (Promega, A1360) and then sequenced. Primer sequences are provided in Table S2.

Uterine stromal cells were isolated largely as previously described (Afonso et al., 2002). In brief, uteri from E3.5 mice were slit longitudinally and disaggregated with 2.5% pancreatin (Sigma-Aldrich, P3293) and 0.5% trypsin (type III) (Sigma-Aldrich, T4799) in Hank's balanced salt solution (HBSS). After vortexing, the medium was removed and the remaining tissue was washed twice in HBSS (discarding luminal epithelial cells) and incubated for 30 min at 37°C in HBSS containing 0.007% collagenase (Thermo Fisher Scientific, 17104-019), 0.02% DNase I (Roche, 10104159001) and 0.008% protease (Sigma-Aldrich, P8811) with vortexing every 5 min. Dissociated tissues were then triturated, filtered through a 70 μm cell strainer (BD Biosciences) and spun at 1000 g for 5 min at 4°C. Cell pellets were resuspended in phenol-red free Dulbecco's modified Eagle's medium: Nutrient Mixture F-12 (DMEM/F-12) (Thermo Fisher Scientific, 21041-025) plus 10% foetal bovine serum, 1 mM sodium pyruvate (Thermo Fisher Scientific, 11360-039), 1× antimycotic/antibiotic (Thermo Fisher Scientific, 15240-062) and 50 µM 2-mercaptoethanol (Thermo Fisher Scientific, 31350-010).

RNA-seq libraries were generated using the SureSelect Strand-Specific RNA Library Preparation kit (Agilent, G9691A) or the ScriptSeq V2 RNA-seq Library Preparation Kit (Illumina, SSV21124). Libraries were quality-controlled using the Bioanalyzer 2100 (Agilent) and quantified using the KAPA Library Quantification Kit (KAPA Biosystems, KK4824). Indexed libraries were pooled and sequenced on an Illumina HiSeq2500 sequencer, using a 100 bp single-end protocol. FastQ data were trimmed for adaptors using Trim Galore! v0.4.4, before mapping to the Mus musculus GRCm38 genome assembly using HISAT2 v2.1.0 (uterine stromal cells), or TopHat v2.0.12 (whole uterus and decidua). Data were quantified using the RNA-seq quantitation pipeline in SeqMonk (www.bioinformatics.babraham.ac.uk), and normalised according to total read count (read per million mapped reads, RPM). Differential expression was calculated using DESeq2 with P<0.05 and adjusted for multiple testing correction using the Benjamini–Hochberg method. For stringent difference analysis, an intensity difference filter of ±2 (after log transformation) was applied in SeqMonk. Pgr KO expression data were retrieved from NCBI Gene Expression Omnibus (www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE39920).

Sonicated DNA (1.5 µg) was subjected to methylated DNA immunoprecipitation (MeDIP) as described previously (Senner et al., 2012). Indexed libraries were pooled and sequenced on an Illumina HiSeq2500 sequencer, using a 50 bp paired-end protocol. FastQ data were trimmed and aligned to the Mus musculus GRCm38 genome assembly using Bowtie v2.3.2. DMRs were identified by either the MEDIPS package in R, using a window size of 200 bp and Bonferroni P-value cut-off of 0.01, or using LIMMA statistics with p_Bonf<0.05 for replicated datasets within SeqMonk. To detect meaningful enrichment differences, minimum read count cut-offs were set per probe and sample set: for 2 kb running windows ≥10 reads per probe, and for CGIs ≥5 reads per probe in at least one of the replicate sample sets. Neighbouring windows were merged, and those located in blacklist regions were removed (ENCODE project, accession ENCSR636HFF). The feature distribution of DMRs was determined using ChIPseeker in R. The enrichment of methylated fragments in different genomic features was determined by dividing the number of reads in each feature by the percentage of the genome occupied by that feature, followed by log transformation.

Heatmaps were generated using heatmap.2 in R. Principal component analysis was performed using DESeq2 rlog-normalised RNA-seq data on read counts using the top 500 most variable genes. The principal component analysis results were plotted using prcomp in R. t-SNE plots were generated using SeqMonk (www.bioinformatics.babraham.ac.uk).

Volcano plots and wiggle plots were generated in R using ggplot2 and core R plotting functionality, respectively. Gene Ontology (GO) and gene set enrichment analyses using the MSigDB (software.broadinstitute.org/gsea/msigdb) were performed on genes found to be significantly differentially upregulated or downregulated, against a background list of genes consisting of those with more than ten reads aligned. GO terms with a Bonferroni P-value of <0.05 were found using DAVID (Huang da et al., 2009) and GREAT (McLean et al., 2010). Venn diagrams were plotted using BioVenn (www.cmbi.ru.nl/cdd/biovenn/). Two-way ANOVA and Fisher's exact test were performed using GraphPad Prism and R, respectively.

E3.5 uteri of young and aged females were fixed in 4% paraformaldehyde, processed for routine paraffin histology and embedded in the same paraffin block. Sections (7 µm thick) were used for staining. Sections were deparaffinised and rehydrated according to standard protocols. A denaturation step was performed using 2N HCl for 20-40 min. Sections were blocked in PBS, 0.1% Tween 20 and 0.5% Bovine Serum Albumin, and incubated with an anti-5hmC antibody (Active Motif, 39769; 1:1000). Detection was carried out with the appropriate Alexa Fluor 488-conjugated secondary antibody (1:600; Life Technologies, A-21206). Nuclear counterstaining was performed with 4′,6-diamidine-2′-phenylindole dihydrochloride (DAPI) (MilliporeSigma, 10236276001). Images were captured using an Olympus BX41 Fluorescence Microscope.

Source data for Fig. S3A,B and Fig. 3E are provided in Tables S3 and S4, respectively.

We thank Dr Sarah Burge for expert help with sequencing data analysis, and Kristina Tabbada for Illumina high-throughput sequencing.