The experimental animals used in this study, Hu Sheep and Pishan Red Sheep, were provided by Xiyu Muyangren Agriculture and Animal Husbandry Co, Ltd. in Pishan County, Hotan Prefecture. Experimental ewes were selected based on lambing records under uniform feeding conditions. The cohort comprised Hu Sheep (FecB^BB), heterozygous FecB^B+ Pishan Red Sheep, and wild-type FecB⁺⁺ Pishan Red Sheep, all being healthy, non-related nulliparous ewes aged 7–8 months with consistent body conformation, the number of animals per group was n = 4.

Estrus synchronization was performed on all 18 selected ewes using a single prostaglandin F2α (PGF2α)-based protocol. The luteal phase was initially confirmed through visual assessment of estrus signs. Each ewe received an intramuscular injection of 0.1 mg PGF2α. At 48 h post-injection, estrus behavior was detected using the ram exposure method to determine optimal timing for sample collection.

Four estrous Pishan Red Sheep with FecB^B+ genotype (Group A), four with FecB⁺⁺ genotype (Group B), and four Hu Sheep (Group C) were selected and humanely slaughtered. Ovarian tissue samples were promptly collected, and tissues surrounding the ovaries were removed. The samples were rinsed with phosphate-buffered saline (PBS), transferred to RNase-free cryotubes, and immediately snap-frozen in liquid nitrogen. All animal experimental procedures strictly adhered to the Institutional Animal Care and Use Guidelines of Xinjiang Agricultural University.

The experiment employed the following instrumentation to meet operational requirements across all experimental stages: PCR system (Life Technologies), dry bath incubator (TIANGEN), mini centrifuge (TIANGEN), vortex mixer (Vortex-Genie 2, Scientific Industries), adjustable pipettes (RAININ), magnetic separation rack (Life Technologies), tube rotator (QLinbeier), Qsep-400 automated electrophoresis system (BiOptic Inc.), and Qubit 3.0 Fluorometer (Thermo Fisher Scientific).

Total RNA was extracted from ovine ovarian tissues using TRIzol reagent (Thermo Fisher Scientific, USA). RNA concentration was measured with a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, USA), while RNA integrity was assessed using an Agilent 2,100 Bioanalyzer (Agilent Technologies, USA). Library quality control was performed with the Qsep-400 system (BiOptic Inc., Taiwan, China). Qualified libraries were subsequently sequenced on the NovaSeq X Plus platform (Illumina, USA) with 150 bp paired-end reads. The RNA-Seq data of Pishan Red Sheep and Hu Sheep have been deposited in the Sequence Read Archive (SRA) database of the National Center for Biotechnology Information (NCBI) under BioProject accession number PRJNA1226111.

Following sequencing, the data were analyzed using the bioinformatics pipeline provided by the BMKCloud platform¹ from Biomarker Technologies. The FASTQ data obtained from high-throughput sequencing were analyzed to predict and identify multiple RNA species, followed by quantitative profiling of distinct RNA categories through alignment with established RNA annotation databases. In this study, differentially expressed genes (DEGs) and lncRNAs were screened using stringent thresholds of |Fold Change| ≥ 1.5 and p-value <0.01. Differential comparison groups: A-L vs. B-L, A-L vs. C-L, B-L vs. C-L. Target genes of lncRNAs were predicted using two distinct approaches. The first method focused on cis-regulatory mechanisms, identifying adjacent protein-coding genes located within a 100 kb genomic window upstream/downstream of lncRNA loci as potential targets, based on the positional relationship between lncRNAs and their neighboring genes. The second approach employed trans-regulatory prediction by analyzing inter-sample expression correlation between lncRNAs and mRNAs, identifying potential target genes through co-expression patterns.

The differentially expressed lncRNAs and mRNAs identified were further analyzed to investigate their expression dynamics across different genotypes. Gene expression dynamics were visualized through trend analysis plots to delineate specific expression pattern alterations.

The Gene Ontology (GO) database, established in 2000 by the Gene Ontology Consortium, constitutes a cross-species standardized annotation framework for genes and their products, structured into three orthogonal categories: Biological Process (BP), Molecular Function (MF), and Cellular Component (CC). Following gene annotation, functional categorization of differentially expressed genes (DEGs) was conducted at the second-level Gene Ontology (GO) classification. Enrichment analysis of DEG sets from each comparison group was performed using the hypergeometric test implemented in the clusterProfiler R package (v4.0.5), with results visualized as bar plots. A lower q-value indicates greater statistical significance of the enriched functional pathways, thereby facilitating the functional inference of candidate genes. Pathway enrichment analysis of differentially expressed genes was conducted using the KEGG (Kyoto Encyclopedia of Genes and Genomes) database. KEGG pathways were utilized as ontological units, and the hypergeometric test was applied to identify significantly enriched pathways (adjusted p-value < 0.01) by comparing against the genomic background. Concurrently, KEGG pathway annotation results for both target genes of differentially expressed lncRNAs and mRNAs were systematically categorized according to KEGG pathway classifications, thereby identifying key biochemical metabolic and signal transduction pathways in which these genes are functionally implicated. The present investigation primarily centers on biological pathways associated with reproductive physiology and ovarian functional regulation.

This study implemented an integrated methodology combining systematic literature review and multi-omics database mining to identify functionally pivotal genes associated with ovine reproductive performance. The selection criteria prioritized core constituents within key reproductive pathways, including the estrogen signaling pathway, MAP kinase cascade, cytokine-cytokine receptor interaction, and Hippo signaling pathway. Corresponding nucleotide sequences and spatiotemporal expression profiles were retrieved from the NCBI GenBank and GEO database. Differential expression analysis of mRNA and lncRNA was conducted using DESeq2 and edgeR packages on ovarian tissue samples from FecB⁺⁺-genotype Pishan Red Sheep, FecB^B+-genotype Pishan Red Sheep, and Hu Sheep, based on high-throughput sequencing data to identify genes with statistically significant expression changes.

During the integrated analysis of hub mRNAs and key lncRNAs, significantly differentially expressed mRNAs and lncRNAs were first identified. Subsequently, a cross-correlation heatmap was constructed to visualize expression covariation patterns between these two RNA classes.

To validate the expression levels of key candidate genes identified from screening, we performed qRT-PCR analysis on additional ovine ovarian tissue samples, thereby confirming the reliability of both sequencing and bioinformatics results. Using β-Actin as the internal reference, the relative gene expression levels were calculated for each sample and group. The processed data were organized in Excel and subsequently analyzed for visualization using Origin 2024 software. Primer sequences are provided in Table 1. Primers were designed using Premier 6 (Premier Biosoft, USA) and synthesized by Biomarker Technologies (Beijing, China). cDNA was reverse transcribed from total RNA, followed by RT-PCR amplification. Quantitative PCR was performed under the following conditions: initial denaturation at 95 °C for 3 min; 39 cycles of 95 °C for 10 s (denaturation), 60 °C for 30 s (annealing/extension) with plate reading; followed by a melting curve analysis (60–95 °C, increment of 1 °C per cycle, 4 s hold time).

The detection conclusions were based on Biomarker Beijing’s sequencing sample detection standard requirements and constituted a comprehensive evaluation of the tested samples. The samples with Test Result A meet the quality requirements for library construction, and the total quantity is sufficient for two or more standard library preparations. RNA quantification via NanoDrop demonstrated a concentration range of 98.9–652.8 ng/μL with total yields ranging from 1.48 to 9.79 μg. RNA Integrity Numbers (RIN) were distributed between 6.1 and 8.8, while 28S/18S rRNA ratios spanned 1.0 to 4.03. Detailed mRNA extraction parameters are presented in Table 2.

Strand-specific RNA sequencing libraries were constructed and sequenced on the NovaSeq X Plus platform (Illumina). Following stringent quality control, a total of 204.58 Gb clean data were obtained with Q20 > 97.45% and Q30 > 93.31%, while maintaining low N ratios and stable GC content, indicating high-quality sequencing outputs (detailed metrics in Table 3). Reference genome alignment statistics revealed mapping rates ranging from 81.40 to 91.19% against the Ovis aries reference genome assembly (Ovis_aries.GCF_000298735.2_Oar_v4.0.genome.fa). Complete alignment results are tabulated in Table 4.

The coding potential of putative lncRNAs was systematically assessed using four distinct computational tools: CNCI, CPC2, Pfam Scan, and CPAT. The intersection analysis of these predictions (Figure 1A) identified 34,651 high-confidence lncRNAs with conserved non-coding features. As depicted in Figure 1B, the composition of lncRNA classes was: intergenic lncRNAs (15,611, 45.1%), antisense lncRNAs (3,462, 10%), intronic lncRNAs (15,062, 43.5%), and sense lncRNAs (516, 1.5%). As illustrated in Figure 1C, the middle ring depicts the chromosomal distribution of differentially expressed mRNAs, while the inner ring represents the chromosomal distribution of differentially expressed lncRNAs. Within the rings, red and green denote upregulated and downregulated mRNAs, respectively; yellow and blue indicate upregulated and downregulated lncRNAs, correspondingly. As shown in Figure 1D, the expression levels of lncRNAs were consistently lower than mRNAs overall. Figure 1E illustrates the transcript length distribution, revealing that lncRNAs are primarily distributed between 200 and 400 bp, whereas mRNAs predominantly range from ≥3,000 bp.

Comparative analysis revealed distinct differential expression profiles across the three comparison groups: A-L vs. B-L (Figure 2A), A-L vs. C-L (Figure 2B), and B-L vs. C-L (Figure 2C). Specifically, 92 differentially expressed lncRNAs (56 upregulated, 36 downregulated), 269 lncRNAs (144 upregulated, 125 downregulated), and 466 lncRNAs (217 upregulated, 249 downregulated) were identified in the respective groups. Concurrently, 141 differentially expressed mRNAs (60 upregulated, 81 downregulated), 1,397 mRNAs (685 upregulated, 712 downregulated), and 1,821 mRNAs (830 upregulated, 991 downregulated) were detected in these comparisons. To comprehensively characterize the transcriptional changes in ovarian tissues during estrus, we identified differentially expressed long non-coding RNAs (lncRNAs) and messenger RNAs (mRNAs). The top five upregulated and downregulated lncRNAs and mRNAs are presented in Supplementary Tables 1 and 2, respectively. Venn analysis demonstrated that 1,481 shared differentially expressed mRNAs (Figure 2D) and 698 shared differentially expressed lncRNAs (Figure 2E) overlapped across all three groups.

Expression trend analysis was performed on differentially expressed mRNAs and lncRNAs (Figures 3, 4). In the trend analysis result figures, each panel is arranged in the order of FecB⁺⁺-genotype Pishan Red Sheep (Group B), FecB^B+-genotype (Group A), and Hu Sheep (Group C). Clusters were selected based on the lambing performance patterns of these three types of sheep. The selection criteria were defined as follows: clusters exhibiting either consistently increasing or decreasing expression trends across all three groups, lower expression levels in Group B, but higher levels in Groups A and C, or higher expression levels in Group B, but lower levels in Groups A and C. This screening identified 888 differentially expressed mRNAs and 6,740 differentially expressed lncRNAs conforming to these expression trends.

LncRNAs orchestrate a broad spectrum of critical physiological processes through transcriptional and post-transcriptional regulation of target genes. Leveraging trend-based expression profiles, we performed systematic prediction of both cis-acting and trans-regulatory target genes for these prioritized lncRNAs. Cis-regulatory target genes were predicted by identifying protein-coding genes located within ±100 kb genomic regions flanking the lncRNAs using a Perl script, yielding 845 cis-target pairs. Furthermore, trans-regulatory targets were determined through Pearson correlation analysis of inter-sample co-expression patterns between lncRNAs and mRNAs, identifying 5,895 trans-target interactions.

Functional enrichment analysis was performed on the 888 trending differentially expressed genes (DEGs), revealing significant enrichment of 60 GO terms and 285 KEGG pathways (Figure 5). Significant enrichments were observed in biological processes (Biological Process, BP), including key terms such as “cellular process” and “metabolic process.” Furthermore, KEGG pathway analysis revealed differential enrichment across multiple pathways, including cysteine and methionine metabolism, axon guidance, and Fc gamma receptor (FcγR)-mediated phagocytosis, among others. Notably, reproduction-associated pathways including progesterone-mediated oocyte maturation, estrogen signaling pathway, and gonadotropin-releasing hormone (GnRH) signaling pathway were prominently enriched in the ovarian mRNA enrichment analysis. These findings suggest that genes within these pathways may play pivotal roles in ovarian physiological activities and reproductive regulation during the estrus phase in sheep.

Functional enrichment analysis of the 845 cis-target genes and 5,859 trans-target genes predicted from the 6,704 trend-conforming lncRNAs revealed significant enrichments of 60 and 63 GO terms, and 289 and 330 KEGG pathways, respectively. GO enrichment analysis of cis-target genes (Figure 6) demonstrated significant enrichment in biological processes (BP), including “cellular process,” “single-organism process,” “biological regulation,” and “metabolic process,” indicating potential functional roles of these genes in diverse biological processes within ovarian tissues. KEGG pathway analysis of cis-target genes (Figure 6) revealed differential enrichment across multiple pathways, including “Cytokine-cytokine receptor interaction,” “Hippo signaling pathway,” and “Estrogen signaling pathway,” among others. GO enrichment analysis of trans-acting target genes (Figure 7) revealed that in the biological process category, entries such as “cellular process” and “single-organism process” exhibited relatively higher numbers of enriched lncRNA trans-acting target genes. This observation suggests that during physiological activities in ovarian tissues at the estrus phase, these biological processes may be significantly influenced by genes under trans-regulatory control of corresponding lncRNAs. KEGG pathway analysis (Figure 7) revealed significant enrichment of lncRNA trans-target genes in ovarian tissues for key pathways including “Cytokine-cytokine receptor interaction,” “MAPK signaling pathway,” and “cGMP-PKG signaling pathway.” These findings suggest these pathways may exert central regulatory functions in ovarian physiological processes during the estrous cycle. Furthermore, pathways including the “MAPK signaling pathway” and “cAMP signaling pathway” demonstrated consistent enrichment in both cis- and trans-target gene KEGG analyses, suggesting their potential roles in regulating signal transduction, cellular proliferation, and differentiation within ovarian tissues. The co-enrichment of these pathways further elucidates the multifaceted mechanisms by which lncRNAs orchestrate ovarian functional regulation.

Functional enrichment analysis identified a cohort of significantly expressed mRNAs and lncRNAs in ovarian tissues during the estrous cycle, particularly enriched in key KEGG pathways including “Cytokine-cytokine receptor interaction,” “Hippo signaling pathway,” and “MAPK signaling pathway.”

Transcriptomic profiling at the mRNA level identified critical genes – including GDF9, GRIA4, HOXC9, HOXD3, MAPK8IP3, AMH, ANGPT2, FGF14, MAPK8IP1, MMP9, and BRINP3 demonstrating multifaceted roles in ovarian folliculogenesis, cell cycle modulation, and MAPK-mediated signal transduction cascades. Bar graph analysis of mRNA expression levels (Figure 8) revealed significant differences in gene expression across breeds, with these transcriptional disparities potentially linked to variations in their reproductive performance.

The AMH gene exhibited significantly higher expression levels in FecB⁺⁺-genotype Pishan Red Sheep (Group B) compared to FecB^B+-genotype (Group A) and Hu Sheep (Group C). This transcriptional pattern suggests that AMH may play a more prominent role in follicular quiescence regulation within mono-ovulatory breeds, while its activity in poly-ovulatory breeds is likely modulated through complementary regulatory mechanisms. The ANGPT2 gene displayed the highest expression levels in FecB⁺⁺-genotype Pishan Red Sheep, compared to FecB^B+-genotype and Hu Sheep. This expression pattern may correlate with its dual roles in modulating ovarian angiogenesis and regulating follicular maturation during reproductive cycles. The FGF14 gene exhibited the highest expression levels in FecB⁺⁺-genotype Pishan Red Sheep, whereas significantly lower expression was observed in FecB^B+-genotype and Hu Sheep. This differential expression pattern may be associated with FGF14’s regulatory functions in cellular proliferation and differentiation. The HOXC9 and HOXD3 genes exhibited significantly lower expression levels in FecB⁺⁺-genotype Pishan Red Sheep compared to other genotypes, suggesting potential associations with their regulatory roles in embryonic morphogenesis and cellular differentiation processes. MAPK8IP1 and MMP9 exhibited the highest expression levels in FecB⁺⁺-genotype Pishan Red Sheep, potentially linked to their functional roles in follicular maturation, ovulation dynamics, and extracellular matrix remodeling during ovarian tissue homeostasis.

Critical lncRNAs identified through systematic analysis include MSTRG.61044.1, MSTRG.2677.1, MSTRG.23016.1, MSTRG.27015.1, MSTRG.60286.1, and MSTRG.15154.3, which exhibited conserved expression patterns across genotype groups. The results (Figure 9) demonstrated significant differential expression of these lncRNAs across ovine breeds, a transcriptional divergence potentially associated with their reproductive performance variations.

Specifically, the lncRNA MSTRG.61044.1, which is highly expressed in FecB⁺⁺-genotype Pishan Red Sheep, demonstrates dual regulatory mechanisms: it exerts cis-regulatory control over the adjacent MAPK8IP1 gene, while trans-regulating reproduction-associated genes such as AMH and CCL25. These findings implicate MSTRG.61044.1 as a potential hub regulator orchestrating follicular selection and gonadotropin-responsive signaling cascades. Further investigation revealed that MSTRG.23016.1, also highly expressed in FecB⁺⁺-genotype Pishan Red Sheep, exhibits functional modulation of MMP9 and ANGPT2 transcripts, suggesting its potential involvement in extracellular matrix remodeling and angiogenic niche regulation during ovarian tissue homeostasis. Simultaneously, MSTRG.15154.3 dynamically regulates granulosa cell proliferation and lineage commitment through cis-mediated transcriptional control of ERBB4, establishing a mechanistic link between lncRNA-driven chromatin interactions and follicular developmental competence. In contrast, Hu Sheep-specific highly expressed lncRNAs exhibited distinct functional orientations: MSTRG.2677.1 modulates ovarian stromal cell homeostasis through trans-regulation of HGF and BRINP3, potentially orchestrating extracellular matrix dynamics and stromal-epithelial crosstalk during follicular wave progression. MSTRG.27015.1 and MSTRG.60286.1 were found to target MAPK8IP3 and PPP3CB, respectively. This mechanistic association suggests their potential roles in regulating cell cycle checkpoints and orchestrating oocyte maturation through calcium-dependent signaling cascades. Notably, MAPK signaling pathway-associated genes (MAPK8IP1 and MAPK8IP3) exhibit breed-specific differential regulation by distinct lncRNAs, which may reflect evolutionary divergence in follicular development synchronization mechanisms across ovine breeds.

Figure 10 demonstrates the correlation heatmap between key differentially expressed mRNAs and lncRNAs identified in this study. The heatmap revealed strong correlations between key differentially expressed mRNAs and lncRNAs, with these RNA molecules distinctly clustering into two groups based on their expression correlation patterns. Specifically, lncRNAs MSTRG.23016.1 and MSTRG.61044.1 demonstrated strong positive correlations with key reproductive regulators including AMH, ANGPT2, MMP9, FGF14, and MAPK8IP1, while exhibiting robust negative correlations with the HOXD3 gene. Furthermore, lncRNAs MSTRG.2677.1, MSTRG.27015.1, MSTRG.60286.1, and MSTRG.50603.1 exhibited strong positive correlations with GRIA4, MAPK8IP3, GDF9, HOXC9, and BRINP3 mRNAs. Notably, lncRNA MSTRG.50603.1 also exhibited strong negative correlations with MAPK8IP1, FGF14, and MMP9 mRNAs. Furthermore, MSTRG.15154.3 displayed strong positive correlations with FGF14 and MAPK8IP1, while showing negative correlations with GRIA4, MAPK8IP3, GDF9, HOXC9, and BRINP3. Such correlation-based co-expression analysis has proven instrumental in pinpointing putative functional interactions within RNA regulatory networks, holding significant implications for deciphering transcriptional circuitry and elucidating molecular mechanisms governing ovarian folliculogenesis. These findings provide critical leads for subsequent functional validation studies, facilitating in-depth exploration of the intricate interplay between mRNAs and lncRNAs, as well as their roles in disease pathogenesis and progression.

To assess the reliability of RNA sequencing, we randomly selected 4 mRNAs (MMP9, AMH, ANGPT2, MAPK8IP3) and 1 lncRNA (MSTRG.60286.1) for qRT-PCR validation. The results demonstrated consistent expression patterns across different reproductive stages between qRT-PCR and RNA-Seq data (Figures 11, 12), confirming the robustness of our transcriptomic analysis.