Genome-wide analysis of miRNAs in the ovaries of Jining Grey and Laiwu Black goats to explore the regulation of fecundity

Goat fecundity is important for agriculture and varies depending on the genetic background of the goat. Two excellent domestic breeds in China, the Jining Grey and Laiwu Black goats, have different fecundity and prolificacies. To explore the potential miRNAs that regulate the expression of the genes involved in these prolific differences and to potentially discover new miRNAs, we performed a genome-wide analysis of the miRNAs in the ovaries from these two goats using RNA-Seq technology. Thirty miRNAs were differentially expressed between the Jining Grey and Laiwu Black goats. Gene Ontology and KEGG pathway analyses revealed that the target genes of the differentially expressed miRNAs were significantly enriched in several biological processes and pathways. A protein-protein interaction analysis indicated that the miRNAs and their target genes were related to the reproduction complex regulation network. The differential miRNA expression profiles found in the ovaries between the two distinctive breeds of goats studied here provide a unique resource for addressing fecundity differences in goats.

Length distribution and base preference of the miRNAs. The length of the mapped miRNAs in this study was around 22 nt. Previous studies suggest that the miRNAs with a base U at the 5′ terminus are more conserved than miRNAs with a non-U base at the 5′ terminus 24 . We analyzed the base at the 5′ terminus of the miRNAs identified, and the results showed that U is the most common base in both of the goat breeds.
Differential regulation of the miRNAs between the two breeds of goats. A total of 5254 miRNAs were identified in the two goat samples, of which 603 (314 mature and 289 predicted, Table S1) had at least 10 read counts in the two goat samples. The expression level of these 603 conserved and novel miRNAs was assessed by a differential analysis, a downstream target prediction and biological annotation enrichment analyses. At a false discovery rate (FDR) < 0.05 and an absolute value of log 2 (fold change) > = 1, 30 miRNAs (  8 predicted) were defined as differentially expressed miRNAs between the Laiwu Black and Jining Grey goats (Table 3), of which 10 miRNAs were down-regulated and 20 were up-regulated in the Jining Grey goat as compared with the Laiwu Black goat.
Prediction of the miRNA targets. The programs TargetFinder 25 and miRanda 26 were used to assess the 3′ -UTR regions of the potential target genes. A total of 15,848 gene targets without redundancy were prediced from the 30 differential miRNAs (30,172 pairs) with a free enengy (Δ G) of less than or equal to − 15.0, with 1-7 miR-NAs targeting a single gene or 17-4357 genes targeted by a single miRNA. Of these 15, 848 gene targets, known miRNAs gene targets without redundancy were 13,699, and predicted miRNAs gene targets without redundancy were 6,703 (Table S2).
Gene Ontology and KEGG pathway analysis. The predicted miRNA-targeted genes were further analyzed by Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. The target genes that interact with conserved and predicted miRNAs were analyzed separately, with regard to the enrichment of the GO terms and the KEGG pathway analysis. For the mature miRNA-targeted genes, there were no significantly enriched GO terms and KEGG pathways. In contrast, the predicted miRNA-targeted genes were enriched significantly in 28 biological processes and 2 KEGG pathways (Fig. 1). Among the 28 biological processes, 19 were associated with reproduction, such as steroid biosynthetic process, the BMP signaling pathway, and the transforming growth factor beta receptor signaling pathway (Figs 1A and 2A). Glycosaminoglycans that are localized to the extracellular matrix of bone are thought to play a key role in mediating aspects of bone development. The two enriched KEGG pathways were the TGF-β signaling pathway and steroid biosynthesis (Figs 1B and 2B). Indeed, these functional annotation results demonstrate that this database contains clues that might lead to uncovering potential regulators of fecundity between these 2 breeds of goats.
Interaction networks between miRNAs and the target genes related to reproduction. Given the fact that we are interested in exploring the fecundity differences between these 2 breed of goats, the relationships between the differentially expressed miRNAs and the targeted genes related with reproduction are described  based on the STRING database analysis. Figure 3 demonstrates that the relationships between the miRNAs and the genes BMPR 1B, SMAD1, BMP7, ACVR1, CHRD, BMPR2, BMP88, ACVR18, SMAD4, THBS1, TGFB1, MAPK3 and SMAD3, which are related to reproduction, are complicated. The BMPR1B gene, which is regulated by chi-miR-493-3P alone, is a very important finding. In addition, chi-miR-187 regulates four genes, including  TGFB1, THBS1, ACVR18 and BMP88, and chr12_10768_star regulates three genes, including CHRD, SMAD1 and BMP7. In addition, chi-miR-874-3P regulates three genes, including MAPK3, BMPR2 and CHRD. Therefore, chi-miR-187, chr12_10768_star and chi-miR-874-3P might play an important role in reproductive regulation processes.
qRT-PCR validation of RNA-Seq data. To confirm the RNA-Seq data, 6 differentially expressed miRNAs were examined by qRT-PCR, and these data were consistent with the above RNA-Seq results (Table 4 and Fig. 4). The genes potentially targeted by the 6 miRNAs are involved in diverse cellular activities, such as signal transduction, kinase, motor activity, transport, metabolic process and DNA binding.

Discussion
With the use of RNA-Seq technology, here, we performed a genome-wide analysis on the miRNAs in Jining Grey and Laiwu Black goats and characterized 30 differentially expressed miRNAs. Because these two breeds of goats demonstrate different fecundity, the differentially regulated miRNAs may contribute to this process. Previous studies indicate that several miRNAs are important for the follicular-luteal transition in the ruminant ovary and for fetal gonadal development 27,28 . Homologs of the identified differentially expressed miRNAs function in various cellular activities. For instance, miR-9-3p is an important regulatory factor in the osteoblastic differentiation of mouse iPS cells and also targets β 1 integrin to sensitize claudin-low breast cancer cells to MEK inhibition 29,30 . miR-449a regulates cell proliferation and causes Rb-dependent cell cycle arrest and senescence 31,32 . Mature miR-183 promotes 3T3-L1 adipogenesis through the inhibition of the canonical Wnt/β -catenin signaling pathway by targeting LRP6 33 . miR-542-3p also has many important functions, including the induction of growth arrest through the survival pathway, the suppression of osteoblast cell proliferation and differentiation, the targeting of BMP-7 signaling and the inhibition of bone formation 34,35 . Given the broad roles of the miRNA homologs, the differentially expressed miRNA molecules might also potentially be involved in the regulation of goat growth and development. miRNAs function mainly by pairing with mRNAs, and by an in silico analysis, we predicted the potential target genes for the differentially expressed miRNAs identified in our RNA-Seq analysis. Using the predicted target genes, the Gene Ontology (GO) enrichment analysis revealed several significantly enriched biological processes that were associated with the steroid biosynthetic process, the BMP signaling pathway and the transforming growth factor beta receptor signaling pathway (Fig. 1). The glutathione metabolic process is involved in a variety of protective processes, and we found that, in this process, 15 differentially expressed miRNAs potentially target 16 genes ( Fig. 2A), which include the homologs of miR-9-3p, miR-183, miR-449a and miR-542-3p, which were mentioned in the prior section. In addtion, more favorable cellular processes in the ovaries, such as the steroid biosynthetic process, the BMP signaling pathway, the transforming growth factor beta receptor signaling pathway, contribute to the high prolificacy observed in the Jining Grey goat. Moreover, among the differentially expressed miRNAs that were identified to be involved in the glutathione metabolic process, 15 of the predicted target genes are also involved in the bone development (Fig. 2B). A recent study reported 36 that certain genes in ovaries were involved in cartilage development and bone healing, which is consistent with our findings. The KEGG pathway analysis was consistent with the GO analysis, indicating that the miRNA targeted genes were significantly enriched for glutathione metabolism and the steroid biosynthesis pathways.
Moreover, the relationship between the miRNAs and the genes BMPR1B, SMAD1, BMP7, ACVR1, CHRD, BMPR2, BMP88, ACVR18, SMAD4, THBS1, TGFB1, MAPK3 and SMAD3, which are related with reproduction 37 , is complicated. The BMPR1B gene, which is regulated by chi-miR-493-3P alone, was a very important finding. BMPR1B (FecB), was the first major prolificacy gene identified in sheep. The FecB locus is autosomal with a codominant expression and has an additive effect for ovulation rate when it is associated with a mutation (Q249R) in the BMPR1B gene 38,39 . BMPR1B is a member of the bone morphogenetic protein (BMP) receptor family of transmembrane serine/threonine kinases. BMPs are involved in endochondral bone formation and embryogenesis. BMPR1B is a member of the family of receptors for transforming growth factor (TGF) ligands, such as TGF, activins, BMPs and growth and differentiation factors (GDFs), which play a role in folliculogenesis 40 .
In conclusion, we identified 30 differentially expressed miRNAs in the ovaries of the two goats and most of these had defined homologs that were key miRNAs related to reproduction. This genome-wide miRNA analysis study will serve as a resource for understanding goat proliferacy. The differential miRNAs identified here are predicted to contribute to different prolificacies of the two goat breeds through a number of biological processes and pathways. In particular, Chi-miR-187, chr12_10768_star, and chi-miR-874-3P may play an important role in the reproductive regulation processes. This study provides a unique resource to address important issues related to goat fecundity, which will aid in elucidating the specific mechanisms responsible for differential gene expression in the ovaries between two distinctive breeds of goats. Moreover, these data also highlight the important roles of miRNAs during genetic regulation in different goat breeds.

Materials and Methods
Ethics statement. All of the experiments were performed in accordance with the relevant guidelines and regulations and were approved by the Institutional Animal Care and Use Committee of Institute of Animal Sciences, Chinese Academy of Agricultural Sciences.
Goat sample preparation. All of the samples in this study were from female adult goats between the ages of 1.5 to 2 years old. The age-matched healthy female Jining Grey goats and female Laiwu Black goats (5 per breed) were obtained from the Qingdao Aote Farm (Shandong, China). The goats were treated with intravaginal sponges (Intervet, Mexico) to achieve estrous synchronization as previously described 41 . Pregnant mare serum gonadotropin (PMSG) (Ningbo Sansheng Pharmaceutical Co., LTD., Zhejiang, China) was injected after the removal of the sponge. Estrus was checked by observing the goats' reactions to an adult male goat every 12 h after the sponge was removed. The female goats were determined to be in estrus when they showed standing estrous behavior. The goats were sacrificed 4-5 hours after estrus was determined, and the whole ovaries were excised. The samples were collected to obtain ovulation points on the surface of the ovaries. All of the samples were immediately snap-frozen in liquid nitrogen and stored at − 70 °C for RNA extraction 41 .

Construction of small RNA libraries and sequencing.
Total RNA was extracted from ovaries with better ovulation points on the surfaces from 5 Jining Grey and 5 Laiwu Black goats using Trizol (Life Technologies/ Invitrogen, California, USA) according to the manufacturer's instructions 42 . The extracted RNA from each breed was pooled. The quality and quantity of the RNA samples were assessed on a Bioanalyzer 2100 system using an RNA 6000 Nano kit (Agilent Technologies, Palo Alto, CA). After the total RNA was extracted, two small RNA libraries were prepared as described in the Illumina ® TruSeq ™ Small RNA Sample Preparation protocol. First, the 3′ adaptor and 5′ adaptor were ligated to the RNA by T4 RNA ligase, and then, reverse transcription was carried out with SuperScriptII reverse transcriptase (Invitrogen) to generate the cDNA. The resulting cDNA was then amplified to generate the small RNA library. The DNA size and the purity of the cDNA library were checked using a high sensitivity DNA 1000 kit on a Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA), and the quantification of the cDNA libraries was performed with Qubit ™ dsDNA HS kit on a Qubit ® 2.0 Fluorometer (Life Technologies, CA, USA). The cDNA libraries were then subjected to single-end sequencing on the Illumina Genome AnalyzerIIx (GAIIx) using the proprietary Solexa sequencing-by-synthesis method at the Shanghai Biotechnology Corporation (Shanghai, China) according to the manufacturer's recommended cycling parameters. The image analysis and the base calling were performed with the Illumina built-in SCS2.8/RTA1.8 software.
Workflow of the bioinformatics analysis. The workflow of bioinformatics analysis of the RNA-Seq results is shown in Fig. 5. Briefly, the adaptor sequences and the low quality sequences were removed from the original reads by fastx (fastx_toolkit-0.0.13.2) (http://hannonlab.cshl.edu/fastx_toolkit/), and the quality read data were aligned with the whole reference genome and were compared with multiple databases for the annotation of the ncRNA, miRNA, repeat sequences, and exon/intron sequences. The differential expression of the miRNAs was identified using IDEG6 software 43 . For the unannotated small RNA sequences, a miRNA prediction process was carried out, and then, the novel miRNAs and the targets of the differentially expressed miRNAs were annotated.
Programs used for the bioinformatics analysis. To predict potential novel miRNAs, Mireap software 23 was used to analyze the unannotated reads. The inverted repeats (step loops or hairpin structure) described by Jones-Rhoades and Bartel 44 were searched, and the secondary structure of the inverted repeat was predicted by RNA fold 45 . To predict the targets of the miRNAs, the programs TargetFinder 25 and miRanda 26 were used. The differentially expressed miRNAs were identified using IDEG6 software 43 . A rigorous significance test for the digital gene expression profiling was used 46 . The differentially expressed gene contigs were analyzed by gene ontology (GO) 47 . Three structured controlled vocabularies (ontologies), including cellular component, molecular function, and biological process, were used to create a consistent description of the gene products as previously described 48,49 . Like the GO enrichment, the association of the genes with different pathways was computed using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg) 50,51 . Moreover, based on the String database (http://string-db.org/), the interaction networks of the miRNAs and their target genes, related with reproduction, were predicted to draw relationships between the miRNAs and their target genes. Validation of RNA-Seq data. To validate the RNA-Seq results, the expression of several of the miRNAs was confirmed by quantitative real-time RT-PCR. Briefly, the cDNA for the miRNA was produced using a miScript II Reverse Transcription Kit (Qiagen, US) in a GeneAmp ® PCR System 9700 (Applied Biosystems, USA).
Next, the cDNA was used as the template for the qPCR, which was performed using a SYBR Green I Master (Roche, Swiss). The primer sequences were designed in the laboratory and were synthesized by Generay Biotech (Generay, PRC) and are based on the mRNA sequences obtained from the NCBI database. The gene U6 was used as the internal control, and the expression levels were calculated using the 2 delta delta Ct method (2− Δ Δ Ct).