Integrated Transcriptome Analysis Reveals the Crucial mRNAs and miRNAs Related to Fecundity in the Hypothalamus of Yunshang Black Goats during the Luteal Phase

Simple Summary The hypothalamus plays an important role in goat reproduction. In this research, RNA sequencing technology was used to detect and identify microRNAs (miRNAs) related to goat fertility in the goat hypothalamus. The results showed that several key miRNAs may affect the expression of their target genes in different ways with the high-fecundity goats during the luteal phase (LP-HF) vs. low-fecundity goats during the luteal phase (LP-LF), thereby affecting the key reproductive processes. These results reveal the mechanism on the reproductive performance of different litter sizes during the luteal phase of goats, from the perspective of the hypothalamus. Abstract A normal estrus cycle is essential for the breeding of goats, and the luteal phase accounts for most of the estrus cycle. The corpus luteum (CL) formed during the luteal phase is a transient endocrine gland that is crucial for the reproductive cycle and pregnancy maintenance, and is controlled by many regulatory factors. However, the molecular mechanism of the hypothalamus effect on the reproductive performance of different litter sizes during the luteal phase of goats has not been elucidated. In this study, RNA-sequencing was used to analyze the mRNA and miRNA expression profiles of the hypothalamic tissues with the high-fecundity goats during the luteal phase (LP-HF) and low-fecundity goats during the luteal phase (LP-LF). The RNA-seq results found that there were 1963 differentially expressed genes (DEGs) (890 up-regulated and 1073 down-regulated). The miRNA-seq identified 57 differentially expressed miRNAs (DEMs), including 11 up-regulated and 46 down-regulated, of which 199 DEGs were predicted to be potential target genes of DEMs. Meanwhile, the functional enrichment analysis identified several mRNA-miRNA pairs involved in the regulation of the hypothalamic activity, such as the common target gene MEA1 of novel-miR-972, novel-miR-125 and novel-miR-403, which can play a certain role as a related gene of the reproductive development in the hypothalamic–pituitary–gonadal (HPG) axis and its regulated network, by regulating the androgen secretion. While another target gene ADIPOR2 of the novel-miR-403, is distributed in the hypothalamus and affects the reproductive system through a central role on the HPG axis and a peripheral role in the gonadal tissue. An annotation analysis of the DE miRNA-mRNA pairs identified targets related to biological processes, such as anion binding (GO:0043168) and small molecule binding (GO: 0036094). Subsequently, the KEGG(Kyoto Encyclopedia of Genes and Genomes) pathways were performed to analyze the miRNA-mRNA pairs with negatively correlated miRNAs. We found that the GnRH signaling pathway (ko04912), the estrogen signaling pathway (ko04915), the Fc gamma R-mediated phagocytosis (ko04666), and the IL-17 signaling pathway (ko04657), etc., were directly and indirectly associated with the reproductive process. These targeting interactions may be closely related to the reproductive performance of goats. The results of this study provide a reference for further research on the molecular regulation mechanism for the high fertility in goats.


Introduction
The traits of mammalian fecundity are determined by complex regulatory factors, including genetic material, nutritional level and feeding environment. The high performance of reproduction depends on the coordinated synthesis and release of hormones in the gonadal axis [1]. The hypothalamic-pituitary-gonadal (HPG) axis is a branch of the endocrine system, which controls the secretion of sex hormones and is the basis of the endocrine control for the reproduction in mammals [2]. The hypothalamus, located at the apex of the HPG axis, contains neurons that secrete a gonadotropin-releasing hormone (GnRH), which can produce GnRH signals and regulate the secretion of the downstream hormones, including the follicle-stimulating hormone (FSH) and the luteinizing hormone (LH) [3]. It can be said that reproduction is controlled by the hypothalamus, which is the key area in the brain to initiate reproductive activities. Alterations in the hypothalamus function, also affects the reproductive activity, including the follicular development, ovulation, and oviposition [4]. The corpus luteum is formed from a ruptured follicle during the luteal phase after ovulation [5]. The formation of the corpus luteum is crucial for the steroid biosynthesis and maintenance of early pregnancy [6], while sex steroid hormones, such as estrogen, progesterone, and androgen, play a key role in gender differentiation, reproductive function, and the sexual behavior of mammals [7].
Furthermore, reproduction is a complex process that is regulated by many genes, transcription factors, and regulatory non-coding RNAs. Since the Booroola fecundity gene (FecB) was identified in Booroola sheep [8], a series of studies have been carried out on the major genes controlling the sheep and goat follicle development, estrus, and ovulation by using molecular biology techniques. It has been found that deiodinase type 2 (DIO2) and deiodinase type 3 (DIO3) are involved in the regulation of the seasonal estrus in goats [9]. There were also other genes, retinol-binding protein 4 (RBP4), GnRH and its receptor GnRHR, which are related to the reproduction in goats. In addition, kisspeptin is a polypeptide encoded by the Kiss-1 gene, which is highly expressed in two major neuronal populations of the hypothalamus. It is a recognized reproductive hormone coordinator, and the kisspeptin/GnRH pathway is very important for the reproductive endocrine system [10]. Moreover, different molecules that play an epigenetic role, such as miRNA, were also investigated in this context. MicroRNAs (miRNAs) can bind to complementary mRNAs and influence the post-transcriptional level of genes by targeting or inhibiting the expression of the transcripts [11]. There are many miRNAs acting as mediators to further regulate the activity of the HPG axis [12]. So far, miR-7a, let-7a, and miR-30b had been found to regulate the reproductive activity. The miR-7a can be co-expressed with agouti-related peptide (AgRP) in AgRP neurons [13] and plays an important role in reproduction by acting on the mTOR pathway within the AgRP neurons [14]. One study reported that the deletion of miR-7a led to hypogonadotropic hypogonadism in mouse models [15]. Let-7a and its related protein lin-28b are key components of the neuroendocrine mechanisms that control the timing of the onset of puberty, regulating reproduction by acting on kisspeptin and/or GnRH neurons [16,17]. The miR-30b is expressed in kisspeptin ARC(arcuate nucleus) neurons and regulates the kisspeptin expression by binding to the makorin ring finger protein 3 (MKRN3) gene, thereby inhibiting its activity [18].
In the development of human civilization, goats played important roles in human agricultural and economic development by providing meat, milk, skin, hair, and other products. The reproductive trait is one of the main economic traits of the goat, and the reproductive performance is directly related to the economic benefits of the goat industry. Therefore, it is important to elucidate the interaction mechanism of mRNAs and miRNAs during goat reproduction, to improve the reproduction performance of goats. However, the molecular mechanisms involved in the regulation of the high-fecundity traits in the hypothalamus tissue during the luteal phase of goats, are still poorly understood. The Yunshang black goat is a new black goat breed bred for meat, with a high fertility in China, which has the characteristics of a fast growth, good meat production performance, and excellent meat quality. In this study, the high-fecundity goat breed was used as the research object, and the RNA-seq technology was used to analyze the expression profiles of miRNAs and mRNAs in the hypothalamus during the luteal phase of high-and lowfecundity goats. The key molecular mechanisms and molecular networks of the breed's fecundity were revealed. The purpose of this study was to help our understanding of the molecular mechanism of the prolific traits. The results of this study provided a useful resource for further research on miRNAs and mRNAs in high-fecundity goats and their potential relationships and functional roles in the reproductive regulation.

Ethics Statement
The experimental animals in this study were in compliance with animal welfare standards. The animal study under a permit (No. IAS2019-63) was reviewed and approved by the Institutional Animal Care and Use Ethics Committee of the Institute of Animal Sciences and Chinese Academy of Agricultural Sciences (IAS-CAAS) (Beijing, China) for all of the experimental procedures mentioned. All protocols complied with animal welfare guidelines and minimized animal suffering.

Goat Selection and Sample Collection
Ten Yunshang black nanny goats, aged 3-5 years (The average age of each lowfecundity and high-fecundity group was 4 years), with no significant difference in weight, height, and physical condition, were selected as the research objects from the goat core breeding farm in Kunming, Yunnan Province. In addition, all goats were raised under the same conditions with free access to water and feed. All of the experimental goats were treated for synchronous estrus, and a CIDR (controlled internal drug-release device) suppository (30 mg flurogestone acetate) was placed in the vagina of the goats for 16 days and then removed. On the 11th day after the CIDR removal, it was the luteal phase, and slaughter and hypothalamus sampling were performed. The specific process is shown in Figure 1. The goats were divided into high fecundity goats at luteal phase (LP-HF) group (n = 5, average litter size: 3.40 ± 0.42) and low fecundity goats at luteal phase (LP-LF) group (n = 5, average litter size: 1.80 ± 0.27) according to their kidding number. The collected hypothalamic tissues were washed with normal saline, placed in the cryopreservation tubes, and immediately placed in liquid nitrogen for short-term storage, and then brought back to the laboratory and placed in a −80 • C refrigerator for long-term storage and the subsequent experiments.

Total RNA Extraction
According to the manufacturer's instructions, the total RNA was extracted from the hypothalamus tissue samples of 10 Yunshang black goats with TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA) and used for the RNA sequencing [19]. High-quality RNA is the basis for the successful sequencing. We used the following methods to detect the samples, and the library can be constructed only after the test results meet the requirements: (1) The Nanodrop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA) was used to detect the purity (OD 260/280), concentration, and nucleic acid absorption peak of RNA; (2) The Agilent 2100 RNA Nano 6000 Assay Kit (Agilent Technologies, Palo Alto, CA, USA) was used to accurately detect the integrity of RNA; (3) We used 1% agarose gel electrophoresis to detect whether the RNA sample was degraded or contaminated with genomic DNA.

Library Preparation and Sequencing
To begin, 3 µg of total RNA was taken from each hypothalamic sample for the mRNA library construction. The library construction was performed, according to the instructions for the NEB Next ® Ultra TM Directional RNA Library Prep Kit and Illumina ® (NEB, Ipswich, MA, USA). Then, the library was sequenced on the Illumina Novaseq6000 platform, to generate 150 bp paired-end (PE150) sequence reads. The sequence reads in fastq format were first filtered by the SOAPnuke (v2.1.0) [20] and HISAT2 (v2.1.0) [21] was used to map the clean reads to the reference genome (GCF_001704415.1). Subsequently, we used String Tie (v1.3.5) [22] to assemble the transcripts, and performed the fragments per kilobase per base (FPKM) transformation to obtain the expression levels of the transcript.
Similarly, 3 µg of total RNA was used for miRNA library construction, the library was sequenced on the Illumina Hiseq 2500 platform, resulting in 50 bp single-end (SE50) reads. Using the clean reads of 18-35 nt in length for the subsequent analyses, and small RNA (sRNA) was aligned to the reference genome (GCF_001704415.1) by bowtie (v1.0.1) [23] to identify the known miRNAs. In this process, the miRbase (v22.0) database [24] was used as a reference. Referring to mirEvo (v2.0.0.5) and miRdeep (v2.0.0.5) [25], novel miRNAs were predicted, based on the signature hairpin structures of the miRNA precursors. The expression amount of the known and novel miRNAs in each sample was counted. The expression levels were normalized with TPM (Transcripts per million).

Analysis of the Differential Expression mRNAs and miRNAs
We used DESeq2 (v3.18.1) [26] to analyze the differentially expressed mRNAs in LP-HF and LP-LF. The default settings for the screening threshold are |log 2 FC| > 1 and q-value < 0.05. Meanwhile, based on the TPM expression levels of miRNAs, we also used DESeq2 to identify the DEMs of LP-HF and LP-LF, with |log 2 FC| > 1, p-value < 0.05 as the differential expression.
We input the differential gene expression values of the high-and low-fecundity comparison groups, and use the log10 as the base to process the data. We used the pheatmap package of the R software (v 4.2.1) (https://www.r-project.org, accessed on 1 May 2022) to read the data for drawing the heat map, in order to visually present the global expression changes and clustering relationships of multiple genes in multiple samples.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Analyses
The functional annotation and pathway enrichment analysis of the DEGs and the predicted target genes of the DEMs were performed using GO [27] (Gene Ontology, http: //www.geneontology.org (accessed on 3 November 2022), of GOseq (v2.12) software and KEGG [28] (Kyoto Encyclopedia of Genes and Genomes) of KOBAS (v2.0) software, to visualize the data. The GO analysis mainly includes a cellular component (CC), molecular function (MF) and biological process (BP). The GO items and KEGG pathways significantly enriched by the DEGs, were calculated using hypergeometric distribution, and the GO items and KEGG pathways with a q-value < 0.05 (If q-value < 0.05 was used to enrich the very few GO items and KEGG pathways, then p-value < 0.05 was used for the enrichment) were considered significantly enriched.

Integrated miRNA-mRNA Co-Expression Network Analysis
The MiRanda v3.3a and qTar (https://github.com/zhuqianhua/qTar.git, accessed on 28 April 2020) were used to predict the target genes of the miRNAs. Then, we analyzed the interactions between the miRNAs and mRNAs. Based on the miRNAs function, the mRNAs that were negatively related to the miRNAs were screened out, in order to accurately identify the key association with the reproductive DEMs and DEGs, and the miRNA-mRNA interaction networks were built by using Cytoscape (v3.8.2, http://www.cytoscape.org/, accessed on 1 May 2016).
For the RT-qPCR analysis of the mRNAs, the RT-qPCR with the SYBR Green qPCR Mix (TaKaRa, Dalian, China) was conducted with a Roche Light Cycler ® 480 II system (Roche Applied Science, Mannheim, Germany). For the miRNAs, the RT-qPCR was conducted by miRcute Plus miRNA qPCR Kit (TIANGEN, Beijing, China).
For the specific operation steps, please refer to the supplier's manual. In addition, RPL19 (for mRNA) and U6 small nuclear RNA (for miRNA) were utilized as the reference gene-miRNA to calculate the relative expression level with the method of 2 −∆∆Ct .

Statistical Analysis
We performed the statistical analysis of the RT-qPCR results and graphs using Graph-Pad Prism (v 5.0) software. The statistical data were tested by performing paired t-tests. The results are presented as the means of three replicates, then the gene expression levels were determined and compared with the RNA-seq data.

The cDNA Library Sequencing and mRNA Transcriptome Analysis
Ten cDNA libraries were sequenced by Illumina, and the RNA-seq for the mRNA generated approximately 1111 million clean reads after the data filtering (Table 2). More than 96.15% of the clean reads were mapped to the goat reference genome. The GC(guaninecytosine content) content ranged from 41.3 to 49.0%, and the scores of Q20 and Q30 were above 97.5% and 92.8%, respectively. This indicated that the reliability and quality of the sequencing data were sufficient for further analysis. A total of 72,220 mRNAs were identified after mapping the goat genome (Supplementary Table S1). Regarding the expression levels of the mRNAs, our results showed that genes with a high RNA-seq expression (genes with FPKM > 60) accounted for about 1.19% (Figure 2A), and the gene expression levels (FPKM) of the 10 hypothalamus tissues were similar ( Figure 2B).

Small RNA Library Sequencing and miRNA Transcriptome Analysis
A total of 249.7 million raw reads were obtained by the miRNA transcriptome sequencing in the hypothalamus from the luteal phase of the high and low fecundity Yunshang black goats, and 215.5 million clean reads remained after filtering. The clean reads were compared with the goat reference genome (Table 3), the alignment rate was about 97.52%, the length of the filtered sequences was between 18-35 nt ( Figure 3A), the GC content was between 47.71% and 51.79%, and the average Q20 content was 99.22% and Q30 content was 97.38%, indicating that the quality of the data generated by sequencing was relatively high. In addition, a variety of non-coding RNAs were also identified (Supplementary Table S4), and most of the small RNAs were known miRNAs, accounting for about 73.29% ( Figure 3B). Finally, a total of 1837 miRNAs were identified in our transcriptome data, including 424 known miRNAs and 1413 novel miRNAs for the subsequent analysis (Supplementary Table S3). Among them, the most abundant miRNAs in the hypothalamus from the luteal phase of Yunshang black goat include miR-9, let-7, and miRNA-200 family members.

Differential Expression and Functional Enrichment Analysis of the mRNAs in the LP-HF vs. LP-LF
In this study, the differential expression analysis of the mRNA of the high-and lowfecundity Yunshang black goats was carried out, with |log 2 FC| > 1 and q value < 0.05 as the screening criteria. In the LP-HF vs. LP-LF, a total of 1963 DEGs were screened. Among them, 890 DEGs were up-regulated and 1073 DEGs were down-regulated (Supplementary  Table S2, Figure 4A). The heatmap of the DEGs in the LP-HF vs. the LP-LF comparison group showed that the DEGs had different expression patterns between the LP-HF and LP-LF ( Figure 4C).  Table S7, Figure 5A). The KEGG analysis showed that the DEGs were significantly enriched in 11 KEGG pathways, of which the hedgehog signaling pathway (ko04340) was the most enriched (Supplementary Table S8, Figure 6A). Moreover, we also discovered several pathways involved in the reproductive biology, such as the GnRH signaling pathway (ko04912), steroid hormone biosynthesis (ko00140) and estrogen signaling pathway (ko04915). Of the DEGs, we found several key genes involved in the reproduction process, such as cytochrome P450 family 19 subfamily a member 1 (CYP19A1), neural cell adhesion molecule 1 (NCAM1), and the fibroblast growth factor receptor (FGFRs) family.

Differential Expression and Functional Enrichment Analysis of the miRNAs in the LP-HF vs. LP-LF
A total of 57 DEMs were screened in the LP-HF vs. LP-LF groups. There were 11 up-regulated DEM expressions and 46 down-regulated DEM expressions (Supplementary  Table S5, Figure 4B). Furthermore, a heatmap of the DEMs in the LP-HF vs. the LP-LF comparison group revealed that the DEMs had different expression patterns between the LP-HF and LP-LF ( Figure 4D). In order to better understand the biological functions of the 57 identified DEMs, we predicted the potential target genes of these miRNAs ( Supplementary Table S6 Table S10, Figure 6B) showed that in the LP-HF vs. the LP-LF comparison group, the most significantly enriched pathways were the Apoptosis (ko04210) and the Notch signaling pathway (ko04330). In addition, the Jak-STAT signaling pathway (ko04630), the Cytokinecytokine receptor interaction (ko04060) and the VEGF signaling pathway (ko04370) also had been enriched. Among them, 30 KEGG pathways were identified as significantly associated with the DEM target genes.

miRNA-mRNA Co-Expression Network Analysis
To further analyze the relationship between the miRNAs and mRNAs, we constructed an interacting co-expression network using the potential target genes of the DEMs and DEGs obtained by the RNA-Seq. The potential target genes of the 57 identified DEMs were intersected with 1963 DEGs, obtained from the RNA-seq, resulting in 199 intersecting genes (Supplementary Table S11, Figure 7A). Cytoscape was then used to construct an interaction network of the miRNA-mRNA pairs (Supplementary Table S12, Figure 7B).
The GO enrichment analysis was performed to reveal the biological process terminology of the enrichment (q-value < 0.05), to gain insight into the function of the miRNA-mRNA pairs with the negatively correlated expression between the LP-HF and LP-LF. The GO terms were shown in Table 4, and the most abundant terms in the LP-HF vs. LP-LF were the anion binding (GO:0043168) and the small molecule binding (GO: 0036094). Among them, the MX2 gene was annotated into almost all GO terms, indicating that this gene has an important regulatory role in the process of this study, and further in-depth research is needed. In addition, we also performed the KEGG pathway analysis of the miRNA-mRNA pairs with negatively correlated miRNAs. The enriched KEGG pathways were shown in Table 5. The two top pathways were the Hedgehog signaling pathway (ko04340) and Ubiquitin mediated proteolysis (ko04120). Interestingly, the GnRH signaling pathway (ko04912), the estrogen signaling pathway (ko04915), Fc gamma R-mediated phagocytosis (ko04666), and the IL-17 signaling pathway (ko04657), were directly and indirectly associated with the reproductive process.

Data Validation
To evaluate the accuracy of the sequencing, we randomly selected the mRNAs and miRNAs from the sequencing data for the RT-qPCR validation. We measured gene expression levels and compared them with the RNA-seq data. The results showed that the RNA-seq data and the RT-qPCR data exhibited similar patterns (consistent up-and down-regulation relationships) (Figure 8), which indicated that the data generated from the RNA-seq were reliable.

Discussion
The hypothalamus is a master regulator of reproduction. The GnRH neurons distributed in the hypothalamic tissue can secrete GnRH, and then GnRH can regulate the secretion of the gonadotropins, such as FSH and LH, by the pituitary gland, causing the secretion of related hormones and affecting the reproductive activities [3]. Strict regulation of the central nervous system and the endocrine system is necessary to achieve reproduction [29]. In this study, the whole transcriptome of the hypothalamus from the luteal phase of Yunshang black goats was established. The expression profiles of the miRNA and mRNA in the hypothalamus were revealed by comparing the reproductive performance of the high and low fecundities, and the miRNA-mRNA network was constructed. The miRNA-mRNA interaction network comprehensively investigated the molecular basis of the candidate genes involved in the goat reproductive traits. A total of 72,220 transcripts and 1837 miRNAs were identified in hypothalamic samples of the goats, including 1963 DEGs and 57 DEMs.
The miRNAs are the most widely distributed regulatory factors in animals. Studies have found that many miRNAs are involved in a variety of biological regulation processes in the hypothalamus, such as immune response [30], osmotic regulation [31], and so on. In this study, miR-9 is the most abundant miRNA in the hypothalamus of Yunshang black goats, during the luteal phase, which was consistent with the results of the follicular phase in the previous [32]. It has been reported that miR-9 is mainly expressed in neural tissues in mammals, such as mice and humans [33], and regulates BAF53a (also known as ACTL6A, actin like 6A) [34], stathmin [35] and other genes involved in the proliferation, migration, and other activities of the nerve cells. Some studies have found that in the hypothalamus of rats, the changes of lin28/let-7 may be involved in the initiation of puberty [17]. In mice, the expression pattern of miRNA-200 was associated with the increase of the GnRH expression in the hypothalamus before sexual maturity [36]. In addition, let-7 and miRNA-200 family members have been found to be highly conserved across species, in sequence and function, and also have a period specificity and fecundity specificity [37], which plays an important role in a variety of cell activities, including the cell cycle regulation and apoptosis [38]. In this study, let-7 and miRNA-200 family members were highly expressed in the hypothalamus of both high-fecundity and low-fecundity groups of Yunshang black goats, during the luteal phase.
In addition, 57 DEMs in the high and low fecundity groups of Yunshang black goats were identified. One of the highly conserved miR-10 family has been shown to be involved in the cell proliferation and apoptosis. Yu et al. found that members of the miR-10 family were highly expressed in goose follicles, during the laying and nesting stages [39]. Tu et al. found that both miR-10a and miR-10b can inhibit the expression of key factors in the brainderived neurotrophic factor (BDNF) and the TGF-β signaling pathways, thereby inhibiting the follicular granule cell proliferation and induction of granulosa cell apoptosis [40]. Wei et al. found that miR-10b could target BDNF and inhibit the activity of goat ovarian granulosa cells by inhibiting the expression of BDNF [41]. However, due to the limited research on the functions of miRNAs in the hypothalamus, the deep functions of the miRNAs screened in this study need to be further studied.
In the functional analysis of the hypothalamic DEGs in the LP-HF vs LP-LF, we found that CYP19A1, NCAM1, and the FGFRs family were involved in the reproduction process.
Aromatase, encoded by the CYP19A1 gene, is an estrogen-synthesizing enzyme that catalyzed many reactions related to steroidogenesis and the conversion of androgens to estrogens. A study has shown that CYP19A1 is involved in the follicle development and follicular atresia [42], and is also important in the sheep gonad development, horse testicular cells, and luteal development [43]. Sun et al. showed that CYP19A1 was expressed in the ovary of Guangxi Ma chicken, during the laying period, and can promote the physiological function of the ovary and negatively regulate the hypothalamus during the process of estrus or ovulation [44]. A previous study on the CYP19A1 gene had also focused on cancerrelated diseases in humans [45]. This study selected Yunshang black goat individuals with different reproductive performances during the luteal phase, and found that CYP19A1 was expressed in the hypothalamus. In addition, the expression level in the hypothalamus of goats with a high fecundity was higher than those with a low fecundity (p < 0.05), and the KEGG pathway enrichment analysis showed that the gene was enriched in the steroid hormone biosynthesis (ko00140) and the ovarian steroidogenesis (ko04913), that the pathways were associated with the reproductive function of the hypothalamus, and that was an important candidate gene.
NCAM1 belongs to the immunoglobulin superfamily of the cell adhesion molecules and is widely present in the nerve cells of the central and peripheral nervous systems [46]. NCAM1 is highly expressed in the hypothalamus, which is similar to the results of this study, especially in the glial cells of the GnRH nerve endings [47]. Furthermore, Rubinek et al. demonstrated that the homologous binding of NCAM1 can induce the GH secretion in fetal pituitary cultures [48]. In addition, studies have shown that NCAM1 can mediate the anti-apoptotic, pro-neurogenesis and neuroprotective biological effects of FGFs and GDNF by binding to FGFR, and the glial cell derived neurotrophic factor (GDNF) receptor GFRα [49,50]. The signal pathway enriched by KEGG was the cell adhesion molecules (ko04514). A study has shown that the intercellular communication mediated by the cell adhesion molecule signaling pathway plays an important role in regulating the expression and secretion of hormones, such as GH [51].
The family of FGFRs is composed of four tyrosine kinase receptors (FGFR1-FGFR4) that bind to the fibroblast growth factors (FGFs) to mediate the cell signaling. The binding of FGFs/FGFR can activate the complex downstream signals, such as PLCγ, MAPK, and PI3K-AKT, and play an important role in the process of the cell proliferation, differentiation, survival and migration [52], and specifically involved in promoting angiogenesis [53], injury repair [54], embryonic development [55], and nerve regeneration [56]. The binding of FGF2 secreted by astrocytes to FGFR1 can promote the differentiation and survival of GnRH neurons, and enhance the synthesis of the GnRH protein [57], thereby participating in the regulation of animal reproduction. FGFR2 is a high-affinity receptor for most FGFs. Studies on animals and humans have shown that FGFR2 is expressed on trophoblast membranes [58]. The inhibition of the FGFR2 expression resulted in the reduced trophoblast formation and delayed the trophoblast growth [59]. Thus, the role of FGFR2 in the placental development and fetal growth is important. Cells that express both FGFR1 and FGFR2 can receive more FGF4 signals and exhibit a higher ERK activity [60]. During development, the FGF-ERK signaling pathway plays an important role in the proliferation, differentiation, and apoptosis.
Overall, the genes mentioned above play important roles in the regulation of animal reproduction. The KEGG enrichment analysis found that these DEGs were significantly enriched in the signal transduction, nervous system, endocrine system, metabolism, and other related pathways. These pathways play crucial roles in the biosynthesis of reproductive hormones. We found that many GO terms related to reproductive biology were also enriched, such as the steroid hormone mediated signaling pathway (GO:0043401), the response to the steroid hormone (GO:0048545) and the cellular response to the steroid hormone stimulus (GO:0071383).
The typical molecular mechanism of the miRNA is to exert its regulatory effect through target genes. An integrative analysis of mRNAs and miRNAs by the RNA-Seq was used to construct regulatory networks to further understand the key genes affecting the reproductive processes.
Among the 57 miRNAs and predicted target genes differentially expressed during the luteal phase, the common target gene of the novel-miR-972, novel-miR-125, and novel-miR-403, was the male-enhanced antigen 1 (MEA1). A study has shown that MEA1 plays a certain role as a gene related to the reproductive development in the HPG axis and its regulatory network, by regulating the androgen secretion [61]. The target gene SDCCAG8 of the novel-miR-972 was found in a previous study at our laboratory [62], indicating that SDCCAG8 was a high-fecundity convergent gene in sheep and goats, and can affect the litter size of goats through the TGFβ pathway. Another target gene predicted by the novel-miR-403 was adiponectin receptor 2 (ADIPOR2). A study has shown that ADIPOR2 is widely distributed in many reproductive organs, including the central nervous system, ovaries, and fallopian tubes [63]. There is evidence that adiponectin affects the reproductive system through the central effects on the HPG axis, as well as the peripheral effects on the gonadal tissue [64]. In addition, adiponectin also affects the secretion of oxytocin, which has a profound effect on the HPG axis and the reproductive function [65]. The novel-miR-1125 targets LINGO-1, a protein encoded by a leucine-rich repeat (LRR), neurite outgrowth inhibitory protein (Nogo), and immunoglobulin domain (Ig) constitutes the central nervous system transmembrane protein. In a variety of animal models, the targeted inhibition of LINGO-1 can promote the neuronal survival, axon regeneration, and the oligodendrocyte differentiation [66].
Meanwhile, miR-1271 was involved in important biological functions, such as neuroregulation. The overexpression of miR-1271 can inhibit the cell proliferation and induce apoptosis [67]. In this study, the miRNA was highly expressed in the hypothalamus, speculating that it could regulate the reproductive performance of goats by affecting the proliferation of nerve cells. A previous study of our subject also demonstrated that miR-1271 targeting TXLNA in ovarian tissue, was a potential regulatory pathway affecting the reproductive performance of goats [68]. The interaction network in this study showed that miR-1271 targets PNN interacting serine and the arginine rich protein (PNISR). PNISR is a part of the multi-protein complex in the nucleus and is involved in the processing of pre-mRNA to affect the cell proliferation and differentiation [69]. In addition, another target gene was MX Dynamin like GTPase 2 (MX2), which has been shown to significantly inhibit the cell proliferation by prolongating the cell cycle [70]. MX2 protein has also been reported as a biomarker for the early gestation diagnosis in buffalos [71]. In conclusion, miR-1271 is worthy of further investigating.
The predicted target genes of the DEMs in high-fecundity and low-fecundity Yunshang black goats, during the luteal phase, were in Apoptosis (ko04210), the Notch signaling pathway (ko04330), Fc gamma R-mediated phagocytosis (ko04666), the Cytokine-cytokine receptor interaction (ko04060), the Jak-STAT signaling pathway (ko04630), cell adhesion molecules (ko04514), and other 30 KEGG pathways, which were significantly enriched. At the same time, both the differential genes and miRNA predicted target genes were enriched in two signaling pathways, the mTOR signaling pathway (ko04150) and the ErbB signaling pathway (ko04012). Studies have shown that the pathways and the sexual maturation initiation are closely related [72,73].

Conclusions
In conclusion, we analyzed the miRNA and transcriptome profiles in the hypothalamus of Yunshang black goats, during the luteal phase, to identify the key genes and miRNAs that may be involved in the reproductive process, and constructed the complete mRNA-miRNA interaction of the Yunshang black goat. We identified several DEGs (such as CYP19A1, NCAM1, and FGFRs) and miRNA-mRNA pairs (including MEA1 that was co-expressed with the novel-miR-972, novel-miR-125, and novel-miR-403). These miRNAs and genes may play important roles in the regulation of the goat hypothalamus, during the luteal phase, which is worthy of further study in the reproductive regulation of goats and the multiple birth mechanism, in the future. Further studies should validate the functions of these key genes and the interaction networks at the cellular level of the goat fecundity.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ani12233397/s1, Supplementary Table S1: The overall genes expressed and the fragments per kilobase per million mapped fragment (FPKM) values and the chromosome distribution of the mRNAs identified in high-fecundity goats, during the luteal phase (LP-HF) versus the low-fecundity goats, during the luteal phase (LP-LF). Supplementary Table S2: Differentially expressed genes identified in the high-fecundity goats, during the luteal phase (LP-HF) versus the low-fecundity goats, during the luteal phase (LP-LF). Supplementary Table S3: The overall miRNAs expressed identified in the high-fecundity goats, during the luteal phase (LP-HF) versus the low-fecundity goats, during the luteal phase (LP-LF). Supplementary Table S4: The identification involving the diverse RNAs and the length distribution of the small RNA in the high-fecundity goats, during the luteal phase (LP-HF) versus the low-fecundity goats, during the luteal phase (LP-LF). Supplementary Table S5: Differentially expressed miRNAs identified in the high-fecundity goats, during the luteal phase (LP-HF) versus the low-fecundity goats, during the luteal phase (LP-LF). Supplementary Table S6: The list of potential target genes of the differentially expressed miRNAs in the high-fecundity goats, during the luteal phase (LP-HF) versus the low-fecundity goats, during the luteal phase (LP-LF). Supplementary Table S7: GO enrichment annotation for the mRNAs, in terms of their molecular function (MF), biological process (BP), and cellular component (CC) level in the high-fecundity goats, during the luteal phase (LP-HF) versus the low-fecundity goats, during the luteal phase (LP-LF). Supplementary Table S8: KEGG enrichment annotation for the mRNAs in the high-fecundity goats, during the luteal phase (LP-HF) versus the low-fecundity goats, during the luteal phase (LP-LF). Supplementary Table S9: GO enrichment annotation for the potential target genes of the differentially expressed miRNAs, in terms of their molecular function (MF), biological process (BP), and cellular component (CC) level in the high-fecundity goats, during the luteal phase (LP-HF) versus the low-fecundity goats, during the luteal phase (LP-LF). Supplementary Table S10: KEGG enrichment annotation for the potential target genes of the differentially expressed miRNAs in the high-fecundity goats, during the luteal phase (LP-HF) versus the low-fecundity goats, during the luteal phase (LP-LF). Supplementary Table S11: The intersected gene list between the DEGs and predicted target genes of the DEMs in the high-fecundity goats, during the luteal phase (LP-HF) versus the low-fecundity goats, during the luteal phase (LP-LF). Supplementary Table S12: Construction of the miRNA-mRNA interaction network in the high-fecundity goats, during the luteal phase (LP-HF) versus the low-fecundity goats, during the luteal phase (LP-LF). Supplementary Table S13: Real-time quantitative PCR (RT-qPCR) validation of the differentially expressed mRNAs and miRNAs in the high-fecundity goats, during the luteal phase (LP-HF) versus the low-fecundity goats, during the luteal phase (LP-LF). Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The original contributions presented in the study are included in the article and supplementary material, further inquiries can be directed to the corresponding authors.

Conflicts of Interest:
The authors declare no conflict of interest.