Integrated analyses of miRNA-mRNA expression profiles of ovaries reveal the crucial interaction networks that regulate the prolificacy of goats in the follicular phase

Background Litter size is an important index of mammalian prolificacy and is determined by the ovulation rate. The ovary is a crucial organ for mammalian reproduction and is associated with follicular development, maturation and ovulation. However, prolificacy is influenced by multiple factors, and its molecular regulation in the follicular phase remains unclear. Methods Ten female goats with no significant differences in age and weight were randomly selected and divided into either the high-yielding group (n = 5, HF) or the low-yielding group (n = 5, LF). Ovarian tissues were collected from goats in the follicular phase and used to construct mRNA and miRNA sequencing libraries to analyze transcriptomic variation between high- and low-yield Yunshang black goats. Furthermore, integrated analysis of the differentially expressed (DE) miRNA-mRNA pairs was performed based on their correlation. The STRING database was used to construct a PPI network of the DEGs. RT–qPCR was used to validate the results of the predicted miRNA-mRNA pairs. Luciferase analysis and CCK-8 assay were used to detect the function of the miRNA-mRNA pairs and the proliferation of goat granulosa cells (GCs). Results A total of 43,779 known transcripts, 23,067 novel transcripts, 424 known miRNAs and 656 novel miRNAs were identified by RNA-seq in the ovaries from both groups. Through correlation analysis of the miRNA and mRNA expression profiles, 263 negatively correlated miRNA-mRNA pairs were identified in the LF vs. HF comparison. Annotation analysis of the DE miRNA-mRNA pairs identified targets related to biological processes such as “estrogen receptor binding (GO:0030331)”, “oogenesis (GO:0048477)”, “ovulation cycle process (GO:0022602)” and “ovarian follicle development (GO:0001541)”. Subsequently, five KEGG pathways (oocyte meiosis, progesterone-mediated oocyte maturation, GnRH signaling pathway, Notch signaling pathway and TGF-β signaling pathway) were identified in the interaction network related to follicular development, and a PPI network was also constructed. In the network, we found that CDK12, FAM91A1, PGS1, SERTM1, SPAG5, SYNE1, TMEM14A, WNT4, and CAMK2G were the key nodes, all of which were targets of the DE miRNAs. The PPI analysis showed that there was a clear interaction among the CAMK2G, SERTM1, TMEM14A, CDK12, SYNE1 and WNT4 genes. In addition, dual luciferase reporter and CCK-8 assays confirmed that miR-1271-3p suppressed the proliferation of GCs by inhibiting the expression of TXLNA. Conclusions These results increase the understanding of the molecular mechanisms underlying goat prolificacy. These results also provide a basis for studying interactions between genes and miRNAs, as well as the functions of the pathways in ovarian tissues involved in goat prolificacy in the follicular phase. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-021-08156-2.


Background
In mammals, litter size is a crucial economic trait because it has a notable impact on profitability in the breeding industry. Ovulation rate and embryo survival are directly linked to mammalian litter size [1]. Most goat breeds produce litter sizes of two kids, which limits production in the goat industry, including the production of meat, milk and fur. Ovaries directly mediate ovulation, which has a significant impact on the fecundity of mammals. Traditional breeding methods require too much time to improve litter sizes and have shown low heritability. Many molecular technologies have been used to analyze the mechanisms underlying reproduction in goats. Molecular mechanisms occurring in the ovarian tissues in the follicular phase may contribute to observing differences in litter size [2]. Many genes, transcription factors, ncRNAs and signaling pathways regulate the complex and precise process of reproduction [3]. These factors have been used to construct a network of many interactions in vivo that regulate the ovulation rate of goats.
Interactions between many genes are involved in ovulation rate and litter size. Several studies have shown that bone morphogenetic protein receptor type 1B (BMPR1B), also known as FecB, plays a major role in sheep fertility. Q249R mutations within FecB are highly associated with ovulation rates in many sheep breeds around the world [4]. Another crucial gene, growth differentiation Factor 9 (GDF9), also plays a critical role in early folliculogenesis [5]. In addition, growth hormones and members of the insulin-like growth factor (IGF) system (IGF-I and IGF-II) may regulate follicular development and atresia [6,7], and FER1L4 and SRD5A2 were found to be associated with high fecundity in goats [8]. Second, many interactions between miRNAs and mRNAs influence biological processes in the goat ovary. miRNAs are endogenous noncoding RNAs that regulate gene expression at the posttranscriptional level by binding to the 3′-UTR of target mRNAs. Multiple miRNAs act on folliculogenesis by regulating transcription factors and genes associated with it. For example, miR-26b promoted apoptosis in porcine granulosa cells (GCs) during follicular atresia by targeting the ataxia-temporal (ATM) gene and increasing the number of DNA breaks [9]. These interactions play important regulatory roles in the regulation of biological processes involved in ovulation rates and fertility. For example, miR-133b represses the expression of its target genes forkhead Box L2 (FOXL2) and transgelin 2 (TAGLN2) in GCs, thereby regulating oocyte growth and maturation [10].
However, the regulation of ovulation rates in goats is a highly complex process involving the interaction of multiple miRNAs and mRNAs within intricate signaling pathways rather than a single network of interactions. miR-141 affects the MAPK signaling pathway by regulating DAPK1 expression, which leads to the development of polycystic ovary syndrome (PCOS) [11]. Moreover, both miR-206 and miR-1 expression were found to be downregulated in the large follicles of a variety of goats, which might be associated with oocyte maturation and ovulation [12]. In addition, mitogen-activated protein kinase kinase kinase 3 (MAP 3 K3) interacts with extracellular signal-regulated kinase (ERK) in the p38 MAPK pathway to regulate litter size traits in Berkshire pigs [13]. miR-122 has been shown to increase LHR mRNA levels by modulating LRBP expression through SREBP activation, a phenomenon that is essential for regulating key reproductive processes, such as ovulation and CL function [14]. These studies suggest that reproductive traits are regulated by multiple interaction networks of genes, miRNAs and signaling pathways. It is therefore important to identify the network of interactions between genes, miRNAs and signaling pathways in the regulation of ovulation rate and prolificacy. There is currently a lack of understanding of these networks in goats.
Goats are important agricultural animals, and as with other animals, genetic factors are the main players in ovulation and litter size. For example, goats have a kidding rate heritability of 0.11-0.14, and the ovulation number of the estrus cycle is 2-6 [15]. Therefore, elucidating the mechanisms by which genes, miRNAs and pathways interact in goat reproduction, particularly in relation to the ovarian ovulation rate, is of great significance for improving genetics and productive performance in goats. Several reports have described gene regulation, miRNA identification and pathway regulation in goat prolificacy. For instance, for litter size, the expression of miRNAs associated with ovulation, such as miR-122, miR-99a, miR-141, miR-493 and miR-200a, and the expression of mRNAs, such as CD36, TNFAIP6, CYP11A1, SERPINA5 and PTGFR, have been identified and comparatively analyzed in the ovaries of prolific and nonprolific goats [16][17][18][19]. Moreover, CYP11A1 plays a key role in the regulation of steroid-producing pathways in granulosa cells [20]. However, these studies have been limited to the role of genes, miRNAs or pathways alone in goat prolificacy, and there is a lack of understanding of the interactions of these factors in goat high yield traits. Furthermore, early studies have focused on goat estrus and entire reproductive cycle periods, but an understanding of the changes in the transcriptomes in goat ovaries in the follicular phase remains lacking.
The Yunshang black goat is an excellent breed due to its fecundity and meat production, and it is native to Yunnan Province, China. This breed has many favorable characteristics, especially its tender, tasty meat and excellent fecundity. Yunshang black goats are suitable for cultivation and reproduction. Many studies have been conducted on the molecular mechanisms related to resource conservation, strain selection, industrial production and trait development of this local breed. Therefore, Yunshang black goats might serve as a representative model to study mammalian prolificacy regulation. To reveal the crucial molecular mechanisms and molecular networks underlying prolificacy in this breed, in the present study, we identified the miRNA and mRNA expression profiles in follicular phase ovaries from goats with large and small litter sizes. Important miRNA-mRNA negative correlation interaction networks and pathways associated with prolific traits were identified. Our aim was to help understand the molecular mechanisms involved in the formation of prolific traits in the follicular phase ovaries of goats and to provide a basis for further investigation of the interactions among miR-NAs, mRNAs and signaling pathways associated with prolificacy. Therefore, the Yunshang black goat might provide a typical model for studying prolific traits in mammals. Our study also provides basic data for breed trait mining and resource conservation.

Results
cDNA library sequencing and transcriptome profiles of the ovary RNA-seq was performed on 10 cDNA libraries using the Illumina Novaseq platform. The AMPure XP system (Beckman Coulter, Beverly, MA, United States) was used to purify library fragments, preferentially selecting cDNA fragments of 150 bp in length. Greater than 50 million raw reads were generated from each library. After filtering out the low-quality reads, the clean reads ranged from 15.1 to 24.6 G, with more than 92.9% of the clean reads mapped to the genome (additional files, Table S1). Although the GC content of the clean data ranged from 45.8 to 53.3%, the quality scores of the clean reads were higher than 98.1 and 92.3% for Q20 and Q30, respectively, indicating that the reliability and quality of the sequencing data were adequate for further analysis. Finally, a total of 43,779 known transcripts and 23,067 novel transcripts were identified in the data (shown in additional files, Table S2). To investigate the key mRNAs involved in the formation of reproductive traits in goats, we analyzed the differentially expressed genes (DEGs) in the ovaries of the LF and HF groups of Yunshang black goats. The number of up-and downregulated genes between the HF and LF groups is shown in Fig. 1A, and a list of the DEGs identified in the comparison is given in additional file Table S3. The log2 (fold change) values ranged from − 30 to 27.62186.
Sequencing of small RNA library and miRNA transcriptome profiling analysis in the ovary The ovarian tissues collected from goats in the LF and HF comparison were also used to construct small RNA libraries. A total of 1.084 to 1.416 G clean reads were generated after filtering the low-quality reads, and more than 93.95% of these reads were mapped to the genome (additional files, Table S4). The Q20 and Q30 quality scores were greater than 99.12 and 96.24%, respectively, and the GC content of the clean data ranged from 48.44 to 49.93%. These results showed that the samples were reliable and of high quality. Finally, 424 known miRNAs and 656 novel miRNAs were identified (additional files, Table S5). To investigate the key miRNAs involved in prolificacy, we analyzed the DEMs in the comparison. The number of up-and downregulated miRNAs in the HF vs. LF comparison is shown in Fig. 1B.

Correlation analysis of differentially expressed miRNAs and mRNAs involved in prolificacy
We integrated the miRNA and mRNA data from Yunshang black goat ovaries, predicted the potential target genes of the DEMs, completed the correlation analysis of the DEM-mRNA pairs, and screened DEM-mRNA pairs with negatively correlated expression. A total of 263 miRNA-mRNA pairs were predicted in the LF vs. HF comparison (Table 1, Pages 23-25). Among these pairs, 115 miRNA-mRNA pairs were found to have negatively correlated expression (q value < 0.1). To gain insights into the functions of the miRNA-mRNA pairs with negatively correlated expression between the groups with different litter sizes, we performed GO enrichment analysis to reveal the enriched biological process terms (corrected P values < 0.05). The associations were mainly with cellular components and the binding of molecular functions (Fig. 2). The GO terms related to reproductive traits are also shown in Table 2, and these terms include reproductive process, ovulation cycle, estrogen receptor binding and ovarian follicle development.

miRNA-mRNA interaction network related to reproduction
To identify key molecular players in the prolificacy of goat ovaries, we focused on some DEG-DEM interactions associated with the ovaries of goats with different litter sizes based on the annotations of the negatively correlated miRNA-mRNA pairs at the expression level ( Fig. 3 and Table 1, Pages 23-25). For example, novel_ 438 was a core gene in the network (Fig. 3), and a total of 55 miRNA-mRNA pairs were identified whose targets were closely related to biological processes such as hormone secretion, hormone-mediated signaling pathways and protein binding. The core gene of the SLC16A5 network was targeted by 3 miRNAs, including novel_438, novel_535 and chi-miR-197-5p. Moreover, the network revealed chi-miR-324-3p as an important node, and its targets were associated with oocyte development. The results indicated that WNT4 is targeted by 2 miRNAs, namely, chi-miR-324-3p and chi-miR-2290.

miRNA pathways related to reproduction
To further understand how DEMs work with DEGs in pathways regulating goat prolificacy, we performed KEGG pathway analysis of the miRNA-mRNA pairs with correlated expression. The KEGG terms enriched in pathway categories for the comparison are shown in Table 3 (Pages 25-27). The two top pathways were autophagy-animal (ko04140) and Hippo signaling pathways (ko04390). Interestingly, we found that five signaling pathways, including oocyte meiosis (ko04114), progesterone-mediated oocyte maturation (ko04914), the GnRH signaling pathway (ko04912), the Notch signaling pathway (ko04330) and the TGF-beta signaling pathway (ko04350), were directly and undirectedly associated with the reproductive process. Eight genes, including CAMK2G, NCOR2, PKMYT1, PKMYT1, SMURF1, BMP4, WNT4 and TCF7L2, were involved in the five pathways and might be important candidate genes for further study.   PPI network related to prolificacy trait The PPI network was constructed using the extracted target gene list from the STRING database in Cytoscape (Fig. 4). The PPI network of the DEGs in the LF vs. HF comparison contained 98 protein-protein nodes. CDK12, FAM91A1, PGS1, SERTM1, SPAG5, SYNE1, TMEM14A, WNT4, CAMK2G and other genes were the key core nodes of the networks, and all of them were targets of DE miRNAs.

The effects of predicted miRNA-mRNA pairs on GC proliferation
To illustrate the functions of the predicted miRNA-mRNA pairs, two pairs were selected for dual luciferase reporter assays and CCK-8 assays. The dual luciferase reporter assay results showed that miR-1271 could bind to the 3'UTR of the TXLNA gene (Fig. 7). To further assess the functional role of miR-1271, we performed miRNA overexpression and repression assays. Goat GCs were used to optimize the dose of miR-1271 in the overexpression and repression experiments. We found that overexpression of miR-1271 significantly reduced the rapid cell proliferation rate and that inhibiting miR-1271 expression significantly increased cell proliferation, as measured by CCK-8 assay (Fig. 8).

Discussion
Precise regulation of hormones, genes, noncoding RNAs and other factors, including their involvement in spatiotemporal expression, signal cascades, transcriptional regulation and feedback mechanisms, affects the prolificacy trait [21]. In recent years, many studies have been conducted on gene identification and functional research, miRNA identification and pathway regulation for prolificacy in goats. However, these studies have focused on estrus and entire reproductive cycle periods [19,22], and there is still a lack of understanding of the molecular mechanisms associated with prolificacy in the follicular phase of goat ovaries. Therefore, it is important to identify these molecular mechanisms. In this study, RNA-seq was performed to assess the follicular phase ovaries of two groups (high-and low-yielding) of Yunshang black goats, and then miRNA and mRNA profiles related to prolificacy were constructed. In total, 43,779 known transcripts, 23,067 novel transcripts, 424 known miRNAs and 656 novel miRNAs were identified in the goat ovarian samples (additional files, Table S2). These data provide a valuable resource for analyzing the molecular mechanisms of prolificacy and transcriptional regulation in goat ovaries. The Yunshang black goat is a local breed used for meat production in China. With its tender flesh, unique flavor and high fecundity, this breed has important    [23,24]. At present, there is a lack of understanding of the molecular mechanisms underlying the prolificacy of Yunshang black goat ovaries in the follicular phase. Differentially expressed genes and miRNAs may play important regulatory roles in prolificacy. In this study, 206  DEGs and 19 DEMs were identified as being differentially expressed in this comparison (Table 1). Thirteen genes, including FAM91A1, DENND1A, SYNE1, TJP3, ST3GAL3, EGFLAM, SERTM1, WNT4, CDK12, TMEM164, TMEM14A, SPAG5, BMP4 and PGS1, were directly associated with fertility. Moreover, we screened DE miRNA-mRNA pairs in Yunshang black goat ovaries in the follicular phase and then performed functional enrichment analysis. The KEGG pathway enrichment results showed that in the LF vs. HF comparison, the genes were mostly enriched in pathways associated with fertility. These results indicate that the follicular phase is an important stage for the prolificacy of Yunshang black goat ovaries, which is consistent with findings from previous studies in goats, bovines and humans [17,25,26]. Therefore, the miRNA and mRNA expression profiles obtained in this study reflect the molecular regulation of the follicular phase of Yunshang black goat ovaries.
Prolificacy involves many biological processes. To gain insight into how miRNAs and their target genes regulate prolificacy, we constructed interaction networks between miRNAs and mRNAs related to prolificacy (Fig. 3). Among these interaction networks, WNT4 was predicted to be a target of miRNAs, such as chi-miR-2290 and chi-miR-324-3p. This gene is also involved in biological processes including reproductive processes, oocyte development, oogenesis and oocyte differentiation. WNT4 plays a regulatory role in ovarian follicular development and female fertility [27]. miR-324-3p has also been suggested to regulate ovarian cancer [28]. miR-2290 is a novel miRNA that was first reported in virus-infected bovines [29]. In our study, miR-2290 was found to be an important regulator of the WNT4 gene in the Wnt signaling pathway. In addition, among the networks, the most complex miRNA-mRNA interaction network was based on miR-324-3p, novel_438 and novel_535 (Fig. 3). miR-1271 targeted TXLNA, and high TXLNA expression could be a favorable biomarker and therapeutic target for pancreatic adenocarcinoma patients [30]. In our study, TXLNA expression was significantly lower in the high-yielding goats than in the low-yielding goats, and the dual luciferase reporter assay showed that its expression was regulated by miR-1271. The CCK-8 assay showed that overexpression of miR-1271 could inhibit the proliferation of GCs and that inhibition of expression promoted the proliferation of GCs. We deduced that miR-1271-TXLNA was a potential regulatory pathway of goat prolificacy. Furthermore, the network is related to ovarian follicular development and oogenesis. The interaction network consists of eight core genes, including ARHGAP17, ALS2CL, ANO2, ATP2A3, ATXN2L, DERL3, MINK1 and SLC16A5, and 9 miRNAs, including miR-197-5p and miR-133a-5p. Some studies have confirmed that ARHG AP17 is a member of the Rho family that is associated with follicular growth and granulosa cell maturation [31]. SLC16A5 belongs to solute carrier proteins, which are included in the top 100 most highly expressed genes in opossum during pregnancy [32]. MINK1 is a member of the germinal center kinases and is known to regulate cytoskeletal organization and oncogene-induced cell senescence, which is essential for cytokinesis [33]. ATP2A3, a member of the SERCA family, is an intracellular calcium pump and translocates cytosolic Ca 2+ ions to the sarcoplasmic reticulum lumen. Its upregulated status in the uterus suggests that ATP2A3 helps in oviposition [34]. In addition, miR-133a-5p is a target of the circular RNA scinderin (circSCIN), which affects estrogen secretion [35]. These findings suggest that these miRNA-mRNA interaction networks may play key roles in the regulation of prolificacy in follicular phase ovaries in Yunshang black goats, thus providing new insights for further understanding the regulatory mechanisms of goat prolificacy.
Follicle development in the follicular phase is the main factor affecting the ovulation rate of mammals [36,37]. Therefore, the pathways involved in adhesion, growth, migration and differentiation of follicles might be critical for goat reproduction. In our study, the GnRH signaling pathway, oocyte meiosis, TGF-beta signaling pathway, Notch signaling pathway and Wnt signaling pathway were enriched in the follicular phase ovaries of goats with different litter sizes (Table 3, Pages 27-29). Similar results have been reported in previous studies of other species [38,39]. In fact, the GnRH signaling pathway has an integral role in the process of ovulation, such as initiating pituitary secretion of luteinizing hormone (LH) and follicle-stimulating hormone [40]. At the end of the follicular phase, estradiol switches from suppressing release to inducing a sustained surge of release, which induces ovulation [41]. Oocyte meiosis is considered a specialized form of cell division required for the formation of haploid gametes, and the mutation of meiosisspecific cohesion components in female mammals results in an increased frequency of oocyte aneuploidy and premature ovarian failure [42][43][44]. In addition, the reinitiation of oocyte meiosis is known to play a key role in the follicle phase [45], thus implicating an enrichment of pathways related to follicle development and ovulation rate. Therefore, the GnRH signaling pathway, oocyte meiosis, TGF-beta signaling pathway, Notch signaling pathway and Wnt signaling pathway may be key to further exploring the molecular mechanisms underlying the regulation of goat prolificacy.
Moreover, we also demonstrated interactions among genes related to the prolificacy of Yunshang black goats in the follicular phase, such as the PPI network of CAMK2G and PXN (Fig. 4). CAMK2G and PXN are regulatory factors related to the dynamics of focal adhesion assembly and disassembly [46,47]. Moreover, CAMK2G is linked to the oocyte meiosis and GnRH signaling pathways (Table 3), which are involved in ovulation and high yielding [48,49]. Among the PPI networks of SYNE1, CDK12, and PGS1 (Fig. 4), the ovarian and other cancer-related gene SYNE1 is located 19 kb downstream of ESR1, and it partially contributes to ovarian cancer [50]. In our data, these three interacting genes were differentially expressed in the follicular phase ovaries of Yunshang black goats with different litter sizes, and the same expression trend was observed in ovarian carcinoma [51,52]. Moreover, BMP4, SMURF1, TCF7L2 and WNT4 belong to the TGF-beta signaling pathway and GnRH signaling pathway, and we identified an interaction between BMP4 and SMURF1 and between WNT4 and TCF7L2 (Fig. 4). Studies have shown that BMP4 is synthesized by vaginal tissues, and its expression is negatively regulated by estrogen [53,54]. Thus, the interactions among genes in these pathways affect the ovulation rate, and the interactions of miRNAs and mRNAs in these pathways may provide new insights into the regulation of prolificacy in goats. The above results indicate a role of complex intergene interactions in the prolificacy of Yunshang black goat ovaries.

Conclusions
In summary, we analyzed the miRNA and transcriptome profiles in follicular phase ovaries from different litters of Yunshang black goats and constructed an interaction network of proteins, miRNAs, mRNAs and signaling pathways involved in the prolificacy trait. These results provide a valuable resource for further research into the molecular mechanism underling the prolificacy in the follicular phase of goat ovary and new insights into the interactions among various factors related to the development and regulation of goat oocytes. Further studies should validate the functions of these key interaction networks and their components at the cellular level in goat prolificacy.

Ethics statement
All the experimental procedures described in the present study were approved by the Science Research Department (in charge of animal welfare issues) of the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences (IAS-CAAS) (Beijing, China). Ethical approval was provided by the animal ethics committee of IASC AAS (No. IAS2019-63).

Sample collection, RNA extraction and quality analysis
Female native domestic goats, known as Yunshang black goats, were used in this study. Ten goats with no significant differences in age, weight, or height were selected and divided into either the high-yielding group (n = 5, average litter size 3.00 ± 0.38, HF group) or the lowyielding group (n = 5, litter size 1.32 ± 0.19, LF group) (P < 0.05) according to their litter size records. Additionally, all the goats were fed under the same conditions, with free access to water, on a goat farm in Yunnan Province. Ovarian tissues were collected from goats in the follicular phase, frozen immediately in liquid nitrogen and stored at − 80°C until RNA extraction.
Total RNA for RNA sequencing (RNA-seq) was isolated from 10 ovarian tissue samples (mixed powders of the entire ovary) with the TRIzol reagent (Invitrogen, Carlsbad, CA, United States) according to the manufacturer's protocol and previous study [55]. The NanoDrop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE, United States) was used to evaluated on the purity and concentration of the RNA samples, and standard denaturing agarose gel electrophoresis was used to monitor degradation and contamination. The Agilent Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, CA, United States) with an RNA Nano 6000 Assay Kit was used to assess the integrity of total RNA. mRNA and miRNA library preparation and sequencing A total of 10 cDNA libraries were constructed from the ovarian tissues of Yunshang black goats from both groups (LF and HF). Three micrograms of RNA from each sample were used as input material for cDNA library and miRNA library construction. First, rRNA was removed with an Epicenter Ribo-ZeroTM rRNA Removal Kit (Epicenter, Madison, WI, United States), and rRNA-free residues were cleaned by ethanol precipitation. Second, following the manufacturer's recommendations, rRNA-depleted RNA was used to prepare sequencing libraries using the NEBNext® UltraTM Directional RNA Library Prep Kit and Illumina® (NEB, Ipswich, MA, United States). Finally, the products were purified by the AMPure XP system, and an Agilent Bioanalyzer 2100 system was used to assess the library quality. The Illumina Novaseq platform was used to sequence libraries, and 150 bp paired-end reads were generated.
Raw data in fastq format were first processed by inhouse scripts. After the removal of reads containing adapters, reads containing poly-N, and low-quality reads from the raw data, the Illumina sequencing raw reads were obtained, among which the number of bases with a quality value Q ≤ 20 was > 50%. The Q20, Q30, and GC contents of the clean data were calculated. Based on the clean data, all high-quality downstream areas were analyzed. Reference genome and gene model annotation files were downloaded directly from the genome website. Bowtie v2.2.3 was used to build the reference genome index, and TopHat v2.0.12 was used to align paired-end clean reads to the reference genome. Both known and novel transcripts from the TopHat alignment results were constructed and identified by the Cufflinks v2.1.1 Reference Annotation Based Transcript (RABT) assembly method. The fragments per kilobase per million reads (FPKM) [56] for each gene were calculated based on the length of the gene and read counts mapped to the gene.
Clean high-quality reads with lengths of 18-35 nt were used in the subsequent analyses. The small RNA tags were mapped to reference sequences in Bowtie21, which were used to search for known miRNAs. miRBase22.0 was used as a reference, and the potential miRNAs and secondary structures were obtained with the modified software mirdeep2 and srna-tools-cli. To remove tags originating from protein-coding genes, repeat sequences, rRNAs, tRNAs, snRNAs, and snoRNAs, small RNA tags were mapped to the Repeat Masker or Rfam databases or data from the specified species itself [57]. The novel miRNAs were predicted by the characteristics of the miRNA precursor hairpin structures. The novel miRNAs were predicted with the software miREvo and mirdeep2 by exploring the secondary structures, Dicer cleavage sites and minimum free energy of the small RNA tags unannotated in the former steps [57,58].

Analysis of DE mRNA and miRNA
The expression levels of mRNAs in both libraries constructed from ovaries were estimated based on the FPKM values of the Illumina sequencing data. The fold change of the expression of each mRNA between the HF and LF groups was calculated according to the comparison. Differential expression analysis of both groups was performed with the DESeq R package (1.8.3). Genes with an adjusted P < 0.05 and |FC| > 1 identified by DESeq were considered to be differentially expressed. The DEGseq R package was used to analyze the DEMs according to the normalized transcripts per kilobase per million reads (TPM) values. The q-valued were adjusted by the P values. A q-value < 0.01 and |FC| > 1 were set as the thresholds for significant DEMs by default.

Gene ontology and Kyoto encyclopedia of genes and genomes analyses
The DEGs targeted by miRNAs were annotated and classified by Gene Ontology (GO2) with the GOSeq (Re-lease2.12) software and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis with the KOBAS (v2.0) software to visualize the data [59][60][61]. Only GO terms or KEGG pathways with corrected P values (t-test) < 0.05 were considered to indicate significant enrichment.

Counter expression miRNA-mRNA analysis
The miRNA and transcriptome profiles of both groups were constructed using ovarian samples from Yunshang black goats. To analyze the interactions between miR-NAs and mRNAs, miRNA target genes were predicted by miRanda with psRobot_tar in psRobot [62,63]. Then, miRNA-mRNA interactions were calculated according to miRNA and transcriptome expression profile data, and negatively correlated miRNA-mRNA pairs were identified using Pearson correlation analysis. Finally, all these analyses were performed using "GeneSymbol" as a unique identifier for all genes/transcripts involved.

Verification of miRNA and mRNA expression profiles with RT-qPCR
For the RT-qPCR analysis, reverse transcription was performed using a PrimeScript™ RT reagent kit (TaKaRa) and miRcute Plus miRNA First-Strand cDNA Kit (TIANGEN, Beijing, China) for mRNA and miRNA according to the manufacturers' instructions and previous study [55], respectively. RT-qPCR was performed by a RocheLight Cycler®480 II system (Roche Applied Science, Mannheim, Germany) with SYBR Green qPCR Mix Kit (TaKaRa, Dalian, China) and miRcute Plus miRNA qPCR Kit (TIANGEN, Beijing, China) for mRNAs and miRNAs, respectively. RT-qPCR analysis of mRNA and miRNA expression was conducted by the following procedure: for mRNA, initial denaturation at 95°C for 5 min, followed by 40 cycles of denaturation at 95°C for 5 s and annealing at 60°C for 30 s; for miRNAs, initial denaturation at 95°C for 15 min, followed by 40 cycles of denaturation at 94°C for 20 s and annealing at 60°C for 34 s. The data were analyzed with the 2 -ΔΔCt method. The goat PRL19 gene and U6 were used as reference genes for the normalization of the target gene data. The sequences of the RT-qPCR primers are listed in Supplementary Table S6.

Protein-protein interaction (PPI) network analysis of DEGs
Known and predicted protein-protein interaction analyses of the DEGs were performed using the STRING database (https://string-db.org/). Because the goat database was not provided in STRING, the highly homologous species sheep was used as an alternative to analyze the correlation of DEGs. PPI analysis of DEGs was based on the STRING database4 (Organism: Ovis aries). The networks were constructed by extracting the target gene list from the database. Otherwise, the target gene sequences were aligned to the selected reference protein sequences by Blastx (v2.2.28), and then the networks were constructed according to the known interactions of the selected reference species.
Validation the function of miRNA-mRNA pairs Vector construction The 3′ UTR of TXLNA containing the predicted target site, including wild type (WT) and mutant type (MUT), was cloned into the pmiR-RB-Report vector (isolated using Xho I and Not I (Takara, Dalian, China)). The recombinant vectors were named pmiR-RB-TXLNA-WT and pmiR-RB-TXLNA-MUT. The insert sequences, including the wild-type/mutant sequences, miRNA mimic and mimic-NC, were synthesized by RiboBio company (Guangzhou, China) (the information is shown in Supplemental Fig. S1).
Cell culture, transfection and dual luciferase assays 293 T cells, a human renal epithelial cell line, were used to validate the miRNA targets. The cells were seeded into 24-well plates. Cotransfection with 200 ng mRNA-3'UTR-WT target or mRNA-3'UTR-MUT target and 10 μL of miRNA mimic or mimic-NC was performed with Lipofectamine 2000 (Invitrogen, USA). Then, the luciferase activities were measured using the Dual Luciferase Reporter Assay System (Promega, WI, USA) at 48 h post-transfection. The assays were performed in triplicate.

Cell counting Kit-8 (CCK-8) assay in granulosa cells (GCs)
Goat GCs were obtained from our laboratory. GCs were inoculated into 96-well plates with approximately 100 μL of cell suspension per well in 3 replicates. The cells were incubated for 2-4 h at 37°C in an incubator, and after cell appositioning, subsequent experiments were conducted. The cell proliferation assay was performed with the CCK-8 method (Beyotime, Beijing, China). Goat GCs were cultured in 96-well plates and assessed by adding 10 μL CCK-8 per well following the protocol. After the addition of the CCK-8 solution, the absorbance at 450 nm was measured at 0 h, 6 h, 12 h, 24 h, 48 h and 72 h to assess cell proliferation.

Statistical analysis
The GraphPad Prism (version 5.0) software (San Diego, CA, United States) was used to statistical analyses of the RT-qPCR results and graphs. The statistical significance of the data was tested by performing paired t-tests [64]. The results are presented as the means ± SEMs of three replicates, and statistical significance was represented by P values (*P < 0.05; **P < 0.01).