Skip to main content
Download PDF

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

BMC Genomics volume 22, Article number: 812 (2021) Cite this article

Abstract

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.

Peer Review reports

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 miRNAs, 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.

Fig. 1
figure 1

Differentially expressed genes (DEGs) and miRNAs (DEMs) in goat ovaries. A Volcano diagram of DEGs in the comparison. B Volcano diagram of DEMs in the comparison

Full size image

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.

Table 1 DE mRNA-miRNA pairs in the ovaries of Yunshang black goats
Full size table
Fig. 2
figure 2

The top 20 GO terms enriched in the comparison

Full size image
Table 2 GO terms related to reproduction traits
Full size table

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.

Fig. 3
figure 3

The DE miRNA-mRNA interaction network in the comparison. Red color: upregulated; Green color: downregulated

Full size image

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.

Table 3 The KEGG pathways enriched in the comparison
Full size table

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.

Fig. 4
figure 4

The protein–protein interaction (PPI) network of differentially expressed proteins in the comparison

Full size image

Verification of miRNA and mRNA expression profiles with RT–qPCR

The frequently occurring genes in the miRNA-mRNA pairs, including CLK3, RPS24, TFDP2, PPP1R10, EXOC7, DENND1A, TXLNA, MS4A1, CD19 and WNT4, were selected for RT–qPCR validation. The gene expression levels were determined and compared with the RNA-seq data, which showed similar patterns (Fig. 5), suggesting that the RNA-seq data were credible. Furthermore, seven miRNA-mRNA connections were negatively correlated, including miR-197-5p-CLK3, miR-133a-5p-RPS24, novel_438-PPP1R10, miR-129-3p-EXOC7, miR-1271-3p-TXLNA and miR-2290-WNT4, indicating that they have potential interactivity (Fig. 6).

Fig. 5
figure 5

Real-time quantitative PCR (RT–qPCR) validation of differentially expressed mRNAs

Full size image
Fig. 6
figure 6

Real-time quantitative PCR validation of the differentially expressed miRNAs and their corresponding target mRNAs

Full size image

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).

Fig. 7
figure 7

Verification of the predicted miRNA-mRNA pairs using a luciferase assay. *P < 0.05

Full size image
Fig. 8
figure 8

Overexpression and inhibition of miR-1271 in GCs stimulated cell proliferation, as determined by the CCK-8 method. *P < 0.05; **P < 0.01

Full size image

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 development and utilization value. Therefore, it is crucial to understand the molecular mechanisms underlying prolificacy in Yunshang black goats for their development and utilization. Accompanied by the increase in FSH production, the follicular phase includes follicle recruitment and dominant follicle development. During this period, multiple follicles develop and ovulate [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 ARHGAP17 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 Ca2+ 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 meiosis-specific 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.

Materials and methods

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 IASCAAS (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 low-yielding 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 in-house 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 (Release2.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 miRNAs 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).

Availability of data and materials

All data generated and analyzed during this study are included in this published article and its supplementary information files. The raw data for this study can be found in BioProject ID PRJNA731513 on NCBI. The URL is https://www.ncbi.nlm.nih.gov/sra/PRJNA731513, accessed date: 23 May 2021.

Abbreviations

GO:

Gene Ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

PPI:

Protein–protein interactions

DE:

Differentially express

GC:

Granulose cellular

PCOS:

Polycystic ovary syndrome

BMPR1B:

Bone morphogenetic protein receptor type 1B

GDF9:

Growth differentiation Factor 9

IGF:

Insulin-like growth factor

ATM:

Ataxia telangiectasia mutated

FOXL2:

Forkhead Box L2

TAGLN2:

Transgelin 2

MAP 3 K3:

Mitogen-activated protein kinase kinase kinase 3

ERK:

Extracellular signal-regulated kinase

DEGs:

Differentially expressed genes

DEMs:

Differentially expressed miRNAs

LH:

Luteinizing hormone

RABT:

Reference Annotation Based Transcript

FPKM:

Fragments per kilobase per million reads

FC:

Fold change

TPM:

Transcripts per kilobase per million reads

References

  1. Haresign W. The physiological basis for variation in ovulation rate and litter size in sheep: a review. Livest Prod Sci. 1985;13(1):3–20. https://doi.org/10.1016/0301-6226(85)90075-2.

    Article  Google Scholar 

  2. Stephens SM, Moley KH. Follicular origins of modern reproductive endocrinology. Am J Physiol Endocrinol Metab. 2009;297(6):E1235–6. https://doi.org/10.1152/ajpendo.00575.2009.

    Article  CAS  PubMed  Google Scholar 

  3. Pan Z, Zhang J, Lin F, Ma X, Wang X, Liu H. Expression profiles of key candidate genes involved in steroidogenesis during follicular atresia in the pig ovary. Mol Biol Rep. 2012;39(12):10823–32. https://doi.org/10.1007/s11033-012-1976-2.

    Article  CAS  PubMed  Google Scholar 

  4. Chu MX, Liu ZH, Jiao CL, He YQ, Fang L, Ye SC, et al. Mutations in BMPR-IB and BMP-15 genes are associated with litter size in small tailed Han sheep (Ovis aries). J Anim Sci. 2007;85(3):598–603. https://doi.org/10.2527/jas.2006-324.

    Article  CAS  PubMed  Google Scholar 

  5. Pan ZY, Di R, Tang QQ, Jin HH, Chu MX, Huang DW, et al. Tissue-specific mRNA expression profiles of GDF9, BMP15, and BMPR1B genes in prolific and non-prolific goat breeds. Czech J Anim Sci. 2015;60:452–8. https://doi.org/10.17221/8525-CJAS.

    Article  CAS  Google Scholar 

  6. Zi XD, Mu XK, Wang Y. Variation in sequences and mRNA expression levels of growth hormone (GH), insulin-like growth factor I (IGF-I) and II (IGF-II) genes between prolific Lezhi black goat and non-prolific Tibetan goat (Capra hircus). Gen Comp Endocrinol. 2013;187:1–5. https://doi.org/10.1016/j.ygcen.2013.03.023.

    Article  CAS  PubMed  Google Scholar 

  7. Zhang C, Liu Y, Huang K, Zeng W, Xu D, Wen Q, et al. The association of two single nucleotide polymorphisms (SNPs) in growth hormone (GH) gene with litter size and superovulation response in goat-breeds. Genet Mol Biol. 2011;34(1):49–55. https://doi.org/10.1590/S1415-47572010005000110.

    Article  PubMed  PubMed Central  Google Scholar 

  8. Ling YH, Quan Q, Xiang H, Zhu L, Chu MX, Zhang XR, et al. Expression profiles of differentially expressed genes affecting fecundity in goat ovarian tissues. Genet Mol Res. 2015;14(4):18743–52. https://doi.org/10.4238/2015.December.28.23.

    Article  CAS  PubMed  Google Scholar 

  9. Lin F, Li R, Pan ZX, Zhou B, Yu DB, Wang XG, et al. miR-26b promotes granulosa cell apoptosis by targeting ATM during follicular atresia in porcine ovary. PLoS ONE. 2012;7:e38640.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Xiao G, Xia C, Yang J, Liu J, Du H, Kang X, et al. MiR-133b regulates the expression of the actin protein TAGLN2 during oocyte growth and maturation: a potential target for infertility therapy. PLoS One. 2014;9(6):e100751. https://doi.org/10.1371/journal.pone.0100751.

    Article  PubMed  PubMed Central  Google Scholar 

  11. Li D, Xu D, Xu Y, Chen L, Li C, Dai X, et al. MicroRNA-141-3p targets DAPK1 and inhibits apoptosis in rat ovarian granulosa cells. Cell Biochem Funct. 2017;35(4):197–201. https://doi.org/10.1002/cbf.3248.

    Article  CAS  PubMed  Google Scholar 

  12. Sheng N, Xu YZ, Xi QH, Jiang HY, Wang CY, Zhang Y, et al. Overexpression of KIF2A is suppressed by miR-206 and associated with poor prognosis in ovarian cancer. Cell Physiol Biochem. 2018;50(3):810–22. https://doi.org/10.1159/000494467.

    Article  CAS  PubMed  Google Scholar 

  13. Hwang JH, An SM, Yu GE, Park DH, Kang DG, Kim TW, et al. Association of single-nucleotide polymorphisms in NAT9 and MAP 3K3 genes with litter size traits in Berkshire pigs. Arch Anim Breed. 2018;61(4):379–86. https://doi.org/10.5194/aab-61-379-2018.

    Article  PubMed  PubMed Central  Google Scholar 

  14. Menon B, Gulappa T, Menon KM. Molecular regulation of LHCGR expression by miR-122 during follicle growth in the rat ovary. Mol Cell Endocrinol. 2017;442:81–9. https://doi.org/10.1016/j.mce.2016.12.002.

    Article  CAS  PubMed  Google Scholar 

  15. Bagnicka E, Wallin E, LUkaszewicz M, Tormod A. Heritability for reproduction traits in polish and Norwegian populations of dairy goat. Small Rumin Res. 2007;68(3):256–62. https://doi.org/10.1016/j.smallrumres.2005.09.033.

    Article  Google Scholar 

  16. Zi XD, Lu JY, Ma L. Identification and comparative analysis of the ovarian microRNAs of prolific and non-prolific goats during the follicular phase using high-throughput sequencing. Sci Rep. 2017;7(1):1921. https://doi.org/10.1038/s41598-017-02225-x.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Ling YH, Ren CH, Guo XF, Xu LN, Huang YF, Luo JC, et al. Identification and characterization of microRNAs in the ovaries of multiple and uniparous goats (Capra hircus) during follicular phase. BMC Genomics. 2014;15(1):339. https://doi.org/10.1186/1471-2164-15-339.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. An X, Song Y, Hou J, Li G, Zhao H, Wang J, et al. Identification and profiling of microRNAs in the ovaries of polytocous and monotocous goats during estrus. Theriogenology. 2016;85(4):769–80. https://doi.org/10.1016/j.theriogenology.2015.09.056.

    Article  CAS  PubMed  Google Scholar 

  19. Zou X, Lu TT, Zhao ZF, Liu GB, Lian ZQ, Guo YQ, et al. Comprehensive analysis of mRNAs and miRNAs in the ovarian follicles of uniparous and multiple goats at estrus phase. BMC Genomics. 2020;21(1):267. https://doi.org/10.1186/s12864-020-6671-4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Bakhshalizadeh S, Amidi F, Shirazi R, Shabani NM. Vitamin D3 regulates steroidogenesis in granulosa cells through AMP-activated protein kinase (AMPK) activation in a mouse model of polycystic ovary syndrome. Cell Biochem Funct. 2018;36(4):183–93. https://doi.org/10.1002/cbf.3330.

    Article  CAS  PubMed  Google Scholar 

  21. Ahn H, Kim KW, Kim HJ, Cho S, Kim H. Differential evolution between monotocous and polytocous species. Asian-Australas J Anim Sci. 2014;27(4):464–70. https://doi.org/10.5713/ajas.2013.13696.

    Article  PubMed  PubMed Central  Google Scholar 

  22. La YF, Tang JS, Di R, Wang XY, Liu QY, Zhang LP, et al. Differential expression of circular RNAs in polytocous and monotocous uterus during the reproductive cycle of sheep. Animals (Basel). 2019;9:797.

    Article  Google Scholar 

  23. Messinis IE, Templeton AA. The importance of follicle-stimulating hormone increase for folliculogenesis. Hum Reprod. 1990;5(2):153–6. https://doi.org/10.1093/oxfordjournals.humrep.a137060.

    Article  CAS  PubMed  Google Scholar 

  24. Ganguly A, Meur SK, Ganguly I. Changes in circulatory FSH of Barbari goats following treatment with high molecular weight inhibin isolated from buffalo follicular fluid. Res Vet Sci. 2013;95(2):374–80. https://doi.org/10.1016/j.rvsc.2013.03.013.

    Article  CAS  PubMed  Google Scholar 

  25. Samuel G, Dessie SW, Ijaz A, Sudeep S, Md Munir H, Michael H, et al. MicroRNA expression profile in bovine granulosa cells of preovulatory dominant and subordinate follicles during the late follicular phase of the estrous cycle. PLoS One. 2015;10(5):e0125912. https://doi.org/10.1371/journal.pone.0125912.

    Article  CAS  Google Scholar 

  26. Jurczak A, Szkup M, Grzywacz A, Safranow K, Grochans E. The relationship between AMH and AMHR2 polymorphisms and the follicular phase in late reproductive stage women. Int J Environ Res Public Health. 2016;13(2):185. https://doi.org/10.3390/ijerph13020185.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Boyer A, Lapointe É, Zheng XF, Cowan RG, Li HG, Quirk SM, et al. WNT4 is required for normal ovarian follicle development and female fertility. FASEB J. 2010;24(8):3010–25. https://doi.org/10.1096/fj.09-145789.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Liu QY, Jiang XX, Tian HN, Guo HL, Guo H, Guo Y. Long non-coding RNA OIP5-AS1 plays an oncogenic role in ovarian cancer through targeting miR-324-3p/NFIB axis. Eur Rev Med Pharmacol Sci. 2020;24(13):7266–75. https://doi.org/10.26355/eurrev_202007_21881.

    Article  PubMed  Google Scholar 

  29. Glazov EA, Kongsuwan K, Assavalapsakul W, Horwood PF, Mitter N, Mahony TJ. Repertoire of bovine miRNA and miRNA-like small regulatory RNAs expressed upon viral infection. PLoS One. 2009;4(7):e6349. https://doi.org/10.1371/journal.pone.0006349.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Lv SY, Zhang GS, Xie LX, Yan ZY, Wang Q, Li YQ, et al. High TXLNA expression predicts favourable outcome for pancreatic adenocarcinoma patients. Biomed Res Int. 2020;2020:2585862.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Agulnik A, Gossett J, Carrillo AK, Kang GL, Ray MR. Abnormal vital signs predict critical deterioration in hospitalized pediatric hematology-oncology and post-hematopoietic cell transplant patients. Front Oncol. 2020;10:354. https://doi.org/10.3389/fonc.2020.00354.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Hansen VL, Schilkey FD, Miller RD. Transcriptomic changes associated with pregnancy in a marsupial, the gray short-tailed opossum monodelphis domestica. PLoS One. 2016;11(9):e0161608. https://doi.org/10.1371/journal.pone.0161608.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Hyodo T, Ito S, Hasegawa H, Asano E, Maeda M, Urano T, et al. Misshapen-like kinase 1 (MINK1) is a novel component of Striatin-interacting phosphatase and kinase (STRIPAK) and is required for the completion of cytokinesis. J Biol Chem. 2012;287(30):25019–29. https://doi.org/10.1074/jbc.M112.372342.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Sah N, Lee Kuehu D, Khadka VS, Deng YP, Peplowska K, Jha R, et al. RNA sequencing-based analysis of the laying hen uterus revealed the novel genes and biological pathways involved in the eggshell biomineralization. Sci Rep. 2018;8(1):16853. https://doi.org/10.1038/s41598-018-35203-y.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Liang GM, Yan JY, Guo J, Tang ZL. Identification of ovarian circular RNAs and differential expression analysis between MeiShan and Large White pigs. Animals (Basel). 2020;10:1114.

    Article  Google Scholar 

  36. Foxcroft GR, Hunter MG. Basic physiology of follicular maturation in the pig. J Reprod Fertil Suppl. 1985;33:1–19.

    CAS  PubMed  Google Scholar 

  37. Clutter AC. Genetic selection for lifetime reproductive performance. Society of Reproduction and Fertility supplement. 2009;66:293–302.

  38. Vanorny DA, Prasasya RD, Chalpe AJ, Kilen SM, Mayo KE. Notch signaling regulates ovarian follicle formation and coordinates follicular growth. Mol Endocrinol. 2014;28(4):499–511. https://doi.org/10.1210/me.2013-1288.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Zhu B, Pardeshi L, Chen YY, Ge W. Transcriptomic analysis for differentially expressed genes in ovarian follicle activation in the zebrafish. Front Endocrinol (Lausanne). 2018;9:593.

    Article  Google Scholar 

  40. Filicori M, Santoro N, Merriam GR, Crowley WFJ. Characterization of the physiological pattern of episodic gonadotropin secretion throughout the human menstrual cycle. J Clin Endocrinol Metab. 1986;62(6):1136–44. https://doi.org/10.1210/jcem-62-6-1136.

    Article  CAS  PubMed  Google Scholar 

  41. Sarkar DK, Chiappa SA, Fink G, Sherwood NM. Gonadotropin-releasing hormone surge in pro-oestrous rats. Nature. 1976;264(5585):461–3. https://doi.org/10.1038/264461a0.

    Article  CAS  PubMed  Google Scholar 

  42. Herrán Y, Gutiérrez-Caballero C, Sánchez-Martín M, Hernández T, Viera A, Barbero JL, et al. The cohesin subunit RAD21L functions in meiotic synapsis and exhibits sexual dimorphism in fertility. EMBO J. 2011;30(15):3091–105. https://doi.org/10.1038/emboj.2011.222.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Hodges CA, Revenkova E, Jessberger R, Hassold TJ, Hunt PA. SMC1beta-deficient female mice provide evidence that cohesins are a missing link in age-related nondisjunction. Nat Genet. 2005;37(12):1351–5. https://doi.org/10.1038/ng1672.

    Article  CAS  PubMed  Google Scholar 

  44. Murdoch B, Owen N, Stevense M, Smith H, Nagaoka S, Hassold T, et al. Altered cohesin gene dosage affects mammalian meiotic chromosome structure and behavior. PLoS Genet. 2013;9(2):e1003241. https://doi.org/10.1371/journal.pgen.1003241.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Jaffe LA, Egbert JR. Regulation of mammalian oocyte meiosis by intercellular communication within the ovarian follicle. Annu Rev Physiol. 2017;79(1):237–60. https://doi.org/10.1146/annurev-physiol-022516-034102.

    Article  CAS  PubMed  Google Scholar 

  46. Brown MC, Turner CE. Paxillin: adapting to change. Physiol Rev. 2004;84(4):1315–39. https://doi.org/10.1152/physrev.00002.2004.

    Article  CAS  PubMed  Google Scholar 

  47. Berrier AL, Yamada KM. Cell-matrix adhesion. J Cell Physiol. 2007;213(3):565–73. https://doi.org/10.1002/jcp.21237.

    Article  CAS  PubMed  Google Scholar 

  48. Viveiros MM. The Ovary || Regulation of Mammalian Oocyte Maturation. 2019;165–80.

  49. Bliss SP, Navratil AM, Xie J, Roberson MS. GnRH signaling, the gonadotrope and endocrine control of fertility. Front Neuroendocrinol. 2010;31(3):322–40. https://doi.org/10.1016/j.yfrne.2010.04.002.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Doherty JA, Rossing MA, Cushing-Haugen KL, Chen C, Van Den Berg DJ, Wu AH, et al. ESR1/SYNE1 polymorphism and invasive epithelial ovarian cancer risk: an ovarian cancer association consortium study. Cancer Epidemiol Biomark Prev. 2010;19(1):245–50. https://doi.org/10.1158/1055-9965.EPI-09-0729.

    Article  CAS  Google Scholar 

  51. Ekumi KM, Paculova H, Lenasi T, Pospichalova V, Bösken CA, Rybarikova J, et al. Ovarian carcinoma CDK12 mutations misregulate expression of DNA repair genes via deficient formation and function of the Cdk12/CycK complex. Nucleic Acids Res. 2015;43(5):2575–89. https://doi.org/10.1093/nar/gkv101.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Kawasaki K, Kuge O, Yamakawa Y, Nishijima M. Purification of phosphatidylglycerophosphate synthase from Chinese hamster ovary cells. Biochem J. 2001;354(1):9–15. https://doi.org/10.1042/bj3540009.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Cai Y. Revisiting old vaginal topics: conversion of the Müllerian vagina and origin of the “sinus” vagina. Int J Dev Biol. 2009;53(7):925–34. https://doi.org/10.1387/ijdb.082846yc.

    Article  CAS  PubMed  Google Scholar 

  54. Chen M, Hsu I, Wolfe A, Radovick S, Huang K, Yu S, et al. Defects of prostate development and reproductive system in the estrogen receptor-alpha null male mice. Endocrinology. 2009;150(1):251–9. https://doi.org/10.1210/en.2008-0044.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Zhang Z, Tang J, Di R, Liu Q, Wang X, Gan S, et al. Integrated hypothalamic transcriptome profiling reveals the reproductive roles of mRNAs and miRNAs in sheep. Front Genet. 2020;10:1296. https://doi.org/10.3389/fgene.2019.01296.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Trapnell C, Salzberg SL. How to map billions of short reads onto genomes. Nat Biotechnol. 2009;27(5):455–7. https://doi.org/10.1038/nbt0509-455.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Friedlander MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012;40(1):37–52. https://doi.org/10.1093/nar/gkr688.

    Article  CAS  PubMed  Google Scholar 

  58. Wen M, Shen Y, Shi S, Tang T. miREvo: an integrative microRNA evolutionary analysis platform for next-generation sequencing experiments. BMC Bioinformatics. 2012;13(1):140. https://doi.org/10.1186/1471-2105-13-140.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11(2):R14. https://doi.org/10.1186/gb-2010-11-2-r14.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Mao X, Cai T, Olyarchuk JG, Wei L. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics. 2005;21(19):3787–93. https://doi.org/10.1093/bioinformatics/bti430.

    Article  CAS  PubMed  Google Scholar 

  61. Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30. https://doi.org/10.1093/nar/28.1.27.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Enright AJ, John B, Gaul U, Tuschl T, Sander C, Marks DS. MicroRNA targets in Drosophila. Genome Biol. 2003;5(1):R1. https://doi.org/10.1186/gb-2003-5-1-r1.

    Article  PubMed  PubMed Central  Google Scholar 

  63. Wu HJ, Ma YK, Chen T, Wang M, Wang XJ. PsRobot: a web based plant small RNA meta-analysis toolbox. Nucleic Acids Res. 2012;40(W1):W22–8. https://doi.org/10.1093/nar/gks554.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Goes CP, Vieceli FM, De La Cruz SM, Simões-Costa M, Yan CYI. Scratch2, a snail superfamily member, is regulated by miR-125b. Front Cell Dev Biol. 2020;8:769. https://doi.org/10.3389/fcell.2020.00769.

    Article  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgments

Not applicable.

Funding

This work was financially supported by the Agricultural Science and Technology Innovation Program of China (CAAS-ZDRW202106 and ASTIP-IAS13), the China Agriculture Research System of MOF and MARA (CARS-38), the Science and Technology Innovation Talents Program of Yunnan Province (2018HC017), and the Basic Research Foundation Key Project of Yunnan Province (202001AS070002).

Author information

Authors and Affiliations

  1. Key Laboratory of Animal Genetics, Breeding and Reproduction of Ministry of Agriculture and Rural Affairs, Institute of Animal Science, Chinese Academy of Agricultural Sciences, No. 2 Yuanmingyuan West Rd, Beijing, 100193, China

    Yufang Liu, Zuyang Zhou, Xiaoyun He, Lin Tao & Mingxing Chu

  2. College of Life Sciences and Food Engineering, Hebei University of Engineering, Handan, 056021, China

    Yufang Liu & Zuyang Zhou

  3. Yunnan Animal Science and Veterinary Institute, Kunming, 650224, China

    Yanting Jiang, Rong Lan & Qionghua Hong

Authors
  1. Yufang Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Zuyang Zhou
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Xiaoyun He
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Lin Tao
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Yanting Jiang
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Rong Lan
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Qionghua Hong
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Mingxing Chu
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

Conceptualization, YFL and MXC; Methodology, YFL, ZYZ, LT and MXC; Validation: YFL and ZYZ; Formal Analysis: YFL and MXC; Investigation: YFL, LT, ZYZ, XYH, YTJ, RL, QHH and MXC; Resources, YFL, LT, ZYZ and MXC; Data Curation, YFL, LT, ZYZ and MXC; Writing-Original Draft Preparation, YFL and MXC; Supervision, MXC; Project Administration, YFL, QHH and MXC; Funding Acquisition, QHH and MXC. All authors have read and agreed to the published version of the manuscript.

Corresponding authors

Correspondence to Qionghua Hong or Mingxing Chu.

Ethics declarations

Ethics approval and consent to participate

All procedures performed in studies involving animals were in accordance with the ethical standards of the institutional and national research committee. Ethical approval was provided by the animal ethics committee of the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences (No. IAS2019–63). The manuscript adheres to the ARRIVE guidelines for the reporting of animal experiments, and the study was carried out in compliance with the ARRIVE guidelines.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests or economic interests in any of the companies cited here.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

The information of RNA-seq data.

Additional file 2: Table S2.

All of the DE mRNAs in the comparison.

Additional file 3: Table S3.

The information of DE mRNAs in the comparison.

Additional file 4: Table S4.

The information of miRNA sequencing data.

Additional file 5: Table S5.

All of the expression miRNAs in the comparison.

Additional file 6: Table S6.

Primer sequences used for qPCR in this study.

Additional file 7: Fig. S1.

The sequence of miRNA-mRNA pair for Dual Luciferase Reporter Assay.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Liu, Y., Zhou, Z., He, X. et al. Integrated analyses of miRNA-mRNA expression profiles of ovaries reveal the crucial interaction networks that regulate the prolificacy of goats in the follicular phase. BMC Genomics 22, 812 (2021). https://doi.org/10.1186/s12864-021-08156-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-021-08156-2

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us