Increased transcriptome variation and localised DNA methylation changes in oocytes from aged mice revealed by parallel single‐cell analysis

Abstract Advancing maternal age causes a progressive reduction in fertility. The decline in developmental competence of the oocyte with age is likely to be a consequence of multiple contributory factors. Loss of epigenetic quality of the oocyte could impair early developmental events or programme adverse outcomes in offspring that manifest only later in life. Here, we undertake joint profiling of the transcriptome and DNA methylome of individual oocytes from reproductively young and old mice undergoing natural ovulation. We find reduced complexity as well as increased variance in the transcriptome of oocytes from aged females. This transcriptome heterogeneity is reflected in the identification of discrete sub‐populations. Oocytes with a transcriptome characteristic of immature chromatin configuration (NSN) clustered into two groups: one with reduced developmental competence, as indicated by lower expression of maternal effect genes, and one with a young‐like transcriptome. Oocytes from older females had on average reduced CpG methylation, but the characteristic bimodal methylation landscape of the oocyte was preserved. Germline differentially methylated regions of imprinted genes were appropriately methylated irrespective of age. For the majority of differentially expressed transcripts, the absence of correlated methylation changes suggests a post‐transcriptional basis for most age‐related effects on the transcriptome. However, we did find differences in gene body methylation at which there were corresponding changes in gene expression, indicating age‐related effects on transcription that translate into methylation differences. Interestingly, oocytes varied in expression and methylation of these genes, which could contribute to variable competence of oocytes or penetrance of maternal age‐related phenotypes in offspring.


| INTRODUC TI ON
Over the past 50 years, the proportion of women delaying childbirth in Western societies has steadily increased (Schmidt et al., 2012).
Age is the most important factor predicting female fertility potential, with the consequences of advanced maternal age ranging from karyotypic defects in oocytes, leading to aneuploidies and early miscarriage, to defects in placentation arising from compromised decidualisation (Hassold & Chiu, 1985;Lopes et al., 2009;Woods et al., 2017). Within the ovary, age affects both the number and quality of oocytes available for fertilisation (Vollenhoven & Hunt, 2018).
Female mammals are born with a fixed pool of oocytes to draw from during their reproductive lifespan. At birth, oocytes are arrested in the ovary at meiotic prophase I, where they will remain quiescent for months or years, depending on the species. Only a small subset will ever fully mature and ovulate, as the vast majority will undergo atresia before reaching ovulation. During the long quiescent phase, molecular changes in the oocyte result in a decline in oocyte developmental competence with advancing maternal age. Amongst the key molecular changes found to occur in oocytes from aged females are mitochondrial dysfunction, telomere shortening, cohesin dysfunction resulting in spindle instability and reduced stringency in spindle assembly checkpoints (Reviewed in Cimadomo et al., 2018;Vollenhoven & Hunt, 2018). In recent years, it has been proposed that transcriptional and epigenetic changes, especially DNA methylation, may contribute to compromised oocyte quality with advanced maternal age (Cimadomo et al., 2018;Ge et al., 2015;Marshall & Rivera, 2018). The changing environment to which an oocyte is exposed in the ageing ovary, such as altered hormone levels, and changes in energy and one-carbon metabolism, could affect gene expression and epigenetic processes that could contribute to the lower developmental capacity of oocytes with age (Ge et al., 2015).
Effects on epigenetic quality of the oocyte are also important because they can potentially be inherited into the embryo, manifesting as later developmental outcomes in offspring.
Gene expression in the oocyte is unique in that it comprises many oocyte-specific transcripts that not only regulate processes in the oocyte, but also in early embryo development. The growing oocyte accumulates and stores transcripts until it undergoes transcriptional arrest prior to ovulation. At this stage, transcription ceases until after fertilisation when it resumes during zygotic genome activation (Reviewed in Svoboda et al., 2015;Winata & Korzh, 2018). This means that all processes occurring during oocyte maturation, fertilisation and zygotic genome activation are regulated by the previously stored maternal transcripts and proteins, and changes in transcription or transcript storage capacity during ageing could therefore influence the developmental potential of the oocyte. Whole-transcriptome analysis using expression microarrays in both human and mouse metaphase-II (MII) oocytes has shown that typically a few hundred transcripts have altered abundance in the ageing oocyte and has identified genes involved in cell-cycle regulation, spindle assembly, oxidative stress and DNA damage as being frequently affected (Barragán et al., 2017;Grøndahl et al., 2010;Hamatani et al., 2004;Pan et al., 2008;Steuerwald et al., 2007). More recent studies using single-cell RNA-seq in human MII oocytes have found similar results, although the number of differentially abundant transcripts detected varied widely between studies (Barone et al., 2020;J.-J. Zhang et al., 2020).
Besides providing maternal factors for oocyte competence, transcription in the oocyte is also required for determining where de novo DNA methylation occurs in the genome. In the female germline, DNA methylation is globally reprogrammed during the course of gametogenesis. DNA methylation is almost completely erased in the primordial germ cells during embryonic development (Gkountela et al., 2015;Guibert et al., 2012;F. Guo et al., 2015; H. Guo et al., 2017;Seisenberger et al., 2012), after which the genome remains essentially unmethylated until de novo methylation initiates in the later stages of oocyte growth in the postnatal ovary (Gahurova et al., 2017;Hiura et al., 2006). De novo methylation during oocyte growth is an intricately regulated process that could be influenced both by internal and external factors, and changes in oocytes with age could influence the ultimate DNA methylation pattern established. For example, changes in nutrient supply and one-carbon metabolism could reduce the amount of substrate available for DNA methyltransferases (DNMTs) (Ge et al., 2015). In mouse oocytes, methylation is catalysed by DNMT3A together with its co-factor DNMT3L (Bourc'his, 2001; Kaneda et al., 2004) and is targeted to regions of active transcription, resulting in a bimodal distribution, unique to the oocyte, of broad unmethylated domains and methylated regions overlapping transcribed genes (Kobayashi et al., 2012;Veselovska et al., 2015).
This link between transcription and DNA methylation predicts that transcriptional changes in oocytes with age could affect DNA methylation by changing recruitment of DNMTs to chromatin and/ or by altering expression levels of DNMTs and related proteins.
Decreased expression of DNMT3A, DNMT3L as well as the maintenance methyltransferase DNMT1, has been reported in ageing mouse oocytes, both at the transcript and protein level (Hamatani et al., 2004;Yue et al., 2012). Whether these changes in DNMT expression are sufficient to affect DNA methylation levels has so far only been assessed on a global level using immunofluorescence (IF): Yue et al. found that the decrease in DNMT protein abundance correlated with a small but significant decline in global DNA methylation levels in MII oocytes and preimplantation embryos (Yue et al., 2012).
Successful epigenetic reprogramming is essential for the developmental capacity of the oocyte, as genomic imprinting is established in this phase. Genomic imprinting, defined as parent-of-origin-dependent gene expression in offspring conferred by differential methylation inherited from the gametes, is the clearest demonstration of long-term intergenerational epigenetic inheritance (Ferguson-Smith, 2011). DNA methylation established in oocytes also contributes to proper development of the trophoblast lineage (Branco et al., 2016). Imprinting has been reported to be stable in embryos and placentas from ageing mice, indicating correct establishment in oocytes and maintenance in the early embryo (Lopes et al., 2009). Consistent with this earlier report, recent analysis of the germline differentially methylated regions (gDMRs) that control the imprinted genes Snrpn, Kcnq1ot1 and H19 has shown correct levels of DNA methylation in oocytes irrespective of advanced maternal age, as well as normal methylation maintenance during preimplantation development (Kindsfather et al., 2019). On the other hand, as we recently showed, there are localised and coherent differences in DNA methylation in oocytes collected from prepubertal female compared with young adult mice, indicating that the DNA methylation landscape of the oocyte is not a completely fixed entity but is modulated during the life-course (Saenz- De-Juano et al., 2019). However, until now, no quantitative genome-wide assessment of methylation at single-base resolution in oocytes from aged females has been done.
Since ovaries from reproductively old mice contain fewer follicles, the number of oocytes in which to study gene expression and DNA methylation at a genome-wide level is severely limited. Such limitations have been overcome by recent advances in single-cell sequencing technologies . Here, using the single-cell M&T-seq approach (Angermueller et al., 2016), which allows the physical separation of RNA and DNA from single cells, we set out to generate RNAsequencing (scRNA-seq) and whole-genome bisulphite sequencing (scBS-seq) libraries from the same oocytes from young (12 weeks) and old (>40 weeks) mice. The profiled cells were fully grown germinal vesicle (GV) oocytes, which are already transcriptionally silent and have acquired most of the DNA methylation of mature oocytes including genomic imprinting marks. In accordance with previous studies, we identified transcript abundance changes in oocytes from aged compared to young mice affecting a relatively small number of genes. However, because of the single-cell resolution, we were uniquely able to identify increased molecular heterogeneity amongst older oocytes. While the majority of aged oocytes showed features suggestive of reduced developmental competence, a few retained a 'young-like' profile. Similar to the transcriptome, the methylome was globally similar irrespective of age, although specific changes in DNA methylation between oocytes from aged and young females were identified. Furthermore, because of the parallel profiling, we were able to link some transcript changes to methylation changes in the same oocyte, indicating that loss of gene expression fidelity in the aged oocyte can result in altered DNA methylation.

| Global patterns of gene expression in oocytes are associated with age and chromatin configuration
In order to investigate the effect of maternal age on global gene expression and DNA methylation in oocytes, we applied scM&T-seq (Angermueller et al., 2016) for joint transcriptome and methylome profiling from the same individual oocytes. Fully grown GV oocytes were isolated from the ovaries of three reproductively young, virgin female mice (C57BL/6Babr, 12 weeks old) and from seven reproductively old, virgin females (44-54 weeks old), an age at which females of this strain are approaching the end of their reproductive lifespan (Woods et al., 2017). scRNA-seq libraries that passed quality control filtering were obtained from 42 young and 45 aged oocytes. On average, each oocyte contained transcripts from 19,400 genes (Table   S1), from which 11,578 were shared across all oocytes. However, a lower number of expressed genes was observed in the aged group (difference in means=1,147 genes; Wilcoxon test, p = 4.2 × 10 −5 ; Figure 1a), indicating a reduced complexity in the transcriptome of GV oocytes from aged females.
To further explore the effect of age on global patterns of gene expression, we used principal component analysis (PCA) as implemented by the denoisePCA method in Scran (Lun et al., 2016). This programme partitions the total variance into its technical and biological components assuming that random noise is the major contributor to the variance of most genes in a scRNA-seq data set, but is uncorrelated across genes and should only be captured by later principal components (PCs) (Lun et al., 2016). In our data set, approximately 25% of the variance, which explained the first five PCs, was attributed to biological components. Age was identified as an important contributor to global variation in transcript abundance as it was associated with PCs 1, 2, 4 and 5, with PC2 being the most strongly associated with age (p = 1.01 × 10 −12 ) ( Figure 1b). PC1 as the main axis of variation was explained in part by age, although not as strongly as PC2, which suggested another important source of variation together with age. By comparing the loadings on PC1 between the age groups, we observed that aged oocytes showed a multimodal distribution compared to the unimodal distribution displayed by the young oocytes (Figure 1c), which suggested that the unknown source of variation was driving heterogeneity to a greater extent in the aged group of oocytes.
A key feature in determining maturation and competence of GV oocytes, and that might contribute a source of variation in oocytes at this stage, is chromatin condensation configuration: the transcriptionally active 'non-surrounded nucleolus' (NSN) configuration or the more mature and transcriptionally inactive 'surrounded nucleolus' (SN) configuration. This chromatin arrangement dictates developmental competence as fertilised oocytes with a NSN configuration mostly arrest at the two-cell stage (Monti et al., 2013). To assign chromatin configuration states in our data set, we classified the 87 oocytes according to the level of expression of genes reported to show at least a two-fold overexpression in SN oocytes compared to NSN oocytes (Ma et al., 2013). Twenty oocytes were found to express these transcripts at a higher level and were classified transcrip- To assess whether further comparisons between age groups would be confounded by chromatin configuration, we tested for a difference in proportions in assigned chromatin configurations between the two age groups. This assessment was important as it has been reported that the SN:NSN ratio amongst GV oocytes increases with age (Zuccotti et al., 1995). Out of the 20 SN oocytes, twelve corresponded to the young group and eight to the aged one.
The proportion of transcriptionally-like SN oocytes was not significantly different between groups (Fisher's exact test, p = 0.46), indicating that our populations are comparable for this parameter, and any observed age effects would not be confounded by oocyte developmental stage. However, we also observed that the number of SN oocytes was low compared to the numbers reported in the literature, especially in our aged group in which more than 80% of oocytes are expected to be SN. It is likely that the classification obtained using gene expression as a proxy does not reflect the actual chromatin state, but instead suggests that most aged oocytes expected to be SN express an immature NSN-like transcriptome. showing the enrichment of maternal effect genes when testing for differences across the four main clusters of oocytes. Each horizontal bar represents one maternal effect gene. The position of the bar along they-axis represents its ranking across all expressed genes tested for differential expression (1 being the most significant)

| Unsupervised clustering reveals both age-and chromatin configuration-driven oocyte populations with differential expression of maternal effect genes
Dimensionality reduction by t-distributed Stochastic Neighbor Embedding (t-SNE) revealed four well-separated clusters of oocytes in our data set ( Figure 1e). Clusters 1 and 2 comprised only NSN aged oocytes; cluster 3 mainly comprised young NSN oocytes plus a small number of SN oocytes; and cluster 4 purely SN oocytes regardless of age, which were closer to young NSN oocytes (cluster 3) than to aged NSN oocytes (clusters 1 and 2). This spatial relationship suggested a trajectory from an immature NSN transcriptome (cluster 1) to a mature SN one (cluster 4). Following this assumption, differential expression was tested using cluster number as a continuous variable to identify transcripts that change in abundance from old transcriptionally-like NSN oocytes to young transcriptionally-like NSN oocytes and, lastly, to SN oocytes (both young and aged). In total, 6,464 transcripts were identified as differentially abundant (Bonferroni adjusted p < 0.05; Table S2). Interestingly, the top ranked transcripts in this analysis were enriched for maternal effect genes (e.g., Padi6, Nlrp5, Gdf9, Tcl1, Dnmt1, Uhrf1, Bmp15 and Dppa3) ( Figure 1f and Table S2). To exclude that the observed differences in transcript abundance of maternal effect genes were an effect solely of chromatin configuration, the analysis was repeated using only the subset of oocytes assigned as NSN and the enrichment for maternal effect genes was also observed ( Figure S1).

| Aged oocytes exhibit greater gene expression variability
Increased heterogeneity in gene expression with age was suggested by the different distributions of the two age groups along the PC1 axis ( Figure 1c). To assess the increased variability at the gene level rather than at global patterns, and to identify those genes becoming variable with age, we tested for differential variability in gene expression between young and aged oocytes.
A total of 2,811 genes showed higher expression dispersion in aged oocytes, while only 74 decreased in variability compared to the young group, which suggests a more stochastic gene expression profile with increased age (Figure 2a and Table S3). We

| A fraction of aged oocytes expresses a younglike transcriptome
In addition to differences in variability, we also tested for differences in mean expression levels. We identified differences between young and aged oocytes at 560 transcripts (FDR <0.05; Figure 2c and  Interestingly, all of the predicted SN aged oocytes were also defined as young-like oocytes, although the young-like group also contained some aged oocytes assigned as NSN (Figure 2e).

| Reduced global DNA methylation but preserved methylation landscape in aged oocytes
A tight relationship between gene expression and DNA methylation has been identified in the oocyte, in which 85-90% of DNA methylation can be attributed to transcriptional events (Veselovska et al., 2015). In order to detect possible coordinate changes between gene expression and DNA methylation, the same GV oocytes from which scRNA-seq libraries were obtained were processed for scBSseq, from which 30 oocytes from the young group and 32 oocytes from the old group yielded libraries that passed quality control filtering (Table S6) itself. Oocytes are unusual in accumulating higher levels of non-CpG methylation (Shirane et al., 2013;Tomizawa et al., 2011). In contrast to the effect observed for CpG methylation, oocytes from aged females had increased non-CpG methylation levels (Wilcoxon test; p = 0.028; Figure 3f).
Amongst the functionally most significant sites of DNA methylation in the oocyte are the gDMRs that control imprinting.
Maternal gDMRs need to acquire methylation fully in oocytes to ensure imprinted monoallelic expression of imprinted genes in the embryo after fertilisation, and most imprinted genes are controlled by gDMRs methylated in oocytes (Schulz et al., 2010). We found that the 19 maternal gDMRs we assessed were highly methylated in oocytes irrespective of maternal age; in addition, three paternal gDMRs and two secondary DMRs were unmethylated similarly in both age groups (Figure 4a). Moreover, there was no evidence of increased variability in gDMR methylation in oocytes from aged females when analysed individually ( Figure S3).

| Age-associated differentially methylated regions
In addition to a global reduction in CpG methylation detected between the oocytes from young and old females, we wished to ascertain whether there were discrete or recurrent age-associated methylation differences. To do so, and to mitigate against sub-genome coverage of single oocyte data sets, we created randomised pseudobulk groups of 10 single oocyte data sets from the young and aged groups. Further, we generated 100 different combinations/permutations of the pseudobulk groups in order to limit the rate of false-positive signals that could arise from the composition of the groups (e.g., differences in coverage). Signals that appeared as significant at a false discovery rate of 5% in at least 95% of the combinations were carried forward as age-associated DMRs. By this approach, we identi-  (Table S7).

| Coordinate changes between gene expression and DNA methylation at a subset of genes
Because de novo methylation in oocytes is strongly associated with gene transcription, we sought to identify the overlap between and DEGs and DMRs to explore the relationship between gene expression and DNA methylation in relation to ageing in the oocyte. In total, 33 of the 166 DMRs overlapped with 23 different genes.
For these genes, a positive correlation between the difference in methylation and fold change in gene expression was observed (Figure 4d), indicating that in these cases the DNA methylation differences could be the result of transcriptional differences between young and old oocytes, given that transcription precedes gene body DNA methylation in oocytes as shown by previous studies (Kobayashi et al., 2012;Veselovska et al., 2015). To explore this in more detail, we estimated the correlation coefficient    (Figure 4f). At these genes, with one exception, we also found that non-CpG methylation correlated positively with the level of CpG methylation ( Figure S5), consistent with the notion that non-CpG methylation is also dependent on the activity of DNMT3A and obeys similar genomic targeting principles as CpG methylation (Shirane et al., 2013).

| DISCUSS ION
In this study, we have undertaken a parallel, single-cell assess- females. It is only through the application of single-cell profiling methods that such age-related heterogeneity can be revealed.
The culmination of oocyte growth is the transition from the transcriptionally active NSN to the inactive SN configuration, in which state the GV oocyte remains until GV breakdown (GVBD) and the resumption of meiosis. We did not determine the chromatin organisation of oocytes on collection, but inferred this by reference to existing RNA-seq data sets (Ma et al., 2013), thereby assigning oocytes as transcriptionally NSN-like and SN-like. NSN and SN oocytes have different developmental potential (Zuccotti et al., 1998), and the associated transcriptome differences are a reflection of this competence (Ma et al., 2013;Monti et al., 2013). We scored a surprisingly high proportion of oocytes from older females as NSN from their transcriptomes at an age when the ovary should contain few NSN oocytes (Zuccotti et al., 1995). On the other hand, it has been reported that more than a quarter of GV oocytes in old female mice cannot be classified as either NSN or SN, but display anomalous chromatin configurations (Manosalva & González, 2010). Together, these observations suggest that aged oocytes may have undergone the NSN-SN transition, but imperfectly and without the full transcriptome maturation that ensures developmental competence.
Consistent with this notion, we find that transcripts distinguishing SN oocytes from the old transcriptionally NSN-like oocytes were enriched in maternal effect genes that are critical for viability of the early embryo.
DNA methylation globally was similar between oocytes from young and aged females, although there was a slight reduction in total CpG methylation but increase in non-CpG methylation in the aged group. In addition, the characteristic DNA methylation landscape of alternating hyper-and hypomethylated domains was preserved. Consistent with our results, a previous study reported a slight reduction in 5-methylcytosine in MII oocytes from aged compared with young females as assessed by immunofluorescence (Yue et al., 2012). IF signal is likely to report predominantly on repetitive DNA sequences such as pericentromeric heterochromatin, a genomic fraction that BS-seq is less able to quantify owing to read mapping and methylation calling ambiguities of repetitive sequences. De novo methylation is progressive over oocyte growth, so a reduction in total CpG methylation could indicate slightly attenuated activity or fidelity of the de novo methylation complex in oocytes from aged females: quite substantial reductions in DNMT3A and DNMT3L protein levels in aged oocytes have been reported (Yue et al., 2012).
Alternatively, there could be elevated turnover of DNA methylation with age, although active demethylation is not known to occur in oocytes. A surprising observation was the reduced CpG methylation in oocytes assigned as SN-like.
Much of the DNA methylation the oocyte accumulates during its growth is of unknown functional significance, as oocytes deficient in DNMT3L or DNMT3A and therefore lacking most DNA methylation are able to be ovulated and fertilised with embryo progression until early post-implantation stages (Bourc'his, 2001;Kaneda et al., 2004). Functionally, the most important methylation in oocytes occurs at gDMRs of imprinted genes. Our results show that de novo methylation of imprinted gDMRs in aged oocytes is not impaired.
This accords with a recent study of the Snrpn and Kcnq1ot1 gDMRs in mouse oocytes that reported normal methylation acquisition irrespective of age, as well as no effect on maintenance of methylation to the blastocyst stage following fertilisation of oocytes from aged females (Kindsfather et al., 2019). Methylation in embryos and placentas at embryonic day 10.5 of aged mice was previously found to be normal at both maternal and paternal gDMRs (Snrpn, U2af1-rs1, Peg1, Igf2r, H19) (Lopes et al. 2009); methylation of the SNPRN gDMR is also reported to be unaffected by age in bovine oocytes (Mattern et al., 2016).
By using a very stringent approach that accounts for the relatively sparse and random coverage of scBS-seq data sets, we were able to identify 166 high-confidence age-related DMRs, the majority of which showed less methylation in oocytes from aged females.
DNA methylation in the oocyte is tightly linked to transcription events (Kobayashi et al., 2012;Veselovska et al., 2015), and for 33 of these DMRs we find concordant changes in transcript and methyla-  (Wang et al., 2004;Zhang et al., 2007).
Reduced transcription of Dync1i1 in aged oocytes could therefore be associated with impaired competence. Ppp3ca encodes calcineurin, which plays important roles in meiosis maturation in oocytes from many species (Zhang et al., 2019) and may contribute to regulation of GVBD in mouse oocytes . There were also changes in genes encoding activities potentially involved in RNA processing, such as the nuclear cyclophilin component PPIL3 (Rajiv & Davis, 2018), but whose relevance in oocytes has not yet been explored. It is possible that the heterogeneity in expression and methylation of these and other genes could contribute to variable developmental competence of the oocyte or predispose to variable phenotypes in offspring.
The mRNA mixture was reverse transcribed by incubation for 60 min at 42 °C followed by 30 min at 50 °C and 10 min at 60 °C.
The cDNA was amplified by PCR by adding 11 μl of 2 × KAPA HiFi HotStart ReadyMix and 1 μl ISPCR primer (2 μM; (Picelli et al., 2013(Picelli et al., , 2014). Samples were incubated in a thermocycler at 98 °C for 3 min, followed by 18 cycles of 98 °C for 15 s, 67 °C for 20 s, 72 °C for 6 min and finally 72 °C for 5 min. The amplified product was purified using Ampure XP beads with a 1:1 ratio and eluted into 20 μl of water.
In parallel, the lysate containing the DNA was purified using a 0.8:1 volumetric ratio of Ampure XP Beads (Beckman Coulter) and eluted into 10 μl of water to use for scBS-seq library preparation.

| Single-cell RNA-sequencing
Libraries were prepared from 100 to 400 pg of cDNA using the Nextera XT Kit (Illumina), per the manufacturer's instructions but with one-fifth volumes. All 88 single-cell RNA-seq libraries were pooled together and sequenced on the Illumina NextSeq platform to an average depth of 4.2 million reads (Table S1), using single-end 100 bp read-length settings.

| Single-cell BS-sequencing
Single-cell BS-seq (scBS-seq) libraries were prepared as previously described (Clark et al., 2017;Smallwood et al., 2014) (Kim et al., 2015) in single-end mode using a high penalty score for soft-clipping (--sp 1000,1000) and --dta for alignment reporting. Primary alignments were used for gene expression quantification over features of the mouse oocyte transcriptome (Veselovska et al., 2015). One oocyte was identified as an outlier in a distance matrix (cutoff =0.95-quantile) and was excluded from downstream analyses.

| scBS-seq data processing
A total of 80 scBS-seq libraries (39 young and 41 aged oocytes) were processed for analysis. The first 6 bp containing the N portion of the random primers, adapters and bases called with poor quality (Phred score <20) were removed using Trim Galore version 0.4.4 (Krueger, 2017) in single-end mode. Reads with a minimum length of 20 bp after trimming were retained for downstream processing.
Bismark version 0.18.2 (Krueger & Andrews, 2011) was used for read alignment, deduplication and methylation calling. Alignment to GRCm38 was performed in single-end non-directional mode followed by deduplication (deduplicate_bismark) and methylation calling (bismark_methylation_extractor). scBS-seq libraries having a mapping efficiency <10%, fewer than 500,000 CpGs covered, or greater than 50% methylated CpGs were discarded. An additional quality control we applied as a means of identifying possible somatic cell (granulosa cell) contamination that might have been missed by visual inspection was to assess the average methylation of X chromosome CpG islands, with the rationale that CpG islands on the inactive X chromosome in female somatic cells are highly methylated, whereas during oocyte development the X chromosomes have undergone reactivation. All oocyte libraries had average X chromosome CpG island methylation rates <10%, compared with sequenced single granulosa cells with X chromosome CpG methylation >20%. In total, 30 young and 32 aged oocytes were kept for downstream analyses. On average, 3,203,801 CpGs per oocyte were interrogated (Table S6).

| Differential expression and differential variability
Expression levels were normalised using size factors, and differential expression was tested using the R package DESeq2 (Love et al., 2014) at a false discovery threshold (FDR) of 5%. Differential variability was tested using the R package BASiCS (Vallejos et al., 2015) with an expected false discovery rate of 10%.

| Differential methylation analysis and cell clustering
In the fully grown GV oocyte, global DNA methylation is around 40% with methylation states shared along broad clusters of unmethylated and methylated CpG sites. Veselovska et al., 2015, identified and termed clusters as hypermethylated (75-100% methylation) and hypomethylated (0-25% methylation) domains with median sizes of 20.9 kbp and 24.9 kbp, respectively. Here, DNA methylation was quantified over these distinctive domains and the interdomains between these two in random groupings of 10 oocytes (pseudobulk groups). Differential methylation was tested by fitting a weighted logistic regression model. To limit the rate of false-positive signals due to the composition of the groups, one hundred different combinations of oocytes belonging to the two groups were computed. Signals called significant at a false discovery rate of 5% in at least 95% of the combinations were called age-associated differentially methylated regions. Single oocytes were clustered into two groups based on their patterns of DNA methylation using Melissa (Kapourani & Sanguinetti, 2019), a Bayesian method for imputation and clustering of single-cell methylomes.

| Gene expression and DNA methylation pairwise correlations
Pearson's correlation was used to identify linear relationships between DNA methylation and gene expression levels. For this analysis, DNA methylation was quantified at gene bodies over single oocytes.

ACK N OWLED G EM ENTS
We should like to thank members of the

EH-P was supported by the National Researcher Training Program of
Colciencias 567-2012.

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests.

AUTH O R CO NTR I B UTI O N S
MH, GK conceived and supervised the study; EH-P generated sequencing data sets with assistance from SJC; JC-F performed bioinformatic analysis with input from HD, CH; JC-F, HD, EH-P, GK wrote the manuscript with input from all other authors.

DATA AVA I L A B I L I T Y S TAT E M E N T
Sequence data that support the findings of this study have been deposited under accession code GSE154370 in the Gene Expression Omnibus database at https://www.ncbi.nlm.nih.gov/geo/query/ acc. cgi?acc=GSE15 4370.