Insight Into Pituitary lncRNA and mRNA at Two Estrous Stages in Small Tail Han Sheep With Different FecB Genotypes

The pituitary is a remarkably dynamic organ with roles in hormone (FSH and LH) synthesis and secretion. In animals with the FecB (fecundity Booroola) mutation, the pituitary experiences hormone fluctuations during the follicular–luteal transition, which is implicated in the expression and regulation of many genes and regulators. Long non-coding RNAs (lncRNAs) are a novel type of regulatory factors for the reproductive process. Nevertheless, the expression patterns of lncRNAs and their roles in FecB-mediated follicular development and ovulation remain obscure. Thus, we profiled the pituitary transcriptome during the follicular (F, 45 h after evacuation vaginal sponges) and luteal (L, 216 h after evacuation vaginal sponges) phases in FecB-mutant homozygous (BB) and wild-type (WW) Small Tail Han sheep. We identified 78 differentially expressed genes (DEGs) and 41 differentially expressed lncRNAs (DELs) between BB_F and BB_L, 32 DEGs and 26 DELs between BB_F and WW_F, 16 DEGs and 29 DELs between BB_L and WW_L, and 50 DEGs and 18 DELs between WW_F and WW_L. The results of real-time quantitative PCR (RT-qPCR) correlated well with the transcriptome data. In both the follicular and luteal phases, DEGs (GRID2, glutamate ionotropic receptor delta type subunit 2; ST14, ST14 transmembrane serine protease matriptase) were enriched in hormone synthesis, secretion, and action. MSTRG.47470 and MSTRG.101530 were the trans-regulated elements of ID1 (inhibitor of DNA binding 3, HLH protein) and the DEG ID3 (inhibitor of DNA binding 3, HLH protein), and EEF2 (eukaryotic translation elongation factor 2), respectively; these factors might be involved in melatonin and peptide hormone secretion. In the FecB-mediated follicular phase, MSTRG.125392 targeted seizure-related 6 homolog like (SEZ6L), and MSTRG.125394 and MSTRG.83276 targeted the DEG KCNQ3 (potassium voltage-gated channel subfamily Q member 3) in cis, while MSTRG.55861 targeted FKBP4 (FKBP prolyl isomerase 4) in trans. In the FecB-mediated luteal phase, LOC105613905, MSTRG.81536, and MSTRG.150434 modulated TGFB1, SMAD3, OXT, respectively, in trans. We postulated that the FecB mutation in pituitary tissue elevated the expression of certain genes associated with pituitary development and hormone secretion. Furthermore, this study provides new insights into how the pituitary regulates follicular development and ovulation, illustrated by the effect of the FecB mutation.


INTRODUCTION
The FecB (fecundity Booroola) mutation was first detected in high-fecundity Australian Booroola Merino sheep (1). Mulsant et al. (2) Souza et al. (3), and Wilson et al. (4) revealed that the effect of FecB mutation was due to the substitution of a base (A to G at 746 bp, A746G) in the BMPR1B (bone morphogenetic protein receptor 1B) coding region on ovine chromosome 6, resulting in the corresponding substitution of one amino acid (glutamine to arginine, Q249R). Through an impairment of the activity of the BMPR1B receptor kinase, FecB mutation weakened the inhibitory effect of the BMPR1B gene on steroid synthesis in granulosa cells, leading to the maturation and ovulation of ovulatory follicles at significantly small diameters and fewer granulosa cells in homozygous FecB sheep than those in non-FecB sheep (2,3,(5)(6)(7). In our previous study, litter size was significantly correlated with the FecB genotype in highprolificacy Chinese Small Tail Han sheep (8). The gonadotropic hormones follicle-stimulating hormone (FSH) and luteinizing hormone (LH) are the dominant hormones in the coordinate regulation of follicle development and antral follicle ovulation (9). A previous study predicted that FecB might impact the pituitary gland through the regulation of hormone release during one estrus cycle in ewes (10). At a particular physiological time during the follicular phase, the plasma concentration of FSH in FecB homozygous ewes was higher than that in non-FecB ewes (11)(12)(13). The follicular-luteal transition is a complex and elaborate process that involves pituitary hormones and several cell-specific gene products that are under strict developmental regulation. The identification of the key regulators of these developmental processes could be beneficial for understanding the mechanisms of how FecB regulates the follicular-luteal transition at the molecular level and for screening candidate protein-coding genes with impacts on reproductive ability in ewes.
Previous transcriptomics studies have focused solely on protein-coding genes (14). Nevertheless, transcriptomic analyses and chromatin signatures showed a plethora of noncoding transcripts produced by eukaryotic cells (15,16). Of these, long non-coding RNAs (lncRNAs) are a heterogeneous group of non-protein-coding transcripts of greater than 200 nucleotides in length with unknown coding potential. They include circular RNAs, long intergenic RNAs, pseudogenes, and antisense RNAs (17). Accumulating evidence has implicated lncRNAs in a wide array of cellular processes, such as transcriptional regulation, cellular pluripotency and reprogramming, and differentiation (reviewed in 18,19). In terms of transcriptional regulation, lncRNAs act in cis or trans through the recruitment of RNAs and proteins. Cis-acting lncRNAs activate, repress, or modulate gene expression in a manner dependent on the location of their own transcription sites (20). Unfortunately, it has been difficult to accurately identify lncRNA sequences from transcribed sequences. Transacting lncRNAs affect nuclear structure and organization and interact with and modulate the activity of RNAs and proteins (21). The cytoplasmic localization of lncRNAs and the biological activities they participate in suggest that some lncRNAs can regulate the secretion of reproductive hormones.
Pituitary gonadotropins form the center of the endocrine axis involved in the reproductive process. At puberty, hypothalamic gonadotropin-releasing hormone (GnRH) stimulates the pituitary tissue to synthesize FSH and LH. Both the FSH and LH receptors are G protein-coupled receptors and have cascade amplification effects on the estrous or menstrual cycle. Gonadotropins are produced at demanded levels and specific times in response to multiple hormonal signals; thus, regulatory elements (lncRNAs) are expected to play key roles in both the stage specificity and hormonal responsiveness of their expression. Several lncRNAs have been proposed to affect the signal transduction of FSH and LH. In rats, the pituitary intergenic lncRNA m433s1 promoted the expression level of FSHB and the secretion of FSH by functioning as a competing endogenous RNA of miR-433 and reducing its inhibition of FSHB messenger RNA (mRNA) (22). In humans, silencing the pituitary lncRNA UCA1 significantly suppressed the secretion of prolactin secretion and the growth of pituitary cancer cells (23). In ram pituitary cells, short-interfering RNA knockdown of lncRNA-TCONS_00066406 reduced the mRNA levels of FSHB and LHB (24). FecB mutation has been reported to cause putative transcriptional changes in the pituitary. Pituitary MSTRG.259847.2 was upregulated in FecB homozygous ewes compared with FecB heterozygous ewes. A subsequently constructed lncRNA-mRNA interaction network highlighted Smad2, and experiments at the cell level showed that MSTRG.259847.2 can inhibit Smad2 and LHB (25). However, FecB mutation reduced the expression levels of some genes in the pituitary, including BMP2, TGFB2, TGFB3, ID1, and ID3 (26). Studies identifying lncRNAs associated with the follicular-luteal transition have been performed mainly on ovarian tissues, and investigation of lncRNAs at the pituitary level in ewes with different FecB status and during different stages of estrus has been rarely reported.
Small Tail Han sheep are an indigenous Chinese highfecundity sheep breed with year-round estrus and a high ovulation rate (27). The estrus cycle consists of the follicular and luteal phases, accompanied by changes in the concentrations of hormones such as FSH and LH. Our previous study observed that FecB-mutant ewes regulate the expression levels of the BMPR1B, BMP15, and GDF9 genes, which are relevant to hormone and signal transduction (28)-that is, FecB mutation might affect the follicular-luteal transition to further influence sheep fertility. In this study, we collected pituitary tissues to analyze how the pituitary regulates follicular development, maturation, and ovulation in the context of FecB mutation. Our data provided a useful resource for further studying the functional roles of these genes and lncRNAs in the follicular development and ovulation of BB genotype ewes.

Ethical Consideration
All experimental ewes were housed and fed by the Science Research Department of the Institute of Animal Sciences, Chinese Academy of Agriculture Sciences (IAS-CAAS), and ethical approval was received from the Animal Ethics Committee of the IAS-CAAS (no. IAS 2019-49).

Experimental Design and Pituitary Collection
Based on TaqMan genotyping (29), six FecB-mutant homozygous (BB, litter size = 2.2 ± 0.6) and six FecB wild-type (WW, litter size = 1.0 ± 0.0) Small Tail Han ewes were selected from a nucleus herd (Shandong, China). The selected ewes were approximately 3 years old and were similar in weight ( Table 1). All non-pregnant animals had free access to food and water under natural temperature and lighting conditions. We divided genotype BB and WW ewes into a follicular phase group and a luteal phase group. Estrus synchronization was performed in September, within the usual breeding period (from August to November) described for this breed in Tianjin (117.2 E, 39.13 N). Firstly, ewes were treated with vaginal sponges (300 mg progesterone; InterAg Co., Ltd., Hamilton, New Zealand), followed by injection of vitamin AD to protect the endometrium. The vaginal sponges were evacuated after 12 days, and the removal time was set as 0 h. In accordance with previous reproductive trait analyses, six ewes (three BB and three WW ewes) were euthanized at 45 h (follicular phase) and six ewes (three BB and three WW ewes) at 216 h (luteal phase) (30,31). Pituitary samples were dissected immediately after euthanasia, frozen in liquid nitrogen, and stored at −80°C for further analysis.

RNA Isolation, Library Construction, and Sequencing
Total RNA was extracted from the pituitary using the TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA). The purity, concentration, and integrity of RNA were checked with a Nano Photometer ® spectrophotometer (IMPLEN, Westlake Village, CA, USA), Qubit ® RNA Assay kits (Thermo Fisher Scientific, Waltham, MA, USA), and RNA Nano 6000 Assay, respectively. The RNA integrity number (RIN) value of all the samples was greater than 7.
Three micrograms of total RNA was used as the starting amount for the construction of an lncRNA library. Ribo-Zero ™ Gold Kits (Epicenter, Madison, WI, USA) were applied to remove ribosomal RNA (rRNA) from the total RNA. Following the instructions for the NEB Next ® Ultra Directional RNA Library Prep Kit for Illumina ® (NEB, Ipswich, MA, USA), 12 RNA sequencing libraries were constructed for paired-end sequencing. Pooled libraries were sequenced on the Illumina HiSeq X platform (Illumina, San Diego, USA) using a chain-specific library construction strategy and the number and structure of the mRNA, known lncRNA, and novel lncRNA transcripts. All the experimental sequencing data were generated by Annoroad Gene Technology Co., Ltd. (Beijing, China).

Data Quality Control and Assembly
To ensure the quality of further analytical data, raw reads were filtered using in-house Perl scripts (Annoroad Gene Technology Co., Ltd, Beijing, China) to obtain high-quality reads. The filtering step involved deleting reads with low quality (Phred quality score <5%), adapter contamination, matches to rRNA sequence, or a rate of ambiguous bases higher than 5% (32). The Phred quality score refers to the rate of sequencing errors for a given base; for example, Q30 indicates that the base sequencing error rate was 0.1%. We aligned paired-end clean reads to the oar4.0 sheep reference genome (https://www.ncbi.nlm.nih.gov/assembly/GCF_ 000298735.2) using HISAT2 (v.2.0.5) (33) with the parameters "-rna-strandness RF" and "-dta -t -p 4". Only the uniquely mapped reads were assembled, and the expression levels were predicted using String Tie (v.1.3.2d) (33) with the parameter "-G ref.gtf -rf -1". Thereafter, we calculated the number and ratio of the uniquely mapped reads within the three gene functional elements: exons, introns, and intergenic elements. Homogeneity analysis was subsequently performed to guarantee that the sequencing results did not impact further transcriptome analysis. Thus, mRNA and known lncRNA transcripts were identified.

Identification of Potential lncRNA Candidates
To decrease the false-positive rate, the screening of candidate lncRNAs is divided into two steps: one is based on the position of the coding reads [intronic lncRNA, long intergenic non-coding RNA (lincRNA), and antisense lncRNA] and the other is to count the encoding potential (34). The following were applied: 1) Transcripts with lengths longer than 200 bp, exon numbers greater than 2, and read coverage greater than 5 were retained. 2) The GffCompare program (v.0.10.1) was used to annotate the assembled transcripts and discard known mRNAs and other non-coding RNAs (rRNAs, tRNAs, snoRNAs, etc.) (35).
3) Transcripts with class_code "i" (potential intronic lncRNA), "x" (potential antisense lncRNA), "u" (potential intergenic lncRNA), "j" (potential novel isoform with more than one splice junction of a reference transcript), or "o" (genetic exonic overlapping with a reference transcript) were retained (36). 4) The Coding-Non-Coding Index (CNCI), Coding Potential Calculator (CPC), the Pfam database, and the Coding Potential Assessment Tool (CPAT) were utilized to predict the coding potential of preliminary candidate lncRNAs. The results predicted by the four aforementioned tools were intersected to obtain potential lncRNA candidates. The CNCI was determined by profiling adjoining nucleotide triplets to distinguish coding/ non-coding transcripts independent of known annotations (37), with the parameter "-score 0 -length 199 -exon_num 2". CPC can access the protein-coding potential of each transcript according to six biologically meaningful sequence features (38

Differential Expression Analysis of lncRNAs and mRNAs and Clustering
To verify the rationality of the samples and the reliability of the experiment, the correlation of the expression levels between samples was examined using Pearson's correlation methods. Gene expression levels are generally measured on the basis of the amount of mRNA transcribed from a gene. Based on the reference GTF file and the HISAT BAM files, the HTSeq Python package (v.0.6.1) was used to calculate read counts, with the parameters "-I gene_id -f bam -s" and "reverse -a 10 -q". Because three biological replicates were performed for each experimental group, the DEseq2 package (v.1.28.1) was used for the analysis of differential expression to calculate the log 2 (fold change) and pvalue based on the negative binomial distribution model (41). Meanwhile, log2(fold change) > 1 and p adj < 0.05 were regarded as the essential criteria to identity the mRNAs and lncRNAs that were significantly differentially expressed in BB_F vs.

Target Gene Prediction for Differentially Expressed lncRNAs and Network Construction
The relationships between lncRNAs and potential proteincoding genes were predicted on the basis of their distances and expression correlations, including the cis-and trans-acting relationships. Protein-coding genes adjacent to the lncRNAs (100 kb upstream and downstream) were identified as potential cis-target genes (43). Potential trans-target genes were predicted using Spearman's correlation coefficients between the expressions of mRNAs and lncRNAs (r 2 ≥ 0.95) (43). According to the cis-or trans-acting relationship between the differentially expressed lncRNAs (DELs) and target genes, we visualized the lncRNA-mRNA networks with Cytoscape (v.3.8.0), which was run with layout = "attribute circle layout".

Functional Annotation of Candidate Genes
To reveal the potential roles of DELs, we performed Gene Ontology (GO) category and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses on their target genes using the clusterProfiler package (v.3.16.0) with the false discovery rate (FDR) pAdjustMethod. GO analysis (http:// geneontology.org) provides structured and quantifiable knowledge concerning the functions of genes and gene products, and the molecular function ontology represents the activity of the gene products, especially focusing on transcription regulator activities (44). KEGG systematically analyzes gene functions, linking genomic information with functional information (45). Furthermore, KEGG analysis can predict the protein interaction networks involved in various cellular processes. A GO term or a KEGG pathway with corrected p < 0.05 was considered significantly enriched. We mainly focused on the pathways related to reproduction and pituitary function.

RT-qPCR Validation
Six mRNAs and six lncRNAs were randomly selected for validation via real-time quantitative PCR (RT-qPCR). b-actin was used as a control. The primers were designed with Primer Premier 5 and synthesized by Qine Ke Biotech (Beijing, China) (Supplementary Table S6). For RT-qPCR analysis, 1 mg of total RNA was reverse transcribed to complementary DNA (cDNA) following the instructions in the PrimeScript ™ RT reagent Kit with gDNA Eraser (Takara, Beijing, China). The genes and lncRNAs were diluted fourfold before use for RT-qPCR according to the instructions in TB Green ® Premix Ex Taq ™ II (Takara, Beijing, China) and analyzed with QuantStudio ® 3 (ABI, Foster City, CA, USA). All genes and lncRNAs were amplified in triplicate. The data of each gene and lncRNA were analyzed with the 2 −DDCt method (46). The correlation between the RT-qPCR data and the FKPM value from the RNA sequencing (RNA-seq) data was calculated using Pearson's correlation method in R (v.4.0.2).

Overview of the Sequencing Data in Small Tail Han Sheep Pituitary Tissue
The raw RNA-seq data from the Small Tail Han sheep were analyzed for quality control before aligning reads to the reference genome. More than 97.50 Gb of clean read sequences were generated from each library, with a corresponding rate of clean Q30 bases higher than 91.24%, indicating that the filtered sequence was appropriate for mapping. All the high-quality clean reads were mapped to the oar4.0 sheep reference genome using HISAT2. The percentage of the total mapped and uniquely mapped reads in these samples reached 92.66% and 88.27%, respectively. The rate of multi-mapping reads from the 12 libraries was less than 4.62% (Supplementary Table S1).

Identification of lncRNAs and mRNAs in Small Tail Han Sheep Pituitary Tissue
According to their locations, the unique reads from 12 samples were categorized into three types: intergenic, intronic, and exon ( Figure 1A). Approximately 15.17% of the lncRNAs were distributed in intergenic regions, 42.85% were introns overlapping lncRNAs, and 41.98% were exons overlapping lncRNAs. After screening preliminary candidate lncRNAs, the intersecting results of the CNCI, CPC, Pfam, and CPAT tools were used to obtain potential lncRNA candidates, shown in Figure 1B. Furthermore, we identified the essential features of the lncRNAs and contrasted them with those of mRNAs. The genome distribution of most mRNAs and lncRNAs was randomly assigned to the 26 autosomes and the X chromosome; 480 mRNAs (2.3%) and 431 lncRNAs (3.8%) not aligned with any chromosomal location, and 10 mRNAs and 1 lncRNA were found to align to a mitochondrial location ( Figure 1C). The length distribution of the lncRNAs was mainly appropriately 2,900-3,000 bp, similar to that of the mRNA transcripts in pituitary tissue ( Figure 1D). Most lncRNAs possessed two exons, significantly fewer than the average number of mRNA transcripts ( Figure 1E). The expression abundance of lncRNAs was lower than that of mRNAs in the 12 samples ( Figure 1F).

Profiling of DE lncRNA and DEGs in Small Tail Han Sheep Pituitary
The correlation coefficient of the gene expression levels between samples was higher than 0.80, indicating that the sample selection was consistent and reliable. Under the criteria of log2 (fold change) > 1 and p adj < 0.05, 78 differentially expressed genes (DEGs; 37 upregulated and 41 downregulated) and 41 DELs (26 upregulated and 15 downregulated) were screened between BB_F and BB_L (Figure 2A and Supplementary Table S2), 32 DEGs (20 upregulated and 12 downregulated) and 26 DELs (15 upregulated and 11 downregulated) were screened between BB_F and WW_F ( Figure 2B and Supplementary Table S2), 16 DEGs (9 upregulated and 7 downregulated) and 29 DELs (9 upregulated and 20 downregulated) were screened between BB_L and WW_L ( Figure 2C and Supplementary Table S2), and 50 DEGs (19 upregulated and 31 downregulated) and 18 DELs (9 upregulated and 9 downregulated) were screened between WW_F and WW_L ( Figure 2D and Supplementary Table  S2). The mRNA and lncRNA expression patterns were similar in both the follicular and luteal phases. Nevertheless, a different expression profile was observed at the follicular-luteal transition, which may indicate that changes in transcript expression accompanied FecB mutation ( Figures 2E, F) Figure 2G).

GO and Pathway Enrichment Analysis of DEGs
The described methods were employed to predict the functions of the DEGs; that is, the enriched GO functional terms and KEGG pathways were considered to be the predicted functional terms and pathways for the protein-coding gene and lncRNA (47). In the comparison of the follicular and luteal phases, GO analysis of the DEGs showed that prolactin secretion, regulation of mitotic cell cycle, and amino acid metabolism were among the top 50 enriched terms (Supplementary Table S3). In addition, 125 pathways were annotated to hormone synthesis, secretion, and action in follicular-luteal transition, including steroid hormone, growth hormone, thyroid hormone, estrogen, and insulin ( Figure 3A and Supplementary Table S3). DEGs between the different FecB genotypes were enriched in pituitary hormone signal transduction processes (in the top 50 enriched terms), such as neurotransmitter metabolic process, choline metabolic process, and regulation of synapse process (Supplementary Table S3). KEGG pathway analysis identified 65 enriched pathways, including neuroactive ligand-receptor interaction, oxidative phosphorylation, dopaminergic synapse, melanogenesis, gap function, Wnt, mTOR, MAPK, PI3K-Akt, and Hippo signaling pathways ( Figure 3B and Supplementary  Table S3).

Screening of Potential Functional LncRNAs Involved in Small Tail Han Sheep Reproduction
To further reveal the potential biological function of pituitary lncRNAs in Small Tail Han sheep, an interaction network of lncRNAs and their corresponding target genes was constructed. In the comparison BB_F vs. BB_L, 12 known lncRNAs corresponded to 86 target genes, and 28 novel lncRNAs corresponded to 189 target genes ( Figure 4A and Supplementary  Table S5). In the comparison WW_F vs. WW_L, two known lncRNAs corresponded to eight target genes, and 15 novel lncRNAs corresponded to 75 target genes ( Figure 4B and Supplementary  Table S5). In the comparison of the follicular and luteal phases, the  Table S4). In addition, 144 pathways were related to circadian entrainment, oocyte meiosis and maturation, and hormone signaling pathways, including thyroid hormone, oxytocin, and growth hormone ( Figure 3C and Supplementary Table S4). These follicular-luteal transition DELs might play similar roles in pituitary hormone synthesis and secretion. In Table 2 In the comparison BB_F vs. WW_F, two known lncRNAs corresponded to four target genes, and 21 novel lncRNAs corresponded to 113 target genes ( Figure 4C and Supplementary  Table S5). In the comparison BB_L vs. WW_L, six known lncRNAs corresponded to 19 target genes, and 21 novel lncRNAs corresponded to 132 target genes ( Figure 4D and Supplementary  Table S5). In the comparison of the FecB genotypes, adenosine receptor binding, postsynaptic cytosol, kinase activity, and neurohypophyseal hormone activity were enriched in the top 50 terms (Supplementary Table S4). In addition, 208 pathways were identified, which included tight junction, TGF-beta, Hippo, Wnt, circadian rhythm, melanogenesis, MAPK, and PI3K-Akt signaling pathways ( Figure 3D and Supplementary Table S4). Enrichment analyses of the DELs and DEGs showed similarities, indicating that DELs might be involved in the process by which FecB mutation had an effect on pituitary function. Two DEGs and five genes targeted by DELs are shown in detail ( Table 2).

Validation of RNA Sequencing Using RT-qPCR
To validate the RNA-seq data, eight genes and eight lncRNAs were examined with RT-qPCR, and the relative gene expression was calculated on the basis of the 2 −DDCt values ( Figure 5). The results indicated that the data on the mRNA and lncRNA expression levels were reliable. Meanwhile, the results of Pearson's correlation analysis of all genes showed a strong positive correlation between the RT-qPCR and RNA-seq data (r 2 > 0.93, p < 0.05).

DISCUSSION
Pituitary gonadotropin is essential for the activation of reproduction through the synthesis and secretion of gonadotropic hormones (FSH and LH) at the required levels. Given the unique function of   the pituitary, i.e., the production of FSH and LH at specific times and appropriate levels, and its response to diverse hormone signals, it is expected that some regulatory elements play dominant roles in this cell-specific expression and hormonal responsiveness. In our study, 8,562 novel lncRNAs identified from the sheep pituitary shared many characteristics with those from other mammals, such as goats, mice, and pigs. In these mammals, most lncRNAs possessed two exons (followed by three exons), and lncRNAs were shorter and lower than mRNAs and were expressed at lower levels, reinforcing the legitimacy of the transcripts identified here (48,49). Fewer lncRNAs were identified here than in the hypothalamus and pituitary gland of sheep, which had 34,293 and 19,672 lncRNAs, respectively (25,50), likely reflecting tissuespecific expression patterns and characteristics and the different sampling times employed. The most lncRNAs in the present study were detected on ovine chromosomes 1-3, which was perhaps unsurprising as these chromosomes are the largest in the ovine genome. This implies a relationship between these chromosomes and pituitary function in the context of the FecB mutant. Similar chromosome locations of lncRNAs were found in previous studies on Hu sheep (chromosomes 1, 2, and 3) (25). In other words, the total distribution of lncRNAs was consistent between Hu sheep and Small Tail Han sheep, indicating that the functions of these lncRNAs might be conserved. Meanwhile, 20,248 mRNAs were identified, including OXTR, which is associated with the neurohypophysis, and hormone genes (FSHB, LHB, PRL, and INHA) and receptors (FSHR, ESR1, and ESR2) associated with the adenohypophysis (51,52).
To understand the regulatory mechanism of the follicularluteal transition, a previous study focused mainly on the microRNA (miRNA)-mRNA interactions in ovarian tissue and revealed that some miRNAs, such as miR-200a, miR-200b, and miR-200c, appeared to be critical regulators of ovarian follicular and luteal development (53,54). However, pituitary hormones (FSH and LH) are required for the transition from the follicle to the corpus luteum, which involves receiving feedback signals from gonadal steroids and activin/inhibin signaling pathways (55,56). In our study, DEGs from the follicular-luteal transition of wild-type ewes were mainly enriched in the G2/M transition of mitotic cell cycle, response to BMP, and Wnt signaling pathway of the GO terms and tight junction, TGF-beta signaling pathway, circadian entrainment, gap junction, and hormone-related signaling pathways of the KEGG pathways. IGFBP3, PLCE1, and GRM1 were annotated to hormone-related signaling pathways, which might belong to the adenohypophysis. Moreover, GRM1 was annotated as a gap junction gene; gap junctions mediate intercellular communication or cytosolic connections to regulate the amount of gonadotropin produced (57). With the FecB mutation, the number of DEGs increased and were annotated to amino acid metabolic process, prolactin secretion, and melanin biosynthetic process of the GO terms and the neuroactive ligand-receptor interaction, amino acid metabolism, hormone biosynthesis, and synapse-related pathways of the KEGG pathways. RLN3, GRM8, GRID2, and VIP were enriched in neuroactive ligand-receptor interaction, and GRM8, KCNQ3, and MAOB were enriched in synapse-related pathways, suggesting that these genes might be part of the neurohypophysis. Among these DEGs, GRID2 and ST14 were downregulated in the follicular phase compared with the luteal phase. A previous study revealed that GRID2 is the receptor for glutamate that controls anterior pituitary hormone secretion and reproduction (58,59); that is, GRID2 may promote LH secretion to be involved in ovulation like glutamate receptor AMPA 1 (60). ST14 reduced cell growth by downregulating the cell cycle-related proteins, and ST14 was detected in the pituitary, indicating that it played an important role in pituitary development (61). Meanwhile, we identified some RAB proteins, among which RAB1A and RAB26 had higher relative expressions. RAB proteins are master regulators of vesicle transport, which is an essential process for secretion (62 (2). It was suggested that, in ewes with the BB genotype, the BMPR1B mutation would advance granulosa cell differentiation and ovulatory follicle maturation, resulting in differential gene expression in the pituitary gland (68). In our study, the BMPR1B gene did not reach a significant difference in the different genotypes. However, estrogen functions in ovarian tissue and is involved in the process by which it downregulates the production and secretion of pituitary gonadotropin via negative feedback to the pituitary (69). Comparison of BB_F vs. WW_F showed that the FecB mutation significantly increased the expression of the WNT11 (Wnt family member 11) gene in the follicular phase. The expression levels of the three DELs (MSTRG.125392, MSTRG.125394, and MSTRG.83276) and the two target genes (SEZ6L and KCNQ3) in the pituitary of BB genotype ewes were significantly higher than those of WW genotype ewes, and the expression trends of lncRNAs and their corresponding cis-target genes were consistent. Moreover, cisregulatory lncRNAs showed enhancer-like activity and promoted the expression of neighboring genes. The DEG SEZ6L was cisregulated by MSTRG.125392 and MSTRG.125394 and has been reported to influence dendritic morphology and synapse numbers (70). Previous conservation genomic analysis identified SEZ6L as a candidate gene associated with reproduction and system development in Black Slavonian pigs through the modulation of key oncogenic pathways such as Wnt (71,72). Meanwhile, the DEG KCNQ3, a cis-regulatory element of MSTRG.83276, was important in the regulation of NPY/AgRP neuron excitability, which promoted hormone communication to maintain homeostasis. These results indicated that these three DELs might regulate the transcription of the target genes by promoting the enhancer and the promoter of the transcriptional machinery molecule, similar to the manner in which classical lncRNAs (HOTTIP and HOTAIRM1) promoted the expression of the adjacent target gene HoxA by cis regulation (73). FKBP4, an element that is trans-regulated by MSTRG.55861, is relevant to steroid hormone receptor binding and transportation and to reproductive and neurological diseases (74). Our study showed that FKBP4 was somewhat differentially expressed and might induce steroid hormone synthesis in pituitary tissue. As an important paralog of this gene, FKBP5 (FKBP prolyl isomerase 5) might produce more FSH in the hypothalamus tissue (75); that is, FKBP4 and FKBP5 are closely tied to hormone synthesis and secretion in the hypothalamic-pituitary-ovarian (HPO) axis, and FecB affects their expression trend. In the luteal phase in different FecB genotypes, the target genes TGFB1 (transforming growth factor beta 1), SMAD3 (SMAD family member 3), and OXT (oxytocin/neurophysin I prepropeptide) were trans-elements of the lncRNAs LOC105613905, MSTRG.81536, and MSTRG. 150434, respectively. It was suggested that, with FecB mutation, lncRNAs in the follicular phase might be involved in pituitary development and steroid hormone communication, while lncRNAs in the luteal phase probably participate in the TGFb/SMAD signaling pathway and in oxytocin secretion.
Overall, our results suggested that some genes were differentially expressed in the pituitary tissue between the follicular and luteal phases via the regulation of reproductive hormone synthesis and secretion. During the follicular phase, the expression of some genes in the pituitary tissue was elevated in the FecB mutant in comparison with the wild-type genotype. Whether the change in the gene expression levels affected the corresponding reproductive pathways and processes remains to be further elucidated.

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: https://www.ncbi. nlm.nih.gov/, accession ID: PRJNA782215.