Seeds from wild-type and insertional ago4-6 (SALK_071772), ago9-3 (SAIL_34_G10), rdr6-15 (SAIL_617) mutants (all in a Col-0 background) were disinfected with 40% chlorine for 10 minutes under constant agitation, rinsed three times with sterile distilled water, planted in petri plates containing 0.5% Murashige and Skoog (MS) medium, vernalized for three days at 4°C under dark conditions, and germinated at 21°C in long daylight conditions (16h light/8h dark). Seedlings were transferred to soil (3:1:1 Sunshine substrate: vermiculite: perlite), grown under greenhouse conditions and fertilized with slow release 14-14-14 fertilizer (1.84kg/m³). Homozygous ago9 and rdr6 individuals were identified by PCR genotyping using primer combinations listed in Supplementary Table S1 . In the case of rdr6-15, amplified DNA products were digested with MscI (New England Biolabs, USA) to distinguish wild-type and mutant alleles.

For each biological replicate, 200 gynoecia of precisely 0.5 mm in length – and containing ovules at stages 1 and 2 (Rodríguez-Leal et al., 2015) - were isolated and collected under a Leica M80 stereomicroscope. Stage 9 and 10 flower buds (Schneitz et al., 1995) were isolated using emasculation tweezers (Dumont #5 Tweezers; Electron Microscopy Sciences) and immobilized on a slide with double-sided adhesive tape ( Supplementary Figure S1 ). Gynoecia measuring 0.5 mm were isolated and collected using 1 ml insulin syringes (Terumo, USA), immediately placed 1.5 ml microcentrifuge tubes in liquid nitrogen, and stored at -80°C before DNA extraction.

DNA extraction was performed with cetyltrimethyl ammonium bromide (CTAB) buffer adding a final incubation step with RNase (modified CTAB protocol on the base of Doyle and Doyle, 1987). For bisulfite conversion of unmethylated cytosines to uracil, we used the EpiTect Fast Bisulfite Conversion kit (Qiagen, #59802). Library construction was performed using the Ovation Ultralow Methyl-Seq Library System kit (Nugen, #0336), following two rounds of bisulfite conversion (Walker et al., 2018). Two biological replicates were generated for each genotype. Wild-type replicate #1 (Wt_R1), ago4 replicate #1 (ago4_R1), ago9 replicate #1 (ago9_R1), and ago9 replicate #2 (ago9_R2) libraries were generated from 59 to 120 ng of total DNA that was treated with bisulfite as described above, and subsequently amplified with 10 cycles of PCR. Wild-type replicate #2 (Wt_R2) and ago4 replicate #2 (ago4_R2) libraries were generated from 70 ng of total DNA treated with bisulfite and subsequently amplified with 11 cycles of PCR. Finally, rdr6 replicate #1 (rdr6_R1) and rdr6 replicate #2 (rdr6_R2) were generated from 20 and 29 ng of total DNA, and amplified with 13 and 12 cycles of PCR, respectively. All libraries were sequenced using the single-end Illumina HiSeq mode.

Library adapters, the initial nine nucleotides (nt) from the 5’ end, short sequences of less than 20 nt, and low-quality reads (Phred score less than 20) were discarded from all methylomes using the TrimGalore version 0.4.2 (Krueger F. https://github.com/FelixKrueger/TrimGalore). Resulting reads were mapped to the Arabidopsis TAIR10 genome after in silico converting all sense cytosines to thymine, and all antisense guanines to adenine, allowing one mismatch per sequence and keeping only reads showing a single best alignment when using the Bismark version v0.19.0 (Krueger and Andrews, 2011). Duplicate reads were removed using the Samtools version 1.9, rmdups option (Danecek et al., 2021). Reads that showed a perfect match to the unconverted TAIR10 genome were eliminated by considering they represent DNA residues unaffected by bisulfite reactions (less than 1.1% of all sequences in all cases). In each of the eight resulting methylomes (two replicates per genotype), the total number of methylated and unmethylated cytosines were determined in each context (CG, CHH, and CHG) using the Bismark-methylation-extractor tool (Krueger and Andrews, 2011) in both sense and antisense genomic strands. The single nucleotide methylation frequency was calculated by estimating the number of cytosines over the number of cytosines and thymines at each individual genomic site. Genomic bins of 50 nt were considered to estimate the methylation frequency corresponding to the ratio of the total number of methylated cytosines by the total number of cytosines included in each bin.

Genes and transposable elements (TEs) annotated in TAIR10 were aligned at the 5’ end or the 3’ end, discarding from the analysis either 1500 bp (genes) or 250 bp (TEs), from the end opposite to the one used for alignment (Ibarra et al., 2012). Methylation levels were averaged across 100 bp intervals using the number of cytosines and thymines sequenced at the single nucleotide level for each sequence context in a given methylome. Finally, average methylation levels were calculated over a region encompassing 5 kb upstream and 5 kb downstream of the alignment site. Resulting methylation profiles were compared to available datasets for seedling (Zhang et al., 2013), rosette leaf (Stroud et al., 2014), central cell (Park et al., 2016), male meiocyte (Walker et al., 2018), and ovule (Zhou et al., 2022).

A potential differentially methylated region (DMR) was defined as a genomic segment with a minimum length of 100 nt that contained at least 20 sites corresponding to cytosines, and showed at least 10% (CHH context), 20% (CHG context) or 40% (CG context) methylation difference between wild-type and a specific RdDM mutant. Potential DMRs of less than 100 nt were discarded. The 100 nt window is a standard arbitrary parameter that accounts for frequent limited cytosine coverage in 50 nt genomic windows, increasing the probability of detecting significant methylation differences over a larger window, but maintaining the length within the range of a single nucleosome to prevent uncertain biological interpretations. Two 50 nt windows separated by more than 100 nt were considered to represent independent DMRs. The statistical value of all DMRs was determined by applying a Fisher’s exact test with a p-value of less than 0.001. DMRs were considered to map to a promoter, gene body, or TEs, when at least 50% of the DMR overlapped with either 200 bp upstream of the first exon (promoter), with a region encompassing the first and last exon of an annotated transcriptional unit (gene), or with 50% of a previously annotated TE. DMRs were considered to map to gene elements if 50% of their length overlapped with a gene unit, considering 200 bp upstream of the first exon as part of the gene unit.

We confirmed the reproducibility of biological replicates for all four genotypes by conducting a principal component analysis (PCA) using cytosine sites in a CG context that were covered in all eight methylomes (two biological replicates per genotype), for a total of 4.1x10⁵ sites. As illustrated in Supplementary Figure S2 , components 1 and 2 explain 27% and 25% of the variation between replicates, respectively. All samples had a reproducible tendency two group within their corresponding genotype (wild-type, ago4, ago9, or rdr6), with wild-type and ago9 methylomes showing less variability than methylomes of ago4 and rdr6.

For genes and TEs, we compared the global methylation profile of the pre-meiotic wild-type gynoecium to methylation profiles from other sporophytic organs or gametophytic cells. These included seedlings (Zhang et al., 2013), rosette leaves (Stroud et al., 2014), the central cell of the female gametophyte (Park et al., 2016), and the male meiocyte (Walker et al., 2018). These results are illustrated in Figure 1 . For TEs, wild-type methylation levels in the CG context were higher in the central cell and the male meiocyte than in rosette leaves or seedlings, consistent with the prior knowledge that reproductive cells have higher CG methylation than somatic tissues (Hsieh et al., 2016; Park et al., 2016). Strikingly, wild-type TEs methylation levels in the CG context of pre-meiotic gynoecia were higher than expected, reminiscent of those prevailing in the male meiocyte rather than those of sporophytic organs ( Figure 1 ). Methylation levels of TEs in the CHG context were also similar to gametophytic cells, significantly higher than those of sporophytic organs. And for the CHH context, TE methylation levels in pre-meiotic gynoecia are also similar to those of the male meiocyte, but lower than those found in seedlings, rosette leaves and the central cell. In the case of genes in the CG context, methylation levels were quite similar for all organs and cells analyzed ( Figure 1 ). The overall methylation levels for TEs and genes had little difference from wild type in pre-meiotic gynoecia of homozygous ago4, ago9, and rdr6 individuals, suggesting that they are not majorly involved in the establishment of high methylation levels observed in TEs and gene bodies ( Supplementary Figure S3 ).

To determine if these results reveal a tendency to establish gametophytic methylation levels early during female reproductive organ development, we compared methylation levels of TEs and genes in pre-meiotic gynoecia and ovules, using a recently available dataset generated by Zhou et al. (2022). As illustrated in Supplementary Figure S4 , methylation levels in TEs and genes are similar in both organs, in both cases reminiscent of those prevailing in male meiocytes, confirming that drastic changes in methylation levels for CHH, CG and CHG contexts occur at the sporophytic level, and prior to the alternation of generations within the ovule primordium.

The wild-type pre-meiotic gynoecia is characterized by overall levels of genomic methylation that reach 22.37%, 8.63%, and 1.14% for the CG, CHG, and CHH context, respectively ( Supplementary Table S2 ). As expected, the pre-meiotic gynoecium of mutants affected in the RdDM pathway show global methylation levels lower than those found in wild-type gynoecia. As expected, in all three mutants analyzed, the most drastic reduction was in the CHH context, with 34.44%, 40.7% and 35.17% less overall methylation for ago4, ago9, and rdr6, respectively ( Supplementary Table S2 ).

Methylomes were compared to identify DMRs corresponding to genomic segments with a minimum length of 100 nt, between mutant and wild-type. The total number of DMRs including all those present in three different contexts was 8,417 for ago4, 5,127 for ago9, and 5,832 for rdr6 ( Table 1 ). In all three mutants, at least 90% of DMRs in the CG and CHH context had a length comprised between 100-200 nt, and between 100-250 nt for the CHG context ( Supplementary Table S3 ), representing a size distribution similar to those of other mutants involved in the RdDM pathway (Yang et al., 2016). As expected for all three mutants, the DMRs in the CHH context are distributed across all 5 chromosomes, reflecting the overall genomic activity of the RdDM pathway ( Supplementary Figure S5 and Table S4 ). The proportion of DMRs corresponding to a CHH context was of 19% (ago4), 27.3% (ago9), and 26.5% (rdr6) for a total coverage of 227.05 kb, 188.9 kb, and 204.7 Kb, respectively ( Table 1 ). This proportion is significantly smaller than DMRs found in the CHH context for other RdDM mutants such as drm2, rdr2, or mutants affecting specific subunits of POLIV and POLV, in either seedlings or the male meiocyte (Walker et al., 2018; He et al., 2021), suggesting that functional redundancy is prevalent in the pre-meiotic gynoecium to maintain DNA methylation and genome stability.

When compared to wild-type, DMRs in ago4, ago9, and rdr6 were most frequently hypomethylated in CG and CHG contexts ( Figure 2A ). The same is true for ago4 in the CHH context, but not for DMRs in ago9 and rdr6 for the CHH context, which were preferentially hypermethylated ( Figure 2A ). Mapping of all these DMRs revealed that most of them correspond to discrete genomic regions located in TEs. The proportion of DMRs mapping to TEs corresponds to 54.4%, 65.5%, and 63.1% for ago4, ago9, and rdr6 respectively ( Table 2 and Supplementary Tables S5 , S6 ). As expected, the majority of DMRs map to TEs of the most abundant superfamilies, including DNA/Harbinger, LTR/Gypsy, RC/Helitron, DNA/MuDR and LTR/Gypsy, representing more than 70% of all TEs targeted by DMRs ( Figure 2B ; Supplementary Table S5 ). SINE retrotransposons represented 7.6%, 3.1%, and 1.5% of all hypomethylated TEs in ago4, ago9, and rdr6, respectively: most abundant hypomethylated superfamilies also included RathE1_cons and RathE2_cons in ago4 and ago9; and DNA/Mariner and LRT/copy in rdr6 ( Supplementary Figure S6 and Supplementary Table S5 ). In the case of hypermethylated DMRs, DNA/Harbinger and LTR/Gypsy TEs were among the most represented superfamilies common to all three mutants, whereas DNA/MuDR was found only in ago4 and ago9, and SINE was among the most represented in ago4 and rdr6 ( Supplementary Figure S6 and Supplementary Table S5 ).

Given that both ago9 and rdr6 had a higher number of hypermethylated DMRs in the CHH context and the majority mapped to TEs, we determined if families that contained more than 10% of their members targeted by hypermethylated DMRs were equivalent in ago9 and rdr6. Among 54 families overall targeted by DMRs in both genes (25 in ago9 and 29 in rdr6), only six were represented in hypermethylated DMRs in both mutants ( Supplementary Figure S7 ), suggesting that the two corresponding genes target a distinct universe of TEs during female reproductive development. In addition, ago4 hypermethylated TE in ago9 and rdr6, with only one of 67 families being common to all three mutants ( Supplementary Figure S7 ).

Of all the TEs targeted by DMRs, 59.58% (485/814), 58.62% (483/824), and 61.30% (537/876) were unique to ago4, ago9, and rdr6, respectively ( Figure 3A and Supplementary Table S6 ). If only hypomethylated TEs are considered, 78.62% (526/669), 46.43% (78/168) and 41.01% (73/178) were unique to ago4, ago9, and rdr6, respectively, confirming the wider role that ago4 appears to play in methylation of TEs during gynoecium development ( Figure 3B ). If only ago9 and rdr6 were compared, only 20.5% (139/678) of DMRs map to the same TEs, confirming that targeted families of TEs are mostly distinct among the corresponding genes. This result is consistent with the fact that the families of hypermethylated TEs in both mutants were different ( Supplementary Figure S7 ). Only 7.19% (10/139) of hypermethylated TEs are shared between ago4, ago9, and rdr6 ( Figure 3C ), confirming the limited redundancy about their function during gynoecium development.

In all the three mutants a small fraction of DMRs map to gene elements, including their promoter region: 12.7% (204) in ago4, 12.4% (174) in ago9, and 10.3% (159) in rdr6. The total number of genes targeted by these DMRs are illustrated in Figure 3D and Supplementary Table S7 . Only 22 genes had methylation patterns altered in all three mutants, opening possibilities for their functional investigation ( Figures 3D-F ). The identity of these 22 genes is presented in Supplementary Table S8 , using annotation and description from TAIR10 (Berardini et al., 2015), complemented with UNIPROT information (UniProt Consortium, 2021) when necessary. Eight out of 22 genes were hypomethylated in all three mutant backgrounds: ago4, ago9 and rdr6 (At1g11785; At2g28400; At3g17750; At3g46470; At4g03380; At4g18010; At4g28397; At5g46295; Supplementary Table S8 ). These genes encode for two transmembrane proteins, the tyrosine kinase DYRKP-1, an RRM/RBD/RNP motif RNA-binding protein, a non-specific lipid-transfer-like protein, a protein with posphatase activity and two proteins of unknown function. A single gene encoding for a kinase-like protein (At3g23650) was hypermethylated in ago4, ago9 and rdr6 backgrounds. Of the remaining 13 genes, two were hypomethylated in ago4 and ago9, but hypermethylated in rdr6 (At1g52990 and At2g44175, encoding for a thioredoxin and acyl-CoA N-acyltransferase, respectively; Supplementary Table S8 ). Two are hypomethylated in ago4 and rdr6, but hypermethylated in ago9 (At1g18130 and At4g22810: encoding for a biotin synthetase and a protein AT-Hook nuclear localized protein, respectively; Supplementary Table S8 ). A single gene is hypermethylated in ago4 and ago9, but hypomethylated in rdr6 (At1g34410, encoding for the auxin response factor ARF21), or hypermethylated in ago4 and rdr6, but hypomethylated in ago9 (the hypothetical protein At4g03830; Supplementary Table S8 ). Four are hypomethylated in ago4 but hypermethylated in ago9 and rdr6 (At2g18480; At3g02610; At4g15440; and At5g17600; encoding for the probable polyol transporter PLT3, the acyl-acyl carrier protein desaturase AAD2, the cytochrome CYP74B2, and the Arabidopsis toxic yeast protein ATL52, respectively; Supplementary Table S8 ). Finally, there are three genes that contain antagonistic but independent DMRs within their promoter or coding sequence: the hypothetical protein At2g39160 and the UMAMIT45 transporter (At3g28100) are hypomethylated in ago4 and rdr6 and contain two antagonistic RDMs in ago9; whereas the putative transmembrane protein DUF239 is hypomethylated in ago4 but contains two antagonistic RDMs in ago9 and rdr6 ( Supplementary Table S8 ). Our results suggest that these genes are regulated by the RdDM pathway and could play unforeseen roles during female reproductive development.

To determine if these 22 genes were also included as RdDM targets in other organs or cells, we analyzed a collection of 9,993 DMR loci identified in leaves, plantlets and roots of drm2 and rdr2 mutants (Walker et al., 2018); and to a collection of 9,051 DMR loci identified in male meiocytes and pollen cells of drm1 drm2 double mutants (Walker et al., 2018). As described in Supplementary Table S9 , 15 out of these 22 genes (68.2%) are also included in RdDM-dependent DMRs of vegetative tissues. In addition, 18 out of 22 (81.8%) are included in RdDM-dependent DMRs of reproductive cells of the male lineage, suggesting that most of these genes are also regulated by RdDM during male gametogenesis ( Supplementary Table S9 ). Only three out of 22 were not identified in RdDM -dependent DMRs of vegetative organs or male reproductive cells (At5g46295, At4g03830, At3g23650), suggesting that their regulation could be specific to female reproductive development ( Supplementary Table S9 ). These results suggest that the RdDM-dependent regulation of these 22 genes is not restricted to reproductive tissues or specifically dependent on the activity of AGO4, AGO9, and RDR6, but rather shared by multiple genes involved in the RdDM pathway throughout Arabidopsis development.