Effect of Concentrate Supplementation on the Expression Profile of miRNA in the Ovaries of Yak during Non-Breeding Season

Simple Summary Yak (Bos grunniens) is an important and remarkable livestock species that survives in the challenging environment of the Qinghai–Tibetan Plateau. However, its growth rate is slower and reproductive ability is generally lower than cattle. This may be due to the yak living in high altitudes all year round where in the whole year, grasses are only available in July, August, and September (warm season), and from November to the next year of May (cold season), there is a scarcity of pastures. So, the reproductive efficiency of yak is very low. Meanwhile, it has been reported that enhanced nutrition improves the reproductive efficiency of animals. Therefore, this study aimed to explore the effects of supplemental nutrition on the growth traits and reproductive performance of yaks in the cold season. In addition, miRNAs related to yak reproductive traits were screened by miRNA sequencing technology. This research might be helpful for improving the reproductive potential of yak during the non-breeding season. Abstract Yak (Bos grunniens) is an important and remarkable livestock species that survives in the challenging environment of the Qinghai–Tibetan Plateau. However, its growth rate is slower and reproductive potential is generally lower than cattle. Meanwhile, it has been reported that enhanced nutrition improves the reproductive efficiency of animals. The purpose of this study was to investigate the effect of concentrate supplementation on the miRNA expression profile in the ovaries of yak during the non-breeding season. The study displayed that non-breeding season supplementation significantly improved growth performance, serum biochemical indicators, and reproductive hormone concentrations in yaks. In this study, we also examined the differential expression analysis of miRNA in the ovaries of yak during non-breeding seasons using Illumina Hiseq sequencing technology. As a result, 51 differentially expressed miRNAs were found in the experimental group (CS) and control group (CON). Gene Ontology (go) and Kyoto Genome Encyclopedia (KEGG) analysis of target genes showed that beta-alanine metabolism; tryptophan metabolism; sphingolipid metabolism; alanine, aspartate and glutamate metabolism; and the inositol phosphate metabolism pathway attracted our attention. Based on qRT-PCR, seven miRNAs were assessed to verify the accuracy of the library database. We predicted and identified potential miRNA target genes, including LEP, KLF7, VEGFA, GNAQ, GTAT6, and CCND2. miRNA and corresponding target genes may regulate yaks’ seasonal reproduction through their nutritional status. This study will provide an experimental basis for improving the reproductive efficiency of yaks by supplementation in the non-breeding season.


Sample Collection
At the beginning and end of the cold season supplementary feed experiment, blood was collected from the jugular vein of 20 yaks, and serum was isolated for blood biochemical and reproductive hormone detection. Before the supplementary feed experiment, the selected animals were all about 6 years old and about the same weight. In addition, 3 yaks were selected from the CS and 3 yaks were selected from the CON and executed painlessly. Ovarian tissue samples were collected from 6 yaks for sequencing of miRNA, with 3 yaks were selected from the CS and 3 yaks from the CON. After ovarian sample collection of each animal, they were kept in liquid nitrogen for transportation and finally stored at −80°C.

Blood Biochemistry and Hormonal Concentration
The concentrations of CP (crude protein), EE (crude fat), CF (crude fiber), ash (crude ash), P (phosphorus), and Ca (calcium) in the diet were quantified by feed analysis and feed quality testing techniques. On the 0th and 120th day of the trial period, each yak was weighed in the morning prior to feeding, and the the ADG (average daily weight gain) calculated as follows: Average Daily Gain (ADG) = (Final Weight − Initial Weight)/Days. Albumin (ALB), total protein (TP), serum blood glucose (GLU), triglycerides (TRIG), and cholesterol (CHOL) were measured using a fully automatic bioanalysis machine. Serum follicle-stimulating hormone (FSH), gonadotropin-releasing hormone (GnRH), estrogen (E2), progesterone (PROG), and luteinizing hormone (LH) concentrations were established using ELISA kits and calibrated against standard curves using a micro plate reader.

Extraction of RNA and Sequencing of Small RNA
Total RNA was isolated from the ovaries of concentrate supplement group (CS) yaks and from control group (CON) yaks, by using Trizol reagent (Takara Biotechnology, Dalian, China). Afterward, low-molecular-weight RNA was isolated by using 1% agarose gel electrophoresis (PAGE) and enriched for RNA molecules in the range of 18-30 nt. The product was ligated to the 3 and 5 ends of the optimal temperature with a proprietary adaptor and amplified by RT qPCR. Later, sufficient product was obtained by PCR and tested by using the Agilent 2100 bioanalyzer and the ABI Step One Plus Real-time PCR System and sequencing using Illumina HiSeq (Biomarker Technologies, Beijing, China), which was tested in tissues collected from both (CS and CON) groups.

Analysis through Nioinformatics for Sequencing of Small RNA
The filtered sRNA sequences were mapped to the yak genome map by SOAP or bowtie software, and the expression of these sRNAs and their dissemination on the genome were analyzed and then Animals 2020, 10, 1640 4 of 17 passed through the GenBank database (http://www.ncbi.nlm.nih.gov/genbank/) and Rfam Database (http://rfam.xfam.org/), which screens and removes related sequences, such as repeating ribosomal RNA (rRNA), small interfer RNA (siRNA), small nucleolar RNA (snoRNA), nuclear small RNA (snRNA), and transport RNA (tRNA). The filtered sequence was then followed by miRBase (version 22.0) (http://www.mirbase.org/) to compare and identify identified miRNAs, and predict new miRNAs with mirdeep software (Versionv2.0.5).

Quantification and Differential Analysis of MicroRNA
The miRNA expression levels in each sample were counted and normalized by the TPM algorithm [20]. The TPM normalized treatment formula is as follows: TPM = Readcount * 1,000,000/over over Mappedreads. DESeq R software package was used for differential expression analysis. In the detection of differentially expressed miRNA, |log2 (FC)| ≥ 1.00 and p-value ≤ 0.05 were used as screening criteria.

Target Gene Prediction, Gene Ontology, and Pathway Analysis
Analysis of the pathways could contribute to understanding the biological functions of genes. The KEGG pathway examination was directed at the prediction results of target genes (herein after referred to as "candidate target genes"). In vivo, different genes are mutually coordinated to make their biological functions. The pathway analysis provides further insight into the biologicalfunctions of the gene. In addition, KEGG is the central public database of pathways, and the paths are categorized by significant enrichment analysis. The paths in KEGG are units, and hypergeometric tests are applied to find pathways that are significantly enriched in candidate target genes compared to the entire reference gene. We used target gene prediction for differentially expressed miRNAs, prediction with target scan and RNAhybrid software, and intersection or union as prediction. The gene ontology (GO) analysis of the potential target genes was based on the terms of the Gene Ontology database (http://www.geneontology.org/) by detecting the gene function-related biology process and putting the genes with similar functions together.

Validation of miRNA Expression by Quantitative Real-Time PCR
From the microRNA library, we randomly selected 7 miRNAs, and U6 was used as the internal reference gene. Afterwards, qRT-PCR was performed for expression analyses of miRNA from the ovaries of yak. The reverse primer was universal primer provided with the Mir-X™ miRNA First Strand Synthesis (Takara, Dalian, China) while the forward primers were designed based on the mature miRNA sequences (Table 2), which were obtained from miRBase (http://www.mirbase.org/), and all reactions were performed in 3 replications and the quantitative expression level of each target gene was calculated by using the threshold cycle 2 −∆∆Ct method [21]. verification results were expressed as mean ± SD. An independent t-test was used for the comparison of groups. Differences were considered significant when p < 0.05.

Effect of Mixed Diet Supplementation during the Cold Season on the Growth Performance of Yaks
The effects of a mixed formula diet during the cold season on the growth of yaks are presented in Table 3. At the start and end of the experiment, each yak was weighed. Weight gain in the CS was significantly higher (p < 0.05) than in the CON after mixed diet supplementation during the cold season.

Effect of Mix Diet Supplementation during the Non-Breeding Season on Serum Biochemical Parameters of the Yaks
The serum biochemical parameters of the yak were also evaluated at the start and end of the study. The concentrations of serum biochemical parameters, such as, GLU, TP, ALB, TRIG, and CHOL, were significantly (p < 0.01) higher in the CS than the CON (Table 4).

Effect of Cold Season Mix Diet Supplementation on Serum Reproductive Hormone of Yaks
The concentrations of reproductive hormones in the serum of yaks were measured at the start and after completion of the study. The results revealed that the concentrations of GnRH, FSH, LH, E2, and P in the serum of the CS were significantly higher than those in the CON (Table 5).

Solexa Sequencing of Ovary Small RNAs
This study was designed to understand the expression patterns of small RNA in the yak ovary during the non-breeding season following different nutrient levels, through small RNA construction library analysis. Therefore, two miRNA libraries were constructed from six yaks' ovary samples, including concentrate supplementation (CS) yaks with an enhanced forage nutrient status and CON yaks with a lower forage nutrient status. The two libraries obtained contained 52,727,701 and 41,208,552 clean reads, respectively (Table S1). Sequence length analysis indicated that the majority of the length was 21-23 nt ( Figure 1). The Bowtie software was used to sequence clean reads with the Repbase database(http://www.girinst.org/repbase/), Silva database (http://www.arb-silva.de/), Rfam database(http://rfam.xfam.org/), and GtRNAdb database(http://lowelab.ucsc.edu/GtRNAdb/), respectively, filtering snRNA, snoRNA, rRNA, and tRNA, and a repeat sequence obtained unannotated reads containing miRNA. The sRNA annotation classification statistics are shown in Table 6 and  Table S2.

Solexa Sequencing of Ovary Small RNAs
This study was designed to understand the expression patterns of small RNA in the yak ovary during the non-breeding season following different nutrient levels, through small RNA construction library analysis. Therefore, two miRNA libraries were constructed from six yaks' ovary samples, including concentrate supplementation (CS) yaks with an enhanced forage nutrient status and CON yaks with a lower forage nutrient status. The two libraries obtained contained 52,727,701 and 41,208,552 clean reads, respectively (Table S1). Sequence length analysis indicated that the majority of the length was 21-23 nt ( Figure 1). The Bowtie software was used to sequence clean reads with the Repbase database(http://www.girinst.org/repbase/), Silva database (http://www.arb-silva.de/), Rfam database(http://rfam.xfam.org/), and GtRNAdb database(http://lowelab.ucsc.edu/GtRNAdb/), respectively, filtering snRNA, snoRNA, rRNA, and tRNA, and a repeat sequence obtained unannotated reads containing miRNA. The sRNA annotation classification statistics are shown in Tables 6 and S2.

Identification of Differentially Expressed miRNAs in the Ovary
Sequence alignment and subsequent exploration were achieved using the specified Bos_mutus.v2.0 (http://www.ncbi.nlm.nih.gov/genome/?term=yak) as the reference genome. The Bowtie software (Version v1.0.0, Johns Hopkins University, Baltimore, MD, USA) was used to sequence the unannotated reads with the reference genome to obtain the positional information on the reference genome, which is mapped reads. For the biological characteristics of miRNA, we used the miRDeep2 software package to compare the reads on the reference genome with the known miRNA precursor sequences in the miRBase database (http://www.mirbase.org/) to identify the expression of known miRNAs. After successive filtering of these data sets, we identified a total of 1074 miRNAs (670 known miRNA and 404 novels miRNA) in the libraries (Table S3). Amongst these miRNAs, 51 differentially expressed miRNAs were found, in which 25 miRNAs were downregulated and 26 miRNAs upregulated in the experimental group vs. the control group ( Figure 2). The MA diagram provides an intuitive view of the overall distribution of the expression levels and differential multiples of the two groups ( Figure 3). Systematic cluster analysis was performed to analyze the similarity of miRNAs with different expressions in the two groups of ovaries. The heat map shows the aggregation of CS and CON due to their similar expression profiles (Figure 4).
Animals 2020, 10, x 7 of 17 the reference genome, which is mapped reads. For the biological characteristics of miRNA, we used the miRDeep2 software package to compare the reads on the reference genome with the known miRNA precursor sequences in the miRBase database (http://www.mirbase.org/) to identify the expression of known miRNAs. After successive filtering of these data sets, we identified a total of 1074 miRNAs (670 known miRNA and 404 novels miRNA) in the libraries (Table S3). Amongst these miRNAs, 51 differentially expressed miRNAs were found, in which 25 miRNAs were downregulated and 26 miRNAs upregulated in the experimental group vs. the control group ( Figure 2). The MA diagram provides an intuitive view of the overall distribution of the expression levels and differential multiples of the two groups ( Figure 3). Systematic cluster analysis was performed to analyze the similarity of miRNAs with different expressions in the two groups of ovaries. The heat map shows the aggregation of CS and CON due to their similar expression profiles ( Figure 4).

Validation of miRNA Expression with qRT-PCR
We performed qRT-PCR to detect the expression levels of the differentially expressed miRNAs and to confirm the reliability of the sequencing data ( Figure 5). For this, we selected seven known

Validation of miRNA Expression with qRT-PCR
We performed qRT-PCR to detect the expression levels of the differentially expressed miRNAs and to confirm the reliability of the sequencing data ( Figure 5). For this, we selected seven known miRNAs confirmed in the CS and CON group ovarian samples of yak. The relative expression levels of the seven selected miRNAs indicated that miR-449a, miR-205, and miR-200a were significantly upregulated in the CS than the CON group (p < 0.05), while miR-483, miR-200b, miR-223, and miR-2431-5p were downregulated in the CS group (p < 0.05). The relative expression levels of the seven differentially expressed miRNAs were consistent with the RNA sequencing results. It showed that the sequencing results were true and reliable.

Validation of miRNA Expression with qRT-PCR
We performed qRT-PCR to detect the expression levels of the differentially expressed miRNAs and to confirm the reliability of the sequencing data ( Figure 5). For this, we selected seven known miRNAs confirmed in the CS and CON group ovarian samples of yak. The relative expression levels of the seven selected miRNAs indicated that miR-449a, miR-205, and miR-200a were significantly upregulated in the CS than the CON group (p < 0.05), while miR-483, miR-200b, miR-223, and miR-2431-5p were downregulated in the CS group (p < 0.05). The relative expression levels of the seven differentially expressed miRNAs were consistent with the RNA sequencing results. It showed that the sequencing results were true and reliable.

miRNA Target Gene Prediction, GO Enrichment, and KEGG Pathway Analysis
Target gene prediction was implemented using miRNA and target scan based on newly predicted miRNAs, known miRNAs, and gene sequence evidence of corresponding species. After sequencing, 670 known miRNAs were found, but 418 known miRNAs were able to predict target genes. The remaining known miRNA did not predict target genes. These known miRNAs predicted a total of 12,485 target genes. In addition, GO enrichment analysis was used on the target genes of differentially expressed miRNAs. The differentially expressed miRNAs were enriched in 117 GO terms (p < 0.01), including 26 molecular functions (MFs), 20 cell component (CC), and 71 biological processes (BP) in the CS vs. the CON group (Table S4). Meanwhile, the KEGG pathway analysis of differentially expressed miRNA target genes help us to understand the biological functions of genes. The differentially expressed miRNAs target genes had no significant enriched signal pathway in CS and CON (p > 0.05). Although these signaling pathways were not significant, the top 20 signal pathways contain some metabolic pathways that may be indirectly related to yak reproduction (Table S5). These signaling pathways include beta-alanine metabolism, tryptophan metabolism; sphingolipid metabolism; alanine, aspartate, and glutamate metabolism; inositol phosphate metabolism, etc. (Figure 6). We also identified some key reproduction-related genes, such as LEP, NOTCH1, KLF7, VEGFA, GNAQ, GTAT6, and CCND2 (Table 7).

miRNA Target Gene Prediction, GO Enrichment, and KEGG Pathway Analysis
Target gene prediction was implemented using miRNA and target scan based on newly predicted miRNAs, known miRNAs, and gene sequence evidence of corresponding species. After sequencing, 670 known miRNAs were found, but 418 known miRNAs were able to predict target genes. The remaining known miRNA did not predict target genes. These known miRNAs predicted a total of 12,485 target genes. In addition, GO enrichment analysis was used on the target genes of differentially expressed miRNAs. The differentially expressed miRNAs were enriched in 117 GO terms (p < 0.01), including 26 molecular functions (MFs), 20 cell component (CC), and 71 biological processes (BP) in the CS vs. the CON group (Table S4). Meanwhile, the KEGG pathway analysis of differentially expressed miRNA target genes help us to understand the biological functions of genes. The differentially expressed miRNAs target genes had no significant enriched signal pathway in CS and CON (p > 0.05). Although these signaling pathways were not significant, the top 20 signal pathways contain some metabolic pathways that may be indirectly related to yak reproduction (Table S5). These signaling pathways include beta-alanine metabolism, tryptophan metabolism; sphingolipid metabolism; alanine, aspartate, and glutamate metabolism; inositol phosphate metabolism, etc. (Figure 6). We also identified some key reproduction-related genes, such as LEP, NOTCH1, KLF7, VEGFA, GNAQ, GTAT6, and CCND2 (Table 7).  were not significant, the top 20 signal pathways contain some metabolic pathways that may be indirectly related to yak reproduction (Table S5). These signaling pathways include beta-alanine metabolism, tryptophan metabolism; sphingolipid metabolism; alanine, aspartate, and glutamate metabolism; inositol phosphate metabolism, etc. (Figure 6). We also identified some key reproduction-related genes, such as LEP, NOTCH1, KLF7, VEGFA, GNAQ, GTAT6, and CCND2 (Table 7).

Discussion
The blood biochemistry index can reflect the change of organism metabolism to a certain extent. The concentration of total protein and albumin in serum is closely related to protein metabolism of the body, which can reflect the protein level in the animal's diet and the degree of digestion and absorption of protein. The results showed that the concentration of total protein and albumin in yak serum in the supplementary group was significantly higher than that in the control group. This indicated that the supplementation of concentrate increased the intake of yak protein, so the synthesis of protein in vivo increased, and the concentration of total serum protein and albumin increased. This is consistent with published research [22][23][24]. Serum triglyceride levels are closely related to fat deposition and energy metabolism in animals, which can be broken down for energy and converted to adipose tissue for storage. After the end of the experiment, the serum triglyceride content of the supplement group was significantly higher than that of the control group. The results showed that the supplementary feed increased the intake of nutrients and provided a substrate for lipid synthesis in yaks, which was beneficial to the fat deposition of the growing yaks. Blood sugar is an important component of the body and an important source of energy. The body needs a lot of sugar every day to provide energy for the normal operation of various tissues and organs to provide power. Cholesterol is an indispensable and important substance for animal tissue cells. It is not only involved in the formation of cell membranes but also the raw material for the synthesis of bile acids, vitamin D, and steroids. The results showed that the serum glucose and cholesterol levels of the supplement group were significantly higher than those of the control group.
For female reproduction, appropriate nutrition is an essential factor at all stages [25]. Under adequate nutrient conditions, follicular growth and the ovulation rate increased magnificently in livestock while these were reduced during malnutrition [26]. Meza-Herrera et al. [27] illustrated that additional supplementation in the diet improves the reproductive potential of animals. Similarly, some other findings displayed that nutritional supplementation influenced the serum biochemical parameters [28][29][30]. The current findings are also in agreement with the above findings and our results revealed that concentrate supplementation enhances the biochemical parameters and concentration of reproductive hormones.
Yak is an endemic species in high-altitude areas, and it grazes all year round. Therefore, yaks have certain wild habits. Compared with dairy cattle, sheep, and other domestic animals, yak sampling and experimental data collection is difficult. Therefore, three yaks were selected in the experimental group and the control group, respectively. Then, the miRNA sequencing analysis was carried out. Ma X et al. [31] studied the transcriptome analysis of the molecular mechanism of yak longissimus dorsi in different development stages. Three biological repeats were selected for follow-up analysis in each group. LAN D et al. [32] studied the comparative transcriptome analysis of yak and yellow cattle ovaries during estrus. Three biological replicates in each group were selected for follow-up analysis. Therefore, three biological repeats in each group were statistically significant.
Yak is an endemic species in high-altitude areas, and it grazes all year round. Therefore, yaks have certain wild habits. Compared with dairy cattle, sheep, and other domestic animals, yak sampling and experimental data collection is difficult. Therefore, three yaks were selected in the experimental group and the control group, respectively. Then, the miRNA sequencing analysis was carried out. Ma X et al. [31] studied the transcriptome analysis of the molecular mechanism of yak longissimus dorsi in different development stages. Three biological repeats were selected for follow-up analysis in each group. LAN D et al. [32] studied the comparative transcriptome analysis of yak and yellow cattle ovaries during estrus. Three biological replicates in each group were selected for follow-up analysis. Therefore, three biological repeats in each group were statistically significant.
Understanding the molecular regulation of miRNAs in the ovaries helps to explore the molecular mechanisms by which yak regulates reproduction during the non-breeding season. Many studies have shown that miRNA plays an essential role in almost all ovarian biological processes, including follicular development, follicular atresia, corpus luteum development, and degeneration [16]. This study used Solexa sequencing technology to compare the small RNA profiles of yak ovaries during the non-breeding season after concentrate supplementation. The objective was to understand the role of sRNA-mediated post-transcriptional regulation in seasonal reproduction. By obtaining miRNAs with significant differences, we found some miRNAs, such as bta-miR-483, bta-miR-449a, bta-miR-223, bta-miR-200b, bta-miR-200a, and bta-miR-2431-5p. Earlier findings illustrated that miR-483 plays a vital role in the development and regulation of livestock reproduction. Besides, another study presented that miR483 is a hormone-sensitive testis miRNA [33]. Gabriella Guelfi et al. [34] found that the reproductive efficiency of buffalo was affected by a long calving interval and late estrus. Therefore, the expression level of mir-200b in buffalo progesterone mature oocytes and pregnancy was evaluated. It was found that there were differences in the expression between pregnant and non-pregnant buffalo. Heng Yang et al. [35] explored the regulation of miRNA in sheep on seasonal reproduction and the nutritional status. It was found that mir-200b was related to the nutritional status and seasonal reproduction of sheep. Meanwhile, Sohel et al. [36] established 27 differentially expressed miRNAs, such as miR-223, between the mature oocytes and non-exosome part of follicular fluid.
In our study, RT-qPCR results showed that miR-483, miR-223, miR-200b, and miR-2431-5p were significantly downregulated in the experimental group as compared to the control group, but miR-449a, miR-205, and miR-200a were upregulated in the experimental group. These results indicate that miR-483, miR-200b, miR-2431-5p, miR-449a, miR-223, and miR-200a may be involved in the reproduction of yak through nutritional changes. In addition, the top 20 pathways, including beta-alanine metabolism; tryptophan metabolism; sphingolipid metabolism; alanine, aspartate, and glutamate metabolism; and inositol phosphate metabolism, have attracted our attention. Through KEGG pathway and software prediction, we screened the target genes of differentially expressed miRNA. Due to different nutritional levels, differentially expressed mir-483, mir-200b, mir-200a, mir-2431-5p, and mir-223 and their corresponding target genes, such as LEP, VEGFA, GNAQ, GTAT6, CCND2, and KLF7, may have potential regulatory effects on yak reproductive traits (Table 8). Leptin may act as a metabolic signal, transmitting nutrients to the reproductive endocrine system. The hypothalamic-pituitary-gonadal axis (HPG) is activated when the body is suitable, which plays an important role in female animal estrus, fertilized egg implantation, early embryo development, pregnancy maintenance, fetal growth and development, and other reproductive activities [37,38]. Previous studies presented that leptin (LEP) controls the reproductive process at the periphery of the hypothalamic-pituitary and reproductive tissues (e.g., ovaries). Additionally, the leptin is expressed in the ovary, indicating that leptin and its receptors might have autocrine and paracrine effects [39]. Leptin may affect the process of steroidogenesis either indirectly or directly by modulating the actions of metabolic hormones, for example, IGF-I GH and insulin [40]. Meanwhile, another study displayed that leptin was associated with a shorter interval of estrus and a higher concentration of leptin was associated with it, demonstrating the relationship between leptin and estrus expression [41]. Additionally, Meikle et al. [42] displayed that in primiparous lean cows, the ovarian cycle started late, with longer intervals from birth to first conception, and the endocrine signal that most likely represents the negative energy balance of the reproductive axis is leptin. In addition, another study illustrated that leptin has different regulatory effects on the gene expression of oocytes and cumulus cells. Additionally, leptin promotes oocyte maturation and development through mechanisms independent of cumulus cells and dependent mechanisms [43]. In our study, mir-483 was upregulated in the CS group, and LEP was the key target gene predicted by mir-483. These results suggest that mir-483 targeting LEP may be involved in the reproductive process of yak.
FSH and LH are the most famous endocrine regulatory factors, which mainly promote the development and growth of follicles. However, the normal occurrence of this process also needs the participation of other factors, including vascular endothelial growth factor A (VEGFA) [44]. VEGFA is a member of the VEGF family of cysteine knot growth factors. Its main function is to stimulate the mitosis of endothelial cells. In the ovary, VEGFA was expressed in both granulosa cells and follicles as ovulation approached, and the expression of VEGFA in granulosa cells was significantly increased [45]. Meanwhile, some evidence shows that VEGFA plays an important role in follicular survival and inhibition of granulosa cell apoptosis [46]. The researchers injected VEGFA into the ovaries to promote the development of follicles [47]. The results of Gao et al. showed that the expression of VEGFA reached a peak in tertiary-generation follicles during the follicular development stage, mainly distributed in the capsule layer and the inner granular layer. During follicular development, VEGFA mRNA was mainly expressed in the inner layer of granules [48]. GNAQ is a member of G protein-coupled receptor superfamily. As an important G protein-coupled receptor, GNAQ promotes glucose-induced insulin release and participates in the process of nutritional metabolism [49,50]. In addition, GNAQ can also affect the ovary and increase the size of ovarian follicles during ovulation. GNAQ affects follicular growth and ovarian atresia by inhibiting PI3K-Akt signaling pathway [51,52]. Heng Yang et al. found that the effect of mir-200b differentially expressed in sheep and target gene GNAQ had a regulatory effect on the seasonal reproduction and nutritional status of sheep [35]. In our study, mir-200b expression was downregulated in the CS group, and VEGFA and GNAQ were key target genes predicted by mir-200b. These results suggest that mir-200b targeting VEGFA and GNAQ may be involved in yak reproductive process.
GATA6 is a member of the GATA family. It plays an important role in the proliferation, differentiation, and apoptosis of ovarian cells. In the mouse ovary, GATA6 mRNA was detected in granulosa cells and strongly expressed in the corpus luteum but not in membrane cells [53]. In the human ovary, GATA6 mRNA was located in the granulosa cells and membrane cells of pre sinus and pre sinus follicles [54]. GATA-6 protein was found in granulosa cells of the pig ovary [55]. Cyclins D (CCNDs) include three subtypes: CCND1/CCND 2, and CCND3. In the ovary, CCND 2 is mainly expressed in granulosa cells. As one of the downstream factors of follicle-stimulating hormone (FSH), CCND 2 promotes the proliferation of granulosa cells during follicular development [56]. In different stages of bovine follicular development, the expression level of CCND 2 in granulosa cells of the active estrogen stage and preovulatory follicle was high, which was closely related to the proliferation of granulosa cells [57]. Through positive regulation of CCND 2 expression, it can promote the cell cycle process and accelerate the proliferation of granulosa cells [58]. In our study, mir-200a was upregulated and mir-2431-5p was downregulated in the CS group. GTAT6 and CCND2 are key target genes predicted by mir-200a and mir-2431-5p. These results suggest that mir-200a and mir-2431-5p targeting GTAT6 and CCND 2 may participate in the reproductive process of yak at reproductive age.
Besides, previous studies stated that KLF7 plays an essential role in cell proliferation, differentiation, and other processes, and participates in physiological processes, such as embryo development and glucose metabolism [59,60]. Whereas, genetic correlation analysis indicated that KLF7 is a candidate gene for human obesity, which inhibits the differentiation of adipose precursor cells and downregulates many genes involved in adipose differentiation [61]. Overexpression of KLF7 in already differentiated adipocytes inhibits the expression of adipocytokines, such as adiponectin [62]. In our study, KLF7 is a key target gene predicted by mir-223. The expression of mir-223 was upregulated in the CS group. These results suggest that mir-223 targeting KLF7 may be involved in the yak reproductive process.

Conclusions
In summary, our findings illustrated that non-breeding season supplementation significantly improved growth performance, serum biochemical indicators, and reproductive hormone concentrations in yaks. In addition, we found six differentially expressed miRNAs (mir-483, mir-449a, mir-223, mir-200b, mir-200a, mir-2431-5p). Through target gene prediction, GO enrichment, and KEGG pathway analysis, we identified that target genes corresponding to the potential miRNA may regulate the seasonal reproduction of yaks and be related to nutritional status. This study provides a new reference for the identification of new genetic markers in yaks during the non-breeding period.
Author Contributions: J.X. and X.G. conceptualization and designed the experiments; J.X. analyzed the data; X.G. Methodology and data curation; J.X. and Q.K. wrote the manuscript; P.Y. supervision and approved the final draft of the manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: This study: including experiments design, samples collection, data analysis and manuscript writing were funded by the China Agriculture Research System (CARS-37).

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