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

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.

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

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.

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). Both CNCI and CPC are dependent on a score <0, thus defining lncRNA transcripts as non-coding RNAs. In Pfam, the profile hidden Markov models (HMMs) were searched for protein domains to obtain transcripts with known protein domains, based on the UniProt Knowledgebase (UniProtKB) (39). LncRNAs could be regarded as non-coding RNAs when the E-value ≥0.001. CPAT (v.1.2.1) works on a logistic regression model with four sequence features: Fickett TESTCODE statistic, open reading frame (ORF) size, ORF coverage, and hexamer usage bias (40).

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₂(fold change) and p-value 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. BB_L (different estrus stages under genotype BB ewes, named gene set 1), BB_F vs. WW_F (different FecB genotypes ewes in the follicular phase, named gene set 2), BB_L vs. WW_L (different FecB genotypes ewes in the luteal phase, named gene set 3), and WW_F vs. WW_L (different estrus stages under genotype WW ewes, named gene set 4) individually. The intersection of BB_F vs. BB_L and WW_F vs. WW_L was named gene set 5; genes in this set were considered highly conserved during the follicular–luteal transition and to be potentially affected at the transcriptional levels by the FecB mutation. The intersection of BB_F vs. BB_L and BB_F vs. WW_F was named gene set 6; these genes were implicated in the process of follicular development and ovulation modulation by FecB. StringTie (v1.3.2d) was adopted to count the fragments per kilobase of transcripts per million mapped reads (FPKM). Subsequently, systematic clustering analysis based on the log₂(FPKM) value of each mRNA and lncRNA was performed with pheatmap (v.1.0.2) to analyze the similarities and relationships between the different libraries (42). The analyses included Pearson’s correlation and Euclidean distance methods.

The relationships between lncRNAs and potential protein-coding 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 ² ≥ 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”.

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.

Six mRNAs and six lncRNAs were randomly selected for validation via real-time quantitative PCR (RT-qPCR). β-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 μg 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^−ΔΔCt 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).

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

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

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 ). Furthermore, we compared the quantity and distribution of the DEGs and DELs in each of the four pairwise comparisons. The DEGs and DELs from the intersection of BB_F vs. BB_L and WW_F vs. WW_L were involved in follicular–luteal transition and conserved across the different genotypes, including GRID2, LOC105610376, MSTRG.163804, and MSTRG.81824. The intersection of BB_F vs. BB_L and BB_F vs. WW_F contained vital genes and lncRNAs, including LRRC7, INSM2, SEZ6L, BMX, KCNQ3, GRM8, LDB3, CDH12, LOC101109936, MSTRG.125394, MSTRG.125388, MSTRG.125392, MSTRG.83276, and MSTRG.81824, which are involved in the process by which FecB mutations affect follicular development. The intersection of BB_F vs. BB_L and BB_L vs. WW_L included LOC105610376, MSTRG.26398, MSTRG.127921, and MSTRG.134082, which are involved in ovulation in FecB mutant ewes ( Figure 2G ).

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

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 top 50 enriched terms of the target genes included tight junction, steroid hormone binding, and transcription elongation factor complex ( Supplementary 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 , four DEGs and six genes targeted by DELs are shown in detail. The targets of lncRNAs MSTRG.125392, MSTRG.125394, and MSTRG.83276 were cis-regulated, while others were trans-regulated. The expression patterns of most lncRNAs were similar to those of their corresponding target genes during the follicular–luteal transition.

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

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^−ΔΔCt 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 ² > 0.93, p < 0.05).