A group of 21 Sunite ewes (35–40 kg, 3 y, clinically normal, and non-pregnant) were selected from Urat Middle Banner, Bayan Nur City, Inner Mongolia Autonomous Region, China and housed in a farm at Tianjin Institute of Animal Sciences, Tianjin (39°N latitude), China. All ewes were fed ad libitum and had free access to water. OVX + E₂ sheep and light control rooms were constructed as previously described (La, et al., 2020; He et al., 2021; Xia, et al., 2021). In brief, estradiol treatment was achieved with an inner diameter of 3.35 mm and an outer diameter of 4.65 mm, packed with 20 mg crystalline 17β-estradiol (Sigma Chemical Co., St. Louis, MO). Implants were inserted into the axillary region and designed to produce the circulated E₂ levels of approximately 3–5 pg/ml for 2 weeks (Smith, et al., 2007). Then, all ewes were maintained in one of three photoperiod-controlled rooms, SP: short photoperiod (8/16 h light–dark), LP: long photoperiod (8/16 h light–dark) and SP–LP: photoperiod from short transfer to long). Along with the photoperiodic treatments, slaughter was performed on SP42 or LP42, SP-LP3 (the third day of SP–LP), SP-LP7, SP-LP15, SP-LP21, and SP-LP42 (n = 3). After slaughtering, the hypothalamus was rapidly removed from the brain, flushed with PBS (pH 7.4), snap-frozen at liquid-nitrogen, and stored at −80°C for subsequent study. Simultaneously, jugular venous blood samples were collected twice weekly (24 times, 12 times for different photoperiods) from three ewes in the SP–LP room. After centrifugation, serum was collected and frozen until assayed for hormones. The hormone preparation and sample collection’s schematic diagrams are shown in Figure 1.

The iodine (¹²⁵I) human follicle-stimulating hormone radioimmunoassay kit and the iodine (¹²⁵I) prolactin radioimmunoassay kit (Beijing North Biotechnology Research Institute Co., Ltd., Beijing, China) were used for FSH and PRL assays by radioimmunoassay (RIA). All samples from different time points were measured in triplicate. The FSH kit had an assay sensitivity of <1.0 mIU/ml, and the PRL kit sensitivity was ≤40 uIU/ml.

Total RNA was isolated from each sample using TRIzol reagent (Invitrogen, Carlsbad, CA, United States), and 1% agarose gels were used to detect the degradation and contamination prior to subsequent. In addition, RNA’s purity, concentration, and integrity were also detected by a NanoPhotometer^® spectrophotometer (IMPLEN, CA, United States), a Qubit^® RNA Assay Kit in Qubit^® 2.0 Fluorometer (Life Technologies, CA, United States), and the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, United States), respectively. Subsequently, libraries were generated using rRNA-depleted RNA with the NEBNext^® Ultra™ Directional RNA Library Prep Kit for Illumina^® (NEB, Ipswich, MA, United States) according to the manufacturer’s recommendations. Finally, the libraries were sequenced on an Illumina HiSeq 4000 platform, and 150 bp paired-end reads were generated.

Clean data were obtained by removing reads-containing adapter, reads-containing poly-N sequence, and low-quality reads from raw data. Simultaneously, the Q20, Q30, and GC contents of the clean data were calculated. All the following analyses were based on the clean data. The Ovis aries reference genome and gene model annotation files were downloaded from the genome website directly (Oar_v4.0). The index of the reference genome was built using bowtie2 (v2.2.8) (Sirén, et al., 2014). Reads were aligned to the reference genome using HISAT2 (v2.0.4). The mapped reads of each sample were assembled by StringTie (v1.3.1) (Pertea, et al., 2016).

To improve the statistical validity of the identification, a series of strict screening conditions were employed to identify lncRNAs. First, we removed transcripts that had fewer than two exons, and transcripts of less than 200 bp in length were discarded. The remaining transcripts were BLAST to known sheep lncRNAs using Cuffcompare, and lncRNAs that overlapped with the transcripts in the database were incorporated into the subsequent analysis. Next, transcripts with an FPKM ≥ 0.5 were retained after calculating the expression of each transcript by Cuffdiff. Finally, CNCI (Sun, et al., 2013), CPC (Kong, et al., 2007), and Pfam (El-Gebali, et al., 2019), three databases were used to identify the candidate lncRNAs. The intersections of the results from each program were defined as the novel lncRNA transcripts, and those identified by all three tools were defined as the final candidate lncRNAs for further analysis.

Cuffdiff (v2.1.1) was used to calculate FPKMs, which represent the gene expression level of both the lncRNAs and coding genes in each sample (Trapnell, et al., 2010). FPKMs were computed by summing the FPKMs of transcripts in each gene group. A criterion of absolute log2 (fold change) > 1 and p ≤ 0.05 between the two groups was used to identify differentially expressed genes in Cuffdiff.

Eight transcripts (four mRNAs and four lncRNAs) were selected for each group to validate the RNA-sequencing results by qRT-PCR. The specific primers are listed in Table 1 cDNA was synthesized from the RNA samples remaining after sequencing with the PrimeScript™ RT reagent Kit (TaKaRa, Dalian, China). qPCR was performed with a SYBR Green assay (TaKaRa, Dalian, China) on a Roche LightCycler 480 (Roche Applied Science, Mannheim, Germany). The qPCR mixture and program has been described in our previous study (Zhang, et al., 2019). All samples were examined in triplicate. Differentially expressed mRNA and lncRNA levels were normalized to β-actin to determine the relative expression using the 2^−ΔΔCt method (Livak and Schmittgen, 2000).

To further elucidate the functions of DE-mRNAs and DE-lncRNAs, we predicted the target genes of the lncRNAs in cis and trans. To understand the functional roles of the target genes of the lncRNAs, we used GOseq R package to implement enrichment analysis. In addition, the differentially expressed protein-coding genes were also analyzed using GO. We employed KOBAS software to detect the enrichment of lncRNA target genes or differentially expressed genes in KEGG pathways. GO terms with -log10 (p-value) < 0.05 and KEGG pathways with p < 0.05 were considered significantly enriched.

To infer the functions of the photoperiod-based DE-lncRNAs and DE-mRNAs in the ovine hypothalamus, a protein–protein interaction (PPI) network between the protein-coding genes was constructed based on information from STRING. Additionally, we constructed a complementary pair network based on mRNA and mRNA as well as between mRNA and lncRNA by using Cytoscape (V3.6.1).

The statistical evaluation of the experimental results was performed with Student’s t-test using SPSS 20.0 statistical software. All data are expressed as the means with standard error (SE). For hormone profiles, RIA data were analyzed by repeated measures ANOVA, and comparisons of hormones between different photoperiods and qPCR validation were performed by one-way ANOVA. p < 0.05 was considered statistically significant.

RIA was used to assess the serum hormone levels of FSH and PRL in three ewes exposed to the changed photoperiod. The FSH data shown in Figure 2A mostly reflected the impact of the photoperiod, and the serum hormone levels differed significantly between SP and LP (p < 0.05) (Figure 2B). For PRL, the concentration increased with the photoperiodic transition from SP to LP (Figure 2C), and the difference also reached significant level (p < 0.01) (Figure 2D).

Twenty-one hypothalamus samples from seven groups (SP42, LP42, SP-LP3, SP-LP7, SP-LP15, SP-LP21, and SP-LP42, n = 3) were used to construct 21 cDNA libraries for sequencing. More than 268 million raw reads were generated for each group. The Q30 of the data was not less than 89%, and the GC contents of each library ranged from 49.34 to 55.81%. In addition, 303.38 Gb clean reads were retained and used in the following analysis after discarding reads with adapters or a poly-N content >10% and other low-quality reads. Approximately, 97% of the clean reads from each library were mapped to the sheep reference genome (Oar_v4.0) (Table 2), and 552,061 transcripts were assembled from these 21 libraries, which were assembled using Scripture and Cufflinks. After a rigorous 5-step screening, 11,965 putative non-coding transcripts were retained (Figure 3A), and the intersection of CNCI, CPC, and Pfam yielded the same result (Figure 3B), identifying 9,645 lincRNAs (80.6%), 2,320 antisense lncRNAs (19.4%), and 0 intronic lncRNAs (0.0%).

In this study, 13,859 lncRNAs and 44,437 mRNAs were identified. Our results indicated that the length and ORFs of lncRNAs were shorter than those of mRNAs (Figure 4A); the expression of lncRNAs was lower than that of mRNAs (Figure 4D); and lncRNAs tended to contain fewer exons (Figure 4B). These results were consistent with those of our previous studies (Zhang, et al., 2019). In addition, some transcripts of uncertain coding potential (TUCP) were expressed in the hypothalamus (Figure 4C).

The criteria of |log2 (fold change) |>1 and p ≤ 0.05 were used to identify differentially expressed transcripts. The summary of DE transcripts is shown in Table 3, and the details of the differences in expression can be found in Supplementary Table S1 and Supplementary Table S2. Importantly, the differences between different phases of the photoperiod change need to be verified prior to further analysis. In this study, four DE-mRNAs and four DE-lncRNAs were selected from each compared group to validate the RNA-seq data, and the expression of these RNAs was analyzed by qPCR (Supplementary Figure S1). qPCR results showed that the expression levels of selected DE-mRNAs (PLCB1, PRKACB, SMAD4, BCL2, ATF4, RASD1, CLOCK, CREB1, IGF2, DI O 2, BHLHE41, LHB, ESR1, GNB2, TSHB, FBXL3, and PRL) from the six compared groups were all significantly different (Supplementary Figure S1A). Similarly, the expression of 19 DE-lncRNAs (LNC_004953, LNC_002738, LNC_008255, LNC_008321, XR_173415.3, XR_173257.3, XR_001043724.2, XR_001045234.2, XR_001023464.2, XR_001020627.1, XR_001034881.2, XR_001433797.1, XR_001433798.1, XR_001024596.2, XR_001434242.1, XR_001027187.1, XR_001039609.1, XR_001023520.1, and XR_001434385.1) also reached the level of significant difference (Supplementary Figure S1B), and the expression patterns of these genes and lncRNAs were consistent with the RNA-seq results.

To understand the functions of DE-mRNAs and DE-lncRNAs in greater depth, DE-mRNAs and DE-lncRNA-colocalized mRNAs were subjected to GO and KEGG pathway analyses. The top five significantly enriched GO terms were in molecular function (MF), cellular component (CC), and biological process (BP),three GO types for DE-mRNA. Compared with SP42, when the photoperiod changes from SP-LP3 to SP-LP42, blinding, RNA binding, and cell cycle items were significantly enriched (Figure 5, data in Supplementary Table S3). For DE-lncRNA, we used their targets to conduct GO enrichment and most of the significantly enriched GO terms were consistent with mRNAs. In addition, most of the significantly enriched GO terms participate in the regulation of biological and cellular processes (Figure 6, data in Supplementary Table S4). KEGG pathway analysis revealed that these differentially expressed mRNA target genes were enriched in several significant items, such as steroid hormone biosynthesis, dopaminergic synapse, and circadian entrainment (Figure 7A, data in Supplementary Table S5), KEGG pathway analyses of differentially expressed lncRNAs were also enriched in pathways mentioned earlier. In addition, insulin secretion, the thyroid hormone signaling pathway, and the Wnt signaling pathway, were also significantly enriched with the photoperiod changes (Figure 7B, data in Supplementary Table S6).

Furthermore, analyses of DE-mRNAs revealed that the key signaling pathways regulating ovine photoperiodic responsiveness and seasonal reproduction, including the thyroid hormone signaling, insulin secretion, Wnt signaling, circadian entrainment, and GnRH signaling pathways were detected after different lengths of photoperiodic treatment (Table 4). For example, compared with LP42, except for PRKACB and PLCB1, all genes were downregulated in the thyroid hormone signaling pathway, and several genes in circadian entrainment (PER1, RASD1, and PRKACB) were significantly changed with the photoperiod from SP to LP. In addition, with the photoperiod change from SP to LP and gradually from SP-LP3 to SP-LP42, eight photoperiod-related genes, namely, PRKACB, AKT3, ATF4, KCNJ5, CLOCK, PLCB1, SMAD4, and MAPK3, were upregulated in SP42, and 11 other photoperiod-related genes (VAMP2, BTRC, PER1, GNB2, MTNR1A, CDC42, ATF4, LHB, TSHB, ESR1, and RASD1) were downregulated. Information on all lncRNAs and their target genes is listed in Supplementary Table S7. In this study, we searched for target protein-coding genes of the differentially expressed lncRNAs to evaluate the potential regulatory functions of lncRNAs, and a total of 130, 75, 132, 101, 145, and 260 target genes of DE-lncRNA were predicted and classified (Supplementary Table S8). Interestingly, we detected some genes related to animal photoperiodic response and estrous seasonality, such as AANAT, TH, CAMK2G, and KCNJ5, which were in near to XR_173257.3, XR_173415.3, XR_001024596.2, and XR_001023464.2, respectively (Table 5).

Several transcripts were selected to construct mRNA–mRNA and lncRNA–mRNA networks impacted by the change of light from SP to LP (Figure 8). For the mRNA–mRNA interaction network, DE-mRNAs are grouped to insulin secretion (VAMP2, PRKACB, PRKCG, and PLCB1), thyroid hormone synthesis (TSHB, CGA, ATF4, and PRKACA), and estrogen secretion (AKT3, ATF4, MAPK3, and ESR1). DE-mRNAs-regulating hypothalamic hormone secretion and photoperiodic responsiveness were also enriched in GnRH (MAPK13, CAMK2D, CGA, CDC42, ATF4, and LHB), Wnt (SMAD4, RHOA, RAC1, and SFRP5) signaling pathways, and circadian entrainment (PRKACB, KCNJ5, PER1, GNB2, MTNR1A, and RASD1) (Figure 8A). Several DE-lncRNAs related to ovine hypothalamic hormone secretion and photoperiodic responsiveness were also selected to construct an lncRNA–mRNA interaction network (Figure 8B). Of these, DE-lncRNAs were highly correlated with steroid hormone biosynthesis, retinol metabolism, and circadian entrainment. The lncRNAs XR_173257.3 targets AANAT, and is potentially involved in the response to photoperiodic change. XR_173415.3 targets TH and IGF2; interestingly, thyroid hormones, encoded by TH, have been proven to play an important role in mammalian seasonal reproduction. In addition, XR_001435315.1 and XR_001023464.2 probably participate in ovine reproductive activity because of their target genes B4GALNT2 (FecL) and KANJ5.

Several lines of evidence indicate that there is a cyclical rhythm in ovine reproduction, and its activation and quiescence can be clearly seen in gonadotrophin secretion in OVX + E₂ ewes (Lincoln, et al., 2006; Menassol, et al., 2012; Dardente, et al., 2019; Guh, et al., 2019). In this study, FSH and PRL were selected as the indicators of the regulation of gonadotropin secretion by photoperiodic changes in OVX + E₂ ewes. For FSH, the concentration was lowest in the anestrous period and significantly increased in the estrus period; Lincoln (2003) found that the concentration of ovine FSH was very low under LP, which is completely consistent with our findings. This means that the release of FSH was also regulated by photoperiod change or season, and a low concentration of FSH may be related to anestrus in ewes. PRL and FSH display opposite responses to the photoperiod, because they are driven by two different neuroendocrine axes (Dardente, 2012), one or several endocrine factors in the PRL axes have been proved to regulate the photoperiodic rhythms (Lincoln, et al., 2006; Dardente, et al., 2019), our findings also demonstrate that the concentration of PRL increases significantly upon photoperiod extension.

In the current study, most omics studies of the hypothalamus have focused on differences in gene expression profiles based on the different estrus periods or some important economic traits with significant differences in animals (Gao, et al., 2018; Hou, et al., 2018; Gley, et al., 2019; Zhang, et al., 2019; Liang, et al., 2021). For example, it was found that the seasonal energy balance of Siberian hamsters was affected by the species-specific responsiveness of hypothalamic T3 (Bao, et al., 2019), and that the hypothalamus was also involved in the circadian of metabolism (Cedernaes, et al., 2019) and daylength response (Sáenz De Miera et al., 2017). Recent studies by David emphasized that photoperiodic transitions in endocrine output are regulated by a multistep signaling cascade within the mammalian neuroendocrine system (Hazlerigg, et al., 2018), and the differences of hypothalamic transcriptome levels can provide more information about photoperiodic response in ewes. In this study, compared with SP42, most DE-mRNAs in the SP-LP3 to SP-LP42 groups compared with the SP groups were associated with the signaling of circadian entrainment, insulin secretion, thyroid hormone, and estrogen (Table 4). As previously reported, SCN lesions lead to the alterations in the estrous cycle, sexual behavior, and tonic and phasic secretion of GnRH and gonadotropins in rodents, and it has been proven that Clock and Clock-related genes are the key factors in mammalian reproductive cycles (Silva and Domínguez, 2019). In our studies, the expression of Clock-related genes, such as MTNR1A, PER1, DI O 2, and ESR1, was changed with the artificial photoperiod from SP42 to LP21 and LP42, and was synchronously downregulated with gonadotropins FSH. Importantly, previous studies have demonstrated that these genes were significantly correlated with reproductive seasonality in goats and sheep (Wood, 2018; He, et al., 2019). Astoundingly, hypothalamus mRNAs that were altered in the SP–LP groups (SP-LP7, SP-LP21, and SP-LP42) compared to the SP42 group were all associated with circadian entrainment (Table 4), which was proven to initiate the maturation of human islets and insulin secretion (Alvarez-Dominguez, et al., 2019). Moreover, recent evidence has shown that insulin plays a vital role in the arrangement of reproductive activities through direct actions on insulin receptors; coincidentally, these receptors are present in several areas that are well-known to play a crucial role in reproduction, such as the ARC, ventromedial hypothalamic nucleus (VMH), and preoptic area (POA) (Caraty, et al., 1998; Brüning, et al., 2000; Cernea, et al., 2016). Several genes that were significantly differentially expressed upon photoperiod changes in our studies participate in insulin secretion (VAMP2, PRKACB, and PLCB1) and circadian entrainment (PRKACB, RASD1, and KCNJ5), and the change of previously mentioned key genes in circadian entrainment triggers insulin secretion; subsequently, the receptors in various hypothalamic functional areas are activated to regulate the ovine photoperiod-induced reproductive activities.

At present, there is vast availability of evidence suggesting that the expression of lncRNAs can show high correlations with the expression of neighboring mRNAs. As stated in the introduction, lncRNA GnRH-E1 RNA targets its nearby gene GnRH1 to contribute to the regulation of GnRH neuronal maturation (Huang, et al., 2016). Furthermore, XLOC_1041225 and XLOC_446331 target neighboring IGFBP5, and regulate the onset of puberty in goats (Gao, et al., 2018). The ovine hypothalamic lncRNAs found in the present study had many similar characteristics with those of other mammals, such as shorter and fewer exons, shorter ORFs, and lower expression levels than protein-coding transcripts (Figure 4). As discussed in 4.3, genes involved in circadian entrainment and insulin secretion participate in photoperiod-induced reproduction in sheep. Several lncRNAs target genes were associated with the signaling of circadian entrainment (KCNJ5 and CAMK2G), prolactin (TH and AKT1), steroid hormone biosynthesis (UGT1A1), and tryptophan metabolism (AANAT) (Table 5 and Figure 6). As previously reported, the AANAT gene is an essential enzyme for the rhythmic synthesis of melatonin (Maronde, et al., 2011), and TH can also be connected with AANAT using dopa decarboxylase (DDC). Combined with the earlier studies, our results indicated that the XR_173257.3-AANAT and XR_173415.3-TH networks may be pivotal in rhythmic melatonin synthesis and secretion and subsequent reproductive seasonality in sheep. Moreover, KCNJ5, also known as a G protein-activated inward rectifier potassium channel 4 (GIRK4), has been confirmed to be associated with circadian rhythmicity in humans (Hsiao, et al., 2019). From our results, XR_001023464.2-KCNJ5 synchronous expression may be involved in the ovine circadian rhythm induced by photoperiod changes. Interestingly, we also detected the differentially expressed lncRNA XR_001435315.1, in which the target gene B4GALNT2 was proven to be associated with the ovine litter size in our earlier study (Guo, et al., 2018). All these findings suggest that lncRNAs may be a regulator involved in the photoperiodic response in sheep.