All sheep were normally reared and maintained in good health. Disturbance to the sheep was kept to the minimum, and suffering or injury was not inflicted on the animals during the experiment. All experiments were performed following the relevant animal experimentation guidelines and regulations in China. This study was approved by the Jilin Academy of Agricultural Sciences with approval No. 542256.

In this study, Small Tail Han sheep (multiparous, marked as Mult) and Xinji fine wool sheep (Chinese Merino, singleton, marked as Sgl) were used as multiparous and singleton bearing ewes, respectively. All the sheep were 5.5 years old ewes with three times lambing record. Only the ewes with more than three lambs per litter were selected for Small Tail Han sheep, and 12 healthy animals were included in each group. The animals were raised on the farm in Jilin Academy of Agricultural Sciences (sheep of the same breed were grouped). They were subjected to drinking and feeding ad libitum. The feeding, management and experiments on the animals were carried out according to the regulations approved by the Experimental Animals Committee of China Agricultural University.

Ewes were checked daily for estrous behavior using vasectomized rams in the estrus season (August 2019, Changchun, China). The first estrous performance time of every ewe was accurately recorded as 0 day. During the estrus period of sheep (from 0 to 17.5 days), the blood of jugular vein for Mult and Sgl was collected using a coagulation promoting tube every morning and evening (3 ml each time, 7:00 a.m. and 19:00 p.m.). The blood tubes were left at room temperature for 30 min and then centrifuged at 1,500 × g for 10 min. Serum samples were collected and stored at −20°C. Enzyme-linked immune sorbent assay (ELISA) was carried out to measure the concentration of follicle-stimulating hormone (FSH), luteinizing hormone (LH), estradiol (E2) and progesterone (P4) in serum using an ELISA kit (BNIBT, Beijing, China), following the instructions of the manufacturer. Briefly, a sandwich ELISA was performed by adding standard dilutions or diluted samples to wells coated with antibodies. After washing, an enzyme-linked polyclonal antibody specific was added to the wells. After washing, an enzyme-substrate was added, and after 30 min, the OD at 450 nm was read using an ELISA plate reader.

After routine ovariectomy at 17.5 days (pre-estrus in the next estrus cycle), the ovaries were collected after the first estrous performance time point. After ovariectomy, the development status of ovarian follicles was observed and recorded. Ovary samples of the sheep were washed with normal saline and snap frozen in liquid nitrogen. The ovary samples were divided into two groups: the Small Tail Han sheep group (Mult) and the Xinji fine wool sheep with large follicles group (Sgl). Each group contained three replicates, each having at least three ovaries from three different sheep. Furthermore, the follicle diameter and follicle number were also counted and recorded with three biological replicates and five technical replicates.

According to the instructions of the manufacturer, total RNA was extracted from ovary tissues of Mult and Sgl using TRIzol reagent (Takara, Shiga, Japan). After RNA purification, ribosomal RNA was removed from total RNA samples using a Ribo-zero rRNA Removal Kit (EPICENTRE, Madison, WI, United States). About 10 μg of total RNA from each mixed group was used to prepare libraries. Sequencing libraries were constructed using the NEB Next Ultra Directional RNA Library Prep Kit for Illumina (NEB, Ipswich, MA, United States). Briefly, the RNA was first converted to cDNA using random hexamer primer and M-MuLV Reverse Transcriptase and then cut into small bands of ∼380 bp using a Covaris DNA ultrasonic interrupter. Adaptors were ligated to the bands, and the libraries were prepared. Then all the libraries were sequenced on an Illumina HiSeq 2,500 platform with 125-base pair-end reads in Novogene Biotechnology Co., Ltd. (Novogene, Tianjin, China). For miRNA, Illumina’s TruSeq small RNA library preparation kit was used to prepare the miRNA library from samples. The sequence data were filtered to obtain clean reads by FastQC (v 0.11.4) with the default parameters. The clean data were assembled and compared with the reference genome of sheep (https://sheephapmap.org/news/OARv2p0) using HISAT2 (Kim et al., 2015). The value of FPKM (expected number of fragments per kb per million reads) of reads in each sample was calculated using Cufflinks (version 2.2.1). The differentially expressed genes (DEGs) by different comparisons were identified using DESeq2 (v1.6.3) (Love et al., 2014). Differential expression analysis was performed on the negative binomial distribution test and Benjamini–Hochberg method. An adjusted p-value (Q-value) <0.05 and |log2 fold change (FC)| >2 was considered as the thresholds for screening DE-mRNA, DE-miRNA, and DE-lncRNA significantly (Qian et al., 2017). Genes identified as log2FC >1 and log2FC <−1 were identified as up and downregulated DEGs, respectively. A heatmap was drawn on the DE-mRNA, DE-miRNA, and DE-lncRNA using the pheatmap R package (Kolde and Kolde, 2015).

Ovary tissues were lysed in sodium deoxycholate (SDC) lysis buffer containing 4% (wt/vol) SDC and 100 mM Tris–HCl (pH 8.5). Proteome preparation was done using StageTip (iST) method (Kulak et al., 2014). Samples were separated using HPLC in a single run (without pre-fractionations) and analyzed with MS. The peptides were separated on a reverse-phase column, 50 cm packed in-house with 1.9-μm C18-Reprosil-AQ Pur reversed-phase beads (Dr Maisch GmbH, Ammerbuch, Germany) for about 120 min (single-run proteome analysis). Eluting peptides were electrosprayed and analyzed using tandem MS on a Q Exactive HF (Thermo Fischer Scientific) using higher-energy collisional dissociation (HCD)-based fragmentation, which was set to alternate between a full scan followed by up to five fragmentation scans. Raw MS data were analyzed with MaxQuant (version 1.5.1.6) (Cox and Mann, 2008) and using the Andromeda engine for database search. The MS/MS spectra were matched against the UniProt database (https://www.uniprot.org), with a false discovery rate (FDR) of <1% at the level of proteins, peptides, and modifications. The search results were then filtered using a cutoff of 1% for the peptide false identification rate. Label-free quantification was used to quantify the proteins. For quantitative changes, a cutoff of ≥1.5 or ≤0.666-fold change and Q-value <0.05 were set for differentially expressed proteins (Wei et al., 2018).

The pIRES2-ZsGreen1 vector sequence was purchased from Takara Bio (Takara, Shiga, Japan), and the structure of the vector is shown in Supplementary Figure S1. The plasmid was prepared by inserting cDNAs from the mRNA (HSD17B1) of interest and lncRNA (MSTRG.28645) into pIRES2-ZsGreen1 plasmid at Xho I and BamH I sites using restriction enzyme (TaKaRa, Otsu, Shiga, Japan) and T4 DNA ligase (NEW ENGLAND BioLabs, Ipswich, MA, United States). After identification and replication, the extracted plasmid was prepared for the transfection of granular cells. According to the instructions of the manufacturer, the transient transfection was achieved by Lipofectamine TM 3000 (Thermo Fisher Scientific, CA, United States) reagent. The siRNA oligonucleotides were synthesized and purified by Genomeditech Inc. (Shanghai, China), which also provided the scrambled negative control. According to the instructions of the manufacturer, granulosa cells (1 × 10⁶ cells/well) were grown to 80% confluence and transfected with 5 nM CD36 siRNA for 6 h using Lipofectamine 2000 (Life Technologies).

The MTT assays were used to determine cell proliferation. The cells were seeded in a 96-well plate with a density of optimized cell number (5,000 cells/well). After 48 h of seeding, 20 µl of MTT (5 mg/ml) was added to the wells. Four hours later, the mixed medium was replaced with 150 µl of dimethyl sulfoxide (Sigma, St. Louis, MO, United States). The 96-well plate was then agitated for 15 min at room temperature. The OD value of each well was measured using a fluorescence microplate reader (Sunrise Remote, Tecan Austria GmbH, Grödig, Austria) at a wavelength of 490 nm. An AnnexinV-APC/7-AAD cell apoptosis detection kit (Nanjing Kaiji Biology, KGA1026) was used to detect cell apoptosis for cytometry analysis. The cells were briefly digested using trypsin and then washed using pre-cooled PBS. The binding buffer suspended cells were added to the flow tube along with 5 µl ANNEXIV-APC and 5 µl 7-AAD. The final solution was shielded from light exposure and left to react at room temperature for 5–15 min. Last, apoptosis was assessed using flow cytometry in each group.

Detailed procedures for qRT-PCR were carried out as previously described by Zou et al. (2016). All samples were performed in triplicate, and the 2^−ΔΔCT calculation method was used to analyze the mRNA expressions after normalization with glyceraldehyde-3-phosphate dehydrogenase or actin housekeeping gene expressions, where necessary. The primers used in qRT-PCR are shown in Supplementary Table S1. qRT-PCR conditions were 95°C for 3 min followed by 40 cycles of 95°C for 15 s, 60°C for 40 s.

Statistical analysis was performed using SPSS statistical software (version 22.0). Statistical analysis was performed by using Student’s t-test to compare two groups and analysis of variance (ANOVA) to compare three or more groups. A value of p < 0.05 was considered statistically significant. The 2^−ΔΔCT method was used to determine the relative expression levels.

The concentration and trend of FSH, LH, E2, and P4 in serum was used to indicate the difference of the main reproductive hormones between the two breeds of sheep. Generally, it was found that the content of FSH in Mult was higher than that in Sgl sheep (Figure 1A). There was no significant difference in the content of LH between the two groups (Figure 1B). At most sampling time points, it was evident that Mult had a higher concentration of E2 than Sgl with the exception of 3, 5, 10, 12, and 13 days (Figure 1C). The content of P4 was found to be higher in Mult than in Sgl during 3–5, 8–11, and 12–17 days (Figure 1D). Furthermore, it was revealed that all these hormones were secreted in pulses. In summary, Mult had a higher average concentration of FSH, E_2, and P₄ than Sgl.

Rising hormone secretion levels may affect folliculogenesis and ovarian development. Therefore, the current study examined folliculogenesis between ovaries in Mult and Sgl. There were differences between ovaries in Mult and Sgl. Before ovulation, there was no significant difference in the number of small follicles (less than 3.5 mm in diameter) between Mult and Sgl groups (Table 1, p > 0.05). The number of large follicles (more than 3.5 mm in diameter) in the Mult group was more than that in the Sgl group (Table 1, p < 0.05). Nonetheless, there were large follicles on the ovaries of both sides. In the Sgl group, the ovary containing a large follicle was found only on one side of the ovary.

Transcriptome analysis was conducted to delineate the genes and pathways difference in Mult and Sgl sheep. A total of 855.71 Mb reads were obtained from the nine ovarian samples. The summary of the RNA-seq reads for three cDNA libraries is listed in Supplementary Table S2. After filtering, a total of 62 canonical genes were identified as differentially expressed genes (DEGs) between Mult and Sgl sheep (Figure 2A). It was found that the DEGs were significantly enriched in steroid hormone biosynthesis (ko00140), ovarian steroidogenesis (ko04913), glycerolipid metabolism (ko00561), and ECM–receptor interaction (ko04512) (Figure 2B).

The details of DEGs and enriched pathways are listed in Supplementary Table S3. All the DEGs, HSD17B1, HSPA1A, INSRR, and NUCB2 genes showed the largest differences in expression between the two groups (Figure 2C). It was evident that the HSD17B1 gene was involved in both the steroid hormone biosynthesis and ovarian steroidogenesis pathways. Nonetheless, all the other genes involved in the two pathways did not show significant differences in their expression. This study hypothesized that the HSD17B1 gene is one of the genes that play an essential role in the biological processes related to sheep ovulation. The expression difference of HSD17B1 between Mult and Sgl sheep may be related to litter size.

In addition, a total of 12 lncRNAs were identified in the current study. Of the 12 lncRNAs, 10 lncRNAs were significantly expressed in both Mult and Sgl (Figure 2D). Furthermore, to screen the significant lncRNAs, the interaction network between mRNAs and lncRNAs was analyzed based on the published criteria (Liao et al., 2011). As shown in Figure 2E, there were 1,118 sets of correlations between lncRNAs and mRNAs. It was found that MSTRG.28645 was significantly related to 11 lncRNAs and 169 mRNAs in the co-expression network. Consequently, MSTRG.28645 should be profiled as a key lncRNA because it had the most links with other genes.

Proteomic profiling revealed the significant differences between Mult and Sgl sheep. Of the approximately 2,100 proteins identified in each group, it was found that there were 2,031 known protein IDs. Principal component analysis (PCA) showed that the samples of the same group were clustered together (Figure 3A). There were 60 DEPs identified between Mult and Sgl sheep, 40 of which were more highly expressed in Sgl, whereas 20 DEPs were more expressed in Mult sheep (Figure 3B; Supplementary Table S4). Six DEPs (POSTN, DARS, VCL, TTN, SRRM1, and ERC1) were revalidated with the parallel reaction monitoring (PRM) method.

The results of this study showed that the protein expression trends detected with the two methods were consistent (Figure 3C). Furthermore, the protein–protein interaction network was constructed by enrichment of DEPs through signaling pathways and biological functions (Figure 3D). It was found that there were 11 proteins (POSTN, ERC1, DARS, PRKG1, MARS, SMS, TNXB, FAM213B, OAS2, SUCLA2, and CYP17A1) that interacted with each other through multiple signal pathways and biological functions. The CYP17A1 played a major role in multiple pathways as a hub protein. Therefore, both CYP17A1 protein and HSD17B1 gene were involved in the regulation of steroid hormone biosynthesis and ovarian steroidogenesis. This confirms that the two pathways play a crucial role in the litter size mechanism in sheep.

Among all the 62 DEGs identified from the canonical genes, the encoded proteins corresponding to 13 genes were detected using label-free quantification. The downregulation trend of HSPA1A, NUCB2, CFH, RBM3, LDHA, and ACTN2 were opposite in transcriptome and proteome, while seven genes (HSD17B1, FKBP5, SERPINE1, KRT8, CRABP1, TOP1, and THY1) were consistent in the expression trend of transcriptome and proteome (Figure 4). Previously, HSD17B1 gene has been shown to participate in the reproductive regulation of mice and rats (Hakkarainen et al., 2015; Yuan et al., 2021). Therefore, in the current study, HSD17B1 was chosen as the key regulator for subsequent verification. Furthermore, the lncRNA MSTRG.28645 was also treated as a key lncRNA for the subsequent validation.

Comprehensive gene expression and gene annotation finally selected HSD17B1 as the candidate gene for the functional verification. As expected, the expression of HSD17B1 was upregulated and downregulated in the overexpressed and siRNA interfered granulosa cells, respectively (Figure 5A). The flow cytometry analysis results showed that the rate of apoptosis in granulosa cells for the HSD17B1 interfered (5.50% ± 0.21%) and overexpressed cells (10.64% ± 0.66%) were significantly lower and higher than in the control group (7.18% ± 0.11%), respectively (Figures 5B–E). In addition, the MTT method was also performed to detect the proliferation rate of sheep granulosa cells after HSD17B1 overexpression and interference treatment.

Results of this study showed that there was no significant difference between interference, overexpression, and the control groups (Figure 5F). The E2, P4, ACT, INH, and FS levels were then detected in HSD17B1 interference and overexpression cells. The results showed that the HSD17B1 gene could promote the secretion of E2, P4, and ACT, but inhibit the secretion of INH and FS (Figures 5G–K). By analyzing the location of MSTRG.28645 in the sheep genome, it was found that MSTRG.28645 was located on chromosome 6. The results of this study revealed that MSTRG.28645 can promote the apoptosis rate after successful overexpression and interference of MSTRG.28645 in granulosa cells (Figure 6A–E). In addition, it was evident that MSTRG.28645 also inhibits the secretion of E2 and P4 (Figures 6G–K). To further verify the role of MSTRG.28645, the expression of genes that were related to follicular development including LHR, FSHR, Erβ, ESR, and AMH were detected as well as the key genes in the TGF-β signaling pathway (TGFβ1 and TGFβ2). Moreover, it was found that SFR1 had a promoting effect on LHR, FSHR, and ERβ in sheep granulosa cells as well as an inhibitory effect on TGFβ1 and TGFβ2 (Figure 6L).