All experimental procedures mentioned in this study were approved by the Science Research Department (responsible for 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 IASCAAS (No. IAS2021-25).

In this study, ten female native domestic goats in a goat farm in Yunnan Province, named Yunshang black goat, were used. The height, weight, and age (2–3 years old) of goats, with no significant differences, were randomly selected and grouped into a high-prolificacy group (n = 5, mean kidding number = 3.4 ± 0.42, H group) and a low-prolificacy group (n = 5, mean kidding number = 1.8 ± 0.27, L group) according to their kidding number records. All animals were treated for estrus synchronization. This was done by treating female goats with progesterone suppositories for 16 days and collecting ovarian tissues at slaughter 11 days after stopping the hormones. The ovarian tissues were collected from goats and then immediately frozen in liquid nitrogen and stored at −80°C until RNA was extracted.

According to the manufacturer’s instruction, the total RNA of 10 ovarian samples was isolated by TRIzol reagent (Invitrogen, Carlsbad, CA, United States) for RNA sequencing (RNA-seq). The NanoDrop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE, United States) was used to assess the purity and concentration of the RNA samples, and the degree of degradation and contamination was monitored by standard denaturing agarose gel electrophoresis. An RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, CA, United States) was used to assess the integrity of RNA.

Ten cDNA libraries of the ovarian tissues from two Yunshang black goat groups with high- and low-prolificacy kidding numbers (H and L group) were constructed. A total of 3 µg of RNA from each sample was used as input material for the miRNA library and cDNA library. First, an Epicentre Ribo-Zero™ rRNA Removal Kit (Epicentre, Madison, WI, United States) was used to remove the ribosomal RNA, and the ethanol precipitation was used to clean up the rRNA-free residues. Secondly, following the manufacturer’s recommendations, the NEBNext^® Ultra™ Directional RNA Library Prep Kit for Illumina^® (NEB, Ipswich, MA, United States) was used in the rRNA-depleted RNA to generate sequencing libraries. Finally, an AMPure XP system was used to purify the products, and an Agilent Bioanalyzer 2100 system was used to assess the quality of the libraries. An Illumina Hiseq 2500 platform was used to sequence the libraries, which generated 150 bp of paired-end reads.

The fastq format raw data were initially processed by in-house scripts. Subsequently, the raw reads of Illumina sequencing were used for the removal of reads containing adapter, ploy-N, and low-quality reads, and the quality value Q ≤ 20 was maintained over 50%. The size of clean data was calculated depending on the Q20, Q30, and GC content. Because all downstream analyses were based on clean data, they must be maintained at a high level of quality. We directly downloaded the goat reference genome and the files of gene model annotation from the genome website. Bowtie v2.2.3 was used to build the index of the reference genome, and TopHat v2.0.12 was used to align the pair of clean reads to the reference genome. The 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. According to the length of the gene and read counts mapped to the gene, the fragment per kilobase per million reads (FPKM) for each gene was calculated (18).

In the subsequent analyses, the clean high-quality reads with lengths of 18–35 nt were used. The small RNA tags were mapped to Bowtie21 reference sequences, which were used to find out the known miRNAs. Using miRBase22.0 as a reference database, the potential miRNAs were obtained by the modified software mirdeep2 and srna-tools-cli, and then their secondary structures were mapped (19). In order to remove tags derived from repetitive sequences, protein-coding genes, tRNAs, rRNAs, snoRNAs, and snRNAs, the tags of small RNA were mapped to the Repeat Masker or Rfam databases or to data from the specified species. The characteristics of the miRNA precursor hairpin structures could be used to predict the novel miRNAs. The secondary structure Dicer cleavage sites and the minimum free energy of unannotated small RNA tags in the previous steps of novel miRNAs were predicted by the software miREvo and mirdeep2 (19, 20).

The mRNA expression levels in the constructed libraries of ovarian tissues were estimated according to the sequencing data of Illumina, based on the values of FPKM. The fold change (FC) for mRNA between the H and L groups was calculated according to comparisons of the combination. The DESeq R package (1.8.3) was used to analyze the DE mRNAs from the two groups. The DE genes were adjusted by q < 0.05 and |log₂FoldChange| > 1 of DESeq. According to the values of normalized transcripts per kilobase per million reads (TPM), the differentially expressed miRNAs (DEMs) were analyzed by the DEGseq R package. The q-values were adjusted by p-values, and the significant DEMs were set by default by the threshold q-values < 0.01 and |log₂FoldChange| > 1.

The miRNA and transcriptome profiles of the two groups with different kidding numbers were constructed by using the ovarian tissues from Yunshang black goats. To analyze the interactions between miRNAs and mRNAs, first, miRanda (21) was used to predict miRNA target genes with psRobot_tar in psRobot (22). Then, according to miRNA expression profiles and transcriptome data, the interactions of miRNA–mRNA were calculated, and Pearson correlation analysis was used to determine the negative correlation of miRNA–mRNA pairs. Finally, the “GeneSymbol” as a unique identifier was created for all genes/transcripts involved in all these analyses.

The Gene Ontology (GO2) with the software GOSeq (Release2.12) (23) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis with the software KOBAS (v2.0) (24) were used to annotate and classify mRNAs and enriched genes among miRNA-targeted mRNAs to visualize data. The enrichment on GO terms or KEGG pathways with corrected p-values (t-test) < 0.05 was indicated as significant.

Using the STRING database (https://string-db.org/), the PPI analysis of DE genes was known and the protein–protein interactions were predicted. The STRING database4 (organism: Capra hicus) was used for PPI analysis of DEGs. The networks were constructed by the extracted target gene list from the database. In addition, the sequences of target genes and the selected reference protein were aligned by Blastx (v2.2.28), and then the known interaction of selected reference species was used to build the networks.

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

The sequencing results were subjected to a BLAST search (NCBI, http://blast.ncbi.nlm.nih.gov/Blast.cgi) to retrieve homologous gene sequences of JAK3. Multiple nucleic acid sequences were compared using DNAMAN software, and a phylogenetic tree was constructed using the Neighbor-Joining method in MEGA7 software.

The total protein of ovary tissues was isolated by the proteinase inhibitor-containing lysis buffer. The 4% SurePAGE gel (GenScript, Nanjing, China) was used to separate the equivalent amounts of protein. The separated proteins were placed onto a PVDF membrane (Pall, Mexico) after electrophoresis, and then blocked with sealing solution (Tiangen, Beijing, China). The blocked membrane was incubated overnight at 4°C with anti-JAK3 (1:500; Santa Cruz, Cambridge, UK), anti-CDK4 (1:1,000), anti-cyclin D2 (1:1,000) (all from Cell Signaling Technology, US), and GAPDH (1:2,000; Proteintech, Chicago, IL, USA). After rinsing 3 times with Tris-buffered saline/Tween, the corresponding HRP-labeled Goat Anti-Rabbit IgG (1:1,000; Beyotime) or HRP-labeled Goat Anti-Mouse IgG (1:1,000; Beyotime) was used to incubate the membranes for 1 h at room temperature. The Enhanced Chemiluminescent Reagent (Beyotime) was used to visualize the protein blots.

Following previously described methods, primary goat GCs were isolated from the follicles of ovarian tissues in the follicular phase (25). The diameter of the GCs in the follicular phase was over 3.5 mm. The isolated cells were seeded in 6-mm plates and maintained in a complete medium [DMEM/F12 (1:1), 10% FBS, and 1% penicillin/streptomycin] as described by Yang et al. (26). When the cell confluence was >90%, the cells were transferred to 10-mm plates for the next experiment.

The region of JAK3 3’ UTR containing the predicted target site, including wild type (WT) and mutant type (MUT), was cloned into a pmiR-RB-Report vector [isolated using XhoI and NotI, respectively (Takara, Dalian, China)]. The recombinant vectors were named pmiR-RB-JAK3-WT/pmiR-RB-JAK3-MUT. The inserted wild/mutant sequence, miRNA mimic, and mimic-NC were synthesized by the RiboBio company (Guangzhou, China). The coding region of JAK3 was cloned into a pcDNA3.1 expression vector, and the recombinant vectors were named pcDNA3.1-JAK3 and pcDNA3.1-NC. siRNA-JAK3, siRNA-NC, and the chi-miR-493-3p mimic and inhibitor were synthesized by Genewiz (Suzhou, China). A human renal epithelial cell line, HEK293T, was used to validate the target of miRNA. Cells were seeded into 24-well plates. Co-transfection with 200 ng of target JAK3-3’UTR-WT or target JAK3-3’UTR-MUT and 10 µl of miRNAs mimic or mimic-NC was performed with Lipofectamine 2000 (Invitrogen, USA). Goat primary GCs were seeded into 6-well plates and transfected with pcDNA3.1-JAK3, siRNA-JAK3, and the chi-miR-493-3p mimic or inhibitor (RiboBio, Guangzhou, China); meanwhile, their NC was used as negative control. Furthermore, after transfection, cells were collected for transfection efficiency analysis, and the proliferation of cells and the expression of genes were examined after 48-h transfection.

The goat GCs with an added growth medium were seeded in 96-well cell culture dishes. When the cell density reached the appropriate level, the cells were transfected. The proliferation of cells after 0, 6, 12, 24, and 48 h of transfection was detected according to the instruction protocol of CCK-8 Kit (Vazyme, Nanjing, China), and the multi-mode micropore detection system (EnSpire, Perkin Elmer, USA) was used to detect the absorbance at 450 nm.

HEK293T cells and GCs were detected by the dual luciferase assay. HEK293T cells and GCs were evenly seeded in a 24-well cell culture dish. When the cells reached the appropriate density, they were co-transfected with pmiR-RB-JAK3-WT (wild type) or pmiR-RB-JAK3-MUT (mutant type) and chi-miR-493-3p mimic or mimic-NC. After 48-h transfection, according to the instructions of the double luciferase detection kit (Vazyme, Nanjing, China), the luciferase activity was detected by the multi-mode micropore detection system (EnSpire, Perkin Elmer, USA).

ELISA was used to detect the levels of E₂ and PROG in GCs with different treatment according to the protocol of the E₂ ELISA Kit and the PROG ELISA Kit (Enzyme-linked Biotechnology Co., Ltd., Shanghai, China). The goat GC culture supernatant was transfected, and pcDNA3.1-JAK3, siRNA-JAK3, and their NC were centrifuged at 2000 × g for 5 min at room temperature. Then, 10 μl of the sample supernatant was collected and added to the bottom of the plate wells along with 40 μl of sample diluents and 100 μl of enzyme-labeling reagents. After sealing the plates with sealing film, they were incubated at 37°C for 60 min. Finally, we followed standard procedures of washing, color development, and termination of the plates. The absorbance at 450 nm was determined by using a Fluorescence/Multi-Detection Microplate Reader (Bio Rad, USA).

The GraphPad Prism (version 5.0) software (San, Diego, CA, United States) was used to analyze the results of RT-qPCR and for graphing. The paired t-tests were used to test the statistical significance of the data. The results were expressed as the means ± SEM of three replicates, and the statistical significance was expressed as p-values < 0.05 (*p < 0.05; **p < 0.01).

The cDNA and small RNA libraries were contrasted by RNA-seq with 10 goat ovaries, and after quality screening, 115 million and 26 million clean read pairs were obtained from all cDNA and small RNA libraries, respectively, of which 96.99% of mRNA and 94.88% of miRNA reads from the samples, on average, were mapped to the goat genome (Capra hircus, ARS1) ( Table 1 ). A series of filter criteria were used to determine mRNA and miRNA candidates. A total of 28,618 transcripts ( Figure 1A ) and 1,038 miRNAs ( Figure 1B ) were screeded in the ovary, and 2,968 differentially expressed mRNAs (DEGs) and 99 differentially expressed miRNAs (DEMs) were identified between the H and L groups. In the H vs. L comparison, there were 1,553 upregulated and 1,415 downregulated DEGs, and 80 upregulated and 19 downregulated DEMs ( Figures 1C, D ). All the identified DEGs and DEMs are presented in Supplementary Table S2 .

The identified DEMs and DEGs from goat ovarian tissues of individuals were integrated to obtain the miRNA–mRNA network associated with goat prolific traits. A total of 671 miRNA–mRNA pairs were predicted in the H vs. L comparison ( Figure 2 and Supplementary Table S3 ), of which the core gene of network DUOX2 is targeted by 5 miRNAs, namely, miR-133b, novel_146, novel_328, novel_353, and novel_59. The major function of miRNAs is to suppress the expression of target genes; therefore, the miRNA–mRNA pairs with negatively correlated expression levels are the focus of research. Among these DEG–DEM pairs, 386 pairs were negatively correlated in expression (q-value < 0.05), which might play key roles in the molecular mechanism of the goat prolific trait. In the network, miR-93-3p was a core miRNA with 73 miRNA–mRNA pairs, 29 of which showed a negative correlation. Another important core miRNA was miR-493-3p, which targets 7 genes including JAK3, MPDZ, and CLEC16A; these genes were closely related to the JAK-STAT signaling pathway, tight junction, and the PI3K-Akt signaling pathway. A previous study has shown that the PI3K-Akt signaling and JAK-STAT signaling pathways were associated with ovarian function, such as the primordial follicle recruitment, the proliferation of GCs, the survival of corpus luteum, and oocyte maturation (27, 28). JAK3 was an important member of the membrane-related intracellular non-receptor tyrosine kinase protein family and participated in the JAK-STAT signaling pathway and PI3K-Akt signaling pathway. According to their potentially preserved functions, two regulators, miR-493-3p and JAK3, were chosen for further analyses in the studies related to goat prolific trait.

In order to further investigate the functions of DEG–DEM pairs between the different fertility performance groups, the GO enrichment analysis was performed to reveal the biological process terms of the important DEGs (corrected q-values < 0.05). In this comparison, the associations were mainly with cytoplasm and nucleoplasm ( Figure 3A ). Then, the GO terms related to the reproduction trait are also displayed in Supplementary Table S4 , including the positive regulation of the receptor signaling pathway by JAK-STAT, estrogen biosynthetic process, estrogen metabolic process, and ovarian follicle development. To further explore the function of DEMs–DEGs in regulation pathways for goat prolific trait, KEGG pathway analysis of miRNA–mRNA pairs was performed ( Figure 3B ). The enrichment of KEGG pathway in the pathway categories of the comparison is shown in Supplementary Table S4 . In the comparison, the neurodegeneration–multiple diseases pathway (chx05022) and the calcium signaling pathway (chx04020) were the top two pathways. Among them, there eight signaling pathways were directly and indirectly associated with fertility trait, namely, the GnRH signaling pathway (chx04912), oocyte meiosis (chx04114), the PI3K-Akt signaling pathway chx04151), steroid biosynthesis (chx00100), ovarian steroidogenesis (chx04913), the JAK-STAT signaling pathway (chx04630), GnRH secretion (chx04929), and the estrogen signaling pathway (chx04915). A total of 21 DEGs were enriched in the eight pathways, including JAK3, ADCY3, ADCY6, IKBKB, STAT3, and other genes shown in Supplementary Table S4 , which might be important candidates for further study.

Subsequently, the list of target genes extracted from the STRING database was used to construct the PPI network ( Figure 3C ). The PPI network on DEGs obtained from the H vs. L comparison contains 741 protein–protein nodes, of which 136 have a combined score of over 0.7. Thereinto, the key core nodes of the network included ITGB1, IRF4, ACTN4, CPSF4, NCOR2, KIT, FLNA, CSF2RA, COPB1, and JAK3, all of which were targets of DE miRNAs (shown in Supplementary Table S5 ).

Ten DEG–DEM pairs from RNA-seq, namely, novel_530-FBXL17, chi-miR-93-3p-ITGA2, chi-miR-133b-RAB3IP, novel_235-ADGRF5, novel_328-FBXO4, novel_353-FOXP1, novel_615-IGF1R, chi-miR-495-3p-RECK, chi-miR-485-5p-TCF7L2, and chi-miR-493-3p-JAK3, were randomly selected for RT-qPCR validation. The results suggested that the expression levels of mRNAs and miRNAs agreed with the data from RNA-seq, indicating that the data from RNA-seq were plausible and that they have a negative correlation, indicating a potential interaction ( Figure 4 ).

The proliferation of GCs is essential for ovulation in mammals. To reveal the molecular mechanism of the goat prolific trait, the function of JAK3 in the goat GCs was furthermore studied. The Western blot results showed that the JAK3 protein expression was significantly higher in the H group than in the L group (p < 0.05) ( Figure 5A ), as opposed to JAK3 mRNA expression. This result showed that the expression of JAK3 was regulated by the miRNAs in the post-transcript level. The phylogenetic tree analysis found that JAK3 was highly conserved in different species and that the goat JAK3 was first clustered with sheep, followed by cattle ( Figure 5B ). Hence, we supposed that the function of JAK3 was similar to other species. To explore the role of JAK3 in goat GC proliferation, JAK3 overexpression and siRNA vectors were constructed. Overexpression of JAK3 induced an increase in mRNA expression of the cell proliferation markers CDK4, cyclin D1, and cyclin D2 after 24 h post-transfection ( Figure 5C ). These results demonstrated that overexpression of JAK3 promotes the proliferation of goat GCs. In contrast, we observed the opposite pattern in the knock-down experiment using si-JAK3 ( Figure 5D ). Western blot results showed that the protein expression of JAK3, cyclin D2, and CDK4 was consistent with mRNA expression ( Figure 5E ). Overall, these results demonstrated that JAK3 promoted the proliferation of goat GCs.

Integrative analysis of DE miRNAs and mRNAs showed that JAK3 was one of the targets of chi-miR-493-3p. To determine the function of the chi-miR-493-3p-JAK3 pathway, the dual-luciferase reporter assays were performed to validate whether there was a direct interaction between chi-miR-493-3p and JAK3. The luciferase activity of the JAK3 wild-type plasmid with the chi-miR-493-3p mimic was significantly lower than that of the JAK3 mutant plasmid (p < 0.01) ( Figure 6A ), supporting the hypothesis that JAK3 is a direct target of chi-miR-493-3p. The expression of JAK3 mRNA and protein in goat GCs was determined after transfection with the chi-miR-493-3p mimic or chi-miR-493-3p inhibitor. The results showed significantly decreased mRNA ( Figure 6B ) and protein ( Figure 6D ) expression of JAK3 in the chi-miR-493-3p mimic transfected in GCs of goat (p < 0.01). In goat GCs transfected with the chi-miR-493-3p inhibitor, both mRNA and protein expression levels of JAK3 were significantly increased (p < 0.01) ( Figures 6B, D ). These results suggested that JAK3 is a direct target of chi-miR-493-3p in goats.

To investigate the function of chi-miR-493-3p in the proliferation of goat GCs, chi-miR-493-3p was overexpressed or inhibited in goat GCs. RT-qPCR results demonstrated that overexpression of chi-miR-493-3p in GCs significantly increased the expression level of chi-miR-493-3p, and reduced the expression of JAK3, as well as decreased the expression levels of proliferation-related genes (CDK4, cyclin D1, and cyclin D2) ( Figure 6B ). Moreover, CCK-8 assays indicated that the proliferative state of these GCs was depressed after transfection ( Figure 6C ). However, the GCs transfected with the chi-miR-493-3p inhibitor showed a significantly decreased level of chi-miR-493-3p mRNA and significantly increased levels of JAK3, CDK4, cyclin D1, and cyclin D2 mRNA ( Figure 6B ). Subsequently, CCK-8 assays indicated that the proliferative state of these cells was promoted after transfection ( Figure 6C ). We also investigated the expression of JAK3, cyclin D2, and CDK4 proteins in GCs after transfected with the chi-miR-493-3p mimic and inhibitor and found that the protein levels of JAK3, cyclin D2, and CDK4 were all decreased in cells transfected with the chi-miR-493-3p mimic, whereas cells transfected with the chi-miR-493-3p inhibitor showed increased levels of JAK3, cyclin D2, and CDK3 proteins ( Figure 6D ). These results indicated that chi-miR-493-3p suppressed goat GC proliferation.

After transfection of pcDNA3.1-JAK3 and siRNA-JAK3 in goat GCs for 24 h, the E₂ and PROG levels in the cell-free supernatants were measured by ELISA to investigate the effect of JAK3 on synthesis of steroid hormones. JAK3 overexpressing cells showed enhanced secretion of E₂ and PROG ( Figure 7A ), while inhibition of JAK3 reduced E₂ and PROG production ( Figure 7B ). These results suggested that chi-miR-493-3p affected steroid hormone synthesis in goat GCs by regulating the expression of JAK3.