Identification of Photoperiod-Induced LncRNAs and mRNAs in Pituitary Pars Tuberalis of Sheep

The pituitary pars tuberalis (PT) is the regulating center of seasonal reproduction, which can sense the melatonin signal and eventually cause downstream changes of GnRH secretion through TSHβ. Recently, lncRNAs have been identified in animal reproductive-related tissues, and they play important roles in reproductive regulation. Therefore, in this study, we expect to identify photoperiod-induced lncRNAs and genes in pituitary PT of sheep by comparison of expression profiles between short photoperiod (SP) and long photoperiod (LP). Through RNA-Seq, a total of 55,472 lncRNAs were identified in pituitary PT of Sunite ewes. The number of differentially expressed (DE) genes and lncRNAs between SP and LP increased gradually with the extension of LP (from LP7 to LP42). The notable LP-induced candidate genes included EYA3, TSHB, SIX1, DCT, VMO1, AREG, SUV39H2, and EZH2, and SP-induced genes involved ENSOARG00000012585, CHGA, FOS, SOCS3, and TH. In enriched pathways for DE genes and lncRNA target genes between SP and LP, the reproduction- and circadian-related pathways were highlighted. In addition, the interactome analysis of lncRNAs and their targets implied that MSTRG.209166 and its trans-target TSHB, MSTRG.288068 and its cis-target SIX1, and ENSOARG00000026131 and its cis-target TH might participate in regulation of seasonal reproduction. Together, these results will help to determine important photoperiod-induced lncRNAs and genes and give us some new insights into the epigenetic regulation of seasonal reproduction in sheep.

Introduction

The reproductive activity of some animal species inhabiting the temperate zone is limited in specific seasons in order to maximize the survival possibility of their offspring. According to the different seasons of breeding, these animals are categorized as long photoperiod (LP) breeders and short photoperiod (SP) breeders (1, 2). Among them, sheep belongs to the SP breeder (3); that is, they are mostly bred in autumn and winter. For example, Sunite sheep in China exhibit obvious seasonal reproductive behavior throughout the year, i.e., estrus from August to March of the next year and anestrus from April to July (2, 4). Seasonal reproduction is an important factor limiting the production effciency of the sheep industry; therefore, research on the molecular basis of seasonal reproduction of sheep is indispensable for possible artificial regulation of this trait in the future. However, the molecular mechanism and regulatory network of seasonal reproduction are not very clear until now.

Pituitary pars tuberalis (PT) plays an important role of transmission center in the seasonal reproduction of animals. The pineal gland first converts external photoperiod signals into biological signals (nocturnal melatonin secretion), which causes photoperiod-induced signals (EYA3, TSHB, etc.) alternation in PT. Thus, the thyrotrophin secreted by the anterior pituitary will increase in long days, which acts on TSH receptor-expressing cells in the adjacent mediobasal hypothalamus, leading to type III thyroid hormone deiodinase (DIO3) switch to type II thyroid hormone deiodinase (DIO2). DIO2 regulates the thyroid hormone in the hypothalamus, controlling the activity of the hypothalamus–pituitary–gonad axis to exhibit the summer phenotypes by direct (GnRH neuron) or indirect (KISS1/RFRP system) ways (5, 6). Besides, the photoperiod-induced gene expression in PT is mainly affected by photoperiod, independently of the TH status (6). However, molecular changes related to seasonal reproduction in medio-basal hypothalamus are regulated not only by photoperiod but also by TH. Therefore, PT is the best ideal tissue for seasonal reproduction analyses. So far, several important genes involved in the regulation of seasonal reproduction have been discovered in PT. For example, LP-induced genes (e.g., EYA3 and TSHβ) and an SP-induced gene (CHGA) were revealed in PT of sheep (3, 7, 8). Their expression levels had a remarkable change in different photoperiods. In recent years, long non-coding RNAs (lncRNAs) are considered as key regulators of gene expression because they play crucial roles in transcriptional and post-transcriptional regulation. In sheep, many lncRNAs have been identified in reproduction-related tissues, such as ovary, uterus, and germ cells (9–12). Moreover, lncRNAs play important roles in many aspects of sheep reproductive regulation, such as fecundity (12, 13), gonadal development (14), and sex hormone response (15). However, the functions of lncRNAs in animal PT on seasonal reproduction are unknown. Therefore, one purpose of this research is to seek the new photoperiod-induced genes and lncRNAs in PT of sheep by transcriptome sequencing. Another objective is to explore the potential relationship between genes and lncRNAs by target prediction and joint analysis of their co-expression in PT of sheep. These results will give us some new insights into the epigenetic regulation of seasonal reproduction in sheep.

Methods

Ethical Statement

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).

Animal and Tissue Acquirement

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.

FIGURE 1

Figure 1. The location of PT tissue in Sunite sheep.

RNA Extraction, Library Construction, and Sequencing

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).

Reference Genome Mapping and Transcriptome Assembly

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.

LncRNA Identification and Differentially Expression Analysis

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.

GO and KEGG Pathway Enrichment Analysis of Differentially Expressed Genes and lncRNAs

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.

Construction of Integral lncRNA–mRNA Interaction Networks

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.

RNA-seq Data Validation

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.

Results

Summary of Transcriptome Sequencing Data

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).

TABLE 1

Table 1. Summary of raw reads after quality control and mapping to the reference genome.

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).

FIGURE 2

Figure 2. Expression characteristics of lncRNAs and genes in the PT of ewes. (A) Venn diagram for screening lncRNAs by four software (CNCI, CPC, CAPT, and PFAM). The sum of the numbers in each large oval represents the total number of non-coding transcripts from each software, and the overlapping parts of the oval represent the non-coding transcripts common to the software. (B) The length distribution of lncRNAs and mRNAs. (C) The expression level of lncRNAs and mRNAs. (D) The exon number distribution of lncRNAs and mRNAs.

Photoperiod-Induced Factor Analysis Through Differentially Expressed Genes and lncRNAs

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.

FIGURE 3

Figure 3. The volcano plots of differentially expressed genes between long photoperiod (LP) and short photoperiod (SP) in PT of Sunite sheep. (A) The results of LP7 vs. SP21; (B) the results of LP21 vs. SP21; (C) the results of LP42 vs. SP21.

TABLE 2

Table 2. The top 10 differentially expressed mRNAs in three comparisons between LP and SP.

GO Annotation and KEGG Enrichment Analysis of DE Genes and lncRNAs

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).

FIGURE 4

Figure 4. Top enriched Gene Ontology (GO) terms of differentially expressed lncRNAs and genes between short photoperiod and long photoperiod in PT of Sunite sheep. (A) The top GO terms of target genes of differentially expressed LncRNAs; (B) the top GO terms of differentially expressed mRNAs. The GO terms in blue are from SP21 vs. LP7; those in orange are from SP21 vs. LP21, and those in green are from SP21 vs. LP42. (C) A simplified network of related statistically significant GO terms (p < 0.01) using the Cytoscape add-on ClueGO (34, 35). The network of overlap DE genes among SP21 vs. LP7, LP21, and LP42 are shown. The filled colored circles (nodes) represent each statistically significant parent GO term. The lines (edges) between the nodes show that there are overlapping genes within each term. The colored ovals group these parent GO terms into more generic functional descriptions. The red boxes represent the seasonal reproduction-related terms.

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.

FIGURE 5

Figure 5. KEGG pathway analysis for target genes of differentially expressed lncRNAs between short photoperiod and long photoperiod in PT of Sunite sheep. (A) Top 20 enrichment pathways in SP21 vs. LP7. (B) Top 20 enrichment pathways in SP21 vs. LP21. (C) Top 20 enrichment pathways in SP21 vs. LP42. Rich_factor is defined as the amount of differentially expressed genes enriched in the pathway/amount of all genes in the background gene set. The red boxes represent the seasonal reproduction-related terms.

FIGURE 6

Figure 6. KEGG pathway analysis of differentially expressed mRNAs between short photoperiod and long photoperiod in PT of Sunite sheep. (A) Top 20 enrichment pathways in SP21 vs. LP7. (B) Top 20 enrichment pathways in SP21 vs. LP21. (C) Top 20 enrichment pathways in SP21 vs. LP42. Rich_factor is defined as the amount of differentially expressed genes enriched in the pathway/amount of all genes in the background gene set. The red boxes represent the seasonal reproduction-related terms.

FIGURE 7

Figure 7. A simplified network of related statistically KEGG pathways using the Cytoscape add-on ClueGO (34, 35). The network of overlap DE genes among SP21 vs. LP7, LP21, and LP42 are shown. The filled colored circles (nodes) represent each statistically parent KEGG pathway. The lines (edges) between the nodes show that there are overlapping genes within each pathway. The colored ovals group these parent KEGG pathways into more generic functional descriptions.

LncRNA–mRNA Network Construction

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).

FIGURE 8

Figure 8. The interaction networks of lncRNAs and their corresponding target gene. (A) The networks for SP21 vs. LP7; (B) the networks for SP21 vs. LP21; (C) the networks for SP21 vs. LP42. The dashed and solid lines represent trans- and cis-regulation functions, respectively. Red and green represent up- and downregulation, respectively. Circles and triangles represent mRNAs and lncRNAs, respectively.

Gene Expression Validation

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).

FIGURE 9

Figure 9. Validation of RNA-Sequencing (RNA-seq) data using RT-qPCR. (A) RNA-Seq and RT-qPCR results of three selected differentially expressed lncRNAs in PT of Sunite sheep at different photoperiods. (B) RNA-Seq and RT-qPCR results of three selected differentially expressed mRNAs in PT of Sunite sheep at different photoperiods.

Discussion

Seasonal reproduction is a result of adaptation of animal reproductive activities to environmental changes, which is essential for breeding success and survival of future generations. However, the detailed molecular mechanism of animal seasonal reproduction has not been fully revealed. So far, several key genes for seasonal reproduction have been found in animals, including EYA3 gene in Japanese quail (38) and sheep (8), the vasoactive intestinal peptide (VIP) gene in Yangzhou goose (39), TSHB in sheep (8), TSHB and pituitary adenylate cyclase activating polypeptide (PACAP) in mice (19, 40), and CHGA and Tachykinin 1 (TAC1) in sheep (3, 7), which are mainly expressed in the pituitary PT. The pituitary PT is an important regulatory center for seasonal reproductive traits. Therefore, it was used as the target tissue in this study and the first objective of this study is to screen new candidate photoperiod-induced genes in PT of sheep.

Combined with the results of gene differential expression, GO, and KEGG pathway enrichment analyses, some important candidate genes were screened out. They included CHGA, FOS, SOCS3, EYA3, TSHB, SIX1, GHRH, DCT, TH, VMO1, AREG, SUV39H2, and EZH2, which were mainly involved in the following pathways: reproduction-related pathways, TNF signaling pathway, circadian entrainment, and tyrosine metabolism. Specifically, LP-induced genes included EYA3, VMO1, TSHB, SIX1, GHRH, DCT, TH, VMO1, AREG, SUV39H2, and EZH2 and SP-induced genes involved ENSOARG00000012585, CHGA, FOS, SOCS3, and TH. Wood et al. (7) reported CHGA as a SP-activator, which had a significantly high expression at SP compared with LP in PT of rams. Our results in PT of ewes also further confirm that CHGA is a SP-induced gene. For FOS, SOCS3, and ENSOARG00000012585, they may be considered as new candidate SP-induced genes in PT, whose function in PT need further analysis. The suprachiasmatic nuclei of the hypothalamus (SCN) are the master circadian clock in mammals. The role of FOS and SOCS3 in SCN has been reported in long-day breeding animals (redheaded bunting, rat, and hamster), which are enriched in the TNF signaling pathway and circadian entrainment. Their expression patterns under SP and LP are opposite to those of short-day breeding animals (41–43). For example, FOS is predominantly expressed in the SCN of redheaded buntings under LP conditions (42). Moreover, our results showed that FOS was highly expressed in pituitary PT of sheep under SP. Photoperiod also modulated SOCS3 gene expression and maintained its expression in a high level in the SCN of hamsters during LP compared with SP (44). Then, SOCS3 conveys seasonal changes into leptin sensitivity in the Siberian hamster (44).

So far, several LP-induced genes have been revealed. EYA3 and TSHB, as LP-induced genes, displayed a marked increase from SP to LP in PT of sheep (3, 7, 45). Similarly, in this study, EYA3 and TSHB were LP-induced and went up at LP compared with SP in PT of ewes, further certifying the findings of previous studies (8, 15, 46). In both sheep and mouse, EYA3 and its partner SIX1 synergistically act as upstream inducers of the TSHB transcription to induce the expression of TSHB (8, 19). Furthermore, TSHB acts locally on TSHR-expressing cells in the adjacent basal hypothalamus, leading to altered expression of DIO2 (5, 6, 47). Then, DIO2 can regulate the secretion of TH in the basal hypothalamus and further modulate the transition of reproductive status through the spatial structure of GnRH neurons (7). When the GnRH neurons were wrapped by ependymal cells, little GnRH were released into pituitary, which caused the decrease of LH and FSH secretion, and eventually sheep entered the state of anestrus (7).

Tyrosine metabolism can mediate the biological effects of many hormones and cytokines on reproduction, immunity, and cell growth (48, 49). In this study, the DE DCT was enriched in the tyrosine metabolism pathway. DCT was found as a novel long-day marker, whose expression levels in tanycytes lining the infra-lateral walls and floor of the ovine third ventricle showed marked increase with the extension of the photoperiod (6, 45). Similarly, in this study, DCT was significantly upregulated during LP compared with SP in PT of sheep. These results suggested that DCT expression was crucial for initiation of anestrus in sheep. Previous research has shown that TH can modulate DCT expression (46). DCT is a part of, or is driven by, the circannual clock in sheep. Therefore, these characteristics—LP induction and circannual changes in expression—place DCT in a list of key genes involved in seasonal timing, which had included TSHB, EYA3, DIO2, and DIO3 (45, 46).

In this study, the expression level of VMO1, AREG, SUV39H2, and EZH2 was also significantly higher during LP than those at SP in PT of sheep. As secreted factors, VMO1 and AREG were found to be acutely induced by LP in MBH of Ile-de-France ewes (6). In polyovular species, the LH-driven signaling can promote oocyte maturation and cumulus expansion by AREG. In addition, murine data indicated that LH binding to LHCGR in mural granulosa cells will up-regulate AREG (50). However, the function of AREG gene in PT is still unknown, and it may be considered as a candidate LP-induced gene to analyze. For two histone methyltransferases genes, SUV39H2 and EZH2, the results about seasonal changes of their expression in the PT are consistent with the hypothesis that epigenetic changes in PT cells are involved in circannual timing (6). Besides, the acute LP responsiveness of most of these markers (TSHB, VMO1, EZH2, and EYA3) is also in accordance with the findings from Dardente et al. (45) and Lomet et al. (6) in Ile-de-France ewes.

In recent years, studies have indicated that lncRNAs play important roles in many aspects of sheep reproductive regulation, such as fecundity (12, 13), gonadal development (14), and sex hormone response (15). As epigenetic regulators, lncRNAs can regulate the expression of reproduction-related genes. Thus, the other objective of this study is to predict lncRNAs that target key genes for seasonal reproduction in PT of sheep. Firstly, our result showed that the sequence length and exon number of mRNAs and lncRNAs in sheep pituitary have different patterns with those in hypothalamus of sheep (3,448 nt and 2.5 exons) (13). This implies that lncRNAs have tissue-specific characteristics. Then, through lncRNA–mRNA network construction, several notable target combinations between lncRNAs and key genes for seasonal reproduction were predicted. They specifically included the trans relationship between MSTRG.209166 and TSHB and cis relationships between MSTRG.235014 and DDC, MSTRG.272793 and KIT, and ENSOARG00000026131 and TH. Importantly, the trans relationship between MSTRG.209166 and TSHB was shared among all of the comparison groups (SP vs. LP), which suggested that MSTRG.209166 might play an important role in seasonal reproduction by regulating the expression of TSHB. TSHB is an important hub in the pathway of seasonal reproduction, so this study provides a new molecular object (MSTRG.209166) for the follow-up epigenetic regulation study of seasonal reproduction in sheep. In addition, the several candidate lncRNAs mentioned above (e.g., MSTRG.235014, MSTRG.272793, and ENSOARG00000026131) are also worthy of in-depth analysis.

Conclusion

In summary, our study provided a genome-wide view of lncRNA and mRNA expression profiling in pituitary PT of sheep during LPs and SPs. Several new candidate photoperiod-induced genes and lncRNAs targeting key genes of seasonal reproduction were predicted in PT of sheep. These results will provide new clues for understanding the molecular regulation of seasonal reproduction in sheep.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI BioProject; PRJNA680667.

Ethics Statement

The animal study was reviewed and approved 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) (No. IAS2018-3,10 April 2018). Written informed consent was obtained from the owners for the participation of their animals in this study.

Author Contributions

QX and MC: design of experiment and analysis for PT of sheep. QL and XH: construction of the OVX+E2 sheep model and light control experiment. RD: Experimental guidance and writing. XZ, JZ, and XG: light control experiment and collection of samples.

Funding

This research was funded by the following bodies: the National Natural Science Foundation of China (31861143012, 31472078, and 31772580), the Earmarked Fund for China Agriculture Research System (CARS-38), the Agricultural Science and Technology Innovation Program of China (ASTIP-IAS13), the China High-level Talents Special Support Plan Scientific and Technological Innovation Leading Talents Program (W02020274), the Tianjin Agricultural Science and Technology Achievements Transformation and Popularization Program (201704020), and the Youth Innovative Research and Experimental Project of Tianjin Academy of Agricultural Sciences (2020013). The APC was funded by the National Natural Science Foundation of China (31861143012). The funding bodies had no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fvets.2021.644474/full#supplementary-material

References

1. Ebling FJP, Foster DL. Photoperiod requirements for puberty differ from those for the onset of the adult breeding season in female sheep. J Reprod Fertil. (1988) 84:283–93. doi: 10.1530/jrf.0.0840283

PubMed Abstract | CrossRef Full Text | Google Scholar

2. La Y, He X, Zhang L, Di R, Chu M. Comprehensive analysis of differentially expressed profiles of mRNA, lncRNA, and circRNA in the uterus of seasonal reproduction sheep. Genes. (2020) 11:301. doi: 10.3390/genes11030301

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Dupré SM, Miedzinska K, Duval CV, Yu L, Goodman RL, Lincoln GA, et al. Identification of Eya3 and TAC1 as long-day signals in the sheep pituitary. Curr Biol. (2010) 20:829–35. doi: 10.1016/j.cub.2010.02.066

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Li X Y, He X Y, Liu Q Y, Wang X Y, Guo X F, Xia Q, et al. Expression pattern analysis of TAC1 and PRLR genes in different reproductive states of sheep. Acta Vet Zootech Sin. (2018) 49:253–62.

Google Scholar

5. Hanon EA, Lincoln GA, Fustin J-M, Dardente H, Masson-Pévet M, Morgan PJ, et al. Ancestral TSH mechanism signals summer in a photoperiodic mammal. Curr Biol. (2008) 18:1147–52. doi: 10.1016/j.cub.2008.06.076

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Lomet D, Juliette C, Chesneau D, Dubois E, Hazlerigg D, Dardente H. The impact of thyroid hormone in seasonal breeding has a restricted transcriptional signature. Cell Mol Life Sci. (2018) 75:905–19. doi: 10.1007/s00018-017-2667-x

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Wood SH, Christian HC, Miedzinska K, Saer BRC, Johnson M, Paton B, et al. Binary switching of calendar cells in the pituitary defines the phase of the circannual cycle in mammals. Curr Biol. (2015) 25:2651–62. doi: 10.1016/j.cub.2015.09.014

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Dardente H, Wyse CA, Birnie MJ, Dupré SM, Loudon ASI, Lincoln GA, et al. A molecular switch for photoperiod responsiveness in mammals. Curr biol. (2010) 20:2193–8. doi: 10.1016/j.cub.2010.10.048

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Miao X, Luo Q, Zhao H, Qin X. Ovarian transcriptomic study reveals the differential regulation of miRNAs and lncRNAs related to fecundity in different sheep. Sci Rep. (2016) 6:35299. doi: 10.1038/srep35299

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Miao X, Luo Q, Zhao H, Qin X. Co-expression analysis and identification of fecundity-related long non-coding RNAs in sheep ovaries. Sci Rep. (2016) 6:39398. doi: 10.1038/srep39398

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Feng X, Li F, Wang F, Zhang G, Pang J, Ren C, et al. Genome-wide differential expression profiling of mRNAs and lncRNAs associated with prolificacy in Hu sheep. Biosci Rep. (2018) 38:BSR20171350. doi: 10.1042/BSR20171350

PubMed Abstract | CrossRef Full Text | Google Scholar

12. He X, Tao L, Zhong Y, Di R, Xia Q, Wang X, et al. Photoperiod induced the pituitary differential regulation of lncRNAs and mRNAs related to reproduction in sheep. PeerJ. (2021) 9:e10953. doi: 10.7717/peerj.10953

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Zhang ZB, Tang JS, Di R, Liu QY, Wang XY. Comparative transcriptomics reveal key sheep (Ovis aries) hypothalamus LncRNAs that affect reproduction. Animals. (2019) 9:152. doi: 10.3390/ani9040152

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Mulvey BB, Olcese U, Cabrera JR, Horabin JI An interactive network of long non-coding RNAs facilitates the Drosophila sex determination decision. Biochim Biophys Acta. (2014) 1839:773–84. doi: 10.1016/j.bbagrm.2014.06.007

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Li W, Notani D, Ma Q, Tanasa B, Nunez E, Chen AY, et al. Functional roles of enhancer RNAs for oestrogen-dependent transcriptional activation. Nature. (2013) 498:516–20. doi: 10.1038/nature12210

CrossRef Full Text | Google Scholar

16. Smith JT, Clay CM, Caraty A, Clarke IJ. KiSS-1 messenger ribonucleic acid expression in the hypothalamus of the ewe is regulated by sex steroids and season. Endocrinology. (2007) 148:1150–7. doi: 10.1210/en.2006-1435

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Karsch FJ, Bittman EL, Foster DL, Goodman TL, Legan SJ, Robinson JE. Neuroendocrine basis of seasonal reproduction. Recent Prog Horm Res. (1984) 40:185–232. doi: 10.1016/B978-0-12-571140-1.50010-4

CrossRef Full Text | Google Scholar

18. Legan SJ, Karsch FJ, Foster DL. The endocrin control of seasonal reproductive function in the ewe: a marked change in response to the negative feedback action of estradiol on luteinizing hormone secretion. Endocrinology. (1977) 101:818–24. doi: 10.1210/endo-101-3-818

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Masumoto KH, Ukai-Tadenuma M, Kasukawa T, Nagano M, Uno KD, Tsujino K, et al. Acute induction of Eya3 by late-night light stimulation triggers TSHβ expression in photoperiodism. Curr Biol. (2010) 20:2199–06. doi: 10.1016/j.cub.2010.11.038

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Wood S, Loudon A. Clocks for all seasons: unwinding the roles and mechanisms of circadian and interval timers in the hypothalamus and pituitary. J Endocrinol. (2014) 222:R39–59. doi: 10.1530/JOE-14-0141

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Tsujino K, Narumi R, Masumoto KH, Susaki EA, Shinohara Y, Abe T, et al. Establishment of TSH β real-time monitoring system in mammalian photoperiodism. Genes Cells. (2013) 18:575–88. doi: 10.1111/gtc.12063

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc. (2016) 11:1650–67. doi: 10.1038/nprot.2016.095

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. (2015) 33:290–5. doi: 10.1038/nbt.3122

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Sun L, Luo HT, Bu DC, Zhao GG, Yu K. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. (2013) 41:e166. doi: 10.1093/nar/gkt646

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Kong L, Zhang Y, Ye Z Q, Liu X Q, Zhao S Q, Wei L, et al. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. (2007) 35:W345–9. doi: 10.1093/nar/gkm391

PubMed Abstract | CrossRef Full Text | Google Scholar

26. El-Gebali S, Mistry J, Bateman A, Eddy S R, Luciani A, Potter S C, et al. The Pfam protein families database in 2019. Nucleic Acids Res. (2019) 47:D427–32. doi: 10.1093/nar/gky995

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Wang L, Jung PH, Surendra D, Wang S, Jean-Pierre K, Li W. CPAT: Coding-Potential Assessment Tool using an alignment-free logistic regression model. Nucleic Acids Res. (2013) 41:e74. doi: 10.1093/nar/gkt006

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Trapnell C, Williams BA, Pertea G, Mortazavi A, Pachter L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. (2010) 28:511–5. doi: 10.1038/nbt.1621

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Young M D, Wakefield M J, Smyth G K, Oshlack A. Gene ontology analysis for rna-seq: accounting for selection bias. Genome Biol. (2010) 11:R14. doi: 10.1186/gb-2010-11-2-r14

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Fatica A, Bozzoni I. Long non-coding RNAs: new players in cell differentiation and development. Nat Rev Genet. (2014) 15:7–21. doi: 10.1038/nrg3606

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Shannon P, Markiel A, Ozier O, Baliga N S, Wang J T, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. (2003) 13:2498–504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-delta delta C(T)) method. Methods. (2001) 25:402–8. doi: 10.1006/meth.2001.1262

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Livak KJ. Analyzing real-time PCR data by the comparative c(t) method. Nat Protoc. (2008) 3:1101–8. doi: 10.1038/nprot.2008.73

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, et al. ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics. (2009) 25:1091–3. doi: 10.1093/bioinformatics/btp101

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Bindea G, Galon J, Mlecnik B. CluePedia Cytoscape plugin: pathway insights using integrated experimental and in silico data. Bioinformatics. (2013) 29:661–3. doi: 10.1093/bioinformatics/btt019

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Lincoln GA, Andersson H, Hazlerigg D. Clock genes and the long-term regulation of prolactin secretion: evidence for a photoperiod/circannual timer in the pars tuberalis. J Neuroendocrinol. (2003) 15:390–7. doi: 10.1046/j.1365-2826.2003.00990.x

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Dardente H, Migaud M. Thyroid hormone and hypothalamic stem cells in seasonal functions. Vitam Horm. (2021) 116:91–131. doi: 10.1016/bs.vh.2021.02.005

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Nakao N, Ono H, Yamamura T, Anraku T, Takagi T, Higashi K, et al. Thyrotrophin in the pars tuberalis triggers photoperiodic response. Nature. (2008) 452:317–22. doi: 10.1038/nature06738

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Zhu HX, Chen Z, Shao XB, Yu JN, Wei CK, Dai ZC, et al. Reproductiveaxis gene regulation during photostimulation and photorefractoriness in Yangzhou goose ganders. Front Zool. (2017) 14:11. doi: 10.1186/s12983-017-0200-6

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Hannibal J, Georg B, Fahrenkrug J. PAC1- and VPAC2 receptors in light regulated behavior and physiology: studies in single and double mutant mice. PLos ONE. (2017) 12:e0188166. doi: 10.1371/journal.pone.0188166

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Veen DRvd, Pol-Meijer MMTvd, Jansen K, Smeets M, Zee EAVD, Gerkema MP. Circadian rhythms of c-FOS expression in the suprachiasmatic nuclei of the common vole (Microtus arvalis). Chronobiol In. (2008) 25:481–99. doi: 10.1080/07420520802254403

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Majumdar G, Yadav G, Rani S, Kumar V. A photoperiodic molecular response in migratory redheaded bunting exposed to a single long day. Gen Comp Endocrinol. (2014) 204:104–13. doi: 10.1016/j.ygcen.2014.04.013

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Jennings KJ, Chasles M, Cho H, Mikkelsen J, Bentley G, Keller M, et al. The preoptic area and the RFamide-related peptide neuronal system gate seasonal changes in chemosensory processing. Integr Comp Biol. (2017) 57:1055–65. doi: 10.1093/icb/icx099

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Alexander T, Claire E, Moar KM, Logie TJ, Adam CL, Mercer JG, et al. Photoperiodic regulation of leptin sensitivity in the Siberian hamster, Phodopus sungorus, is reflected in arcuate nucleus SOCS-3 (suppressor of cytokine signaling) gene expression. Endocrinology. (2004) 145:1185–93. doi: 10.1210/en.2003-1382

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Dardente H, Lomet D, Chesneau D, Pellicer-Rubio MT, Hazlerigg D. Discontinuity in the molecular neuroendocrine response to increasing daylengths in Ile-de-France ewes: is transient Dio2 induction a key feature of circannual timing? J Neuroendocrinol. (2019) 31:e12775. doi: 10.1111/jne.12775

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Dardente H, Lomet D. Photoperiod and thyroid hormone regulate expression of L-Dopachrome tautomerase (Dct), a melanocyte stem-cell marker, in tanycytes of the ovine hypothalamus. J Neuroendocrinol. (2018) 30:e12640. doi: 10.1111/jne.12640

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Trivedi AK, Sur S, Sharma A, Taufique ST, Kumar V. Temperature alters the hypothalamic transcription of photoperiod responsive genes in induction of seasonal response in migratory redheaded buntings. Mol Cell Endocrino. (2019) 493:110454. doi: 10.1016/j.mce.2019.110454

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Rui H, Xu J, Mehta S, Fang H, Williams J, Dong F, et al. Activation of the Jak2-Stat5 signaling pathway in Nb2 lymphoma cells by an anti-apoptotic agent, aurintricarboxylic acid. J Biol Chem. (1998) 273:28–32. doi: 10.1074/jbc.273.1.28

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Zhao Y, Kan FWK. Human OVGP1 enhances tyrosine phosphorylation of proteins in the fibrous sheath involving AKAP3 and increases sperm-zona binding. J Assist Reprod Genet. (2019) 36:1363–77. doi: 10.1007/s10815-019-01502-0

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Park JY. EGF-like growth factors as mediators of LH action in the ovulatory follicle. Science. (2004) 303:682–4. doi: 10.1126/science.1092463

PubMed Abstract | CrossRef Full Text | Google Scholar