We previously reported the details of the procedures for handling experimental animals and sample collection (Yao et al., 2021), which were approved by the Institutional Animal Ethics Committee of the Nanjing Agricultural University (SYXK 2011-0036). In brief, we firstly chose twenty non-pregnant ewes (2, 3 years old) based on their litter size numbers of three records, categorizing them as HP ewes (litter size = 3, n = 4) and LP ewes (litter size = 1, n = 16). Meanwhile, BMPR1B polymorphism genotyping of those ewes was detected as described previously (Yao et al., 2021). Finally, nine ewes were selected and divided into HPBB (n = 3), LPBB (n = 3) and LPB+ (n = 3) groups. All ewes were slaughtered during estrus and the ovarian samples were immediately collected and stored at −80°C for WGBS and RNA-seq.

Genomic DNA and total RNA were extracted from ovarian tissue of the nine sheep using a Genomic DNA kit (Cat.#DP304-03,Tiangen, Beijing, China) and TRIzol reagent (Cat.# 15596-026; Invitrogen, Carlsbad, CA), respectively. The concentration and purity of isolated DNA and RNA were determined using NanDrop spectrophotometer (NanoDrop Technologies, Wilmington, DE, United States) and agarose gel electrophoresis. The preparation of nine genome-wide DNA methylation and transcriptome libraries and sequencing (WGBS and RNA-seq), respectively, were performed by Biomarker Technologies Corporation (Beijing, China). Notably, we previously reported the DNA methylation (Zhang et al., 2017) and mRNAs (Feng et al., 2018) profiles in ovaries of only HPBB and LPBB Hu sheep, and using the same described methods, we analyzed the WGBS and RNA-Seq data in this study, with some modifications. Briefly, the raw data were first filtered to remove low-quality reads, and the clean data were then aligned with the reference genome of sheep (Ovis aries v4.0).

Differential methylated regions (DMRs) were defined by the presence of at least three methylation sites in the region, and the difference in methylation levels was >0.2 (>0.3 for CG type) with P-value (Fisher’s) < 0.05. Differentially expressed (DE) mRNAs were identified with the false discovery rate (FDR) < 0.05 and absolute value of log₂ (fold change, FC) >1. Primarily, DMGs were detected through mapping the DMRs to genes based on their genomic location. In this study, we defined the genomic region from −3000 bp to the transcription start site (TSS) as the promoter region and from the TSS to the transcriptional termination site (TTS) as the gene body region. To visualize the overlapping gene sets, we generated venn diagrams using Venn diagram online tool (http://bioinformatics.psb.ugent.be/webtools/Venn/).

The correlation between gene expression and methylation was calculated using only differentially expressed genes (DEGs) and DMR-related genes (DMGs). For the global methylation around TSS regions (±2000 bp) and gene expression correlation analysis, the genes was divided into four groups according to their expression level: hightest, lowest, medium-high and medium-low. T visualize the relationship between methylation and gene expression, heatmaps were generated using the R package (heatmap).

Gene Ontology (GO) enrichment analyses of DMGs and DEGs were implemented by the Wallenius non-central hyper-geometric distribution in the GOseq R packages (Young et al., 2010). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of DMGs and DEGs were evaluated using KOBAS software (Mao et al., 2005).

The interaction network of DMGs and DEGs associated with female reproduction was analyzed using the STRING database (http://string-db.org/) and visualized with Cytoscape software (V3.4.0).

For BSP analysis, genomic DNA of ovaries from each groups was modified with sodium bisulfite using an EZ DNA Methylation-Direct Kit (Cat.#D5020; Zymo Research, Irvine, CA, United States). The PCR products were purified using a gel extraction kit (Cat.#DC301-01; Vazyme, Nanjing, China), then ligated and cloned into the pMD19-T vector (Cat.#6013; Takara, Osaka, Japan). Subsequently, ten positive clones of each sample were randomly selected for DNA sequencing. Data were analyzed and visualized using BIQ Analyzer software. BSP primers are listed in Supplementary Table S1 and all operations were conducted according to the manufacturer’s instructions.

Hu sheep granulosa cells (GCs) were isolated from healthy follicles (2–5 mm) and cultured as our previously described (Yao et al., 2021). Briefly, GCs were seeded into different plates (6-well: 1 × 10⁶ cells/well) in culture medium (Dulbecco’s modified Eagle’s medium/nutrient mixture F-12 supplemented with 10% fetal bovine serum, 2 mM L-glutamine, 100 IU/ml penicillin, and 100 μg/ml streptomycin) at 37°C and 5% CO₂. When the cultured GCs attained 60–70% confluence, the culture medium was replaced with new medium containing various concentrations (0, 1, 2, 5, 10, and 20 μM) of 5-Aza-deoxycytidine (5-Aza; Cat.#a3656; Sigma-Aldrich, United States), as a DNA methyltransferase inhibitor, and incubated for 48 h.

ITGB2 (integrin β2 subunit) knockdown was achieved using siRNAs, which were synthesized by GenePharma (Shanghai, China) and the sequences listed in Supplementary Table S2. After cultured GCs to 60–70% confluence, siRNA-NC and siRNA-ITGB2 were transfected into GCs with or without 10 μM 5-Aza for 48 h. Subsequently, all treated cells were collected for further analysis. All reagents used in cell culture were purchased from Life Technologies (Pleasanton CA, United States).

Immunohistochemistry was performed following our previously described method (Yao et al., 2017). Rabbit anti-LAPTM4B (lysosome-associated protein transmembrane-4 beta; Cat.#18895 -1-AP; ProteinTech, Rosemont, IL, United States) were used as primary antibody. For negative control, the primary antibody was replaced with Tris-buffered saline. Digital images were examined using a light microscope (Nikon, Tokyo, Japan).

For qRT-PCR analysis, cDNA was synthesized using reverse transcription reagent kit with gDNA wiper (Cat.#R323-01; Vazyme, Nanjing, China). Subsequently, qRT-PCR was performed on an ABI 7500 Real-Time PCR System (Applied BioSystems, Carlsbad, CA, United States) using SYBR Green Master Mix (Cat.#Q711-02; Vazyme, Nanjing, China). Relative mRNA expression levels were quantified using 2^–ΔΔCT method, with Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) as the internal control. Primers were designed using the Primer 5.0 software and listed in Supplementary Table S3.

Data are presented as means ± SEM based on three independent experiments, with GraphPad Prism v7.0 (GraphPad Software, CA, United States). Statistical analyses were performed by the one-way analysis of variance followed by a post hoc test in the SPSS25.0 software package (Chicago, IL, United States). P-values < 0.05 were considered statistically significant.

To determine the effects of DNA methylation on the phenotypic variation of prolificacy, WGBS analysis was performed in the ovaries of HPBB, LPBB and LPB + groups of Hu sheep. After data filtering, each sample had approximately 210 million clean reads. Of mapped reads, 74.26, 73.94, and 66.65% in the HPBB, LPBB and LPB + groups, respectively, were used for subsequent analysis, with an average of 3.58% of methylated cytosine sites (Supplementary Table S4).

Global DNA methylation profiles indicated three contexts, CG, CHH, and CHG (where H is A, C, or T), existed in the ovaries of Hu sheep. The genomic methylation proportions of CG, CHG, and CHH contexts was respectively 72.23, 0.50 and 0.50% in the LPB + group, 69.60, 0.63 and 0.63% in the LPBB group, and 70.77, 0.70 and 0.66% in the HPBB group (Figure 1A). Among the three methylated contexts (CG, CHG and CHH), the methylated CG context represented 95.55, 95.20 and 95.06% in the LPB+, LPBB and HPBB groups, respectively (Figure 1B). In addition, the methylated level of CG context was mainly between 90 and 100%, whereas the methylated levels of CHG and CHH contexts apparently exhibited a more uniform distribution range from 10 to 30% (Figure 1C).

The Model-Based Analysis of Bisulfite Sequencing package was used to identify DMRs and compare the DNA methylation profiles between the different groups. Overall, 10,644 DMRs (6,281 hyper-methylated and 4,363 hypo-methylated; 10,634 CG, 1 CHG and 9 CHH), 9,594 DMRs (4,359 hyper-methylated and 5,235 hypo-methylated; 9,583 CG, 2 CHG and 9 CHH) and 12,214 DMRs (4,562 hyper-methylated and 7,652 hypo-methylated; 12,199 CG, 1 CHG and 14 CHH) were identified in the LPBB vs. HPBB, LPB + vs. HPBB and LPB + vs. LPBB groups, respectively (Figure 2A and Supplementary Tables S5–S7). Moreover, the distribution of DMRs were mainly located at intergenic and intron regions (Figures 2B–D).

We identified 7,933 DMGs in the ovaries of the different comparison groups, of which 1,944 DMGs were common to all comparison groups (Figure 3A and Supplementary Table S8). Of these DMGs, 4,721 (3,393 hyper-methylated and 2,670 hypo-methylated), 4,426 (2,694 hyper-methylated and 2,951 hypo-methylated) and 5,152 (2,813 hyper-methylated and 3,812 hypo-methylated) were differentially methylated in the LPBB vs. HPBB, LPB + vs. HPBB and LPB + vs. LPBB groups, respectively (Figure 3B). Moreover, DMGs that exhibited DMRs in their promoter and gene body are shown in Figures 3C–E. The LPBB vs. HPBB, LPB + vs. HPBB and LPB + vs. LPBB groups contained 487, 404 and 528 genes, respectively, including both hyper- and hypo-methylated DMRs in their promoter and gene body.

There were 87 (19 up-regulated and 68 down-regulated), 1,121 (662 up-regulated and 459 down-regulated) and 2,375 (1,563 up-regulated and 812 down-regulated) DE mRNAs identified in the LPBB vs. HPBB, LPB + vs. HPBB and LPB + vs. LPBB groups, respectively (Figures 4A–C). Hierarchical clustering of the DE mRNAs (Figures 4D–F) revealed the expression patterns of the individuals for each comparison.

To confirm the reliability of the WGBS and RNA-seq data, four regions and seven genes were randomly selected for BSP and qRT-PCR, respectively. The results were consistent with the WGBS and RNA-seq data, suggesting that the WGBS and RNA-seq data were reliable for further study (Figure 5).

As previously mentioned, the majority of methylated cytosines were of the CG type; thus, we focused on DMGs of methylated CG for functional enrichment analysis. GO and KEGG analyses were performed to evaluate the functions of DMGs and DEGs in the ovaries of Hu sheep with different prolificacies. Across all comparisons, the DMGs were significantly enriched in the TGF-β and Wnt signaling pathways (Supplementary Figure S1). Furthermore, we selected the female reproduction associated DMGs (Supplementary Tables S9–S11) and DEGs (Supplementary Tables S12–S14) for functional enrichment and the interaction networks construction for each comparison (Figure 6). Interestingly, most of DMGs from all the comparison groups were enriched in the ovarian follicle development/rupture, BMP signaling pathway and ovulation GO terms, as well as the Wnt and TGF-β signaling pathways among all comparison groups (Supplementary Figure S2). Specifically, INHBA, TGFBR2 and SMAD7 genes of the TGF-β pathway and SFRP1, FZD1 and MAP3K7 genes of the Wnt pathway were differentially methylated among all comparison groups.

Similarly, the female reproduction associated DEGs from the LPBB vs. HPBB group were enriched in the hormone biosynthetic process and embryo implantation in GO terms (Supplementary Figure S3A), ae well as the ovarian steroidogenesis, PI3K-Akt and Rap1 signaling pathways (Supplementary Figure S3B). Meanwhile, female reproduction associated DEGs from both the LPB + vs. HPBB and LPB + vs. LPBB groups were enriched in the female gonad and ovarian follicle development GO terms (Supplementary Figures S3C,E), as well as the ovarian steroidogenesis, PI3K-Akt, TGF-β and Wnt signaling pathways (Supplementary Figures S3D,F).

As shown in Figures 7A–C, a significant negative correlation was observed between DNA methylation level around the TSS and gene expression. In the LPBB vs. HPBB, LPB + vs. HPBB and LPB + vs. LPBB groups, we found that 13, 195 and 530 of the DMGs, respectively, associated with DEGs (Figures 7D–F and Supplementary Tables S15–S17). Moreover, the heatmap was constructed to visualize the relationship between DMGs and DEGs in the different comparison groups (Figures 7G–I). Many DMGs contained more than one DMR. For example, five DMRs in NDST4 were hyper-methylated and one DMR in NDST4 was hypo-methylated in the HPBB group compared to that in the LPBB group. In the LPBB vs. HPBB group, we identified 10 hyper-methylated genes with down-regulated expression levels in the HPBB group, including genes related to the Hippo (ITGB2), PI3K-Akt (SYK) and Rap1 (ITGB2) signaling pathways, while three genes were hypo-methylated and up-regulated in the HPBB group compared to that in the LPBB group (Figure 7G and Supplementary Table S15). In the LPB + vs. HPBB group, 77 genes were hyper-methylated and down-regulated in the HPBB group, and these genes were related to the regulation of the TGF-β (ID2), estrogen (LOC101114987) and oocyte meiosis (LOC101112318 and PTTG1) signaling pathways. In contrast, 49 genes were hyper-methylated and down-regulated in the LPB + group, and these genes were related to the TGF-β (TGFBR2), PI3K-Akt (COL4A1, LOC101115805 and COL24A1), FoxO (TGFBR2 and FOXO3) and insulin (CBLB, PDE3A and RIMS2) signaling pathways (Figure 7H and Supplementary Table S16). In the LPB + vs. LPBB group, 153 genes were hyper-methylated and down-regulated in the LPBB group, and these genes were related to the PI3K-Akt (PRKAA2, FGF10, FGF12 and LOC101103187), TGF-β (ID2) and ovarian steroidogenesis (FSHR) pathways. We found that 201 genes were hypo-methylated and up-regulated in the LPBB group compared to that in the LPB + group, and these genes were related to the PI3K-Akt (COL24A1, TLR4, SYK, ITGAV, ITGA7, VWF, RELN, FLT4, KIT, SGK1, THBS2, COL4A1 and LOC101115805), TGF-β (CDKN2B and TGFBR2) and ovarian steroidogenesis (STAR and PRKX) pathways (Figure 7I and Supplementary Table S17). In addition, some genes exhibited coinciding methylation and expression patterns in three comparison groups (Figures 7G–I and Supplementary Tables S15–S17). Interestingly, we found that ITGB2 and LAPTM4B genes, which were hyper-methylated (DMRs in their promoter) and down-regulated in LPBB vs. HPBB group and both LPB + vs. HPBB and LPB + vs. LPBB group, respectively, to confirm the negative correlation between DNA methylation and gene expression.

ITGB2 mRNA expression level in the LPBB group was significantly higher than that in the LPB+ and HPBB groups; with the expression being higher in the HPBB group than in the LPB + group (Figure 8A). LAPTM4B mRNA expression level in the LPB + group was significantly higher than that in the LPBB and HPBB groups (Figure 8B). These results were consistent with the RNA-seq data. Meanwhile, LAPTM4B protein was predominantly localized in the GCs of the antral follicle (Figure 8C).

Next, we evaluated the effect of 5-Aza on ITGB2 and LAPTM4B expression in cultured GCs, and the results showed that ITGB2 and LAPTM4B mRNA expression were significantly higher in the 5-Aza treatment group (Figure 8D). Moreover, following the treatment of GCs with 5-Aza, the methylation level of the LAPTM4B and ITGB2 promoter decreased compared to that of the control cells (Figures 8E,F). Additionally, the mRNA expression of ITGB2 was significantly decreased using siRNA-specific ITGB2 (Figure 8G), subsequently, we demonstrated that the treatment of 5-Aza inhibited ITGB2-specific siRNA-reduced ITGB2 expression in GCs (Figure 8H).