Summary of raw sequence reads
After removing low-quality sequences, a total of 288,342,450, 250,073,062, 289,224,844 and 277,834,922 clean reads with greater than 91.91% of Q30 were obtained in MM_F, MM_L, ww_F and ww_L, respectively. Approximately 86% to 93% of the reads were successfully mapped to the Ovis aries reference genome (Table 1).
Differential expression analysis of lncRNAs and mRNAs
A total of 21,282 lncRNAs (including 1797 known lncRNAs and 19,485 novel lncRNAs) and 43,674 mRNAs were identified from four groups (MM_F, MM_L, ww_F and ww_L) (Supplementary material 1A, 1B, 2). Overall, 10,785 intronic_lncRNAs, 7091 lincRNAs and 1609 antisense_lncRNAs were screened in the novel lncRNAs (Fig. 1A). Four comparison groups were set based on their genotypes and estrous cycle, MM_F_P vs MM_L_P, MM_F_P vs ww_F_P, MM_L_P vs ww_L_P, and ww_F_P vs ww_L_P. For MM_F_P vs MM_L_P, 17 lncRNAs and 414 mRNAs were upregulated, 22 lncRNAs and 350 mRNAs were downregulated (Fig. 1B, Supplementary material 3A, 4A). For MM_F_P vs ww_F_P, 11 lncRNAs and 122 mRNAs were upregulated, 29 lncRNAs and 116 mRNAs were downregulated (Fig. 1C, Supplementary material 3B, 4B). For MM_L_P vs ww_L_P, 12 lncRNAs and 86 mRNAs were upregulated, 18 lncRNAs and 154 mRNAs were downregulated (Fig. 1D, Supplementary material 3C, 4C). For ww_F_P vs ww_L_P, 64 lncRNAs and 208 mRNAs were upregulated, 32 lncRNAs and 388 mRNAs were downregulated (Fig. 1E, Supplementary material 3D, 4D). All DE lncRNAs (P <0.05) and mRNAs (P <0.05) were statistically significant.
Venn diagram visually showed the numbers of common and unique DE lncRNA_targets and mRNAs among four comparison groups, as shown in Fig. 2A-2D. In addition, distribution of these DE lncRNAs and mRNAs on chromosomes showed they were located on NC_019459.2, NC_019460.2, NC_019458.2 with greater proportion (Fig. S1-S8), and reliable for their exon size and ORF length mostly within 1000 bp (Fig. S9).
GO analysis of the biological function of DE lncRNAs and mRNAs
GO annotation enrichment was used to describe functions of the DE lncRNAs and mRNAs involved in cellular components, molecular function and biological processes, as shown in Fig. 3. Between MM_F_P and MM_L_P, targeted genes for DE lncRNAs were most enriched, and the terms were related to regulation of trans-membrane transport, antigen processing and presentation, immune system process. DE mRNAs were most enriched, the meaningful terms were related to the regulation of C-terminal protein methylation, C-terminal protein amino acid modification, post-translation protein modification, cellular macromolecular complex assembly and cellular metabolic process (Fig. 3A, Supplementary material 5A, 6A).
Between MM_F_P and ww_F_P, targeted genes for DE lncRNAs were enriched, the terms were related to regulation of protein complex assembly and biogenesis, protein complex subunit organization, spindle assembly involved in mitosis process. DE mRNAs were most enriched, the meaningful terms were related to regulation of secondary metabolic and biosynthetic process, viral protein processing, single-organism biosynthetic process (Fig. 3B, Supplementary material 5B, 6B).
Between MM_L_P and ww_L_P, targeted genes for DE lncRNAs were enriched, the terms were mainly related to regulation of single organism signaling, signal transduction, cellular response to stimulus and cellular communication. DE mRNAs were enriched, the meaningful terms were related to regulation of RNA methylation, metabolic process, organic substance metabolic process (Fig. 3C, Supplementary material 5C, 6C).
Between ww_F_P and ww_L_P, targeted genes for DE lncRNAs were enriched, the terms were related to regulation of immune response, glycerolipid metabolic process, cellular response to abiotic stimulus. DE mRNAs were enriched, the terms were related to regulation of nucleosome and chromatin assembly, nucleosome organization process (Fig. 3D, Supplementary material 5D, 6D).
KEGG pathway analysis
KEGG is a primary public pathway database to understand potential function of DE genes. The top 20 pathways were showed in Fig. 4-7. Between MM_F_P and MM_L_P, DE lncRNA targeted mRNAs were associated with pathways such as cell adhesion molecules (CAMs), glutathione (GSH) metabolism and bile secretion pathway (Fig. 4A, Supplementary material 7A). DE mRNAs were enriched in RNA transport, protein processing in endoplasmic reticulum and carbon metabolism pathway (Fig. 4B, Supplementary material 8A).
Between MM_F_P and ww_F_P, DE lncRNA targeted mRNAs were associated with pathways such as phosphatidylinositol signaling system, TNF signaling and p53 signaling pathway (Fig. 5A, Supplementary material 7B). With regard to DE mRNAs, which were enriched in 2-oxocarboxylic acid metabolism, RNA transport, endocrine and other factor-regulated calcium reabsorption pathways (Fig. 5B, Supplementary material 8B).
Between MM_L_P and ww_L_P, DE lncRNA targeted mRNAs were associated with pathways such as olfactory transduction, gap junction and thyroid hormone signaling pathway (Fig. 6A, Supplementary material 7C). With regard to DE mRNAs, which were enriched in ubiquitin mediated proteolysis, vasopressin-regulated water reabsorption, non-homologous end-joining and cell cycle (Fig. 6B, Supplementary material 8C).
Between ww_F_P and ww_L_P, DE lncRNA targeted mRNAs were associated with pathways such as cell adhesion molecules (CAMs), GSH metabolism and tight junction pathway (Fig. 7A, Supplementary material 7D). DE mRNAs were enriched in spliceosome, notch signal pathway, RNA polymerase and adherens junction, ras signaling pathway (Fig. 7B, Supplementary material 8D).
Hence, we acquired DE mRNAs closely related to reproductive signal pathways on the whole from above four comparison groups (Table S1).
Interaction Analysis of DE lncRNAs-mRNAs and function prediction
To better understand the relationship between lncRNA and mRNA, we constructed network of co-expression of DE lncRNAs and DE target mRNAs, after screening the overlaps between target mRNAs and DE mRNAs in each comparison group, which indicated regulation of lncRNA and mRNA in reproduction (|Pearson correlation| >0.95). Between MM_F_P and MM_L_P, a total of 5 DE lncRNAs and 9 DE mRNAs were involved in the network, and it consists of 9 edges (Fig. 8A, Supplementary material 9A). Between MM_F_P and ww_F_P, a total of 10 DE lncRNAs and 14 DE mRNAs were involved in the network, and it consists of 18 edges (Fig. 8B, Supplementary material 9B). Between MM_L_P and ww_L_P, a total of 6 DE lncRNAs and 10 DE mRNAs were involved in the network, and it consists of 10 edges (Fig. 8C, Supplementary material 9C). Between ww_F_P and ww_L_P, a total of 30 DE lncRNAs and 12 DE mRNAs were involved in the network, and it consists of 47 edges (Fig. 9, Supplementary material 9D).
Based on analysis of co-expression, we screened DE lncRNAs and the DE target mRNAs that closely related to reproductive pathways in different reproductive cycles and genotypes sheep. In MM sheep, related pathways were enriched with 4 DE lncRNAs (XLOC_466330, XLOC_391199, XLOC_503926, XLOC_517836) and 4 DE targets (RRM2B, GSTK1, STMN1, RAG2) (Table 2). In ww sheep, related pathways were enriched with 6 DE lncRNAs (XLOC_532771, XLOC_347557, XLOC_339502, XLOC_187711, XLOC_028449, 105604037) and 7 DE targets (GPX2, LOC101111397, RRM2B, GPX1, GSTK1, MGST1, DLG4) (Table 3). Additionally, related pathways were enriched by 7 DE lncRNAs (XLOC_448033, XLOC_252740, XLOC_241702, XLOC_079038, XLOC_078000, XLOC_065274, XLOC_009682) and 9 DE targets (DCT, PLCB4, PIK3CG, S1PR1, BRCA1, OSMR, PDGFD, RRM2B, CHEK1) in two groups of sheep (MM vs ww) at follicular phase (Table 4). And they were also enriched by 3 DE lncRNAs (XLOC_283279, XLOC_187695, XLOC_023278) and 11 DE targets (PRKACB, PRKAA1, PPP2R2A, PLCB4, NOS3, NCOA2, MAP2K6, MAP2K1, LOC101121082, LOC101111988, CAMKK2) in two groups of sheep (MM vs ww) at luteal phase (Table 5).