Screening of Differentially Expressed Genes and miRNAs in Hypothalamus and Pituitary Gland of Sheep under Different Photoperiods

The reproduction of sheep is affected by many factors such as light, nutrition and genetics. The Hypothalamic-pituitary-gonadal (HPG) axis is an important pathway for sheep reproduction, and changes in HPG axis-related gene expression can affect sheep reproduction. In this study, a model of bilateral ovarian removal and estrogen supplementation (OVX + E2) was applied to screen differentially expressed genes and miRNAs under different photoperiods using whole transcriptome sequencing and reveal the regulatory effects of the photoperiod on the upstream tissues of the HPG axis in sheep. Whole transcriptome sequencing was performed in ewe hypothalamus (HYP) and distal pituitary (PD) tissues under short photoperiod 21st day (SP21) and long photoperiod 21st day (LP21). Compared to the short photoperiod, a total of 1813 differential genes (up-regulation 966 and down-regulation 847) and 145 differential miRNAs (up-regulation 73 and down-regulation 72) were identified in the hypothalamus of long photoperiod group. Similarly, 2492 differential genes (up-regulation 1829 and down-regulation 663) and 59 differential miRNAs (up-regulation 49 and down-regulation 10) were identified in the pituitary of long photoperiod group. Subsequently, GO and KEGG enrichment analysis revealed that the differential genes and target genes of differential miRNA were enriched in GnRH, Wnt, ErbB and circadian rhythm pathways associated with reproduction. Combined with sequence complementation and gene expression correlation analysis, several miRNA-mRNA target combinations (e.g., LHB regulated by novel-414) were obtained. Taken together, these results will help to understand the regulatory effect of the photoperiod on the upstream tissues of HPG in sheep.


Introduction
Animals can adapt themselves to their survival environment by regulating themselves, for example, females will choose to breed in specific seasons to ensure that they give birth in the most appropriate season, which will facilitate lactation and early pup development, and a stable physiological mechanism has evolved over the years [1]. Light perception is one of the important mechanisms that trigger the response of organisms and many physiological processes in mammals are influenced by the photoperiod [2]; a long photoperiod induces circadian transcription factor BMAL2 in pars tuberalis (PT) and triggers summer biological activity through the EYA3/TSH pathway. In contrast, the continuous prolongation of a short photoperiod leads to melatonin secretion inducing circadian repressors, including DEC1, which inhibit BMAL2 and EYA3/TSH pathways, into winter biological activity, and ewes with good estrous traits as described above, estrogen silicone tubes were implanted subcutaneously, and the sheep were penned into artificially controlled light sheds after the hormone concentration was stabilized at 7.23 ± 2.50 pg/mL. OVX + E 2 sheep and light control rooms were constructed as previously described [27]. After 21 days of short photoperiod (8/16 h light-dark), six were euthanized (intravenous pentobarbital 100 mg/kg). After slaughtering, hypothalamic and distal pituitary tissues were taken and were labeled as SP21−HYP and SP21−PD, respectively. Six other Small-tailed Sheep were transferred to a long photoperiod (16/8 h light-dark) sheep barn and after 21 days of rearing, tissues were slaughtered and collected and labeled LP21−HYP and LP21−PD. All collected tissue was snap-frozen at liquid nitrogen, and then brought back to the laboratory and stored in a −80 • C refrigerator for long-term storage.

RNA Extraction, Library Construction, Sequencing and Raw Data Processing
Total RNA from each sample was isolated using the TRIzol reagent (Invitrogen, Carlsbad, CA, USA). then the Nanodrop 2000 instrument was used to detect the concentration and purity of the extracted RNA, and the integrity of RNA was detected by agarose gel electrophoresis. After the samples were tested and passed, cDNA and sRNA library construction were performed according to the previous method [29]. To finally obtain valid data, we filtered the raw sequencing data to obtain clean reads. Clean reads were aligned to the sheep reference genome using bowtie (Oar_rambouillet_v1.0).

Differential Expression Analysis of mRNAs
Quantitative analysis was performed using StringTie software to obtain the number of Reads per sample compared to each transcript, which was converted to FPKM (Fragments Per Kilobase per Million bases). Quantitative statistical analysis of transcript expression was performed to filter out transcripts with significantly different expression levels in the samples. We used DESeq2 to analyze the significance of transcript expression differences. The differentially expressed transcripts were screened from two levels of fold change (FC) and corrected significance level (padj/q-value). In addition, |log 2 FC| > 1 and q-value < 0.05 were considered DEGs.

Differential Expression Analysis of miRNAs and Target Gene Prediction
Identification of known miRNAs was performed using bowtie software by comparing sRNAs to miRNA sequences from the specified range of miRBase (v22) database, and prediction of new miRNAs was performed using mirEvo and miRDeep software and expression information was obtained. The expression of miRNAs represented by TPM was obtained by comparing to the reads of miRNA precursor and mature body regions. After calculating the expression, differential analysis was performed using DESeq2, and the screening threshold for differential miRNAs was set to |log 2 FC|> 1, p value < 0.05.
The previously obtained known miRNAs and new miRNAs were analyzed for miRNA target gene prediction. The miRanda and qTar software were selected for target gene prediction, and the intersection of the two software was taken as the final result.

Functional Annotation and Enrichment Analysis of Target Genes of DE miRNAs and DE mRNAs
The DEGs and the predicted target genes of DEMs were analyzed by Gene Ontology (GO) (referred to as GO, http://www.geneontology.org/, accessed on 6 January 2022), which includes biological process, molecular function and cellular component, and Kyoto Encyclopedia of Genes and Genomes (KEGG). GO items or KEGG pathways with a hypergeometric p-value < 0.05 were those that were significantly enriched.

Construction of Integral miRNA-mRNA Interaction Networks
To accurately identify the key association with reproductive DEMs and DEGs, based on miRNA function, the mRNAs that were negatively related to miRNAs were screened out, and the miRNA-mRNA interaction networks were built by using Cytoscape (v3.8.2, http://www.cytoscape.org/, accessed on 1 April 2022).

Quantitative PCR Validation
To verify the accuracy of sequencing data, 7mRNAs including hypothalamus (DIO2, SIX1, RASD1, KCNH3) and distal pituitary (OXT, THUMPD1, LOC101114319) were randomly selected. The first primers were synthesized by Bioengineering (Shanghai) Co., Ltd. For the qPCR analysis of mRNAs, reverse transcription was performed using the PrimeScriptTM RT reagent kit (TaKaRa, Dalian, China). Furthermore, qPCR with the SYBR Green qPCR Mix (TaKaRa, Dalian, China) was conducted with a RocheLight Cycler480 ® II system (Roche Applied Science, Mannheim, Germany) as follows: initial denaturation at 95 • C for 5 min, followed by 40 cycles of denaturation at 95 • C for 5 s, then annealing at 60 • C for 30 s. The sequences of qPCR primers are listed in Table 1. Eight miRNAs were randomly selected, including hypothalamus (oar-miR-200a, novel_10, novel_101, novel_4, oar-miR-143) and distal pituitary (novel_142, oar-miR-181a, oar-miR-323b). For the qPCR analysis of miRNAs, reverse transcription was performed using miRcute Plus miRNAs First-Strand cDNA Kit (TIANGEN, Beijing, China). Then, qPCR was conducted miRcute Plus miRNA qPCR Kit (TIANGEN, Beijing, China) using a RocheLight Cycler480 ® II system (Roche AppliedScience, Mannheim, Germany) in the following procedure: initial denaturation at 95 • C for 15 min, followed by 40 cycles of denaturation at 94 • C for 20 s, then annealing at 60 • C for 34 s. The primers specific to miRNAs were synthesized by Beijing Tianyi Huiyuan Biotechnology Co., Ltd. (Beijing, China). The primer design information is shown in Table 2. In addition, RPL19 (for mRNA) and U6 small nuclear RNA (for miRNA) were utilized as reference gene-miRNA to calculate the relative expression level with the method of 2 −∆∆Ct .

Summary of Sequencing Data for mRNA and miRNA
RNA-seq for mRNA generated about 673 million (hypothalamus) and 640 million (distal pituitary) clean reads after data filtering, and more than 94% of the clean reads were located in the genome. The GC content of clean reads was between 47.97 and 53.37%, and the clean reads quality scores of Q20 and Q30 were above 96.29% and 89.87%, respectively, indicating that the reliability and quality of the sequencing data were sufficient for further analysis (Supplementary Table S1). Regarding the expression level of mRNAs, details are provided in Supplementary Table S2, Figure 1A,C. In addition, the chromosomal distribution of mRNA on both tissues showed that genes were mostly present on chromosomes 1, 2 and 3, occupying 27.53% and 27.46%, respectively, (Supplementary Table S2, Figure 1B    RNA-seq for miRNAs generated about 140 million (hypothalamus) and 132 million (distal pituitary) clean reads after data filtering, and in the range of 92.99-98.65% of the clean reads, which were located in the genome (Supplementary Table S1), with the length between 21-24nt (Supplementary Table S3, Figure 2B,D). The average Q20 content was 98% (hypothalamus) and 96% (distal pituitary) (Supplementary Table S1). In addition, a variety of non-coding RNAs (ncRNAs) were also identified, including transfer RNAs (tRNAs), snR-NAs and miRNAs (Supplementary Table S3, Figure 2A

Identification of Differential Expressed Genes and miRNAs
A total of 1813 DEGs were identified from SP21−HYP vs. LP21−HYP compari Compared to the short photoperiods,966 genes were up-regulated, and 847 genes w down-regulated in the long photoperiods (Supplementary Table S4, Figure 3A). Simil

Identification of Differential Expressed Genes and miRNAs
A total of 1813 DEGs were identified from SP21−HYP vs. LP21−HYP comparison. Compared to the short photoperiods,966 genes were up-regulated, and 847 genes were down-regulated in the long photoperiods (Supplementary Table S4, Figure 3A). Similarly, a total of 2492 DEGs were identified from SP21−PD vs. LP21−PD comparison. Among these DEGs, 1829 genes were up-regulated, and 663 genes were down-regulated in the long photoperiods (Supplementary Table S5, Figure 3C). The 10 mRNAs with the most significant difference are shown in Table 2.
A total of 145 DEMs were detected among the 1051 identified miRNAs. Compared to the short photoperiods, 73 were up-regulated and 72 were down-regulated in the long photoperiods (Supplementary Table S4, Figure 3B). A total of 59 DEMs were detected among the 1051 identified miRNAs, of which 49 were up-regulated and 10 were downregulated (Supplementary Table S5, Figure 3D). The 10 mRNAs with most significant difference was shown in Table 3, the 10 miRNAs with most significant difference was shown in Table 4, Sequence details of miRNAs are available at Supplementary Table S5.

Functional Enrichment Analysis of the DEGs
GO and KEGG analyses were conducted for the differentially expressed mRNAs. In GO terms, the most enriched terms in SP21−HYP vs. LP21−HYP are actin cytoskeleton organization and protein phosphorylation (Supplementary Table S6, Figure 4A). KEGG analysis showed that MAPK signaling pathway, GnRH signaling pathway, and insulin signaling pathway were significantly enriched (Supplementary Table S6, Figure 4C). Table S4, Figure 3B). A total of 59 DEMs were detected among the 1051 identified miRNAs, of which 49 were up-regulated and 10 were downregulated (Supplementary Table S5, Figure 3D). The 10 mRNAs with most significant difference was shown in Table 3, the 10 miRNAs with most significant difference was shown in Table 4, Sequence details of miRNAs are available at Supplementary Table S5.    Figure 4C). In GO terms, the most enriched terms in SP21−PD vs. LP21−PD are modulation of chemical synaptic transmission and regulation of trans-synaptic signaling (Supplementary Table S7, Figure 4C). All of these pathways are related to the basic activities of life, such as circadian entrainment, oxytocin signaling pathway and GnRH signaling pathway, and are associated with reproduction (Supplementary Table S7, Figure 4D).   Table S7, Figure 4C). All of these pathways are related to the basic activities of life, such as circadian entrainment, oxytocin signaling pathway and GnRH signaling pathway, and are associated with reproduction (Supplementary Table S7, Figure 4D).

miRNA-mRNA Interaction Network
To further analyze the relationship between miRNAs and mRNAs, we constructed a genomic network of interactions using potential target genes of DEMs with DEGs obtained by RNA-Seq. Crossing the identified potential target genes of DEMs with the DEGs obtained from RNA-Seq, there were 468 hypothalamic overlapping genes and 426 pituitary overlapping genes ( Figure 5). Then, Cytoscape was used to construct a co-regulatory network of DE miRNA-mRNA pairs with multiple targeting relationships ( Figure 6, Supplementary Tables S10 and S11).
To further determine the role of genes and miRNAs in the regulation of the gonadal axis by light conversion. GO and KEGG analysis was undertaken for the above overlapping genes (i.e., simultaneous differential expression of miRNAs and target genes). The two most enriched entries of GO are the regulatory activity of GTPase regulator activity and nucleoside-triphosphatase regulator activity in SP21−HYP vs. LP21−HYP (Supplementary  Table S8, Figure 7A). GO results found that positive regulation of glial cell production was significantly enriched in SP21−PD vs. LP21−PD (Supplementary Table S9, Figure 7B). The pathways of KEGG enrichment of target genes in DEMs and DEGs were similar in both groups (Supplementary Table S8, Figure 7C, Supplementary Table S9, Figure 7D). network of DE miRNA-mRNA pairs with multiple targeting relationships ( Figure 6, Supplementary Tables S10 and S11).
To further determine the role of genes and miRNAs in the regulation of the gonadal axis by light conversion. GO and KEGG analysis was undertaken for the above overlapping genes (i.e., simultaneous differential expression of miRNAs and target genes). The two most enriched entries of GO are the regulatory activity of GTPase regulator activity and nucleoside-triphosphatase regulator activity in SP21−HYP vs. LP21−HYP (Supplementary Table S8, Figure 7A). GO results found that positive regulation of glial cell production was significantly enriched in SP21−PD vs. LP21−PD (Supplementary Table S9, Figure 7B). The pathways of KEGG enrichment of target genes in DEMs and DEGs were similar in both groups (Supplementary Table S8, Figure 7C, Supplementary Table S9, Figure 7D).  tary overlapping genes ( Figure 5). Then, Cytoscape was used to construct a co-regulatory network of DE miRNA-mRNA pairs with multiple targeting relationships ( Figure 6, Supplementary Tables S10 and S11).
To further determine the role of genes and miRNAs in the regulation of the gonadal axis by light conversion. GO and KEGG analysis was undertaken for the above overlapping genes (i.e., simultaneous differential expression of miRNAs and target genes). The two most enriched entries of GO are the regulatory activity of GTPase regulator activity and nucleoside-triphosphatase regulator activity in SP21−HYP vs. LP21−HYP (Supplementary Table S8, Figure 7A). GO results found that positive regulation of glial cell production was significantly enriched in SP21−PD vs. LP21−PD (Supplementary Table S9, Figure 7B). The pathways of KEGG enrichment of target genes in DEMs and DEGs were similar in both groups (Supplementary Table S8, Figure 7C

Data Validation
In order to assess the accuracy of sequencing, mRNAs and miRNAs were selected randomly for qPCR validation. We measured the gene expression level and compared it with the RNA-seq data. The results demonstrated that RNA-seq data and qPCR data showed similar patterns (Figure 8), which indicate the reliability of the data generated from RNA-seq.

Data Validation
In order to assess the accuracy of sequencing, mRNAs and miRNAs were selected randomly for qPCR validation. We measured the gene expression level and compared it with the RNA-seq data. The results demonstrated that RNA-seq data and qPCR data showed similar patterns (Figure 8), which indicate the reliability of the data generated from RNA-seq.

Discussion
In most animals, reproduction is regulated by the hypothalamic-pituitary-gonadal axis [30]. The hypothalamus and pituitary gland are important upstream tissues that regulate HPG in sheep and control a wide range of animal behaviors [31,32]. To better probe the effects of the photoperiod on hypothalamic and pituitary function, we applied the OVX + E2 sheep model, which is often used to study the response of seasonally breeding animals to photoperiodic and reproductive endocrine changes [33][34][35]. Therefore, this study used transcriptomic approaches to mine differential genes and miRNA expression patterns associated with reproduction in sheep with different photoperiods.

Identification and Analysis of mRNAs and miRNAs Data in SP21− HYP vs. LP21− HYP
A total of 145 DEMs were identified for SP21−HYP vs. LP21−HYP, including miR-200 [36]. Family members were identified multiple times in seasonal and non-seasonal estrus ovaries of sheep. Among them, miR-200a is a member of the miR-200 family, which is associated with ovarian development in goats [37]. One study showed that miR-494-3p significantly inhibited apoptosis of endometrial epithelial cells [38], and oar-miR-485-3p can regulate the developmental quality of oocytes [39]. From the above, we can speculate that DEMs target and regulate genes related to reproductive estrus, which may be responsible for the reproductive variation in sheep due to different photoperiodic changes.
A total of 1813 DEGs were identified in SP21−HYP vs. LP21−HYP, including APAF1, CASP3 associated with apoptosis, and DIO2, LHB, and TSHB associated with estrous reproduction in sheep. Some studies have demonstrated that haplotype H2 and diploid H2/H4 of TSHB may be associated with year-round estrus [40]. GO was enriched for many terms associated with intercellular communication, and differential genes expressed may influence changes in reproductive traits in sheep by affecting the expression of intercellular communication-related receptor proteins. KEGG pathways include enriched ErbB signaling pathway; GnRH signaling pathway (PTK2B and CAM2B), gonadotropin secretion (LHB), and prolactin signaling pathway were also significantly enriched; these pathways have been studied and proven to be associated with reproductive activity [27,41].

Discussion
In most animals, reproduction is regulated by the hypothalamic-pituitary-gonadal axis [30]. The hypothalamus and pituitary gland are important upstream tissues that regulate HPG in sheep and control a wide range of animal behaviors [31,32]. To better probe the effects of the photoperiod on hypothalamic and pituitary function, we applied the OVX + E 2 sheep model, which is often used to study the response of seasonally breeding animals to photoperiodic and reproductive endocrine changes [33][34][35]. Therefore, this study used transcriptomic approaches to mine differential genes and miRNA expression patterns associated with reproduction in sheep with different photoperiods.

Identification and Analysis of mRNAs and miRNAs Data in SP21− HYP vs. LP21− HYP
A total of 145 DEMs were identified for SP21−HYP vs. LP21−HYP, including miR-200 [36]. Family members were identified multiple times in seasonal and non-seasonal estrus ovaries of sheep. Among them, miR-200a is a member of the miR-200 family, which is associated with ovarian development in goats [37]. One study showed that miR-494-3p significantly inhibited apoptosis of endometrial epithelial cells [38], and oar-miR-485-3p can regulate the developmental quality of oocytes [39]. From the above, we can speculate that DEMs target and regulate genes related to reproductive estrus, which may be responsible for the reproductive variation in sheep due to different photoperiodic changes.
A total of 1813 DEGs were identified in SP21−HYP vs. LP21−HYP, including APAF1, CASP3 associated with apoptosis, and DIO2, LHB, and TSHB associated with estrous reproduction in sheep. Some studies have demonstrated that haplotype H2 and diploid H2/H4 of TSHB may be associated with year-round estrus [40]. GO was enriched for many terms associated with intercellular communication, and differential genes expressed may influence changes in reproductive traits in sheep by affecting the expression of intercellular communication-related receptor proteins. KEGG pathways include enriched ErbB signaling pathway; GnRH signaling pathway (PTK2B and CAM2B), gonadotropin secretion (LHB), and prolactin signaling pathway were also significantly enriched; these pathways have been studied and proven to be associated with reproductive activity [27,41].

Analysis of miRNA-mRNA Interaction Network in SP21−HPY vs. LP21−HPY
To better understand the function of miRNAs, the SP21−HPY and LP21−HPY interaction negative network was constructed. In the upregulated network, a key oar-miR-200a targeting CACNA2D1 was identified, and it has been shown that miR-200a-3p plays an important role in chicken [42] and goat [37] reproductive regulation, and some other miRNAs target genes that regulate cell proliferation apoptosis-related biological processes [43].
Among the down-regulated networks, the important DEGs and DEMs were selected to construct the reciprocal network, in which DENND1A was the target gene of novel_550, and it was shown that DENND1A gene was expressed higher in the low yield group than in the high yield group in ovarian and showed an opposite expression trend with chi-miR-324-3p. CCK-8 assay showed that chi-miR-324-3p overexpression significantly inhibited the proliferation of GCs and knockdown of chi-miR-324-3p promoted GC proliferation [44], thus affecting ovulation and reproduction in ewes, while novel_550 inhibited DENND1A gene expression in this study. We speculate that novel_550 and chi-miR-324-3p act similarly. It has also been reported that DENND1A is a susceptibility gene for polycystic ovaries, which is essential for ovarian and embryonic development [45][46][47] and can affect the reproductive process. Another target gene of novel_550, EML6 was also shown to be associated with reproduction. The production of high-quality eggs in mouse reproduction requires normal segregation during oocyte meiosis, and EML6 is highly expressed in oocytes and is responsible for the accurate segregation of homologous chromosomes [48]. GH is a target gene of novel_984; GH plays an inductive role in ovulatory reproduction [49], and is elevated during growth hormone secretion during major transitions in reproductive status such as puberty and pregnancy [50]. It is an important regulator of female reproduction, and is involved in gonadal steroidogenesis, gametogenesis and ovulation [51], and alterations to the GH axis can have a reproduction. In growth regulation, it is mainly used with growth hormone (GH) to promote cellular metabolism, promote growth and regulate reproduction [52]. Novel_414 targets LHB, luteinizing hormone (LH), also known as luteinizing hormone, which is released by gonadotropins and promotes blood flow and secretion of progesterone in the ovary, affecting ovulation [53]. SIX1 is a common target gene of novel_618 and novel_424, and long light signaling molecules are mainly mediated by EYA3 and TSHβ. EYA3 can form a complex with thyrotropic embryonic factor (TEF) and SIX1, which synergistically promotes TSHβ transcription upon binding to the TSHβ promoter D element [54]. In summary, in the down-regulated network, the expression of the above DEGs in this study was significantly higher during the shift from short to long photoperiods, and we hypothesize that DEMs such as chi-miR-324-3p, novel_550, and novel_414 lead to alterations in reproductive traits in sheep by negatively regulating and reproduction-related genes.

Identification and Analysis of mRNAs and miRNAs Data in SP21−PD vs. LP21−PD
A total of 2942 DEGs were identified for SP21−PD vs. LP21−PD, including CREB3L1 and OXT. CREB3L1 acts as a downstream effector of TSH to regulate the expression of transport proteins and increase the ability of secretory amplification [55]. Under short photoperiods, the pituitary nodules secrete the thyrotropic hormone, which stimulates T3 in the hypothalamus, which in turn affects the secretion of gonadotropin-releasing hormone and gonadotropin. Studies have shown that OXT and its receptor (OXTR) play a central role in reproduction and metabolism, especially in females, by functioning in an estrogen-dependent manner in rats [56] and are associated with human social and emotional behavior, as well as physical and mental health and disease [57][58][59]. In addition, we performed functional enrichment analysis of target genes of DEGs and DEMs in sheep, respectively. The results of KEGG pathway enrichment in both showed many common pathways in LP and SP comparison, among which the most interesting reproductive pathways included GnRH signaling pathway, Wnt signaling pathway, ErbB signaling pathway, and MAPK signaling pathway. Circadian rhythm-related pathways included circadian and TNF signaling pathways, and these results are consistent with the Kazakh sheep [60], suggesting that reproduction in Small-tailed Sheep is also influenced by the photoperiod. In KEGG analysis, enriched for pathways related to calcium regulation, light influences vitamin D synthesis and thus changes in calcium processing expression. It was shown that the rapid, light-dependent changes in melatonin levels in the pineal organ are apparently regulated by novel calcium signaling pathways [61,62].
We found that among the identified miRNAs, the most highly expressed ones, oar-miR-148a and oar-miR-379-5p, oar-miR-495-3p, oar-miR-143, oar-miR-106b and oar-miR-218a, have been shown to play important roles in regulating GnRH release [63]. In addition, a variety of miRNAs are expressed in the sheep pituitary and at different levels under different light conditions. However, some miRNAs are species-specific, which is also likely to lead to different sensitivity to photoperiods in different animals and different PD secretion of hormones leading to differences in reproductive cyclicity. A total of 6 of the 59 DEMs are known, of which oar-let-7c, oar-let-7a, and oar-let-7b are of the same miRNA family, in animal genomes. The high conservation of let-7 in different animal species suggests that they may play important (and possibly similar) roles in biological processes in various organisms [64,65]. The effect of miR-150 on GC cell apoptosis by targeting the STAR, thereby affecting ovarian function and subsequently hormone secretion related to estrus [66].

Analysis of miRNA-mRNA Interaction Network in SP−PD vs. LP−PD
In the up-regulated network, containing 45 mRNAs and 14 miRNAs negative interactors, of which ZNF787 and CAMK2 are both target genes of novel_146, CAMK2b is one of the most prominent isoforms of CAMK2 [67], and calcium/calmodulin-dependent protein kinase II (CAMK2) is a key player in synaptic plasticity and memory formation. It has been shown that CAMK2b is associated with apoptosis and is involved in the protection of neurons from homocysteine-induced apoptosis through the HIF-1α signaling pathway [68] and that ZNF787 is a neuronal inhibitory molecule [69].
The down-regulated network contained 169 mRNAs and 6 miRNAs. Bok was targeted by novel_142, a pro-apoptotic Bcl-2 protein, and the loss of Bok significantly reduced susceptibility to apoptosis and had an important role in ovarian and follicular maturation and apoptosis in humans [70] and mice [71]. Therefore, we hypothesize that miRNAs are affecting estrus-related hormone activity by inhibiting apoptosis and thus ultimately estrus. Interestingly, novel_414, novel_427 and novel_587 collectively target Lhx1, and Lhx1 deficiency causes severe loss of circadian rhythm and sleeps' light control in mice [72], and in the present study, Lhx1 expression was elevated during long photoperiods, in which we hypothesize that light sensitivity in sheep is reduced through miRNA regulation of Lhx1.
The two pathways that regulate seasonal estrus in sheep, are the hypothalamicpituitary-gonadal axis and KISS1/GPR54. In birds [73] and sheep [74] TSHβ, a molecular signature of long photoperiods, increased in expression in the shift from short to long photoperiods. Under long photoperiods t, high levels of TSHβ bind to TSHR and then drive the conversion of DIO3 to DIO2 through the cAMP signaling pathway, which regulates GnRH neurons, in addition to higher levels of KISS1 expression in the sheep hypothalamus than in the non-breeding season [75]. In the present study, we did a collated comparison of differential genes in the comparison groups of hypothalamus and distal pituitary, and found consistent trends in KISS1 and DIO3 (Supplementary Table S12). In this study, during the short to long photoperiod shift, TSHR expression increased, KISS1 gene expression decreased, DIO2 significantly increased, and DIO3 was decreasing, consistent with the above studies. LHB expression decreased and PRL increased, consistent with the trend of change during the shift from short to long photoperiods in seasonal estrus [60,75], but FSHR expression was slightly decreased, and these changes may also be the reason for the Small Tail Han sheep.

Conclusions
The present study provides, for the first time, the complete miRNA-mRNA network of sheep under photoperiodic regulation from the hypothalamus and pars tuberalis. We identified many DEGs (e.g., LHB, TSHB) and miRNA-mRNA pairs (e.g., LHB is a targeting gene of novel_414) from the RNA-seq data obtained from hypothalamus. We identifiedmany DEGs (e.g., OXT, GH and Lhx1) and miRNA-mRNA pairs (e.g., Lhx1 is a targeting gene of novel_414, novel_427, and novel_587) from the RNA-seq data obtained from pars tuberalis. These data provide significant data for further exploring the reproduction of sheep.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/genes13061091/s1. Table S1: Overview of quality control of mRNA and miRNAs reads. Table S2: Introduction of mRNA. Table S3: Introduction of miRNA. Table S4: Differentially expressed mRNAs in long photoperiods to short photoperiods were identified in hypothalamic and distal pituitary. tissues. Table S5: Differentially expressed miRNAs in long photoperiods to short photoperiods were identified in hypothalamic and pituitary tissues. Table S6: GO and KEGG enrichment annotation for mRNAs in hypothalamus. Table S7: GO and KEGG enrichment annotation for mRNAs in distal pituitary. Table S8: Crossover genes GO and KEGG between predicted target genes of DEGs and DEMs in the hypothalamus. Table S9: Crossover genes GO and KEGG between predicted target genes of DEGs and DEMs in the distal pituitary. Table S10: Construction of miRNA-mRNA Interaction Network in hypothalamus. Table S11: Construction of miRNA-mRNA Interaction Network in distal pituitary. Table S12: Summary of changes in the genes of the hypothalamus and pituitary gland.

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