All experimental procedures were performed following the relevant guideline and regulations set by the Ministry of Agriculture of the People′s Republic of China. Ethical approval on animal survival was given by the animal ethics committee of the Institute of Animal Science, Chinese Academy of Agricultural Sciences (No. IAS 2019-49).

TMT technology was used to identify and characterize the differentially abundant proteins (DAPs) in the ovaries of polytocous and monotocous ewes at the luteal and follicular phases, respectively. Each group of biological samples consisted of ovarian tissues isolated from three STH sheep. The protein extraction, TMT labeling, and mass spectrometry were performed for protein identification, followed by bioinformatics analysis (Supplementary Figure S1).

Experimental STH ewes were selected from the Chenglian STH sheep breeding farm in the southwest region of Shandong Province, P. R. China. They were housed under the same condition, with free access to water and food. Firstly, jugular vein blood of healthy non-pregnant ewes (n = 890) was collected which was used for the identification of the FecB genotypes using TaqMan assay (28), thus three genotypes were obtained. The ewes chosen were approximately aged 2–4 years and weighed 65–75 kg. A total of 12 ewes (six polytocous with FecB^BB genotype and six monotocous with FecB⁺⁺ genotype) with similar age and weight indexes were then selected for our study. Estrus synchronization was performed on the whole selected ewes, where a controlled internal drug-releasing device (CIDR, progesterone 300 mg, Zoetis Australia Pty. Ltd., NSW, Australia) was put into these ewes′ vagina for 12 days, and 5 mL of vitamin AD was intramuscularly injected to protect the vaginal epithelium. Ovulation rate (the number of corpus luteum) was detected by the laparoscopy procedure after 7 days from the CIDR removal (luteal phase), which is a key index determining sheep prolificacy. Then, the estrus synchronization mentioned above was performed again after 45 h from the CIDR removal, and ovaries with follicle diameter ≥ 3.5 mm were obtained (follicular phase) (12, 29). The 12 ewes were divided into four groups, including polytocous ewes in the follicular phase (PF, n = 3), polytocous ewes in the luteal phase (PL, n = 3), monotocous ewes in the follicular phase (MF, n = 3), and monotocous ewes in the luteal phase (ML, n = 3). Fresh entire ovarian tissues from 12 ewes were collected, frozen immediately in liquid nitrogen, and stored at −80°C for subsequent proteomic experiments.

Firstly, according to the method given by Tang et al. (25), the protein was extracted by lysis with SDT lysate (including 4% w/v SDS, 150 mM Tris/HCl (pH8), and 100 mM DTT). Also, 30 mg of each ovarian tissue was added to 900 μL SDT; the lysate was sonicated (80 W, 10 cycles of 10 s with 15 s intervals) and then boiled for 15 min. After centrifugation at 14,000 × g for 40 min, the final supernatant was quantified using a BCA assay (Bio-Rad, CA, United States) and detected by SDS-PAGE (Supplementary Figure S2). Protein was treated with the filter-aided sample preparation (FASP) method for enzymatic hydrolysis (30), and the filtrate was collected. The peptide fraction was qualified at 280 nm, OD280. After trypsin digestion, the peptide samples were dried by vacuum centrifugation, then 100 μg peptides of each sample were resuspended in 0.5 M TEAB. The TMT labeling kit (Thermo Fisher Scientific, Waltham, MA, United States) with the isobaric labels 126, 127, 128, and 129 was used to label the PF, PL, MF, and ML samples, respectively.

The peptide mixture was fractionated using the pierce high pH reversed-phase fractionation kit (Thermo Fisher Scientific, Waltham, MA, United States) according to the manufacturer′s instructions.

The freeze-dried samples were dissolved in solvent A (0.1% formic acid, v/v) and loaded onto an Acclaim PepMap 100 pre-column (reversed-phase C18, Thermo Fisher) at a flow rate of 300 nL/min and then separated on an Acclaim pepMap RSLC analytical column (reversed-phase C18, Thermo Fisher). The gradient was comprised of an increase in solvent B (84% acetonitrile with 0.1% formic acid) as follows: 0–35% for 50 min, 35–100% for 5 min, and 100% for 5 min. LC–MS/MS analysis was performed on a Q Exactive mass spectrometer that was coupled to an Easy-nLC 1,000 UPLC system (Thermo Fisher Scientific, Waltham, MA, United States) for 60 min. Positive ions were detected, and the parameters of the precursor ions scan ranged from 300 to 1,800 m/z. Primary MS resolution was 70,000 at 200 m/z, automatic gain control (AGC) was 3e6, maximum inject time was 10 ms, and dynamic exclusion duration was 40.0 s. Secondary MS resolution was 17,500 at 200 m/z, maximum inject time was 1.60 ms, and normalized collision energy was 30 eV, with an underfill ratio of 0.1%.

The raw MS data were output in RAW format and searched using the database (uniprot_Ovis_aries_27500_20170904.fasta), which was converted to the mascot generic format (MGF) file format using Proteome Discoverer 1.4 (Thermo Fisher Scientific, 2012) and retrieved through the Mascot 2.2 server¹ to identify and quantify proteins. The related parameters were as follows: maximum missed cleavages = 2; fixed modifications = Carbamidomethyl (C), TMT-6 plex (N-term), and TMT-6 plex (K); variable modifications = Oxidation (M) and TMT-4 plex (Y); peptide ion mass tolerance = ± 20 ppm; fragment ion mass tolerance = 0.1 Da; database = uniprot_Ovis_aries_27500_20170904.fasta, database pattern = Decoy; and peptide false discovery rate (FDR) ≤0.01 as the cutoff to obtain highly reliable qualitative results. The protein ratios were calculated as the median of only unique peptides of the protein, all peptide ratios were normalized by the median protein ratio, and the median protein ratio should be one after the normalization. The p-values of trusted proteins calculated with Student′s t-test showed significant differences (fold-change >1.2 or < 0.833; adjusted p-value <0.05) between groups.

In order to validate the proteomic data obtained by TMT analysis, the PRM method (a technique based on MS analysis) was used to determine the abundance levels of selected ovarian function-related DAPs (unique peptides ≥2, fold change (FC) >1.2 or < 0.833) (31). After enzymatic hydrolysis, the peptides were desalted, lyophilized, and redissolved in 0.1% formic acid. The concentration of the peptides was determined by OD280. The peptides’ information suitable for the PRM analysis was imported into the Xcalibur™ software (Thermo Fisher Scientific) for the PRM setting. Each sample contained 1 μg peptides, and 200 fmol standard peptides were added (PRTC: GILFVGSGVSGGEEGAR) for chromatographic separation using the HPLC system. Q-Exactive HF MS (Thermo Fisher Scientific) was used for the PRM/MS analysis for 60 min and positive ions were detected. All PRM scans were performed according to the inclusion list after each first-order full mass spectrum scan. The raw data of LC-PRM/MS were analyzed using Skyline software (MacCoss Lab, University of Washington) (32).

GO enrichment² analysis for identifying protein was performed using Blast2GO software, BLAST tool³ as the protein sequence database, mapping, annotation, and annotation augmentation. KAAS software was used to annotate the KEGG pathway,⁴ mapping DAPs to the KEGG genes database; the mapped proteins were categorized based on KEGG orthology. GO and KEGG pathways were executed using Fisher′s exact test (20). Additionally, protein–protein interaction (PPI) networks of DAPs were performed by STRING⁵ combined with Cytoscape software to attempt to capture the underlying protein interaction information.

Ovarian tissues from two comparison groups were isolated and pooled for protein extraction, TMT labeling, and MS analysis. A total of 167,020 peptide spectrum matching numbers and 4,749 unique peptides were obtained in this study (Supplementary Table S1), which were mapped to 5,237 proteins, wherein 732 were uncharacterized and 4,505 were known functions. A complete list of the proteins and peptides detected in the ovary, with the UniProt accession number, protein name, and the score number of unique peptides and amino acids, molecular weight as well as sequence coverage are shown in Supplementary Table S2. Among these results, most distributions of the protein molecular weights (MWs) ranged from 10 to 200 kDa (Figure 1A), and the isoelectric point (pI) values ranged from 4 to 10 (Figure 1B). More than half the protein sequence coverage was ≤15% (Figure 1C). Most identified proteins contained more than one or even up to 201 unique peptides, wherein 29.4% of the proteins had only one unique peptide (Figure 1D).

DAPs among the two comparison groups (PL vs. ML and PF vs. MF) were classified according to proteins with FC > 1.20 and p < 0.05. There were 49 DAPs (9 upregulated and 40 downregulated) and 44 DAPs (18 upregulated and 26 downregulated) in PL vs. ML and PF vs. MF comparison groups, respectively (Figure 2B; Supplementary Table S3). In the luteal phase, the differentially abundant proteins, such as TIA1 cytotoxic granule associated RNA binding protein-like 1 (TIAL1), nicotinamide phosphoribosyl transferase (NAMPT), and cellular retinoic acid-binding protein 1 (CRABP1) were higher in PL than that in ML, whereas complement component 4 binding protein alpha (C4BPA), interferon-stimulated gene 17 (ISG17), yes-associated protein 1 (YAP1), sex hormone binding globulin (SHBG), and uncharacterized protein (W5P9M5) were lower in PL than in ML. In the follicular phase, the differentially abundant proteins, including TIAL1, uncharacterized protein (W5PG55), inhibin beta-A-subunit (A2ICA4), and crystallin mu (CRYM) were higher in PF than in MF, whereas lysosomal associated membrane protein 1 (LAMP1), collagen type XII alpha 1 chain (COL12A1), and CD74 molecule (CD74) were lower in PF than in MF. Hierarchical clustering analysis of DAPs is shown in Figure 3; the similarity of protein abundance patterns within each group was relatively high, whereas it was relatively low between polytocous and monotocous groups, and the groups were clearly divided, which suggested the rationality of DAPs in this experiment.

To identify the potential roles of DAPs in the two comparison groups, enrichment analyses of GO terms and KEGG pathways were performed. In general, GO analysis contains three significant ontologies of biological functions: biological process (BP), cellular component (CC), and molecular function (MF). Among them, each DAP was assigned more than one term; the 20 enriched GO terms related to reproduction are shown in Figure 4. In the PL vs. ML group, the BP terms included positive regulation of signaling, reproductive process, and sexual reproduction (Figure 4A; Supplementary Table S4), as well as predominantly enriched proteins including TIAL1, NAMPT, CRABP1, C4BPA, ISG17, W5P9M5, and SHBG proteins, while gamete generation, embryo development, reproductive process, and cellular process involved in reproduction in multicellular organisms were enriched in PF vs. MF (Figure 4B), including proteins of TIAL1, W5PG55, LAMP1, and A2ICA4.

In addition, DAPs in PL vs. ML and PF vs. MF comparison groups were enriched in 64 and 42 pathways, respectively (Supplementary Table S5). p < 0.05 was set as the cutoff for further enrichment analysis. In the luteal phase, several pathways related to ovarian function were significantly enriched such as nicotinate and nicotinamide metabolism (NAMPT), hippo signaling (W5P9M5), NOD-like receptor signaling (NAMPT, W5NY14, and W5PH55), and RIG-I-like receptor signaling pathways (ISG17 and DDX58) (Figure 5A). In the follicular phase, pathways related to follicular development were significantly enriched, including TGF-beta signaling and signaling pathways regulating the pluripotency of stem cells (A2ICA4) (Figure 5B). In general, GO and KEGG enrichment analyses provided crucial clues on the molecular mechanisms underlying ewes′ prolificacy.

Using the STRING database and Cytoscape software (v.3.9.1), we constructed a PPI network for the DAPs. STRING data is shown in Supplementary Table S6, and the betweenness was selected during CytoNCA analysis. Our results showed that a few DAPs were involved in the PPI, and hub proteins were found. Importantly, YAP1 and NAMPT interacted with other DAPs in PL vs. ML (Figure 6A), while COL12A1, LAMP1, and CD74 interacted with other DAPs in PF vs. MF (Figure 6B), respectively.

Six DAPs (NAMPT, CRABP1, TIAL1, and ISG17 in PL vs. ML, and W5PG55, TIAL1, and LAMP1 in PF vs. MF) containing at least two unique peptides were selected for verification by PRM quantitative analysis. Among the DAPs, FC of NAMPT, CRABP1, TIAL1, and W5PG55 were above 1.2, while ISG17 and LAMP1 were below 0.833. The abundance trend of PRM was consistent with that of TMT-based results, which indicated the credibility of our proteomics data (Figure 7; Supplementary Table S7).