Whole-genome bisulfite sequencing (WGBS) was performed to examine DNA methylation in the 129/DBA hybrid Zfp57 ^+/− (M⁺Z⁺), Zfp57 ^−/+ (M⁻Z⁺), and Zfp57 ^−/− (M⁻Z⁻) embryos derived from timed mating with Zfp57 ^+/− or Zfp57 ^−/− 129 female mice being crossed with Zfp57 ^+/− (DBA*) male mice mainly on the DBA/2J genetic background as reported in the previous study (Jiang et al., 2021). We found that the two M⁻Z^- mutant embryos had much higher levels of DNA methylation at the Rasgrf1 imprinted region than two other M⁻Z^- mutant embryos (Jiang et al., 2021). We tested these embryo samples to find out what may cause the differences of DNA methylation at this imprinted region when ZFP57 was absent.

We wonder if gender might contribute to the effects of ZFP57 on DNA methylation at the Rasgrf1 ICR as well as other ICRs. This may have resulted in the substantial differences of DNA methylation observed at the Rasgrf1 ICR among four M⁻Z⁻ embryos in the previous WGBS study. There was no significant difference in DNA methylation at all 24 known ICRs comparing female M⁺Z⁺ embryos with the male M⁺Z⁺ embryos, except that DNA methylation was slightly significantly increased at the Nap1I5 ICR but decreased at the Grb10 ICR in the male M⁺Z⁺ embryos compared with the female M⁺Z⁺ embryos (Supplementary Figure S1A). DNA methylation was similar at all ICRs except for the slight increase of DNA methylation at the Peg3 ICR in the comparison of male M⁻Z⁺ embryos with the female ones (Supplementary Figure S1B). DNA methylation was significantly reduced at the AK008011 ICR in the comparison of male M⁻Z⁻ embryos with the female ones (Supplementary Figure S1C). It was dramatically reduced at the Rasgrf1 ICR in the male M⁻Z^- embryos compared with the female ones, although it was not statistically significant, which will be discussed later (Supplementary Figure S1C). Therefore, we examined DNA methylation at the ICRs in individual embryos subjected to WGBS in our previous study (Jiang et al., 2021). As expected, DNA methylation was similarly lost at most ICRs in the female or male M⁻Z⁻ mutant embryos (Figures 1A, 1A’). It was largely retained at five imprinted regions (Peg10, Kcnq1ot1, Gpr1, Slc38a4, and H19) in both female or male M⁻Z^- mutant embryos (Figures 1B, 1B’). Intriguingly, DNA methylation appeared to be more severely lost at the Peg13, Rasgrf1, and AK008011 ICRs in the male M⁻Z^- mutant embryos than in the female M⁻Z^- mutant embryos (Figures 1C, 1C’). This finding can be easily visualized on the methylation IGV plots for three ICRs in these M⁺Z⁺, M⁻Z⁺, and M⁻Z^- embryos (Figure 2). One of two male M⁻Z⁻ embryos displayed more severe loss of methylation at the Peg13 ICR, whereas DNA methylation was more severely lost at the Rasgrf1 and AK008011 ICRs in both male M⁻Z⁻ embryos (Figures 1, 2). We also found that one female M⁻Z⁻ mutant embryo had reduced DNA methylation at the Rasgrf1 ICR, whereas DNA methylation did not appear to be lost at the Rasgrf1 ICR in the other female M⁻Z⁻ mutant embryos (Figures 1C, 2B).

DNA methylation was significantly reduced at the Peg13 ICR in the male M⁻Z⁻ embryos compared with male M⁺Z⁺ and M^-Z⁺ embryos, whereas it was close to being significantly reduced at the Peg13 ICR in the female M⁻Z⁻ embryos compared with female M⁺Z⁺ embryos (Supplementary Figures S1D, S1D’). Indeed, DNA methylation was significantly lost at the Rasgrf1 ICR in the male M⁻Z⁻ embryos compared with male M⁺Z⁺ and M⁻Z⁺ embryos, whereas it was not much different at the Rasgrf1 ICR comparing female M⁻Z⁻ embryos with female M⁺Z⁺ or M⁻Z⁺ embryos (Supplementary Figures S1D, S1D’). DNA methylation at the AK008011 ICR was significantly reduced in both female M⁻Z⁺ and M⁻Z⁻ embryos compared with the female M⁺Z⁺ embryos (Supplementary Figure S1D). It was close to being significantly reduced in both male M⁻Z⁺ and M⁻Z⁻ embryos compared with the male M⁺Z⁺ embryos (Supplementary Figure S1D’). It was even more significantly lost at the AK008011 ICR in the male M⁻Z⁻ embryos than in the female M⁻Z⁻ embryos (Supplementary Figure S1C). We performed similar statistical comparisons for DNA methylation at the Snrpn, Impact, and Cdh15 ICRs that will be discussed for the maternal effect of Zfp57 on their DNA methylation below (Supplementary Figure S1E, S1E’).

To confirm if DNA methylation was, indeed, more severely lost at the Rasgrf1 ICR in the male M⁻Z⁻ mutant embryos than in the female M⁻Z⁻ mutant embryos, we carried out bacterial colony bisulfite sequencing analysis of this Rasgrf1 ICR with another set of M⁻Z⁺ and M⁻Z⁻ embryos that consisted of 3–4 female and male embryos for each genotype (Figure 3). Since DNA methylation levels were similar at the Rasgrf1 ICR in the M⁺Z⁺ and M⁻Z⁺ embryos of the same gender based on WGBS (Figures 1C, 1C’), we think it should be sufficient to just perform bacterial colony bisulfite sequencing analysis of Rasgrf1 ICR to compare M⁻Z⁺ embryos with M⁻Z⁻ embryos in this study. Much higher levels of DNA methylation were obtained from female and male M⁻Z⁺ embryos than with their M⁻Z⁻ counterparts (Figures 3A, 3A’, 3B, 3B’). DNA methylation at the Rasgrf1 ICR was lower in four female M⁻Z⁻ embryos, but it was lowest in three male M⁻Z⁻ embryos (Figures 3B, 3B’). Indeed, DNA methylation was significantly reduced at the Rasgrf1 ICR in both female and male M⁻Z⁻ embryos compared with their M⁻Z⁺ counterparts (Figures 3C, 3C’). It was similar in the female and male M⁻Z⁺ embryos (Figure 3D). However, DNA methylation was significantly reduced at the Rasgrf1 ICR in male M⁻Z⁻ embryos compared with that of female M⁻Z⁻ embryos (Figure 3E). Taken together, DNA methylation was more severely lost at the Rasgrf1 ICR in male M⁻Z⁻ embryos than in female M⁻Z^- embryos (Figure 3, Supplementary Figure S1).

Similarly, we performed bacterial colony bisulfite sequencing analysis of the AK008011 ICR in a set of M⁺Z⁺, M⁻Z⁺, and M⁻Z⁻ embryos that consisted of 3–4 female and male embryos for each genotype (Figure 4). DNA methylation was significantly reduced at the AK008011 ICR in the male M⁻Z⁻ embryos compared with male M⁺Z⁺ or M⁻Z⁺ embryos, whereas it was close to being significantly reduced in the female M⁻Z⁺ or M⁻Z⁻ embryos compared with female M⁺Z⁺ embryos (Figures 4A, 4A’, 4B, 4B’, 4C, 4C’, 4D, 4D’). It was also significantly decreased in male M⁻Z⁺ embryos in comparison to male M⁺Z⁺ embryos (Figure 4D’). There was no significant difference in DNA methylation at the AK008011 ICR comparing female M⁺Z⁺ or M⁻Z⁺ embryos with their male counterparts, respectively (Figures 4A, 4A’, 4B, 4B’, 4E). Interestingly, DNA methylation was more significantly reduced at the AK008011 ICR in male M⁻Z⁻ embryos than in female M⁻Z⁻ embryos (Figure 4E). These results suggest that DNA methylation at the AK008011 ICR was more significantly lost in the male M⁻Z⁻ embryos than in female M⁻Z⁻ embryos (Supplementary Figure S2). Therefore, Zfp57 exhibited sexually dimorphic effect on DNA methylation at the AK008011 ICR in M⁻Z⁻ embryos lacking both maternal and zygotic Zfp57.

We examined DNA methylation at few other known ICRs. As expected, DNA methylation was similarly lost at the ICRs of Inpp5f, Zac1, and IG-DMR in the female and male M⁻Z^- mutant embryos based on the methylation IGV plots (Supplementary Figure S2). Therefore, it seemed that loss of ZFP57 caused more severe loss of DNA methylation imprint at three ICRs in the male M⁻Z⁻ embryos than in the female M⁻Z^- embryos, with statistical significance observed at the Rasgrf1 and AK008011 ICRs.

DNA methylation appeared to be partially lost at a few ICRs in female M⁻Z⁺ or male M⁻Z⁺ embryos (Supplementary Figures 5A, 5A’). This caused partial allelic expression switch of many imprinted genes at the Snrpn imprinted region as well as Impact in male M⁻Z⁺ embryos, which is described in more detail below (Figures 5B, 5B’). Indeed, DNA methylation was partially lost at the Snrpn and Cdh15 ICRs in female M⁻Z⁺ embryos, although it was almost completely lost at these two ICRs in female M⁻Z⁻ embryos (Figure 5A). One female M⁻Z⁺ embryo had more severe loss of DNA methylation at the Snrpn ICR than the other female M⁻Z⁺ embryo. DNA methylation was similarly partially lost at the Impact and AK008011 ICRs in female M⁻Z⁺ embryos as well as in female M⁻Z⁻ embryos (Figure 5A). Similar results were obtained in bacterial colony bisulfite sequencing analysis for the AK008011 ICR (Figures 4A–4C). Indeed, these observations are confirmed by statistical analyses (Figure 4D, Supplementary Figures S1D–S1E). DNA methylation was significantly reduced at these four ICRs in the female M⁻Z⁺ embryos compared with those of the female M⁺Z⁺ embryos (Supplementary Figures S1D–S1E). Furthermore, it was more significantly or close to significantly reduced at the Snrpn and Cdh15 ICRs, but not at the Impact and AK008011 ICRs, while comparing the female M⁻Z⁻ embryos with the female M⁻Z⁺ embryos (Supplementary Figures S1D–S1E). These results suggest that both maternal and zygotic Zfp57 are necessary but partially redundant for the maintenance of DNA methylation imprint at the Snrpn and Cdh15 ICRs in female embryos. However, only maternal Zfp57, but not zygotic Zfp57, is required for maintaining DNA methylation imprint at the Impact and AK008011 ICRs in female embryos.

DNA methylation was also partially lost at the Snrpn ICR in both male M⁻Z⁺ embryos (Figure 5A’). It was partially lost at the Cdh15 ICR in one but not in the other male M⁻Z⁺ embryo (Figure 5A’). DNA methylation was partially lost at the Impact ICR in male M⁻Z⁻ embryos, but it was largely intact at the Impact ICR in male M⁻Z⁺ embryos (Figure 5A’). DNA methylation was partially lost at the AK008011 ICR in one male M⁻Z⁺ embryo, which was similar to loss of DNA methylation at the AK008011 ICR in both male M⁻Z⁻ embryos (Figure 5A’). Unfortunately, there were no sequence reads at the AK008011 ICR in the other male M⁻Z⁺ embryo and, therefore, we could not determine if DNA methylation was similarly lost at the AK008011 ICR in that male M⁻Z⁺ embryo based on WGBS. These results were also confirmed by statistical analyses (Supplementary Figures S1D’–S1E’). DNA methylation was significantly reduced at the AK008011 ICR in male M⁻Z⁺ embryos compared with male M⁺Z⁺ embryos in bacterial colony bisulfite sequencing analysis, and it was even more significantly reduced in male M⁻Z⁻ embryos than in male M⁻Z⁺ embryos (Figures 4A’–D’). DNA methylation was close to significantly reduced at the Snrpn ICR in the male M⁻Z⁺ embryos compared with the male M⁺Z⁺ embryos, but even more significantly reduced in the male M⁻Z⁻ embryos than in the male M⁺Z⁺ or M⁻Z⁺ embryos (Supplementary Figure S1E’). DNA methylation was only significantly reduced at the Cdh15 ICR in the male M⁻Z^- embryos versus male M⁺Z⁺ embryos, whereas it was significantly reduced at the Impact ICR in the male M⁻Z⁻ embryos compared with that in the male M⁺Z⁺ or M⁻Z⁺ embryos (Supplementary Figure S1E’).

Actually, DNA methylation was more severely lost at the AK008011 ICR in male M⁻Z⁻ embryos compared with that in female M⁻Z⁻ embryos (Figures 4E, 5A, 5A’, Supplementary Figure S1C). It seems that maternal Zfp57 is required for maintaining DNA methylation imprint at the AK008011 ICR in both female and male embryos, whereas zygotic Zfp57 seemed to contribute to maintenance of DNA methylation at the AK008011 ICR in male but not female embryos. Both maternal and zygotic Zfp57 play partially redundant roles in maintaining DNA methylation at the Impact ICR in male embryos, although only maternal but not zygotic Zfp57 is involved in maintaining DNA methylation at the Impact ICR in female embryos (Figure 5A, 5A’, Supplementary Figure S1C). These are also caused by the gender-specific effect of Zfp57 on ICR DNA methylation.

Loss of DNA methylation at these ICRs in female M⁻Z⁺ and M⁻Z⁻ embryos or male M⁻Z⁺ and M⁻Z⁻ embryos is also clearly seen on the methylation IGV plots (Figures 2, 6). DNA methylation was similarly partially lost at the AK008011 ICR in female and male M⁻Z⁺ embryos as well as in female M⁻Z⁻ embryos, with more severe loss observed in the male M⁻Z⁻ embryos than in the female ones (Figures 2C, 2C’). Partial loss of DNA methylation was observed at the Snrpn ICR in female M⁻Z⁺ and male M⁻Z⁺ embryos on the methylation IGV plots, although DNA methylation was almost completely missing at the Snrpn ICR in female M⁻Z⁻ and male M⁻Z⁻ embryos (Figures 6A, 6A’). Partial loss of DNA methylation was observed at the Impact ICR in male M⁻Z⁻, but not male M⁻Z⁺ embryos, whereas DNA methylation was partially lost at the Impact ICR in female M⁻Z⁺ as well as in female M⁻Z⁻ embryos on the methylation IGV plots (Figure 6B, 6B’). Germline-derived DNA methylation imprint at the Snrpn and Impact ICRs was present on the maternal chromosomes in M⁺Z⁺ embryos as shown in the allelic methylation IGV plots (Supplementary Figure S3, Supplementary Table S1, S2). Indeed, loss of DNA methylation imprint occurred at the maternal Snrpn ICR in M⁻Z⁺ and M⁻Z⁻ embryos based on the allelic methylation IGV plots (Supplementary Figure S3A-S3A′). DNA methylation was also lost at the maternal Impact ICR in the female M⁻Z⁺ and M⁻Z⁻ embryos and in male M⁻Z⁻ embryos based on the allelic methylation IGV plots (Supplementary Figure S3B-S3B′). The allelic methylation results further confirm what has been observed on the methylation IGV plots of the Snrpn and Impact ICRs (Figure 6, Supplementary Figure S3A). Partial loss of DNA methylation was observed at the Cdh15 ICR in both female M⁻Z⁺ embryos and one male M⁻Z⁺ embryo (Figures 6C, 6C’). DNA methylation was almost completely lost at the Cdh15 ICR in female M⁻Z⁻ and male M⁻Z⁻ embryos compared with M⁺Z⁺ embryos on the methylation IGV plots (Figures 6C, 6C’).

Since DNA methylation imprint was partially lost at a few ICRs in M⁻Z⁺ embryos, we wonder if this may have any effect on the expression of the corresponding imprinted genes at these ICRs. RNA-seq analyses were carried out with another set of mouse embryos to examine expression of imprinted genes in the female M⁺Z⁺, M⁻Z⁺, and M⁻Z⁻ 129/DBA hybrid E13.5 embryos and male M⁺Z⁺, M⁻Z⁺, and M⁻Z^- 129/DBA hybrid E13.5 embryos, with at least three embryos in each group. These mouse embryos were different from the ones used for WGBS in the previous study. We used P-score for measuring the proportion of the transcripts expressed from the paternal allele of each imprinted gene (Figure 7, Supplementary Table S3–S4, Materials and Methods). Since there are no exonic SNPs on many imprinted genes, we could not analyze their P-score differences in these embryos. Ftx and Jpx were only two out of 48 known imprinted genes with an exonic SNP that showed significant P-score difference in statistical analysis (∆P-Score≥ 0.1, p < 0.05) comparing six female M⁺Z⁺ embryos with four male ones (Figure 7A). As expected, X-linked Ftx and Jpx, which are involved in X chromosome inactivation in the female embryos, along with Xist and Tsix, exhibited almost biallelic expression in the female M⁺Z⁺ embryos, whereas they were exclusively expressed from the maternal allele on the X chromosome in the male M⁺Z⁺ embryos (Figure 7A) (Collombet et al., 2020; Patrat et al., 2020). Thus, Ftx and Jpx serve as good internal positive controls for our data analyses. They were also the only two imprinted genes that displayed significant differences in the expression levels comparing six female M⁺Z⁺ embryos with four male ones that will be discussed further below (Figure 7B).

The intensity difference heatmaps were generated for 21 imprinted genes with statistically different P-score (∆P-Score≥ 0.1, p < 0.05) by comparing the female M⁻Z⁻ embryos with the female M⁺Z⁺ embryos and for 20 imprinted genes (∆P-Score≥ 0.1, p < 0.05) when the male M⁻Z⁻ embryos were compared with the male M⁺Z⁺ embryos (Table 1, Figure 7C, 7C’). Zim1 at the Peg3 imprinted region was maternally expressed in female and male M⁺Z⁺ embryos or male M⁻Z⁻ embryos, but it became partially biallelic in female M⁻Z⁻ embryos (Supplementary Figure S4A).

Similarly, the intensity difference heatmaps (∆P-Score≥ 0.1, p < 0.05) were generated for seven imprinted genes by comparing three female M⁻Z⁺ embryos with six female M⁺Z⁺ embryos and for 11 imprinted genes when three male M⁻Z⁺ embryos were compared with four male M⁺Z⁺ embryos (Figures 7D, 7D’). Indeed, allelic expression of the imprinted genes at the Snrpn imprinted region was mostly affected in female or male M⁻Z⁺ embryos as well as in female or male M⁻Z⁻ embryos compared with female or male M⁺Z⁺ embryos of the same gender (Figures 5B, 5B’, 7, Supplementary Table S3–S4). Snurf, Snrpn, Ndn, Magel2, Peg12, Mkrn3, and Snord64 at the Snrpn imprinted region were all paternally expressed in female and male M⁺Z⁺ embryos. They became biallelic in female or male M⁻Z⁻ embryos, except that Snord64 was slightly preferentially maternally expressed in male M⁻Z⁻ embryos (Figures 5B, 5B’, 7, Supplementary Table S3–S4). They were largely biallelic in female or male M⁻Z⁺ embryos at variable levels. Consistent with their tissue-specific or species-specific imprinting phenomena, as well as our previous results in mouse embryos, Atp10a and Ube3a at the Snrpn imprinted region were almost biallelic in the female M⁺Z⁺, M⁻Z⁺, and M⁻Z⁻ embryos or male M⁺Z⁺, M⁻Z⁺, and M⁻Z⁻ embryos (Figures 5B, 5B’, Supplementary Table S3-S4) (DuBose et al., 2010; Hsiao et al., 2019; Jiang et al., 2021).

Impact was preferentially paternally expressed in female and male M⁺Z⁺ embryos (Figure 5C-5C’, Figure 7C-7C’, Supplementary Table S3–S4). It became almost biallelic in female or male M⁻Z⁻ embryos and partially biallelic in male M⁻Z⁺ embryos, whereas it remained preferentially paternally expressed in female M⁻Z⁺ embryos (Figure 5C-5C’, Figure 7C-7C’, 7D’). Zdbf2 and Nnat (Peg5) were preferentially paternally expressed in female or male M⁺Z⁺ embryos (Supplementary Figure S4, Supplementary Table S3-S4). Zdbf2 was partially biallelic in female M⁻Z⁺ embryos and almost completely biallelic in male M⁻Z⁺ embryos (Supplementary Figure S4B-S4B'). Nnat (Peg5) remained preferentially paternally expressed in female M⁻Z⁺ embryos, but it was almost completely biallelic in male M⁻Z⁺ embryos (Supplementary Figure S4C-S4C'). Slightly preferential maternal expression was observed for H13 in female and male M⁺Z⁺ embryos or female M⁻Z⁺ embryos (Supplementary Figure S4D-S4D'). H13 was biallelic in male M⁻Z⁺ embryos (Supplementary Figure S4D'). Taken together, loss of maternal Zfp57 caused variable allelic expression switch of Impact, Zdbf2, Nnat (Peg5), and H13 in the male M⁻Z⁺ embryos, whereas it only had a minor effect on allelic expression of Zdbf2 in female M⁻Z⁺ embryos. The gender effects on allelic expression of these imprinted genes will be further discussed below.

We also analyzed expression levels of some imprinted genes based on RNA-seq results. Among 105 analyzed imprinted genes with the TPM values above the threshold (see Materials and Methods), only two imprinted genes (Jpx and Ftx) had statistically significant differences (|log₂FC |> 0.3, p < 0.05) in gene expression levels comparing the female M⁺Z⁺ embryos with the male ones (Figure 7B). Since they are involved in X chromosome inactivation only in the females (Furlan et al., 2018; Collombet et al., 2020; Patrat et al., 2020), it is expected that the expression of Jpx and Ftx was higher in the female M⁺Z⁺ embryos than in the male ones (Figure 7B).

The intensity difference heatmaps were generated for 39 imprinted genes with statistically significant changes (|log₂FC |> 0.5, p < 0.05) in their expression levels when the female M⁻Z⁻ embryos were compared with the female M⁺Z⁺ embryos and for 44 imprinted genes when the male M⁻Z⁻ embryos were compared with the male M⁺Z⁺ embryos (Table 2; Figure 8A, 8A’). Similarly, the intensity difference heatmaps were generated for eight imprinted genes in the comparison of the female M⁻Z⁺ embryos with the female M⁺Z⁺ embryos and for nine imprinted genes when the male M⁻Z⁺ embryos were compared with the male M⁺Z⁺ embryos (Table 3, Figure 8B, 8B’). Eight imprinted genes at the Snrpn imprinted region were differentially expressed either in the comparison of the female M⁻Z⁺ embryos with the female M⁺Z⁺ embryos or in the comparison of the male M⁻Z⁺ embryos with the male M⁺Z⁺ embryos (Figure 8B, 8B’).

Peg3os was significantly increased in male M⁻Z⁺ embryos compared with male M⁺Z⁺ embryos (Figure 8B’). Impact was significantly increased in female M⁻Z⁺ embryos compared with female M⁺Z⁺ embryos (Supplementary Figure S5B-S6B', Supplementary Table S5-S6). There was not much change in the expression levels of Cdh15 in female M⁻Z⁺ and M⁻Z^- embryos or male M⁻Z⁺ and M⁻Z⁻ embryos compared with M⁺Z⁺ embryos (Supplementary Figure S5C-S6C', Supplementary Table S5-S6). Zrsr1 was also significantly increased in male M⁻Z⁺ embryos compared with male M⁺Z⁺ embryos (Supplementary Figure S8A-S8A', see below).

Taken together, loss of maternal Zfp57 caused increased expression of eight imprinted genes at the Snrpn imprinted region in both female and male M⁻Z⁺ embryos. It also resulted in increased expression of Peg3os and Zrsr1 in male M⁻Z⁺ embryos, as well as increased expression of Impact in female M⁻Z⁺ embryos.

Allelic expression switch occurs to some imprinted genes in mouse embryos when Zfp57 is lost (Jiang et al., 2021). Since DNA methylation imprint might be different at a few ICRs in the male M⁻Z⁻ embryos compared with the female ones (Figure 9A), we wonder if allelic expression of the imprinted genes could also be differentially affected in the male M⁻Z⁻ embryos compared with the female M⁻Z⁻ embryos. Unfortunately, no exonic SNP was present in Rasgrf1 in the hybrid 129/DBA embryos, and its expression level was lower than the threshold of our RNA-seq TPM expression analysis (see Materials and Methods). There is no known imprinted gene located at the AK008011 imprinted region. Therefore, we could not determine in this study if there is any gender-specific effect on allelic expression switch and expression level differences for the corresponding imprinted genes at the Rasgrf1 and AK008011 imprinted regions in female or male M⁻Z⁻ mutant embryos compared with the control M⁺Z⁺ embryos. This will be tested in a future study.

Peg13 was paternally expressed in female M⁺Z⁺ and M⁻Z⁺ embryos or male M⁺Z⁺ and M⁻Z⁺ embryos (Supplementary Figure S6A-S6A', Supplementary Table S3-S4). It became partially biallelic in female M⁻Z⁻ and male M⁻Z⁻ embryos (Supplementary Figure S6A-S6A', Supplementary Table S3-S4). By contrast, Kcnk9 at the Peg13 imprinted region was maternally expressed and Trappc9 at the Peg13 imprinted region was biallelic in female M⁺Z⁺, M⁻Z⁺, and M⁻Z⁻ embryos or male M⁺Z⁺, M⁻Z⁺, and M⁻Z⁻ embryos (Supplementary Figure S6A-S6A', Supplementary Table S3-S4). Thus, partial loss of DNA methylation at Peg13 ICR caused Peg13 to be partially biallelic in female M⁻Z^- or male M⁻Z⁻ embryos. However, it had no apparent effect on allelic expression of Kcnk9 and Trappc9 at the Peg13 imprinted region in female M⁻Z⁻ or male M⁻Z⁻ embryos. The residual DNA methylation at the Peg13 ICR may be sufficient to maintain allelic expression of Kcnk9 at the Peg13 imprinted region in female M⁻Z⁻ or male M⁻Z⁻ embryos. Despite this, expression of Kcnk9 was decreased in M⁻Z⁻ embryos compared with M⁺Z⁺ and M⁻Z⁺ embryos, and it was more reduced in male M⁻Z⁻ embryos than in female M⁻Z⁻ embryos, although it was not statistically significant (Supplementary Figure S6B-S6B', Supplementary Figure S6C, Table S5-S6). This is consistent with DNA methylation at the Peg13 ICR in these embryos (Figure 1, Supplementary Figure S1).

Usp29 at the Peg3 imprinted region was paternally expressed in female M⁺Z⁺ and M⁻Z⁺ or male M⁺Z⁺ and M⁻Z⁺ embryos, but it was biallelic in female M⁻Z⁻ or male M⁻Z⁻ embryos (Supplementary Figure S4A, Supplementary Figure S4A', Supplementary Table S3-S4). Zim1 at the Peg3 imprinted region was maternally expressed in female M⁺Z⁺ and M⁻Z⁺ or male M⁺Z⁺ and M⁻Z⁺ embryos (Supplementary Figure S4A, Supplementary Figure S4A', Supplementary Table S3-S4). Although it was still maternally expressed in male M⁻Z⁻ embryos, Zim1 was partially biallelic in female M⁻Z⁻ embryos (Table 1, Supplementary Figure S4A, Supplementary Figure S4A', Supplementary Table S3-S4). Thus, Zim1 appeared to be affected differently in female M⁻Z⁻ or male M⁻Z⁻ embryos lacking ZFP57. Consistent with this, Zim1 is the only imprinted gene that displayed gender-specific effect of ZFP57 with significant P-score difference (∆P-Score≥ 0.1) in the comparisons of the female and male M⁻Z⁻ embryos with their counterpart M⁺Z⁺ embryos (Figures 7C–7C’).

Impact was partially biallelic in male M⁻Z⁺ embryos, although it remained preferentially paternally expressed in female M⁻Z⁺ embryos (Figures 5C–5C’). Zdbf2 was biallelic in male M⁻Z⁺ embryos and partially biallelic in female M⁻Z⁺ embryos (Supplementary Figure S4B-S4B'). Nnat (Peg5) became almost biallelic in male M⁻Z⁺ embryos, whereas it was preferentially paternally expressed in female M⁻Z⁺ embryos (Supplementary Figure S4C-S4C'). H13 was slightly preferentially maternally expressed in female M⁻Z⁺ embryos, but it was biallelic in male M⁻Z⁺ embryos (Supplementary Figure S4D-S4D'). These four imprinted genes displayed gender-specific differences in allelic expression switch upon loss of maternal Zfp57.

We also examined allelic expression of the imprinted genes at the Inpp5f, Zac1, and Dlk1–Dio3 imprinted regions (Supplementary Figure S7, Supplementary Table S3-S4). For the tested imprinted genes at these three imprinted regions, they all became biallelic in female or male M⁻Z⁻ embryos. Their allelic expression did not change upon loss of maternal Zfp57 in either female or male M⁻Z⁺ embryos. Therefore, allelic expression was only affected in M⁻Z^- embryos for the imprinted genes at these three examined imprinted regions, and there was no difference in their allelic expression comparing male embryos with female embryos with or without ZFP57.

Loss of maternal Zfp57 resulted in mis-expression of eight imprinted genes at the Snrpn imprinted region in both female M⁻Z⁺ embryos and male M⁻Z⁺ embryos (Figures 8B, 8B’). Two other imprinted genes (Zrsr1 and Peg3os) were deregulated in the male M⁻Z⁺ embryos compared with male M⁺Z⁺ embryos, whereas the Impact expression level was significantly affected in the female M⁻Z⁺ embryos compared with female M⁺Z⁺ embryos (Figures 8B–8B’, 9B, Supplementary Figure S5B-S5B', Supplementary Figure S8A-S8A', Supplementary Figure S10A-10A').

We also observed some gender-specific differences in the expression levels of some imprinted genes in the female or male M⁻Z^- embryos compared with their counterpart female or male M⁺Z⁺ embryos (Figure 9C, Supplementary Figure S8). Expression levels of Slc22a3, Ascl2, and Th were significantly affected in male M⁻Z⁻ embryos compared with male M⁺Z⁺ embryos but not in female M⁻Z⁻ embryos compared with female M⁺Z⁺ embryos (Figures 8A, 8A’, Supplementary Figure S8B-S8B', Supplementary Figure S8C-S8C'). By contrast, Calcr was only significantly reduced in female M⁻Z⁻ embryos compared with female M⁺Z⁺ embryos (Figures 8A, 8A’, Supplementary Figure S8D-S8D').

Many imprinted genes showed similarly significant differences in their expression in both female and male M⁻Z⁻ embryos compared with the M⁺Z⁺ embryos of the same gender (Table 2, Figure 8A–8A’). Indeed, most imprinted genes at the Dlk1–Dio3 imprinted region appeared to be similarly affected upon loss of ZFP57 comparing M⁻Z⁻ embryos with M⁺Z⁺ embryos of the same gender (Table 2, Figure 8A–8A’, Supplementary Figure S9, Supplementary Table S5-S6). Begain, Dio3, Dio3os, Dlk1, and Rtl1 were reduced in the female M⁻Z⁻ embryos compared with female M⁺Z⁺ embryos, whereas Rian, Meg3, Mirg, AF357359, AF357355, AF357425, and Mir410 displayed increased expression in the female M⁻Z⁻ embryos compared with female M⁺Z⁺ embryos (Table 2, Supplementary Figure S9A). These genes behaved similarly in the male M⁻Z^- embryos compared with male M⁺Z⁺ embryos (Table 2, Supplementary Figure S9A'). Expression of Mir431 and Mir335 was increased in male but not in female M⁻Z⁻ embryos compared with M⁺Z⁺ embryos of the same gender (Table 2, Figures 8A, 8A’, Supplementary Figure S9). Therefore, Mir431 and Mir335 were more affected in the male embryos upon loss of Zfp57. Xlr3b was significantly increased in male M⁻Z⁻ embryos compared with male M⁺Z⁺ embryos but not in the similar comparison of the female embryos (see Supplementary Figure S11B-S11B').

We also examined expression of the imprinted genes at the Peg3, Inpp5f, Zac1, and Peg5 imprinted regions in the female or male embryos (Table 2, Figures 8A, 8A’, Supplementary Figure S10, Supplementary Table S5-S6). For the genes at the Peg3 imprinted region, decreased expression of Zim1 and increased expression of Usp29, Peg3, and Peg3os were similarly observed comparing female M⁻Z⁻ embryos with female M⁺Z⁺ embryos or comparing male M⁻Z⁻ embryos with male M⁺Z⁺ embryos (Table 2, Supplementary Figure S10A-S10A'). Expression of Inpp5f was increased in M⁻Z⁻ embryos compared with M⁺Z⁺ embryos, regardless of the gender (Table 2, Supplementary Figure S10B-S10B'). We observed increased expression of Zac1 (Plagl1) and Hymai and decreased expression of Phactr2 at the Zac1 imprinted region in M⁻Z⁻ embryos compared with M⁺Z⁺ embryos, regardless of whether they are females or males (Table 2, Supplementary Figure S10C-S10C'). Increased Nnat (Peg5) expression and decreased Blcap expression were similarly observed at the Peg5 imprinted region in M⁻Z⁻ embryos compared with M⁺Z⁺ embryos of the same gender (Table 2, Supplementary Figures S10D–S10D’). Decreased expression of Ampd3 was also similarly observed in both female and male M⁻Z⁻ embryos compared with female and male M⁺Z⁺ embryos (Supplementary Figure S11A-S11A'). Therefore, no gender-specific effect was observed for the expression levels of these imprinted genes when Zfp57 was absent in the female or male embryos.

The mice carrying the Zfp57 deleted mutant allele on the 129S6/SvEvTac genetic background were generated in the original study, and they were called the 129 mice in this study (Li et al., 2008). We also performed 12 backcrosses for the 129 mice carrying the Zfp57 deleted mutant allele with the wild-type DBA/2J mice to obtain Zfp57 heterozygous mice mainly on the DBA/2J background that was named Zfp57 ^+/− (DBA*) in a recent study (Jiang et al., 2021). Then, timed pregnancy mating was set up between Zfp57 ^−/− homozygous 129 female mice and Zfp57 ^+/− (DBA*) male mice to obtain 129/DBA hybrid E13.5 embryos. The mouse vaginal plug was checked daily. The female mice were presumed to be potentially pregnant for 0.5 days when their vaginal plug was found around the noon of the day. The pregnant female mice were dissected for 129/DBA hybrid E13.5 embryos used for RNA-seq analysis as well as whole-genome bisulfite sequencing (WGBS) analysis.

To determine the gender of these 129/DBA hybrid E13.5 embryos, SRY genotyping was carried out with the primers SRY-F1 and SRY-R1 for a PCR product of 266 bp present in the male embryos but absent in the female embryos. The sequence for SRY-F1 is 5’- CCA​CTC​CTC​TGT​GAC​ACT, whereas the sequence for SRY-R1 is 5’- GAG​AGC​ATG​GAG​GGC​CAT.

Total RNA samples were obtained from the 129/DBA hybrid E13.5 embryos from the timed pregnancy mating. At least three female embryos of each genotype of M⁺Z⁺, M⁻Z⁺, and M⁻Z^- were used for RNA purification and sent out for RNA-seq analysis. Similarly, total RNA samples were purified from at least three male embryos of each genotype of M⁺Z⁺, M⁻Z⁺, and M⁻Z^- and subjected to RNA-seq analysis. The data analyses were performed as described below.

RNA-seq analysis was performed for the RNA samples obtained from the mouse embryos that were different from the ones used in WGBS of the previous study. At least three embryos were analyzed for each genotype of female M⁺Z⁺, M⁻Z⁺, and M⁻Z⁻ 129/DBA hybrid E13.5 embryos and male M⁺Z⁺, M⁻Z⁺, and M⁻Z⁻ 129/DBA hybrid E13.5 embryos.

First, the quality of the sequence reads obtained from RNA-seq was assessed by FastQC (v0.11.9) downloaded from the website (www.bioinformatics.babraham.ac.uk/projects/fastqc/). Then, the adapter sequences were trimmed from the sequence reads by using Trim_Galore (v0.6.6) (https://github.com/FelixKrueger/TrimGalore). Only good quality sequence reads were aligned to the mouse reference genome (mm9) on the UCSC website using STAR (v2.7.8a), with the setting of “–outFilterMultimapNmax 1 –alignEndsType EndToEnd–outSAMattributes NH HI NM MD” plus other default parameters (Dobin et al., 2013). The sequence reads were mapped to the annotated genes (UCSC mm9. gtf) with the featureCounts function (-O–s 2 option) of Subread (v.2.0.1) (Liao et al., 2014). Out of 148 mouse imprinted genes listed on the website of the geneimprint database: www.geneimprint.com, 139 had been mapped to the mouse reference genome (mm9) and then used for expression analyses in this study.

Transcripts per million (TPM) values were calculated for 139 mapped imprinted genes with more than 10 sequence reads. 105 imprinted genes with the TPM value of more than 1.0 in at least one embryo sample were selected for further analysis of their expression differences across different genotypes or genders. The log₂FC values were used to quantify the fold change (FC) of gene expression based on the ratios of TPM for this imprinted gene in three genotypes of M⁺Z⁺, M⁻Z⁺, and M⁻Z⁻. The differentially expressed genes in the RNA-seq data of different genotypes were first identified by DESeq2 before they were used for constructing the intensity difference heatmaps (Love et al., 2014). The cutoff threshold was set for |log₂FC |> 0.3 (p < 0.05) for the intensity difference heatmap when the expression levels of the imprinted genes were compared between six female M⁺Z⁺ embryos and four male M⁺Z⁺ embryos (Figure 7B). It was set for |log₂FC |> 0.5 (p < 0.05) in other intensity difference heatmaps comparing their expression levels in the M⁺Z⁺, M⁻Z⁺, and M⁻Z⁻ embryos of the same gender (Figure 8).

The SNPs present in the 129/DBA hybrid embryos were determined by using SNP calling with the Genome Analysis Toolkit (GATK) software package (McKenna et al., 2010; Van der Auwera et al., 2013). These SNPs were used for creating N-masked genomes from the reference genome (mm9) using BEDTools with incorporation of “N” at the position of an SNP (Quinlan and Hall, 2010). Mapping of the sequence reads to the N-masked genomes was performed using STAR (v2.5.4a), with the setting of “–outFilterMultimapNmax 1 –alignEndsType EndToEnd–outSAMattributes NH HI NM MD” and other default parameters (Dobin et al., 2013). The mapped sequence reads were assigned to the maternal (129) or paternal (DBA) genome using SNPsplit (v0.4.0), and then a file was generated with the SNPs (Krueger and Andrews, 2016). The allele-specific sequence reads were assigned to the maternal and paternal alleles of the imprinted genes using featureCounts (Liao et al., 2014). Further analyses in this study were limited to the genes with more than 10 counts based on the sum of the mapped reads of maternal and paternal alleles.

P-score was calculated for assessing the proportion of the transcripts derived from the paternal allele of each imprinted gene based on the following equation: P-score = P/(M + P)-0.5. P stands for the number of the paternal (P) allele transcript reads of an imprinted gene, whereas M stands for the number of the maternal (M) allele transcript reads of an imprinted gene. This was adapted from the calculation method for allelic expression of X-linked genes in a recent study (Yu et al., 2021). The imprinted genes with statistically significant allelic expression difference (∆P-Score≥ 0.1, p < 0.05) used for intensity difference heatmap analyses were identified by comparing their P-score values among different genotypes by using the Kruskal–Wallis test (Figure 7). Similar results were obtained by using the binomial test used in a recent study (Yu et al., 2021).

The WGBS data were based on the previous study (Jiang et al., 2021). The gender of these female or male embryos used for WGBS was determined by PCR genotyping of the Sry gene as described above. Please refer to the previous article for the WGBS sequence depth and genome coverage (Jiang et al., 2021). The WGBS sequence reads of these embryo samples were mapped to the abovementioned N-masked reference genome. The CpG sites present in at least three unique reads were used for quantification of DNA methylation of an ICR in the subsequent analyses. The methylated and unmethylated C residues were identified at each CpG site of an ICR. Then, the DNA methylation level was quantified for each CpG site of this ICR based on the number of methylated C residues divided by the total number of C residues for this CpG site. The percentage (%) of DNA methylation for the ICR was the average of DNA methylation levels of all CpG sites at this ICR.

SNPsplit (v0.4.0) was used to separate the allelic DNA methylation reads of the ICRs with the default parameters (Krueger and Andrews, 2016). Only CpG sites covered by at least one unique sequence read were subject to measurement of DNA methylation levels of the maternal or paternal ICR. Allelic DNA methylation was determined for only a few ICRs with SNPs.

Genomic DNA samples purified from the mouse embryos were subjected to bisulfite mutagenesis with the EZ DNA Methylation-Gold™ Kit (Zymo Research). Then, bisulfite PCR was carried out for a Rasgrf1 ICR region with two rounds of nested PCR reactions (Zuo et al., 2012; Takikawa et al., 2013a). The primers used for the first round of bisulfite PCR were Ras-Bis-OF with the sequence of 5’- GTT​ATT​ATT​ATG​TGT​TAT​GTG​TAG​TAA​G and Ras-Bis-OR1 with the sequence of 5’- TAA​TAC​AAC​AAC​AAC​AAT​AAC​AAT​C. The primers used for the second round of nested bisulfite PCR were Ras-Bis-IF with the sequence of 5’- GGT​GTA​GAA​TAT​GGG​GTT​G and Ras-Bis-IR with the sequence of 5’- ATA​CAA​CAA​CAA​CAA​TAA​C. Two specific PCR products of 321 bp and 408 bp were obtained after the second round of nested bisulfite PCR, and they had completely overlapping 321-bp long sequences at one end. Both specific PCR products were cloned into the pUCm-T vector (Sangon, cat# B522213), and the resultant bacterial colonies were sequenced to determine the methylation status of the CpG sites within the 321-bp region of the Rasgrf1 ICR. The obtained sequences were analyzed with the web-based program called QUMA (http://quma.cdb.riken.jp/). Then, DNA methylation status was determined for the CpG sites in the Rasgrf1 ICR as shown in Figure 3.

Bacterial colony bisulfite sequencing was similarly carried out for the AK008011 ICR with the genomic DNA samples derived from 3–4 M⁺Z⁺, M⁻Z⁺, or M⁻Z⁻ E13.5 embryos of each gender that had also been used for RNA-seq analyses in this study (Figures 4, 7). Two rounds of nested PCR reactions were performed for the bisulfite genomic DNA product. The primers used for the first round of bisulfite PCR were AK-Bis-OFn1 with the sequence of 5’- GGT​TTA​GTT​AGG​GAA​AGG​GT and AK-Bis-ORn1 with the sequence of 5’- CAC​ACA​CCT​AAA​TCC​TAA​CAC​T, with the PCR product of 765 bp. The primers used for the second round of nested bisulfite PCR were AK-Bis-OFn2 with the sequence of 5’- GTG​GTT​ATA​TAT​TGT​AGG​GTA​GG and AK-Bis-ORn2 with the sequence of 5’- CCT​ACA​TAA​TTA​AAA​CCT​ACC​TC. The obtained PCR product of 605 bp was cloned into the pUCm-T vector (Sangon, cat# B522213), and the resultant bacterial colonies were sequenced to determine the methylation status of the CpG sites of the AK008011 ICR. The obtained sequences were analyzed with the web-based program called QUMA (http://quma.cdb.riken.jp/). Then, DNA methylation status was determined for the CpG sites in the AK008011 ICR as shown in Figure 4.

Except for statistical analyses used to generate the intensity difference heatmaps, ANOVA (Fisher LSD) was used for statistical analysis of imprinted gene expression (FC value and P-score) differences among the RNA-seq samples of different genotypes of the same gender. And, two-tailed Student’s t-test was used for analyzing the statistical significance when expression of the imprinted genes was compared between the female samples and the male samples within the same genotype.

ANOVA (Fisher LSD) was similarly utilized in analyzing the statistical significance of DNA methylation at the ICRs of different genotypes of the same gender based on the WGBS data. And, two-tailed Student’s t-test was used to compare DNA methylation at the ICRs of the same genotype between the female samples and male samples. Calculation of the p-value for all statistical analyses was carried out by using the open source R software (https://www.R-project.org).