All the animals were authorized by the Science Research Department (in charge of animal welfare issue) of the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences (IAS-CAAS; Beijing, China). In addition, ethical approval of animal survival was given by the animal ethics committee of IAS-CAAS (No. IAS2018-3, April 10, 2018).

Experiments were conducted on 12 adult Sunite ewes (2–3 years old; weight 30–40 kg), which were selected from a farm in Urat Middle Banner (40° 75′ north latitude), Bayan Nur City, Inner Mongolia Autonomous Region, China, and maintained in a farm in the Tianjin Institute of Animal Sciences, Tianjin (39° 13′ north latitude), China. All ewes were raised under the same conditions, with free access to water and feed. Construction of ovariectomized (OVX) and estradiol-implanted sheep model and light control experiment were previously described in detail (2, 12). Ewes were ovariectomized and estradiol-implanted (E2, Sigma Chemical Co., St. Louis, MO, USA) to maintain plasma estradiol levels of 3–5 pg/ml (6, 16) in October, 2016, according to the model developed by Karsch et al. (1984) (17). This OVX + E2 model normalizes the level of circulating E2 (12), which uncovers the well-documented central seasonal shift in the negative feedback action of E2 on gonadotropin secretion (18). After the surgery, the ewes recovered for 30 days before artificial light control. Then, in November 2016, all ewes were brought indoors and submitted to a light program simulating the outdoor photoperiodic condition by time-control switch. Firstly, all ewes were kept in artificial SP with lights on during 10:30–18:30 (SP, 8:16 h light/dark) for 21 days and switched to LP with lights on during 06:30–22:30 (LP, 16:8 h light/dark) for 42 days, with free access to water and food. Ewes were euthanized [intravenous pentobarbital (100 mg/kg)] at ZT4 [4 h after lights on] of SP21, LP7, LP21, and LP42 (8, 19–21). Consistent with the specific location described by Wood et al. (7) and Lomet et al. (6), the PT tissue (Figure 1) of each ewe was immediately collected and stored at −80°C for total RNA extraction.

Pituitary PT tissues were used for RNA extraction with TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instruction. Examination involving the integrality and quality of the isolated RNA was performed via electrophoresis and the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA).

The rRNA was depleted from 3 μg of total RNA using Ribo-Zero TM Gold Kits (Epicentre, Madison, WI, USA). Sequencing libraries of the 12 samples (SP21, n = 3; LP7, n = 3; LP21, n = 3; LP42, n = 3) were generated using NEB Next Ultra Directional RNA LibraryPrep Kit for Illumina (NEB, Ipswich, MA, USA) according to the manufacturer's instructions, and index codes were used to label the sequences of each sample. After cluster generation, the library preparations were sequenced on an Illumina Hiseq platform (Illumina, San Diego, CA, USA). Raw data of the performed RNA-seq have been recorded in the SRA public database (Accession number of BioProject: PRJNA680667).

Raw data in fastq format were processed through in-house perl scripts. In this step, clean reads were obtained by removing reads with adapter contamination, reads that contained poly-N, and low-quality reads from raw data. Simultaneously, the Q20, Q30, and GC contents of the clean data were calculated. All downstream analysis was based on high-quality clean data. HiSAT2 (22) was used to align clean reads of each sample to the sheep reference genome Oar_v3.1. StringTie (23) was used for transcriptome assembly and reconstruction. Thus, known lncRNA and mRNA transcripts were identified, and the position of transcripts was obtained.

Novel lncRNAs with more than two exons and lengths of more than 200 nt were predicted by CNCI (24), CPC (25), PFAM (26), and CPAT (27) software after transcriptome assembly. The fragments per kilobase per million mapped reads [FPKM (28)] values were calculated to represent the expression of genes and lncRNAs. To determine the effect of photoperiod on genes and lncRNAs expression in PT of sheep, the expression of genes and lncRNAs at every time point of LP (LP7, LP 21, or LP 42) was compared with the SP21 group using DESeq, i.e., SP 21 vs. LP 7, SP 21 vs. LP 21, and SP 21 vs. LP 42. In addition, p < 0.05 and |Fold change| >1 was considered standard of differential expression between the SP and LP.

Gene Ontology (GO) enrichment analysis of differentially expressed genes or lncRNA target genes was implemented by the GOseq R package, in which gene length bias was corrected (29). GO classifies functions into three groups: cellular components, molecular functions, and biological processes. The KEGG biological pathways database (http://www.genome.jp) is a central public database for understanding high-level functions and regulatory network research. Enrichment analysis was performed on each pathway in KEGG using a hypergeometric test. According to significant threshold (p-value: 0.05), the genes were screened and enriched for the pathways. Next, the significance of the pathway enrichment analysis was corrected by FDR, and the corrected p-value (q-value) was obtained.

The primary role of lncRNAs, which are a type of non-coding RNA, is to regulate their target genes by cis-regulating nearby protein-coding genes and trans-regulating distal protein-coding genes. Here, protein coding genes with a distance of <100 Kb were assumed to be the cis-target genes, and Pearson correlation coefficients with the lncRNAs of >0.95 were assumed to represent the trans-target genes (30).

To further reveal the potential roles of lncRNAs that are involved in modulating the reproductive process, integral interaction networks containing lncRNAs and their corresponding target genes were built using Cytoscape software (31) in each comparison (SP 21 vs. LP 7, SP 21 vs. LP 21, and SP 21 vs. LP 42), which include cis- and trans-forms of regulation.

EYA3, TSHB, CHGA, MSTRG.209166, MSTRG.22781, and MSTRG.42035 were selected to validate the accuracy of RNA sequencing via the reverse-transcription quantitative polymerase chain reaction (RT-qPCR). GAPDH was used as an internal reference to normalize target gene expression. All primers used in the RT-qPCR are shown in Supplementary Table 1. cDNA was used to perform RT-qPCR after reverse transcription from total RNA to cDNA. The RT-qPCR reaction conditions were as follows: 95°C for 15 min, followed by 40 cycles of 95°C for 10 s and 60°C for 30 s. The data obtained from RT-qPCR reaction were then calculated using the 2^−ΔΔCt method (32, 33) and processed by SPSS 19.0 with a one-way analysis of variance. The results are presented as means ± standard deviation. Furthermore, p < 0.05 was regarded as statistically significant.

To identify differentially expressed lncRNAs and genes between SP and LP, RNA libraries of different photoperiods were constructed. After removing low-quality sequences, a total of 1,394,924,578 clean reads with >92.52% of Q30 were obtained after sequencing all 12 libraries. Approximately 92–95% of the reads were successfully aligned to the Ovis aries reference genome (Oar_v3.1) (Table 1).

Overall, a total of 55,472 lncRNAs were identified in PT of 12 ewes using four programs (CNCI, CPC, PFAM, and CAPT) (Figure 2A). The length distribution of lncRNAs was consistent with that of protein-coding gene (Figure 2B). Of them, the transcripts of lncRNAs and mRNAs with lengths of more than 3,000 bp accounted for the majority, and its number in lncRNAs was significantly greater than that in mRNAs (Figure 2B). However, transcript levels of lncRNAs were lower than those of mRNAs in the PT of Sunite ewes (Figure 2C). Most of the lncRNAs have only two or three exons, whereas mRNAs contain a wide range of exons from 2 to 30 (Figure 2D).

The number of differentially expressed (DE) genes between SP and LP increased gradually with the extension of LP (from LP7 to LP42). Specifically, there were 499 DE mRNAs in SP21 vs. LP7, in which FOS, SOCS3, EYA3, TSHB, SIX1, DCT, VMO1, AREG, and EZH2 were related to reproduction (Figure 3A). In SP21 vs. LP21, 625 DE mRNAs were detected, including CHGA, FOS, SOCS3, GHRH, TH, EYA3, TSHB, DCT, VMO1, SUV39H2, and EZH2 that were associated with reproduction (Figure 3B). In addition, 874 DE mRNAs were screened between the SP21 and LP42 group. Of them, the genes related to reproduction included CHGA, FOS, SOCS3, GHRH, TH, EYA3, TSHB, SIX1, VMO1, AREG, and EZH2 (Figure 3C). The top 10 DE mRNAs in each contrast group are shown in Table 2. Among these genes, there are specific differentially expressed genes (FOS, TMEM200A, MLIP, PMP2, PRX, and SPTBN5 and seven novel genes) that are different from previous reports (6, 7). Besides, for the top genes, it is worth paying attention to those genes that appear in the above contrast groups simultaneously, and they are the most likely photoperiodic-induced factors in the regulation of seasonal reproduction. To sum up, the notable LP-induced genes included EYA3, TSHB, VMO1, DCT, and EZH2. SP-induced genes involved ENSOARG00000012585, CHGA, FOS, SOCS3, and TH. For lncRNAs, the number of DE lncRNAs was also increasing with the extension of the LP. Specifically, 1,887, 2,038, and 2,614 DE lncRNAs were detected in the SP21 vs. LP7 group (Supplementary Table 5), SP21 vs. LP21 group (Supplementary Table 6), and SP21 vs. LP42 group (Supplementary Table 7), respectively. These DE lncRNAs may be involved in the regulation of the expression of important genes for seasonal reproductive trait. Besides, our results had some common overlap DE mRNAs with the reports of Wood et al. (2015) and Lomet et al. (2017) (6, 7). The detailed information is shown in Supplementary Table 8.

GO and KEGG analyses were conducted for the differentially expressed mRNAs and target genes of differentially expressed lncRNAs. The top 10 enriched GO terms for lncRNAs and mRNAs of SP21 vs. LP7 (blue), SP21 vs. LP21 (orange), and SP21 vs. LP42 (green) are shown in Figures 4A,B (data in Supplementary Tables 9–11 for lncRNAs and Supplementary Tables 12–14 for mRNAs). The following GO terms are noteworthy, including neuron terms (such as neuron part, axon, synapse part, and response to stimulus), receptor terms (such as cell surface receptor signaling pathway and receptor complex), regulation of apoptosis process, regulation of cell death, regulation of cell proliferation, and metabolic process (such as regulation of metabolic process and positive regulation of phosphorus metabolic process). In order to further focus on the function of differentially expressed genes between SP and LP, we constructed a simplified network of related statistically significant GO terms (Figure 4C) for differentially expressed genes overlapping between SP and three of LP points using the Cytoscape add-on ClueGO (34, 35). The results emphasized enrichment for cell proliferation and differentiation, cell response and pathway, cell death and regulation, and the cytokine process, which might play an important role in the regulation of seasonal breeding in sheep at the level of the PT/MBH (36, 37).

The top 20 enriched KEGG pathways for target genes of differentially expressed lncRNAs and pathways for differentially expressed mRNAs are shown in Figures 5, 6, respectively. The most interesting pathways included reproduction-related pathways (Prolactin signaling pathway, ErbB signaling pathway, Wnt signaling pathway, MAPK signaling pathway, and GnRH signaling pathway in Figure 5; ErbB signaling pathway and MAPK signaling pathway in Figure 6), circadian related pathways (phototransduction-fly in Figure 5; circadian entrainment in Figure 6), and neuroactive ligand–receptor interaction, TNF signaling pathway, tyrosine metabolism, and neuroactive ligand–receptor interaction pathway (detailed data in Supplementary Tables 15–17 for lncRNAs and Supplementary Tables 18–20 for mRNAs). Here, we also combined the overlap DE genes from three comparisons and assigned KEGG terms to create a simplified network of related statistically significant KEGG terms (34, 35) (Figure 7). Figure 7 also showed that many differentially expressed genes between SP and LP are enriched in pathways related to reproduction, such as signaling (prolactin signaling pathway, cAMP signaling pathway, and JAK-STAT signaling pathway), neuron, neurotransmitters, and cell process.

The differentially expressed lncRNAs and their target genes were selected to construct the lncRNA–mRNA network with co-expression information. In the network of SP21 vs. LP7, 23 differentially expressed lncRNAs and 21 target genes were included, among which the trans relationship between MSTRG.209166 and TSHB was hinted (Figure 8A, Supplementary Table 21). In SP21 vs. LP21, 29 differentially expressed lncRNAs and 26 target genes were selected to construct network. Notably, MSTRG.209166 was also found to trans-regulate TSHB and MSTRG.235014 cis-regulated DDC, because the two genes were important for seasonal reproduction (Figure 8B, Supplementary Table 22). In SP21 vs. LP42, 102 differentially expressed lncRNAs and 91 target genes formed a more sophisticated network. For several key genes of seasonal reproduction, several target combinations between lncRNAs and them were highlighted. For example, they implied that MSTRG.209166 trans-regulated TSHB, MSTRG.288068 cis-regulated SIX1, MSTRG.272793 cis-regulated KIT, and ENSOARG00000026131 cis-regulated TH (Figure 8C, Supplementary Table 23).

A total of six genes, namely, three mRNAs (EYA3, TSHB, and CHGA) related to seasonal reproduction and three random lncRNAs (MSTRG.290436, MSTRG.22781, and MSTRG.17707), were selected for RT-qPCR verification. The results indicated that there is a similar expression pattern between RNA-Seq and RT-qPCR data (Figure 9).