Integrative Proteomics and Transcriptomics Profiles of the Oviduct Reveal the Prolificacy-Related Candidate Biomarkers of Goats (Capra hircus) in Estrous Periods

The oviduct is a dynamic reproductive organ for mammalian reproduction and is required for gamete storage, maturation, fertilization, and early embryonic development, and it directly affects fecundity. However, the molecular regulation of prolificacy occurring in estrous periods remain poorly understood. This study aims to gain a better understanding of the genes involved in regulating goat fecundity in the proteome and transcriptome levels of the oviducts. Twenty female Yunshang black goats (between 2 and 3 years old, weight 52.22 ± 0.43 kg) were divided into high- and low-fecundity groups in the follicular (FH and FL, five individuals per group) and luteal (LH and LL, five individuals per group) phases, respectively. The DIA-based high-resolution mass spectrometry (MS) method was used to quantify proteins in twenty oviducts. A total of 5409 proteins were quantified, and Weighted gene co-expression network analysis (WGCNA) determined that the tan module was highly associated with the high-fecundity trait in the luteal phase, and identified NUP107, ANXA11, COX2, AKP13, and ITF140 as hub proteins. Subsequently, 98 and 167 differentially abundant proteins (DAPs) were identified in the FH vs. FL and LH vs. LL comparison groups, respectively. Parallel reaction monitoring (PRM) was used to validate the results of the proteomics data, and the hub proteins were analyzed with Western blot (WB). In addition, biological adhesion and transporter activity processes were associated with oviductal function, and several proteins that play roles in oviductal communication with gametes or embryos were identified, including CAMSAP3, ITGAM, SYVN1, EMG1, ND5, RING1, CBS, PES1, ELP3, SEC24C, SPP1, and HSPA8. Correlation analysis of proteomics and transcriptomic revealed that the DAPs and differentially expressed genes (DEGs) are commonly involved in the metabolic processes at the follicular phase; they may prepare the oviductal microenvironment for gamete reception; and the MAP kinase activity, estrogen receptor binding, and angiotensin receptor binding terms were enriched in the luteal phase, which may be actively involved in reproductive processes. By generating the proteome data of the oviduct at two critical phases and integrating transcriptome analysis, we uncovered novel aspects of oviductal gene regulation of fecundity and provided a reference for other mammals.


Introduction
Goats (Capra hircus) have been raised for a long time to satisfy the human requirements for meat, and with the rapid development of the China economy, goat meat consumption has been gradually increasing [1]. High fecundity in goats is a desirable and genetically correlated trait in mutton production. Many factors that affect fertility, of which the maternal environment is an important cause [2]. Previous work has started to show that maternal interactions with gametes and embryos are fundamental to pregnancy and influence number in goats by combining proteomic and transcriptomic techniques, which revealed the molecular functions of the oviduct under different estrous cycles, as well as the key components regulating fecundity in goats.

Overview of Quantitative Proteomics Analysis of Goat Oviducts
The median intra-group coefficient of variation (CV) distribution for the four-quality control (QC) replicates was around 20% ( Figure 1A). The Pearson correlation coefficient at the protein level between every two QC samples indicated that the replicates of the four QC samples were highly correlated, with a Pearson correlation coefficient of approximately 0.9 ( Figure 1B). Principal component analysis (PCA) showed that samples from four different groups and QC samples could be separated, and samples from the same stage clustered well ( Figure 1C). The DIA-based quality control is shown in Supplementary Figure S2. The false-discovery rate (FDR) method was applied to conduct multiple testing corrections. A q value of 0.01 was set as the threshold, which was equivalent to FDR 0.01. In this study, the C-score was 1.024 with a q value of <0.01, indicating that the qualitative results are highly reliable ( Figure 1D). Furthermore, we counted the protein expression abundance of the twenty samples, and the results showed that the protein expression was relatively consistent among the five biological replicates in different groups ( Figure 1E). These results confirmed the high degree of reproducibility of MS data acquired by the DIA method.
In this study, each sample was quantified, and we selected proteins at constant levels found in more than 50% of samples for follow-up statistics ( Figure 1F, Supplementary Table S1). In total, 81,147 peptides and 8536 proteins were quantified in the data-dependent acquisition (DDA) spectral library ( Figure 1G). The DIA data with the DDA spectra library were analyzed by Spectronaut Pulsar X; about 76,197 peptides and 5409 proteins were further quantified ( Figure 1G, Supplementary Data S1). The quantitative heat map of the total DIAidentified proteins from all samples is shown in Supplementary Figure S2D. In addition, we statistically counted the overlap of proteins between different groups, and 4732 proteins overlapped in the four groups ( Figure 1H).

Oviduct Protein Co-Expression Network Analysis Identifies Module-Trait Relationships and Points to Core Mediators of Fecundity
In total, 5076 proteins were screened from the 5409 to reveal the expression patterns of proteins inside various WGCNA modules following the determination of a soft threshold at R 2 = 0.9 (Supplementary Data S2). Thereafter, 45 modules were identified through the average linkage hierarchical clustering in the four groups and are visualized in Figure 2A (Supplementary Data S2). The co-expressed protein heat map revealed that proteins in the same model had a high biological relevance ( Figure S3A). To assess the significance of the modules with prolificacy traits, we correlated the modules with fecundity hallmarks of goats in the different estrus cycles ( Figure S3B; Supplementary Data S2). In the follicular phase, the blue module (475 proteins, cor = −0.45, and p = 0.04) negatively correlated with the trait of high fecundity (FH); light-yellow, brown, and magenta modules (69,385, and 119 proteins, respectively; cor = 0.5, and p = 0.02) showed the highest correlation with low fecundity (FL). In the luteal phase, the tan module (96 proteins, cor = 0.76, and p = 9 × 10 −5 and steel-blue module (49 proteins, cor = 0.55, and p = 0.01) showed the strongest correlation with the trait of high (LH) and low (LL) fecundity, respectively. These modules suggest that the proteins in the modules were significantly associated with the corresponding trait. Of note, the tan module had a significantly positive correlation with the high-fecundity trait in the luteal phase and was selected for further analysis (p < 0.05) ( Figure 2B, Supplementary Figure S3B). The scatter plot demonstrated that the significance and module membership of the correlated proteins in the tan modules were significantly associated with LH, as shown in Figure 2C.   We built a protein-protein interaction (PPI) network using Cytoscape software (V.3.9.0) for the tan module with several hub proteins with the most connectivity ( Figure  2D, Supplementary Data S2). Some of the top hub proteins in this module are essential for cell migration; germ-cell and embryo development; reproductive hormone secretion; and oviduct function, including NUP107, ANXA11, TEX264, DPAGT1, MED20, COX2, AKAP13, IFT140, and Q58LV0 (somatotropin), which may create an optimal We built a protein-protein interaction (PPI) network using Cytoscape software (V.3.9.0) for the tan module with several hub proteins with the most connectivity ( Figure 2D, Supplementary Data S2). Some of the top hub proteins in this module are essential for cell migration; germ-cell and embryo development; reproductive hormone secretion; and oviduct function, including NUP107, ANXA11, TEX264, DPAGT1, MED20, COX2, AKAP13, IFT140, and Q58LV0 (somatotropin), which may create an optimal intrafallopian environment for early embryo development. In addition, the gene ontology (GO) databases were used to analyze the functional of tan module proteins based on three categories, including biological process (GO BP), cellular component (GO CC), and molecular function (GO MF) ( Figure 2E, Supplementary Data S3). In the GO BP category, reproduction and biological adhesion were associated with oviduct functions. Binding, catalytic activity, and the molecular function regulator were the three groups that were most prevalent in GO MF. The most prevalent GO CC categories were related to cells, cell parts, organelles, and membranes. Mitogen-activated protein kinase (MAPK14), POC1 centriolar protein A (POC1A), angiotensin-converting enzyme (ACE), nuclear pore complex protein (NUP107), tektin 4 (TEKT4), F-box and WD repeat domain containing 8 (FBXW8), and NAD(P) dependent steroid dehydrogenase-like (NSDHL) were enriched in the reproductive process, reproduction, and development process, which may be significantly associated with the fecundity of the goat. As for the KEGG enrichment analysis, just the N-Glycan biosynthesis pathway was significantly enrichened (p < 0.05) ( Figure 2F, Supplementary Data S3), and GlcNAc-1-P transferase (DPAGT1) and dolichol-phosphate mannosyltransferase subunit 1 (DPM1) participated in this pathway.

Identification of Differentially Abundant Proteins
A total of 4975 and 4851 proteins were identified in the FH vs. FL and LH vs. LL comparisons, respectively (Table S1), while 98 and 167 DAPs were identified in each comparison group. In the FH vs. FL comparison group, 15 of these DAPs showed upregulation, 83 showed downregulation ( Figure 3A), and 52 showed upregulation; 115 showed downregulation in the LF vs. LL comparison group ( Figure 3B). DAPs were further selected to build a heatmap ( Figure 3C,D), which revealed the expression patterns in the follicular and luteal phases of goat oviducts. Moreover, the similarity of data patterns within each group is high, while it is low between groups, which effectively differentiates the groups. The Venn diagram of the DAPs showed that just one downregulated DAP was co-expressed between the two phases ( Figure 4A). These results suggested that the follicular phase of the oviduct and luteal phase have different protein profiles in response to the prolificacy trait of goats.
There are a few key DAPs that enriched in the follicular phase (Supplementary Data S1), demonstrating the important roles the oviduct plays. For example, proteins involved in the energy processes (ATP synthase protein 8, ATP8; ATPase H+ transporting V1 subunit C1, ATP6V1C1; ATP synthase membrane subunit f, ATP5MF; mannosidase alpha class 2A member 1, MAN2A1), the cellular adhesion and transporters (adhesion G protein-coupled receptor E5, ADGRE5; actin filament associated protein 1 like 2, AFAP1L2; mediator of cell motility 1, MEMO1), and the calcium regulatory (calmodulin-regulated spectrin-associated protein family member 3, CAMSAP3; calcium/calmodulin-dependent serine protein kinase, CASK). We also found that reproductive-related proteins, such as TAB1 (TGF-beta activated kinase 1 (MAP3K7) binding protein 1), were highly enriched in the FH vs. FL comparison group. For the luteal phase, such as cytoskeleton (cytoskeleton-associated protein 5, CKAP5; SWI/SNF related, matrix associated, actin dependent regulator of chromatin, subfamily a, member 1, SMARCA1) and metabolic transporters (sphingolipid transporter 1, SPNS1; solute carrier family 50 members 1, SLC50A1; actin-related protein 2/3 complex subunit 3, ARPC3; transmembrane p24 trafficking protein family member 8, TEMD8; trafficking protein particle complex subunit 5, TRAPPC5; target of myb1 like 2 membrane trafficking protein, TOM1L2; kinesin family member 3A, KIF3A). Furthermore, cytochrome b (CYTB), dedicator of cytokinesis 6 (DOCK6), cyclin-dependent kinase 6 (CDK6), mitogen-activated protein kinase kinase 5 (MAP2K5), cilia and flagella associated protein 69 (CFAP69), myosin heavy chain 14 (MYH14), epithelial splicing regulatory protein 2 (ESRP2), activating signal cointegrator 1 complex subunit 3 (ASCC3), and remodeling and spacing factor 1 (RSF1), which may be involved in embryo and gamete development, embryo-maternal communi- . The X-axis represents the fold change in protein expression of different comparison groups. The Y axis represents the statistical significance of proteins. The blue dots represent the proteins that were not significantly downregulated, the red dots represent proteins that were significantly upregulated, and the gray dots represent proteins that were not significantly different. Cluster analysis of DAPs in the group of FH vs. FL (C) and LH vs. LL (D). The hierarchical clustering results are represented as a tree heat map, with the ordinate representing significantly differentially abundant proteins and the abscissa representing sample information. Red represents significantly upregulated proteins, blue represents significantly downregulated proteins, and gray represents no quantitative information for proteins. . The X-axis represents the fold change in protein expression of different comparison groups. The Y axis represents the statistical significance of proteins. The blue dots represent the proteins that were not significantly downregulated, the red dots represent proteins that were significantly upregulated, and the gray dots represent proteins that were not significantly different. Cluster analysis of DAPs in the group of FH vs. FL (C) and LH vs. LL (D). The hierarchical clustering results are represented as a tree heat map, with the ordinate representing significantly differentially abundant proteins and the abscissa representing sample information. Red represents significantly upregulated proteins, blue represents significantly downregulated proteins, and gray represents no quantitative information for proteins.  There are a few key DAPs that enriched in the follicular phase (Supplementary Data S1), demonstrating the important roles the oviduct plays. For example, proteins involved in the energy processes (ATP synthase protein 8, ATP8; ATPase H+ transporting V1 subunit C1, ATP6V1C1; ATP synthase membrane subunit f, ATP5MF; mannosidase alpha class 2A member 1, MAN2A1), the cellular adhesion and transporters (adhesion G protein-coupled receptor E5, ADGRE5; actin filament associated protein 1 like 2, AFAP1L2; mediator of cell motility 1, MEMO1), and the calcium regulatory (calmodulin-regulated spectrin-associated protein family member 3, CAMSAP3; calcium/calmodulin-dependent serine protein kinase, CASK). We also found that reproductive-related proteins, such as TAB1 (TGF-beta activated kinase 1 (MAP3K7) binding protein 1), were highly enriched in

Bioinformatic Analysis of DAPs
The subcellular localization, domain enrichment, GO annotation, and Kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis were used to understand the functions of the DAPs. We performed subcellular fractionation analyses, which revealed that 98 and 167 DAPs were classified into six subcellular localization categories in the FH vs. FL and LH vs. LL groups, respectively. In the FH vs. FL comparison groups, a total of 53 DAPs were localized to the nuclear (54.08%), followed by the cytoplasmic (30.61%) ( Figure 4C). Regarding the LH vs. LL comparison, the subcellular localizations were like those in the FH vs. FL comparison of DAPs. Domain enrichment revealed that DAPs between FH and FL mainly included the LSM domain, protein phosphatase 2C, and some other domains (p < 0.01) ( Figure 4D). Meanwhile, the most enrichened domains in the LH vs. LL were the PH domain, RhoGAP domain (p < 0.01), etc. ( Figure 4E). It is important to emphasize that in sheep, changes in PH may be associated with periodic changes in mucin secretion in the oviduct [31].
GO enrichment analysis of all DAPs revealed that they were mostly localized in the cell part, and a large number of DAPs were mainly involved in molecular function with binding ability, catalytic activity, and molecular function regulation (Supplementary Data S3). Among the annotated terms of the FH vs. FL comparison group ( Figure 5A), in the group within the GO BP category, a few DAPs were localized to biological adhesion, locomotion, cell proliferation, and growth, while many DAPs were assigned to respond to the negative regulation of biological processes. At the GO CC category level, 16 major processes were discovered from the DAPs, including the cell part, cell, and organelle. Binding, catalytic activity, and molecular function regulators were predominant in GO MF. Meanwhile, the main enriched GO categories in the case of LH vs. LL ( Figure 5B) were essentially the same as in the FH vs. FL comparison groups. Noticeably, spindling 1 (SPIN1), cystathionine beta-synthase (CBS), and cilia and flagella-associated protein 69 (CFAP69) were enriched in the BP category of the reproductive process (GO:0022414) and reproduction (GO:0000003), similar to the LH vs. LL comparison groups.
As for the KEGG pathway analysis in the FH vs. FL, DAPs were enriched in protein processing in the endoplasmic reticulum, necroptosis, mitophagy-animal, and NF-kappa B signaling pathway (p < 0.01) ( Figure 5C, Supplementary Data S3). In terms of LH vs. LL comparison, the pathways with the highest number of DAPs were the PI3K-Akt signaling pathway, protein processing in the endoplasmic reticulum, and focal adhesion (Supplementary Data S3). Furthermore, aldosterone synthesis and secretion, one carbon pool by folate, and the renin-angiotensin system were significantly enriched (p < 0.05) in KEGG terms belonging to the endocrine system ( Figure 5D). Interestingly, the N-glycan biosynthesis pathway was also enriched in the tan module obtained by WGCNA analysis, reminding us that it is a specific function in the oviduct and needs further investigation. As noted above, the GO categories and KEGG pathways provided further clues to the molecular mechanisms of the oviduct underlying goat prolificacy.

Protein-Protein Interaction (PPI) Network Analysis for Goat Prolificacy Trait
In the FH vs. FL comparison, an interaction network was constructed ( Figure 6A, Supplementary Data S3), which can be divided into three clusters. It is shown that there were significantly enriched interactions among 41 DAPs, and SNRPG, SSB, SNRPD1, SNRNP70, TRMT61A, RRS1, SEC63, SYVN1, EMG1, RRS1, ERLEC1, ND5, DAP3, and RING1 as the hubs in the network. Meanwhile, in the LH vs. LL comparison, the interaction network was composed of 96 DAPs and divided into five clusters, among which PES1, GART, APP, TSR1, CBS, ELP3, SSP1, CAPZB, CYTB, CKAP5, ARMC8, SEC24C, RPL21, HSPA8, and STK36 were the center of the interaction network of each cluster ( Figure 6B). In addition, trifunctional purine biosynthetic protein adenosine-3 (GART), heat shock protein family A (Hsp70) member 8 (HSPA8), F-actin-capping protein subunit beta (CAPZB), SEC24 homolog C, COPII coat complex component (SEC24C), and cystathionine beta-synthase (CBS) were also shown to be key proteins.  As for the KEGG pathway analysis in the FH vs. FL, DAPs were enriched in protein processing in the endoplasmic reticulum, necroptosis, mitophagy-animal, and NF-kappa B signaling pathway (p < 0.01) ( Figure 5C, Supplementary Data S3). In terms of LH vs. LL comparison, the pathways with the highest number of DAPs were the PI3K-Akt signaling pathway, protein processing in the endoplasmic reticulum, and focal adhesion (Supplementary Data S3). Furthermore, aldosterone synthesis and secretion, one carbon pool by

Western Blot Analysis of Hub Proteins
Given that differently characterized proteins were observed in the oviducts in the estrous cycles, we set out to identify specific protein expression signatures that could be used as future specific biomarkers affecting reproduction in the oviducts of goats. The expression of SYVN1, EMG1, HSPA8, and PES1 were subsequently validated in oviduct tissues by Western blotting and GAPDH as the reference protein. Coincidentally, the expression patterns of the four proteins in the Western blot experiments were highly consistent with the quantitative proteomics through the significant analysis ( Figure 6C-F). molecular mechanisms of the oviduct underlying goat prolificacy.

Combined Analysis of Proteome and Transcriptome (DAPs/DEGs) Responding to Prolificacy Trait of Yunshang Black Goat in Oviducts
Comparing the two phases, our transcriptomics and proteomics results showed that many more transcripts could be identified than proteins in the oviduct. In total, 62,322 and 62,591 genes (Supplementary Data S5) were detected from an RNA-seq analysis in goat oviduct in the follicular and luteal phases after stringent quality checks and data analysis, respectively, whereas the number of detected proteins was 4975 and 4851 (Supplementary Data S1). An integrated analysis of proteomic and transcriptomic profiling screened 33,199  We focused GO analysis on the differentially expressed genes identified in the each proteomic and transcriptomic in the follicular and luteal phases, respectively ( Figure 7C,E), of which were more significant at the mRNA level than the protein level ( Figure 7D,F). These findings suggest that fecundity-related changes at the transcriptional level in the oviduct are more prevalent. As shown in Figure 7C,E, the DEGs/DEPs shared a common GO annotation pattern. For instance, proteolysis, intrinsic apoptotic signaling pathway, regulation of macromolecule metabolic process, regulation of the metabolic process, and regulation of cellular metabolic process were enriched in the follicular phase ( Figure 7C). In the luteal phase, the combined enrichment of DEGs and DAPs was mainly shown in the regulation of MAP kinase activity, Rab GTPase binding, estrogen receptor binding, angiotensin receptor binding, and RSF complex ( Figure 7E). These results suggest that many of the altered biological processes in the oviduct are driven, at least partially, by a combination of transcription and translation, and may be in response to goat fertility.

RT-qPCR Validation of Candidate DEGs
To further confirm the differences in gene expression at the transcriptome level, we selected 4 genes from the correlated DEGs in the FH vs. FL comparison to perform RT-qPCR validation with three replicates of each sample ( Figure 8A), including JAK1, RBM3, SRSF5, and MATR3. Simultaneously, four genes were also detected in the LH vs. LL comparison group ( Figure 8C), including RPL21, GPT2, SPOCK1, and APOD. The mean expression ratio was comparable between the RT-qPCR and RNA-Seq datasets. Linear regression analysis between RT-qPCR and RNA-seq showed an overall correlation coefficient (R) = 0.61 and 0.96 ( Figure 8B,D) for FH vs. FL and LH vs. LL, respectively, which indicates the good correlation of these eight genes were consistent in transcriptomics databases.  parison group (Figure 8C), including RPL21, GPT2, SPOCK1, and APOD. The mean expression ratio was comparable between the RT-qPCR and RNA-Seq datasets. Linear regression analysis between RT-qPCR and RNA-seq showed an overall correlation coefficient (R) = 0.61 and 0.96 ( Figure 8B,D) for FH vs. FL and LH vs. LL, respectively, which indicates the good correlation of these eight genes were consistent in transcriptomics databases.

Discussion
Mammalian reproduction is a complex process with coordinated regulation at multiple reproductive organs and molecular levels. In this study, mechanisms involved in the prolificacy trait of goats induced by oviduct with estrous periods were characterized using the DIA-based proteomic method and bioinformatic analyses. Results revealed that a deep proteome map of the oviduct in Yunshang black goats. We demonstrate that the known effects of oviducts, such as DAPs indicative of biological adhesion, gamete and embryo transport, and embryonic development, are present at the follicular or luteal phase. In addition, we also used a combined proteomic and transcriptomic profiling approach to identify coordinated-yet-specific changes in the oviduct during estrous periods, revealing the goat oviduct's molecular mechanisms. Like previous studies that have been performed on other organisms [21,[32][33][34], we utilized combinatorial functional mapping of goat oviducts.
Weighted gene co-expression network analysis (WGCNA) has the potential to reconstruct significant biological pathways from proteomics data by grouping all proteins into modules, regardless of fold-changes or statistical significance, while remaining independent of existing databases of functional annotation enrichment and protein-protein interactions [35]. The protein co-expression network shows that the module related to fecundity in the luteal phases is led by the hub proteins, including NUP107, ANXA11, MED20, COX2, AKAP13, and IFT140. Among these proteins, a functional study revealed that nucleoporin 107 (NUP107) is phosphorylated on its N-terminal sites in vivo during mitosis and can be efficiently dephosphorylated by trimeric protein phosphatase 2A-B55 upon mitotic exit, suggesting that NUP107 is involved in cell cycle progression regulation [36]. As a member of the annexins family, annexin A11 (ANXA11) belongs to Ca 2+ -regulated phospholipid-binding proteins widely expressed in eukaryotic cells, and are always implicated in mediating vesicle trafficking, cell division, differentiation, and cell growth [37][38][39]. Extracellular vesicles (EVs), biological nanoparticles that mediates cell communication and are present in oviductal fluid, play a crucial role in the communication between gameteoviductal and embryo-oviductal cells [40]. The organs of the female reproductive system exhibit basal cyclooxygenase 2 (COX2) expression [41]. At both the mRNA and protein levels, COX2 expression positively correlated with estrogen receptor (ER) expression [42]. According to early reports, female mice lacking COX2 exhibit decreased ovulation and fertilized eggs [43]. Additionally, COX2 modulates the generation of reactive oxygen species [44] and inflammation and inhibits apoptosis [45]. A-kinase anchoring proteins (AKAPs) are a family of proteins with signaling platform functions; A-kinase anchoring protein 13 (AKAP13) is a member of the AKAP protein family that contains the C-terminal nuclear receptor interacting domain (NRID), which can bind estrogen receptors (ERs) and progesterone receptors (PRs) [46] and increase their ligand-dependent activity of those. Therefore, AKAP13 is usually considered one of the candidate genes for regulating female fertility. Intraflagellar transporter protein 140 (IFT140) is a member of the intraflagellar transport (IFT) protein family and is a conserved mechanism essential for the assembly and maintenance of most eukaryotic cilia and flagella. Ciliated and secretory cells are the two main cell types that constitute the oviductal epithelium [47], whereas embryo transport is realized by the complex interaction of smooth muscle contraction, cilia beat, and fluid flow in the oviduct [47,48]. Therefore, IFT140 might display a key role in maintaining the activity of the oviduct cilia. Interestingly, the significant enrichment in N-Glycan biosynthesis confirms the possibility of cell adhesion, migration, and other physiological processes [49], these are essential for maternal-embryo communication and the proper transport of embryos in the oviduct function [50]. Together, these results reflect the multifaceted control of prolificacy trait of Yunshang black goats' oviducts, which are necessary to control the fertility of female mammals and meri further study.
In this study, we quantified 98 and 167 DAPs in the oviduct tissues derived from follicular and luteal phases of the estrus cycle, respectively, which represents a major contribution to the understanding of oviduct function in the Capra hircus species. For DAPs in the follicular and luteal phases, we focused on biological adhesion and transporter activity terms. The process of biological adhesion is indispensable for intercellular communication, signal transduction, proliferation, and apoptosis [51] that occur at follicular and luteal phases in the oviduct. As the site of the origin of mammalian life, the oviducts are involved in the first reproductive activities, many of which involve the formation of tight junctions between gametes, embryos, and reproductive organs highly dependent on the ability to adhere [50]. For example, oocytes form tight junctions with cumulus cells (CCs) in cumulus-oocyte complexes (COCs), allowing the flow of stimulating factors required for oocyte maturation [52]. After the sperm enters the oviduct, countercurrent movement occurs, and some sperm migrate directly to the ampulla and combine with the oocytes to complete fertilization, while most of them adhere to the oviductal epithelial cells (OECs), waiting for capacitation during ovulation [53]. At the same time, this process is crucial for successful fertilization, which can ensure that a few specific sperm reach the ampulla isthmus junction where fertilization takes place [54]. Obviously, this is not specific to one function of the oviduct at the follicular or luteal phase. The effect of the transport process on reproduction has been documented in the literature. Lopez-Ubeda et al. recently found that the main pathways affected by insemination were related to molecular transport, protein trafficking, and cell-to-cell signaling [55]. Thus, proteins involved in biological adhesion and transporter activity may play roles in reproduction and are associated with the kidding number of goats.
During the follicular phase, compared with the FL group, proteins CAMSAP3, GLUL, SYVN1, ITGAM, PEX19, PSTPIP2, and ACSS3 were significantly upregulated in the FH group. At the same time, two of these proteins, CAMSAP3 [56] and ITGAM [57], participated in the ciliary transport function. As described above, the ciliary function is essential: through it, the oviduct completes gamete and embryo transport, and it is directly involved in reproductive regulation. Synoviolin 1 (SYVN1) was shown to promote the migration, invasion, and proliferation of liver cancer cells and always regulates endoplasmic reticulum stress, blood vessel growth, oxidative stress, and cell apoptosis [58], and therefore may also promote gamete maturation and migration, thereby negatively regulating goat fecundity. In addition, some reproduction-related proteins showed downregulation trends, such as EMG1, ND5, and RING1. Essential for mitotic growth 1 (EMG1) is a highly conserved nucleolar protein identified to be widely expressed during early mouse embryonic development, and the loss of EMG1 function in mice prevents embryonic development before the blastoderm stage [59]. ND5 is the mitochondrial-encoded respiratory chain I (NADH dehydrogenase) subunit 5, the deletion of which leads to respiratory chain disorders. Further, mitochondrial dysfunction may be responsible for female infertility because it affects not only sperm and oocyte quality but also the fertilization process and early embryonic development [60]. Ring finger protein 1 (RING1) has been suggested to play roles in regulating germ-cell-specific gene expression both in the male [61] and female germline [62], which occupy the central locations in the PPI network and act as hubs that interact with other DAPs.
Notably, oocytes enter the oviduct after ovulation, and in the luteal phase, changes in the expression of genes that play a role in cell growth and development occur in the oviduct, thus preparing for early embryo development and successful pregnancy [63]. Hydrogen sulfide (H 2 S) is a potent vasodilator and angiogenic factor with potent angiogenic activity, and cystathionine β-synthase (CBS) upregulation increases endogenous H 2 S production, which plays a role in uterine vasodilation and is involved in pregnancy [64,65]. The trophic role of the oviduct is widely supported during early mitotic active embryonic development [66], where the vasculature acts as a nutrient delivery system and therefore provides nutrients to the embryo and gametes in the oviduct. In mammalian cells, PES1 is associated with more than 50% of telomerase activity, and data suggest that PES1 is a key cofactor in telomerase assembly and can regulate telomerase activity and positively regulate cellular senescence and accretion [67]. In addition, we focused attention on the proteins in the luteal phase that were associated with embryonic development. ELP3 and SEC24C, two downregulated proteins have been previously described to have a direct embryo support function [68,69] and may be important mediators of maternal-embryonic communication for embryo development. Of note, among downregulated proteins, two proteins, SPP1 and HSPA8, have a known role in reproductive processes. Secreted phosphoprotein 1 (SPP1) usually binds to integrins to mediate cell-cell and cell-extracellular matrix communication and promote cell adhesion, migration, and differentiation; it is expressed in the mouse uterus and placenta, and it has been reported to provide nutrients for embryonic or fetal development by increasing angiogenesis within the decidua [70]. The conserved nature of heat shock 70 kDa protein 8 (HSPA8) between mammalian species suggests that this protein may represent a common biological mechanism for maintaining sperm survival in the oviduct [71]. These observations list the possible roles of the oviduct protein within the fecundity and merit further study.
In recent years, proteomics and transcriptomics methods have identified abundant proteins and genes, laying the groundwork for more precise and detailed descriptions of molecular processes and the elucidation of complex physiological processes and their genetic regulation [72,73]. Previous work has started uncovering the amounts of protein and mRNA expression in the same sample to better understand probable molecular mechanisms and identify reliable biomarkers of specific biological processes [24,74]. In this study, we utilized transcriptomic data generated in the oviduct from the same tissue sample, which we previously published [75,76], to complement the proteomic analysis. Comparing the two phases, our transcriptomics and proteomics results showed that many more transcripts can be identified than proteins in the oviduct. This situation is commonly reported, revealing the inconsistent or lagging relationship between transcription and translation levels [77]. In addition, a correlation analysis showed a poor negative correlation between the proteome and transcriptome, indicating an inconsistency between genes at the transcriptional level and protein abundance. The reason for the above can be explained by the fact that in the process of protein translation from mRNA, factors such as mRNA splicing, protein overturning, and protein degradation may lead to a poor correlation between mRNA expression patterns and their respective proteins [78]. Our results indicated that the mechanisms affecting fecundity in the luteal phase were more complex than in the follicular phase, as indicated by the differentially expressed genes and differentially abundant proteins in the oviduct. Meanwhile, we examined the functional annotations of the DEGs/DAPs, and the result indicated that some DEGs and DAPs shared a common GO annotation pattern. Metabolic pathways are extensively enriched during the follicular phase in preparation for the maintenance of a stable oviductal microenvironment and the reception of gametes. In the luteal phase, the regulation of MAP kinase activity, estrogen receptor binding, and angiotensin receptor binding were mainly enriched for DAPs and DEGs. MAP kinase is extensively involved in reproductive processes [79]; we suggested that the MAP kinase activity in the oviduct is one of the main reasons for the prolificacy trait of goats. Moreover, the estrogen receptor (ER) binding term might be one of the main functions of the oviduct during embryo transport [80]. The functional classification of the transcriptome and proteome can improve our understanding of the molecular physiology of the mammalian oviduct.

Ethics Statement
In this research, animal samples were approved by the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences, Beijing, China. Ethical oversight was given by the Animal Ethics Committee of the IAS-CAAS (No. IAS2019-63).

Animals and Sample Acquisition
We obtained tissue specimens of Yunshang black goats, known as a new breed of domestic goat in China with high fecundity trait, which were raised in Yixingheng Animal Husbandry Technology Co., Ltd. Tuanjie Township Base in Kunming City (Yunnan, China). Twenty healthy female Yunshang black goats (ages ranged between 2 and 3 years old) with an average weight of 52.22 ± 0.43 kg were selected and divided into high-fecundity group (H: n = 10, kidding number = 3.3 ± 0.11), and low-fecundity group (L: n = 10, kidding number = 1.75 ± 0.08), according to their kidding number records. All Yunshang black goats were treated with simultaneous estrus based on the progesterone vaginal suppository during the nonpregnant period. Briefly, a controlled internal drug release (CIDR, Inter Ag Co., Ltd. Hamilton, New Zealand) device (progesterone 300 mg) was inserted into the vagina for 16 days to synchronize estrus in April 2020, then the vaginal sponge was removed, and the removal time was set to 0 h. Then, ten goats (high fecundity, n = 5; FL: low fecundity, n = 5) were euthanized (intravenous pentobarbital 100 mg/kg) within 48 h and were considered to be in the follicular phase (FH and FL) and another ten goats (high fecundity, n = 5; low fecundity, n = 5) were euthanized (intravenous pentobarbital 100 mg/kg) within 168 h, which corresponds to the luteal phase (LH and LL). All fresh oviduct tissues were immediately collected in the follicular and luteal phases, respectively, and were washed by ice-cold phosphate-buffered saline (PBS) to remove contaminates, rapidly frozen with liquid nitrogen, and stored in a −80 • C refrigerator until protein extractions.

Protein Preparation and Fractionation for Data-Dependent Acquisition (DDA) Library Generation
Twenty oviduct tissues were first homogenized with a MP FastPrep-24 homogenizer (24 × 2, 6.0 M/S, 60 s, twice), and then SDT buffer (4% SDS, 100 mM DTT, 150 mM Tris-HCl pH 8.0) was added. Further, the lysates were sonicated and boiled for 15 min. After centrifuging at 14,000× g for 40 min, the supernatant was quantified with the BCA Protein Assay Kit (Bio-Rad, Hercules, CA, USA), according to the manufacturer's protocol. Then, the protein samples were stored at −80 • C. An equal aliquot of the 20 samples in this experiment were pooled into one sample for generating a DDA spectral library and quality control. Protein digestion was performed according to the FASP procedure described by Wisniewski [81]. Briefly, an equal quantity of proteins in each sample was incorporated into SDT buffer (4% SDS, 100 mM DTT, 150 mM Tris-HCl pH 8.0). The detergent, DTT, and other low-molecular-weight components were removed using UA buffer (8 M Urea, 150 mM Tris-HCl pH 8.0) by repeated ultrafiltration (Microcon units, 10 kD). Then 100 µL iodoacetamide (100 mM IAA in UA buffer) was added to block reduced cysteine residues and the samples were incubated for 30 min in a dark environment. The filters were washed with 100 µL UA buffer three times and then 100 µL 25 mM NH 4 HCO 3 buffer twice. Finally, the protein suspensions were digested with 4 µg trypsin (Promega, Madison, WI, USA) in 40 µL 25 mM NH 4 HCO 3 buffer overnight at 4 • C, and the resulting peptides were collected as a filtrate. The peptides of each sample were desalted on C18 Cartridges (Empore TM SPE Cartridges C18 (standard density), bed I.D. 7 mm, volume 3 mL, Sigma, St. Louis, MI, USA), concentrated by vacuum centrifugation and reconstituted in 40 µL of 0.1% (v/v) formic acid. The peptide content was estimated by UV light spectral density at 280 nm using an extinction coefficient of 1.1 of 0.1% (g/L) solution that was calculated based on the frequency of tryptophan and tyrosine in vertebrate proteins.
Digested pool peptides were then fractionated to 10 fractions using Thermo Scientific TM Pierce TM High pH Reversed-Phase Peptide Fractionation Kit [82]. Each fraction was concentrated by vacuum centrifugation and reconstituted in 15µL of 0.1% (v/v) formic acid. Collected peptides were desalted on C18 Cartridges (Empore TM SPE Cartridges C18 (standard density), bed I.D. 7 mm, volume 3 mL, Sigma) and reconstituted in 40µL of 0.1% (v/v) formic acid. The iRT-Kits (Biognosys, Schlieren, Switzerland) were added to correct the relative retention time differences between runs with volume proportion 1:3 for iRT standard peptides versus sample peptides.

Data-Dependent Acquisition (DDA) Mass Spectrometry Assay
All fractions for DDA library generation were injected on a Thermo Scientific Q Exactive HF X mass spectrometer connected to an Easy nLC 1200 chromatography system (Thermo Scientific). The peptide was first loaded onto an EASY-Spray TM C18 Trap column (Thermo Scientific, P/N 164946, 3 µm, 75 µm × 2 cm), then separated on an EASY-SprayTM C18 LC Analytical Column (Thermo Scientific, ES802, 2 µm, 75 µm × 25 cm) with a linear gradient of buffer B (84% acetonitrile and 0.1% Formic acid) at a flow rate of 250 nl/min over 90 min. MS detection method was positive ion, the scan range was 300-1800 m/z, resolution for MS1 scan was 60,000 at 200 m/z, target of AGC (Automatic gain control) was 3 e6, maximum IT was 25 ms, dynamic exclusion was 30.0 s. Each full MS-SIM scan followed 20 ddMS2 scans. Resolution for MS2 scan was 15,000, AGC target was 5 e4, maximum IT was 25 ms and normalized collision energy was 30 eV.

Data-Independent Acquisition (DIA) Mass Spectrometry Assay
The peptides of each sample were analyzed by LC-MS/MS operating in the dataindependent acquisition (DIA) mode by Shanghai Applied Protein Technology Co., Ltd. Each DIA cycle contained one full MS-SIM scan, and 30 DIA scans covered a mass range of 350-1800 m/z with the following settings: SIM full scan resolution was 120,000 at 200 m/z; AGC 3e6; maximum IT 50 ms; profile mode; DIA scans were set at a resolution of 15,000; AGC target 3e6; Max IT auto; normalized collision energy was 30 eV. Runtime was 90 min with a linear gradient of buffer B (80% acetonitrile and 0.1% Formic acid) at a flow rate of 250 nl/min. QC samples (pooled samples from equal aliquots of each sample in the experiment) were injected in the DIA mode at the beginning of the MS study and after every 6 injections throughout the experiment, which was used to monitor the MS performance.

Raw Mass Spectrometry Data Analysis
For DDA library data, the FASTA sequence database was searched with Spectronaut software (Spectronaut Pulsar X TM _12.0.20491.4, Biognosys). The database was downloaded at website: http://www.uniprot.org.iRT, accessed on 20 October 2021. The following peptides sequence was also added (>Biognosys|iRT-Kit| Sequence_fusionLGGNEQVTRYILA-GVENSKGTFIIDPGGVIRGTFIIDPAAVIRGAGSSEPVTGLDAKTPVISGGPYEYRVEATF-GVDESNAKTPVITGAPYEYRDGLDAASYYAPVRADVTPADFSEWSKLFLQFGAQGSPFLK). The parameters were set as follows: enzyme was trypsin, max missed cleavages were 2, the fixed modification was carbamidomethyl (C), and the dynamic modification was oxidation (M) and acetyl (Protein N-term). All reported data were based on 99% confidence for protein identification as determined by the false discovery ratrunsee (FDR = N(decoy) × 2/(N(decoy) + N(target)) ≤ 1%. The spectral library was constructed by importing the original raw files and DDA search results into Spectronaut Pulsar X TM_12.0.20491.4 (Biognosys).
DIA data were analyzed with Spectronaut TM 14.4.200727.47784, searching the above constructed spectral library. Main software parameters were set as follows: retention time prediction type was dynamic iRT, interference on MS2 level correction was enabled, and cross-run normalization was enabled. All results were filtered based on a Q value cutoff of 0.01 (equivalent to FDR < 1%). After normalized to total peak intensity, the processed data were uploaded before importing into SIMCA-P (version 14.1, Umetrics, Umea, Sweden), which was subjected to multivariate data analysis. We performed the Pareto-scaled principal component analysis (PCA) of quantified proteins to visualize the relationship between the oviduct samples obtained from different fecundity goats using R packages.

Weighted Gene Co-expression Network Analysis
We applied the weighted gene co-expression network analysis (WGCNA) method in R package (Version 1.69) [83] to identify distinct protein modules among all identified proteins, with the following parameters: corType = pearson; minModuleSize = 30; merge-CutHeight = 0.1. A weighted protein co-expression network was generated using the log 2 protein abundance sample matrix. The R software package gplots was used to perform a hierarchical clustering analysis and identify the sample outliers. After performing the necessary soft threshold power calculations, the standard scale-free network was set up. The closer that correlation value was to 1 and the correlation test p value of <0.05, meaning the module was the important one in determining the trait, the more relevant the models or networks created using WGCNA were to the external sample features. Next, tissue-specific indicators were found using correlation functional networks. Proteins were grouped hierarchically to create a dynamic tree-cut technique to identify the modules. The protein information in each module was then extracted using WGCNA, which was utilized to build the protein co-expression modules. The sequences of the proteins were locally searched using the NCBI BLAST+ client software (NCBI-blast-2.2.28+-win32.exe) and InterProScan to find homolog sequences, then GO annotation, KEGG enrichment, and co-expressed protein network analyses were performed to analyze, identify, and interpret diverse biological functions and hub proteins in the functionally related modules.

Data Analysis and Identification of Differentially Abundant Proteins
Differential proteome analysis between high-and low-fecundity goats' oviducts of each phase were conducted. Proteins with a fold change (FC) >1.5 or <0.67 and p value of <0.05 in comparable groups (FH vs. FL or LH vs. LL) were considered differentially abundant proteins (DAPs). Volcano plots were generated using -log 10 (p value) as the y-axis and ±log 2 (|FC|) as the x-axis. Cluster software (R version 3.6.3) and Java Treeview software were used to perform hierarchical clustering analysis. Euclidean distance algorithm for similarity measure and average linkage clustering algorithm (clustering uses the centroids of the observations) for clustering were selected when performing hierarchical clustering.

Functional Annotation of Differentially Abundant Proteins
The functions of differentially abundant proteins were annotated to further explore their impact on biological processes. CELLO (http://cello.life.nctu.edu.tw/, accessed on 11 April 2022), which is a multi-class SVM classification system, was used to predict protein subcellular localization. Protein sequences were searched using the InterProScan software [82] to identify protein domain signatures from the InterPro member database Pfam. Domain enrichment analysis of differentially expressed proteins was performed using Fisher (Fisher's Exact Test), considering the whole quantified proteins as the background dataset. Benjamini-Hochberg correction for multiple testing was applied to adjust derived p values.
Gene ontology (GO) terms were mapped, and sequences were annotated using the software program Blast2GO [84] (https://www.blast2go.com/, accessed on 5 May 2022). The GO annotation results were plotted by R scripts. Next, following annotation steps, the studied proteins were blasted against the online Kyoto encyclopedia of genes and genomes (KEGG) database (http://geneontology.org/, accessed on 5 May 2022) to retrieve their KEGG identifications and were subsequently mapped to pathways in KEGG. The protein-protein interaction (PPI) information of the DAPs was retrieved from the STRING online tool (http://string-db.org/, accessed on 10 May 2022, confidence score ≥ 0.4 and classified through k-means). Furthermore, the degree of each protein was calculated to evaluate the importance of the protein in the PPI network.

Transcriptome Sequencing
High-throughput RNA sequencing (RNA-seq) was conducted following our previously published studies [75,76] and using the same oviduct samples as in protein extraction. Briefly, total RNAs from the 20 oviduct tissues in follicular and luteal phases were extracted by TRIzol reagent (Invitrogen, Carlsbad, CA, USA), respectively. The RNA concentration and purity were measured using a NanoDropTM 2000 (Thermo Scientific TM , Wilmington, DE, USA) instrument, and RNA integrity was evaluated using an Agilent 2100 System (Agilent Technologies, Santa Clara, CA, USA). After quality assessment, these RNA samples were sent to Wuhan Frasergen Bioinformatics Co., Ltd. (Wuhan, China) for cDNA library construction and then sequenced on the Illumina NovaSeq platform. Immediately, HISAT2 (v2.1.0) [85] was used to map the clean reads of each sample to the reference genome. Only the uniquely mapped reads were assembled, and the expression levels were predicted using String Tie software (v.1.3.5) [86]. Next, the FPKM (fragments per kilobase of exon model per million mapped fragments) for each gene was obtained. DESeq2 [87] was used to calculate fold change and p value based on the normalized counts. The threshold of adjusted p value < 0.05 and fold change >1.5 or <0.67 were used to screen significantly differentially expressed genes (DEGs).

Correlation Analyses of Transcriptome and Proteome Profiles
We performed an in-depth analysis of 20 oviduct tissues from goats to understand the oviduct function and facilitate the discovery of prolificacy-related biomarkers, and the sketch map of correlation analyses of the transcriptome and proteome is shown in Supplementary Figure S1. A joint analysis of the proteome and transcriptome will help us grasp gene expression regulation in general. If an mRNA and its corresponding protein were expressed at the same phase, they were considered correlated in the high-and low-fecundity comparison groups, and if the expression levels of both an mRNA and its corresponding protein are significantly different at the same phase, they were defined as differentially expressed associated transcripts (DEG; threshold of fold change >1.5 or <0.67 and p < 0.05) and proteins (DAP; threshold of fold change > 1.5 or <0.67 and p < 0.05), respectively. The expression matrix of all genes/proteins was converted into log 2 format. Two omics databases for the comparison groups of FH vs. FL and LH vs. LL were calculated by Spearman correlation coefficients, respectively. In addition, if the expression levels of both a gene and its corresponding protein were significantly different at the same phase, they were defined as differentially expressed associated transcripts and proteins, respectively. GO annotation and KEGG pathway analysis were then conducted for DEGs and DAPs.

Validation of Proteomic Data with Parallel Reaction Monitoring (PRM)
Six proteins' expression levels selected from DIA quantitative proteomics analysis were validated using LC-PRMMS at Shanghai Applied Protein Technology Co., Ltd. (Shanghai, China). Peptides were prepared according to the DIA protocol, and the internal standard reference was a Peptide Retention Time Calibration Mixture (thermo) stable isotope peptide spiked in each sample. Prior to reversed-phase chromatography on an Easy nLC-1200 system, tryptic peptides were loaded on C18 stage tips for desalting (Thermo Scientific). One-hour liquid chromatography gradients with acetonitrile ranging from 5 to 35% in 45 min were used. A PRM analysis was carried out using a Q Exactive HF mass spectrometer (Thermo Scientific). Experimentally, using unique peptides of high intensity and confidence for each target protein, methods optimized for collision energy, charge state, and retention times for the most significantly regulated peptides were generated. The mass spectrometer was set to positive ion mode with the following parameters: a resolution of 70,000 (at 200 m/z), automatic gain control (ACG) target values of 3.0 × 10 −6 , and a maximum ion injection time of 200 ms. Following full MS scans, 20 PRM scans were performed at 35,000 resolution (at m/z 200), with AGC 3.0 × 10 −6 and a maximum injection time of 200 ms. A 1.5 Th window was used to isolate the targeted peptides. Ion activation/dissociation was carried out in a higher energy dissociation (HCD) collision cell at a normalized collision energy of 27. Skyline (MacCoss Lab, University of Washington) [88] was used to analyze the raw data, and signal intensities for individual peptide sequences for each of the significantly altered proteins were quantified relative to each sample and normalized to a standard reference.

Measuring mRNA Levels Using Real-Time Quantitative PCR
Four co-expressed DE mRNAs were randomly screened from per group and quantified. The primers of mRNAs were designed by Primer Premier 6 software and RPL19 was used as the reference gene, primer information is shown in Table 1. All primers were synthesized by Sangon Biotech (Shanghai, China). Three replicates of each sample were analyzed, and standard curves were established using the Roche Light Cycler ® 480II system (Roche Applied Science, Mannheim, Germany). Relative expression level of mRNA was calculated using the comparative cycle threshold (2 -∆∆Ct ) method [89] dependent on the t-test.

Statistical Analysis
The normality of experiment data was investigated with the Shapiro-Wilk normality test prior to all statistical analyses. According to the results of the normality investigation, the difference in the expression level in the FH vs. FL and LH vs. LL comparison groups was analyzed by a Student's t-test. Statistical analyses were calculated using SPSS v.25.0, and GraphPad Prism v.9.3.1 was performed to visualize the data. Spearman rank correlation analysis was performed to assess the correlation between proteins and mRNAs' abundance in the oviduct. Each experiment was performed with at least 3 independent biological replicates; all data were expressed as means ± standard error of the mean (SEM). P values of < 0.05 indicate a statistical significance and are presented as * p < 0.05, ** p < 0.01, and *** p < 0.001.

Conclusions
Our analysis provides comprehensive insight into the underlying mechanisms by which oviduct function causes differences in the fecundity of goats. The WGCNA revealed that the oviduct function during the luteal phase might play an important role in the prolificacy trait. Meanwhile, the stage-specific changes of proteins were identified, which expanded our understanding of the oviductal maternal-gamete communication factors in regulating kidding numbers. Additionally, we utilized the linkage between the proteome and transcriptome to suggest the main process for oviductal executive functions. The genes at both levels are closely associated with the same stage-specific functions, such as metabolic processes, MAP kinase activity, estrogen receptor binding, and angiotensin receptor binding, highlighting the importance of combining different gene-expression measurements. In future works, candidate biomarkers need to be added or overexpressed in the in vitro co-culture of oviductal cells and gametes/embryos to validate the function of these genes.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ijms232314888/s1, Figure S1: Experimental design and work flow for the integrative proteomics and transcriptomics comparison of high-and low-fecundity goat oviducts from Yunushang black goat in the follicular and luteal phases; Figure S2: Quality control (QC) and quantitative heat map of DIA; Figure S3: Weighted correlation network analysis (WGCNA); Figure S4: Combined Analysis of Proteome and Transcriptome; Table S1: Peptide quality statistics; Table DataS1: The identification and quantitative results of DIA; Table DataS2: WGCNA