Skip to main content

Pineal gland transcriptomic profiling reveals the differential regulation of lncRNA and mRNA related to prolificacy in STH sheep with two FecB genotypes

Abstract

Background

Long noncoding RNA (lncRNA) has been identified as important regulator in hypothalamic-pituitary-ovarian axis associated with sheep prolificacy. However, little is known of their expression pattern and potential roles in the pineal gland of sheep. Herein, RNA-Seq was used to detect transcriptome expression pattern in pineal gland between follicular phase (FP) and luteal phase (LP) in FecBBB (MM) and FecB++ (ww) STH sheep, respectively, and differentially expressed (DE) lncRNAs and mRNAs associated with reproduction were identified.

Results

Overall, 135 DE lncRNAs and 1360 DE mRNAs in pineal gland between MM and ww sheep were screened. Wherein, 39 DE lncRNAs and 764 DE mRNAs were identified (FP vs LP) in MM sheep, 96 DE lncRNAs and 596 DE mRNAs were identified (FP vs LP) in ww sheep. Moreover, GO and KEGG enrichment analysis indicated that the targets of DE lncRNAs and DE mRNAs were annotated to multiple biological processes such as phototransduction, circadian rhythm, melanogenesis, GSH metabolism and steroid biosynthesis, which directly or indirectly participate in hormone activities to affect sheep reproductive performance. Additionally, co-expression of lncRNAs-mRNAs and the network construction were performed based on correlation analysis, DE lncRNAs can modulate target genes involved in related pathways to affect sheep fecundity. Specifically, XLOC_466330, XLOC_532771, XLOC_028449 targeting RRM2B and GSTK1, XLOC_391199 targeting STMN1, XLOC_503926 targeting RAG2, XLOC_187711 targeting DLG4 were included.

Conclusion

All of these differential lncRNAs and mRNAs expression profiles in pineal gland provide a novel resource for elucidating regulatory mechanism underlying STH sheep prolificacy.

Background

Reproduction, one of the major factors significantly affecting profitability of sheep production, is a complicated physiological process and determined by the integrated hypothalamic-pituitary-ovarian axis in breeding season [1]. Reproductive traits like litter size directly determine benefit of sheep production, are controlled by poly-gene at the micro level. How to undertake at molecular level to improve reproduction, thereby serving macro production is a hotspot in recent years. BMPRIB, BMP15 [2] and GDF9 [3] are major fecundity genes which significantly influence sheep prolificacy. FecB is a mutation in BMPRIB occurring in base 746 from A to G, one copy of this mutation significantly increases ovulation rate in sheep about 1.5 and two copies by 3.0 [4]. To date, this mutation has been detected in diverse sheep species such as Booroola Merino sheep (Australia) [5], Small Tail Han (STH) and Hu sheep (China) [6]. Wherein STH sheep is a famous native breed with year-round estrus and high fecundity, being officially recognized as one of the polytocous breeds in China. The average litter size and lambing rate of STH sheep are 2.61, 286.5%, respectively [7]. There are three genotypes based on effects of FecB mutation in STH sheep, namely FecBBB (with two-copy FecB mutation), FecBB+ (with one-copy FecB mutation) and FecB++ (with no FecB mutation), which is closely correlated with litter size of ewes [8]. Therefore, this breed can be used as a classic model for study molecular mechanism of FecB gene regulation of reproductive traits in sheep.

Long noncoding RNA (lncRNA) is polymerase II transcript with length longer than 200 nucleotides that lacks the protein coding ability, its expression has high tissue specificity and distributes in cytoplasm or nucleus [9]. LncRNA is proposed to be the largest transcript class in mammalian transcriptome [10], less than 2% of mammalian genome actually code for protein, 70–90% is transcribed in some context as lncRNA, originally thought to be ‘transcriptional noise’ in genome. Subsequently, studies have gradually shown that lncRNA exerts important roles in various biological processes such as cell proliferation, apoptosis and differentiation [11], signal transduction [12], immune regulation [13]. In terms of reproduction, there have many reports on lncRNA. For example, Miao et al. (2017) compared transcripts in ovaries of low fecundity ewes and high fecundity ewes, found that differentially expressed (DE) lncRNA significantly enriched in the oxytocin signaling pathway [14]. Then, Feng et al. (2018) identified 5 lncRNAs and 76 mRNAs in ovaries of Hu sheep with high and low prolificacy, respectively [15]. Yang et al. (2020) analyzed lncRNA and mRNA in male sheep pituitary and found that 5 candidate lncRNAs and their targeted genes enriched in growth and reproduction related pathways [16]. Su et al. (2020) screened differential lncRNA through high-throughput sequencing, concluded that XLOC-2222497 and its target AKR1C1 could interact with progesterone in porcine endometrium for controlling pregnancy maintenance [17]. These studies indicated the presence and role of lncRNA in reproductive tissues. It is known that the sheep pineal gland as an important reproductive-related gland, that is closely related to hormone and signal transduction. However, studies on function of sheep lncRNA in this organ are limited.

In light of this, the study presented herein was focused on analyzing transcriptomics of pineal gland in STH sheep with FecBBB (MM) and FecB++ (ww) genotypes, to determine the DE lncRNAs and genes, and predict their potential function that related to reproduction. Which is essential for better understanding the molecular mechanisms by lncRNAs regulate sheep reproduction with different genotypes, also providing insight for other female mammals.

Results

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.10 to 92.89% of the reads were successfully mapped to the Ovis aries reference genome (Table 1).

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

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, B, 2). Overall, 10,785 intronic lncRNAs, 7091 intergenic lncRNAs (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_FP vs MM_LP, MM_FP vs ww_FP, MM_LP vs ww_LP, and ww_FP vs ww_LP. For MM_FP vs MM_LP, 17 lncRNAs and 414 mRNAs were upregulated, 22 lncRNAs and 350 mRNAs were downregulated (Fig. 1b, Supplementary material 3A, 4A). For MM_FP vs ww_FP, 11 lncRNAs and 122 mRNAs were upregulated, 29 lncRNAs and 116 mRNAs were downregulated (Fig. 1c, Supplementary material 3B, 4B). For MM_LP vs ww_LP, 12 lncRNAs and 86 mRNAs were upregulated, 18 lncRNAs and 154 mRNAs were downregulated (Fig. 1d, Supplementary material 3C, 4C). For ww_FP vs ww_LP, 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.

Fig. 1
figure1

Gene expression characterization. a The classification and proportion of novel lncRNAs. b Histogram representing the numbers of upregulated and downregulated lncRNAs and mRNAs in sheep pineal body between MM_F_P and MM_L_P. c Histogram representing the numbers of upregulated and downregulated lncRNAs and mRNAs in sheep pineal body between MM_F_P and ww_F_P. d Istogram representing the numbers of upregulated and downregulated lncRNAs and mRNAs in sheep pineal body between MM_L_P and ww_L_P. e Histogram representing the numbers of upregulated and downregulated lncRNAs and mRNAs in sheep pineal body between ww_F_P and ww_L_P

Venn diagram visually showed the numbers of common and unique DE lncRNA_targets and mRNAs among four comparison groups, as shown in Fig. 2a-d. In addition, distribution of these DE lncRNAs and mRNAs on chromosomes showed they were located on Chr2 (NC_019459.2), Chr3 (NC_019460.2), Chr1 (NC_019458.2) with greater proportion (Figures S1, S2, S3, S4, S5, S6, S7, S8), and reliable for their exon size and ORF length mostly within 1000 bp (Figure S9).

Fig. 2
figure2

Venn diagram visualization of DE lncRNA_targets and mRNAs among four comparison groups. a Venn diagram representing the overlapping numbers of differentially expressed lncRNA_targets and mRNAs in MM_F_P vs MM_L_P. b Venn diagram representing the overlapping numbers of differentially expressed lncRNA_targets and mRNAs in MM_F_P vs ww_F_P. c Venn diagram representing the overlapping numbers of differentially expressed lncRNA_targets and mRNAs in MM_L_P vs ww_L_P. d Venn diagram representing the overlapping numbers of differentially expressed lncRNA_targets and mRNAs in ww_F_P vs ww_L_P

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_FP and MM_LP, 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).

Fig. 3
figure3

GO analyses of differentially expressed lncRNA targets and mRNAs. a The top 15 enrichment biological processes for differentially expressed lncRNA targets and mRNAs in MM_F_P vs MM_L_P. b The top 15 enrichment biological processes for differentially expressed lncRNA targets and mRNAs in MM_F_P vs ww_F_P. c The top 15 enrichment biological processes for differentially expressed lncRNA targets and mRNAs in MM_L_P vs ww_L_P. d The top 15 enrichment biological processes for differentially expressed lncRNA targets and mRNAs in ww_F_P vs ww_L_P

Between MM_FP and ww_FP, 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_LP and ww_LP, 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_FP and ww_LP, 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 Figs. 4, 5, 6, 7. Between MM_FP and MM_LP, 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).

Fig. 4
figure4

KEGG analyses of differentially expressed genes between MM_F_P and MM_L_P groups. a The top 20 KEGG enrichment pathways for differentially expressed lncRNA targets between MM_F_P and MM_L_P groups. b The top 20 KEGG enrichment pathways for differentially expressed mRNAs between MM_F_P and MM_L_P groups

Fig. 5
figure5

KEGG analyses of differentially expressed genes between MM_F_P and ww_F_P groups. a The top 20 KEGG enrichment pathways for differentially expressed lncRNA targets between MM_F_P and ww_F_P groups. b The top 20 KEGG enrichment pathways for differentially expressed mRNAs between MM_F_P and ww_F_P groups

Fig. 6
figure6

KEGG analyses of differentially expressed genes between MM_L_P and ww_L_P groups. a The top 20 KEGG enrichment pathways for differentially expressed lncRNA targets between MM_L_P and ww_L_P groups. b The top 20 KEGG enrichment pathways for differentially expressed mRNAs between MM_L_P and ww_L_P groups

Fig. 7
figure7

KEGG analyses of differentially expressed genes between ww_F_P and ww_L_P groups. a The top 20 KEGG enrichment pathways for differentially expressed lncRNA targets between ww_F_P and ww_L_P groups. b The top 20 KEGG enrichment pathways for differentially expressed mRNAs between ww_F_P and ww_L_P groups

Between MM_FP and ww_FP, 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_LP and ww_LP, 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_FP and ww_LP, 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_FP and MM_LP, 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_FP and ww_FP, 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_LP and ww_LP, 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_FP and ww_LP, 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).

Fig. 8
figure8

Construction of the DE lncRNAs-target mRNAs co-expression network. a Co-expression of DE lncRNA-mRNA after lncRNA targets coincided with DE mRNAs in MM_F_P vs MM_L_P. b Co-expression of DE lncRNA-mRNA after lncRNA targets coincided with DE mRNAs in MM_F_P vs ww_F_P. c Co-expression of DE lncRNA-mRNA after lncRNA targets coincided with DE mRNAs in MM_L_P vs ww_L_P. Tangerine and green represent upregulated and downregulated, respectively. Octagons and triangles represent lncRNAs and mRNAs, respectively

Fig. 9
figure9

Co-expression of DE lncRNA-mRNA after lncRNA targets coincided with DE mRNAs in ww_F_P vs ww_L_P. Tangerine and green represent upregulated and downregulated, respectively. Octagons and triangles represent lncRNAs and mRNAs, respectively

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, 105,604,037) 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).

Table 2 Summary of co-expression of differential genes closely related to reproductive cycle (follicular phase vs luteal phase) in MM sheep
Table 3 Summary of co-expression of differential genes closely related to reproductive cycle (follicular phase vs luteal phase) in ww sheep
Table 4 Summary of co-expression of differential genes closely related to reproduction in different genotypes (MM vs ww) sheep at follicular phase
Table 5 Summary of co-expression of differential genes closely related to reproduction in different genotypes (MM vs ww) sheep at luteal phase

Discussion

Studies have found that lncRNA is involved in multiple reproductive functions such as spermatogenesis [18], placentation [19], signaling pathway of sex hormone response [20, 21] and gonadgenesis [22]. It is known that the melatonin synthesized in pineal gland is closely related to the estrus cycle [23]. Herein, the study focused on examining expression profiles of pineal gland lncRNAs and mRNAs in sheep with two genotypes at different phases of estrous cycle using RNA-Seq technology. Analysis of relationship between DE lncRNAs and mRNAs by generating a co-expression network. To our knowledge, this is the first genome-wide analysis of pineal gland in sheep with different genotypes, and might provide valuable resource for searching functional lncRNAs associated with sheep prolificacy.

In present study, we screened 21,282 lncRNAs and 43,674 mRNAs. LncRNAs have synergetic relationship with mRNAs as most lncRNAs are located near protein-coding genes [24, 25]. Additionally, wide location of lncRNAs in chromosomes of sheep indicated its pluripotency. Obviously, distribution ratio of lncRNAs and mRNAs on Chr2 (NC_019459.2), Chr3 (NC_019460.2), Chr1 (NC_019458.2) were greater than those on other chromosomes, which could be explained by close relationship between three chromosomes and pineal gland function. The exon size and ORF length of lncRNAs and mRNAs are mostly within 1000 bp. These results showed the potential lncRNAs were reliable in the pineal gland.

Overall, we screened 135 (39 + 96) DE lncRNAs and 1360 (764 + 596) DE mRNAs in pineal gland at follicular and luteal phases between high and low prolificacy STH sheep (WW vs ww). GO annotation and KEGG enrichment analysis of top 20 terms indicated that DE mRNAs were enriched in reproduction-related pathways such as GnRH, cGMP-PKG, thyroid hormone, MAPK, phototransduction, circadian rhythm, steroid biosynthesis, hippo, mTOR and melanogenesis. It is well known that productive cycle of mammals is regulated through association or acting alone of hypothalamic-pituitary-thyroid (HPT) axis and hypothalamic-pituitary-gonadal (HPG) axis [26, 27]. In the HPT axis, thyrotropin-releasing hormone (TRH) produced in hypothalamus stimulates pituitary to secrete thyroid-stimulating hormone (TSH), which promotes TH synthesis in the thyroid gland [26, 28]. In the HPG axis, GnRH in hypothalamus regulates synthesis and secretion of FSH and LH in the anterior pituitary. These two hormones stimulate gonadal estrogen synthesis by binding to their receptors for affecting development and maturation of follicles and the ewes litter size. Estrogen in turn positively or negatively acts GnRH synthesis, and affects FSHβ gene expression, a hormone specific β subunit that is mainly regulated by GnRH [29, 30]. In the process, binding of GnRH to its receptor activates signaling cascades like MAPK, PI3K-Akt [31]. MAPK pathway is essential for cell proliferation and differentiation, survival, death and transformation [32, 33]. PI3K-Akt can interact with mTOR pathway to effectively regulate growth hormone in pituitary [34]. Additionally, pathways as hippo modulates organ size growth by controlling stem cell activity, proliferation and apoptosis. For instance, its’ effect on the development of pituitary progenitor cells [35]. Our results showed that DE genes like AKT3, MYC, PIK3CB, MAP2K2, PLCB1 and TEAD1 related to thyroid hormone, MAPK, cGMP-PKG, hippo, and up regulated, while CTNNB1, YAP1, PIK3CG, TEAD1, CAMK2A, CACNA1D mainly related to hippo, thyroid hormone, cGMP-PKG, AMPK, GnRH, oxytocin, circadian entrainment, and down regulated, which implied the important roles of these genes mainly involved in regulation of hormone-related pathways. It’s worth exploring their function in pineal gland as candidate genes.

Co-expression analysis of differential lncRNA-mRNA and functional prediction of target genes revealed that lncRNA affects sheep fecundity by modulating genes associated with above signaling pathways and biological processes. In FecBBB genotype sheep, XLOC_466330 and the targets (RRM2B, GSTK1) up regulated at follicular phase, which related to GSH metabolism. Whereas XLOC_391199 and the target (STMN1), XLOC_503926, XLOC_517836 and the target (RAG2) up regulated at luteal phase, which mainly enriched in MAPK, FoxO signaling pathways, respectively. In FecB++ genotype sheep, XLOC_347557 and the target (GPX2), XLOC_532771 and the targets (LOC101111397, RRM2B), XLOC_339502 and the target (GPX1), XLOC_028449 and the target (GSTK1) up regulated at follicular phase, which also related to GSH metabolism. Meanwhile, 105,604,037 and the target (MGST1), XLOC_187711 and the target (DLG4) down regulated at the same phase that related to GSH metabolism and hippo signaling. Wherein GSTK1 and RRM2B involved in GSH metabolism at follicular phase, but their targeted regulators lncRNAs were markedly different among two FecB genotypes. RRM2B gene encodes p53R2, and p53R2 is expressed at all phases of cell cycle to ensure ample supply of mitochondrial DNA [36]. GSTK1 gene encodes a member of GSTK superfamily of enzymes that function in cellular mitochondria and peroxisomes detoxification during GSH metabolism [37, 38], a critical pathway protecting cells from free radicals and oxidative damage, could increase intracellular NADPH [39]. With increase of NADPH oxidase, ROS level tend to be low, whereas the level of intracellular ATP enhanced, as well as mitochondrial activity, which promote oocyte maturation [40], and so forth, the other DE genes involved in GSH metabolism were also novel direction of interest for their effects on the downstream reproductive system.

Furthermore, DE target genes like STMN1 is a highly conserved gene that codes for cytoplasmic phosphoproteins, acting role in cell cycle progression, signal transduction and cell migration through diverse intracellular signaling pathways. Studies have found the potential role of STMN1 in regulation of hormone secretion in rodent pituitary and insulinoma cell lines [41]. Over-expression of STMN1 stimulates progesterone production by modulating the promoter activity of Star and Cyp11a1 in mouse granulosa cells [42]. Besides, RAG2 is indispensable for generation of antigen receptor diversity in immune cells [43]. We found STMN1, RAG2 were down regulated at follicular phase in FecBBB sheep, and mainly related to MAPK, FoxO signaling pathways, respectively. DLG4 was down regulated at follicular phase in FecB++ sheep and enriched in hippo signaling term. As known that DLG4 encodes a member of MAGUK family, is widely expressed and playing an essential role in regulation of cellular signal transduction, circadian entrainment [44]. Taken together, the DE lncRNAs identified in this study might cooperate with their target genes and DE genes to regulate pineal gland physiological function, and involved in hormone synthesis for effecting reproductive cycle and final lambing.

Conclusion

In summary, the pineal gland transcriptomic study reveals differential regulation of lncRNAs and mRNAs related to prolificacy in sheep with different FecB genotyping. We screened several sets of target genes of DE lncRNAs and DE genes under reproductive cycle and genotypes, they were annotated to multiple biological processes such as phototransduction, circadian rhythm, melanogenesis, GSH metabolism and steroid biosynthesis, which directly or indirectly participate in hormone activities to affect sheep reproductive performance. Additionally, we predicted potential role of these DE lncRNAs and constructed network of lncRNAs-mRNAs to expand our understanding. All of these differential lncRNAs and mRNAs expression profiles provide a novel resource for elucidating regulatory mechanism underlying STH sheep prolificacy.

Methods

Ethics statement

Experimental animals in this study were authorized by the Science Research Department (in charge of animal welfare issues) of the Institute of Animal Science, Chinese Academy of Agricultural Sciences (IAS-CAAS; Beijing, China). Additionally, ethical approval of animal survival and the sample collection was given by the animal ethics committee of IAS-CAAS (No. IAS2019–49).

Animals preparation

Animals were from a core population of STH sheep in Luxi district of Shandong province, China. We collected jugular vein blood of healthy non-pregnant sheep aged 2–4 years (n = 890), to identify the FecB genotypes using TaqMan probe [45]. Then, 12 sheep (6 MM and 6 ww, respectively) with no significant difference in age, weight, height, body length, chest circumference and tube circumference were selected for this experiment.

Twelve sheep were managed and raised on a farm, with free access to water and feed. All sheep were processed by estrus synchronization with Controlled Internal Drug Releasing device (CIDR, progesterone 300 mg, Inter Ag Co., Ltd., New Zealand) for 12 days. 3 MM and 3 ww ewes were euthanized (Intravenous pentobarbital at 100 mg/kg) on the 50th hour after CIDR removal, pineal tissues were collected (follicular phase, MM_FP and ww_FP, respectively). The other 6 sheep were euthanized (Intravenous pentobarbital at 100 mg/kg) on the 7th day after CIDR removal, and pineal tissues were collected (luteal phase, MM_LP and ww_LP, respectively) [21]. Obtained samples were stored immediately at − 80 °C for the next step.

RNA extraction and detection

Total RNA was extracted from 12 samples using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to manufacturer’s instruction. 1% agarose gel was used to monitor whether isolated RNA was degraded or contaminated. Quality, integrity and concentration of RNA were assessed by NanoPhotometer® spectrophotometer (IMPLEN, CA, USA), RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA) and Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA), respectively. Among them, the ratio of intact RNA with RIN ≥ 7, 28S: 18S ≥ 1.5:1.

Construction of RNA libraries and sequencing

A total amount of 3 μg RNA per sample was used as input material for the RNA sample preparation. Firstly, rRNA was removed by Epicentre Ribo-zero™ rRNA Removal Kit (Epicentre, USA) and rRNA free residue was cleaned up by ethanol precipitation. Subsequently, libraries were generated using the rRNA-depleted RNA by NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, USA) following manufacturer’s recommendation. After the clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia), libraries were then sequenced through an Illumina Hiseq 4000 platform and 150 bp paired-end reads were generated.

Reference genome mapping and transcriptome assembly

Raw reads of fast-q format were firstly processed through in-house perl scripts to obtain clean reads. At the same time, Q20, Q30 and GC content of the clean data were calculated. All downstream analyses were based on the high quality clean reads. HISAT2 v. 2.0.4 was used to align paired-end clean reads of each sample to sheep reference genome Oar v. 4.0 [46]. The mapped reads of each sample were assembled by StringTie v. 1.3.1 [46].

Identification of potential lncRNA candidates

Potential lncRNA candidates were identified using the following workflow. (1) Transcripts with uncertain chain direction were removed by Cuffmerge. (2) Transcripts length > 200 nt with exon number ≥ 2 were selected. (3) Cuffcompare v. 2.1.1. was used to compare different classes of class_code annotated by “i”, “u” and “x” that were retained, which corresponded to intronic, intergenic, and anti-sense transcripts, respectively. (4) Transcripts with FPKM ≥0.5 were obtained after calculating the expression level of each transcript by Cuffdiff v. 2.1.1. (5) Three tools of CNCI v.2.0 [47], CPC v. 2.81 [48] and PFAM v.1.3 [49] were used to predict the protein-coding potential. After that, Pfam was implemented to exclude transcripts that overlapped with any protein-coding genes. Intersection of these transcripts without coding potential was used as the lncRNA data set. Additionally, mRNAs were obtained from the same RNA-seq libraries in this study.

Analysis of differentially expressed genes

The fragments per kilobase of transcript per million reads mapped (FPKM) value was used to estimate the expression levels of transcripts (lncRNAs and mRNAs) [50]. For experiments with three biological replicates, the R package DEseq2 was used to identify differentially expressed transcripts after a negative binomial distribution [51]. LncRNAs and mRNAs with P-value < 0.05 and a fold change (FC) > 2.0 were considered as differentially expression between two groups of data.

Bioinformatics analysis

LncRNA targets could be predicted by the correlation or co-expression of lncRNA and mRNA expression levels between samples. Among them, Pearson correlation coefficient (PCC) was used to analyze the correlation between lncRNA and mRNA among samples, mRNAs with |PCC| ≥0.95 for functional enrichment to predict lncRNAs [52]. Statistical enrichment of DE lncRNA targets and DE mRNAs in GO annotation and KEGG pathway were evaluated, P-value ≤0.05 defined as the significant threshold, significance of the term enrichment analysis was corrected by FDR and corrected P-value (Q-value) was obtained [53].

Construction of co-expression networks

To predict function of DE lncRNAs and their target genes in sheep reproduction, a network based on lncRNAs and mRNAs was bulit using Cytoscape software (v. 3.8.0) [54].

Statistical analysis

All data were assessed as the “means ± SD”. Student’s t-test was performed and P < 0.05 was considered statistically significant.

Availability of data and materials

All data sets used and analyzed during the current study are available: data is available at the Sequence Read Archive (PRJNA679918).

Abbreviations

STH:

Small tailed Han sheep

FP:

Follicular phase

LP:

Luteal phase

MM:

FecBBB genotype

ww:

FecB++ genotype

LncRNA:

Long noncoding RNA

HPT:

Hypothalamic-pituitary-thyroid

HPG:

Hypothalamic-pituitary-gonadal

TSH:

Thyroid-stimulating hormone

FPKM:

The fragments per kilobase of transcript per million reads mapped

CIDR:

Controlled internal drug releasing device

References

  1. 1.

    Bartlewski PM, Baby TE, Giffin JL. Reproductive cycles in sheep. Anim Reprod Sci. 2011;124:259–68.

    CAS  PubMed  Article  Google Scholar 

  2. 2.

    Chu MX, Liu ZH, Jiao CL, He YQ, Fang L, Ye SC, Wang JY. Mutationsin BMPR-IB and BMP-15 genes are associated with litter size in small tailed Han sheep (Ovis aries). J Anim Sci. 2007;85:598–603.

    CAS  PubMed  Article  Google Scholar 

  3. 3.

    Chu MX, Yang J, Feng T, Cao GL, Fang L, Di R, Huang DW, Tang QQ, Ma YH, Li K, Li N. GDF9 as a candidate gene for prolificacy of small tail Han sheep. Mol Biol Rep. 2011;38:5199–204.

    CAS  PubMed  Article  Google Scholar 

  4. 4.

    Davis GH. Fecundity genes in sheep. Anim Reprod Sci. 2004;82-83:247–53.

    CAS  PubMed  Article  Google Scholar 

  5. 5.

    Mulsant P, Lecerf F, Fabre S, Schibler L, Monget P, Lanneluc I, Pisselet C, Riquet J, Monniaux D, Callebaut I, Cribiu E, Thimonier J, Teyssier J, Bodin L, Cognie Y, Chitour N, Elsen JM. Mutation in bone morphogenetic protein receptor-IB is associated with increased ovulation rate in Booroola merino ewes. Proc Natl Acad Sci U S A. 2001;98:5104–9.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  6. 6.

    Davis GH, Balakrishnan L, Ross IK, Wilson T, Galloway SM, Lumsden BM, Hanrahan JP, Mullen M, Mao XZ, Wang GL, Zhao ZS, Zeng YQ, Robinson JJ, Mavrogenis AP, Papachristoforou C, Peter C, Baumung R, Cardyn P, Boujenane I, Cockett NE, Eythorsdottir E, Arranz JJ, Notter DR. Investigation of the Booroola (FecB) and Inverdale (FecX(I)) mutations in 21 prolific breeds and strains of sheep sampled in 13 countries. Anim Reprod Sci. 2006;92(1–2):87–96.

    CAS  PubMed  Article  Google Scholar 

  7. 7.

    Di R, Chu MX, Li YL, Zhang L, Fang L, Feng T, Cao GL, Chen HQ, Li XW. Predictive potential of microsatellite markers on heterosis of fecundity in crossbred sheep. Mol Biol Rep. 2012;39(3):2761–6.

    CAS  PubMed  Article  Google Scholar 

  8. 8.

    Chu MX, Jia LH, Zhang YJ, Jin M, Chen HQ, Fang L, Di R, Cao GL, Feng T, Tang QQ, Ma YH, Li K. Polymorphisms of coding region of BMPR-IB gene and their relationship with litter size in sheep. Mol Biol Rep. 2011;38(6):4071–6.

    CAS  PubMed  Article  Google Scholar 

  9. 9.

    Wilusz JE, Sunwoo H, Spector DL. Long noncoding RNAs: functional surprises from the RNA world. Genes Dev. 2009;23(13):1494–504.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  10. 10.

    Derrien T, Johnson R, Bussotti G, Tanzer A, Djebali S, Tilgner H, Guernec G, Martin D, Merkel A, Knowles DG, Lagarde J, Veeravalli L, Ruan X, Ruan YJ, Lassmann T, Carninci P, Brown JB, Lipovich L, Gonzalez JM, Thomas M, Davis CA, Shiekhattar R, Gingeras TR, Hubbard TJ, Notredame C, Harrow J, Guigo R. The GENCODE v7 catalog of human long noncoding RNAs: analysis of their gene structure, evolution and expression. Genome Res. 2012;22(9):1775–89.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  11. 11.

    Peng WX, Koirala P, Mo YY. LncRNA-mediated regulation of cell signaling in cancer. Oncogene. 2017;36(41):5661–7.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  12. 12.

    Mattick JS. Long noncoding RNAs in cell and developmental biology. Semin Cell Dev Biol. 2011;22(4):327.

    PubMed  Article  Google Scholar 

  13. 13.

    ChenYG SAT, Chang HY. Gene regulation in the immune system by long noncoding RNAs. Nat Immunol. 2017;18(9):962–72.

    Article  CAS  Google Scholar 

  14. 14.

    Miao XY, Luo QM, Zhao HJ, Qin XY. An integrated analysis of miRNAs and methylated genes encoding mRNAs and lncRNAs in sheep breeds with different fecundity. Front Physiol. 2017;8:1049.

    PubMed  PubMed Central  Article  Google Scholar 

  15. 15.

    Feng X, Li FZ, Wang F, Zhang GM, Pang J, Ren CF, Zhang TT, Yang H, Wang ZY, Zhang YL. Genome-wide differential expression profiling of mRNAs and lncRNAs associated with prolificacy in Hu sheep. Biosci Rep. 2018;38(2):BSR20171350.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  16. 16.

    Yang H, Ma JY, Wang ZB, Yao XL, Zhao J, Zhao XY, Wang F, Zhang YL. Genome-wide analysis and function prediction of long noncoding RNAs in sheep pituitary gland associated with sexual maturation. Genes (Basel). 2020;11(3):320.

    CAS  Article  Google Scholar 

  17. 17.

    Su T, Yu HL, Luo G, Wang MX, Zhou CF, Zhang L, Hou B, Zhang C, Liu M, Xu DQ. The interaction of lncRNA XLOC-2222497, AKR1C1, and progesterone in porcine endometrium and pregnancy. Int J Mol Sci. 2020;21(9):3232.

    CAS  PubMed Central  Article  PubMed  Google Scholar 

  18. 18.

    Arun G, Akhade VS, Donakonda S, Rao MR. Mrhl RNA, a long noncoding RNA, negatively regulates wnt signaling through its protein partner ddx5/p68 in mouse spermatogonial cells. Mol Cell Biol. 2012;32(15):3140–52.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  19. 19.

    Gao WL, Liu M, Yang Y, Yang H, Liao Q, Bai Y, Li YX, Li D, Peng C, Wang YL. The imprinted h19 gene regulates human placental trophoblast cell proliferation via encoding mir-675 that targets nodal modulator 1 (nomo1). RNA Biol. 2012;9(7):1002–10.

    CAS  PubMed  Article  Google Scholar 

  20. 20.

    Li W, Notani D, Ma Q, Tanasa B, Nunez E, Chen AY, Merkurjev D, Zhang J, Ohgi K, Song X, Oh S, Kim HS, Glass CK, Rosenfeld MG. Functional roles of enhancer RNAs for oestrogen-dependent transcriptional activation. Nature. 2013;498(7455):516–20.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  21. 21.

    Xia Q, Li QL, Gan SQ, Guo XF, Zhang XS, Zhang JL, Chu MX. Exploring the roles of fecundity-related long non-coding RNAs and mRNAs in the adrenal glands of small-tailed Han sheep. BMC Genet. 2020;21(1):39–50.

    PubMed  PubMed Central  Article  Google Scholar 

  22. 22.

    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(9):773–84.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  23. 23.

    Dalum JV, Melum VJ, Wood SH, Hazlerigg DG. Maternal photoperiodic programming: melatonin and seasonal synchronization before birth. Front Endocrinol (Lausanne). 2020;10:901.

    Article  Google Scholar 

  24. 24.

    Ulitsky I, Shkumatava A, Jan CH, Sive H, Bartel DP. Conserved function of lincRNAs in vertebrate embryonic development despite rapid sequence evolution. Cell. 2011;147:1537–50.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  25. 25.

    Nam JW, Bartel DP. Long noncoding RNAs in C. elegans. Genome Res. 2012;22:2529–40.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  26. 26.

    Nishiwaki-Ohkawa T, Yoshimura T. Molecular basis for regulating seasonal reproduction in vertebrates. J Endocrinol. 2016;229:117–27.

    Article  CAS  Google Scholar 

  27. 27.

    Yoshimura T, Yasuo S, Watanabe M, Iigo M, Yamamura T, Hirunagi K, Ebihara S. Light-induced hormone conversion of T4 to T3 regulates photoperiodic response of gonads in birds. Nature. 2003;426:178–81.

    CAS  PubMed  Article  Google Scholar 

  28. 28.

    Ikegami K, Yoshimura T. The hypothalamic-pituitary-thyroid axis and biological rhythms: the discovery of tsh’s unexpected role using animal models. Best Pract Res Clin Endocrinol Metab. 2017;31:475–85.

    CAS  PubMed  Article  Google Scholar 

  29. 29.

    Jin JM, Yang WX. Molecular regulation of hypothalamus-pituitary-gonads axis in males. Gene. 2014;551:15–25.

    CAS  PubMed  Article  Google Scholar 

  30. 30.

    Choi YS, Lee HJ, Ku CR, Cho YH, Seo MR, Lee YJ, Lee EJ. FoxO1 is a negative regulator of FSHβ gene expression in basal and GnRH-stimulated conditions in female. Endocrinology. 2014;155(6):2277–86.

    PubMed  Article  CAS  Google Scholar 

  31. 31.

    Naor Z. Signaling by G-protein-coupled receptor (GPCR): studies on the GnRH receptor. Front Neuroendocrinol. 2009;30(1):10–29.

    CAS  PubMed  Article  Google Scholar 

  32. 32.

    Torii S, Yamamoto T, Tsuchiya Y, Nishida E. ERK MAP kinase in G1 cell cycle progression and cancer. Cancer Sci. 2006;97:697–702.

    CAS  PubMed  Article  Google Scholar 

  33. 33.

    Kang HM, Jeong CB, Lee YH, Cui YH, Kim DH, Lee MC, Kim HS, Han J, Hwang DS, Lee SJ, Lee JS. Cross-reactivities of mammalian MAPKs antibodies in rotifer and copepod: application in mechanistic studies in aquatic ecotoxicology. Mar Pollut Bull. 2017;124:614–23.

    CAS  PubMed  Article  Google Scholar 

  34. 34.

    Pasquale CD, Gentilin E, Falletta S, Bellio M, Buratto M, Uberti ED, Zatelli MC. PI3K/Akt/mtor pathway involvement in regulating growth hormone secretion in a rat pituitary adenoma cell line. Endocrine. 2018;60(2):308–16.

    PubMed  Article  CAS  Google Scholar 

  35. 35.

    Lodge EJ, Russell JP, Patist AL, Francis-West P, Andoniadou CL. Expression analysis of the hippo cascade indicates a role in pituitary stem cell development. Front Physiol. 2016;7:114.

    PubMed  PubMed Central  Article  Google Scholar 

  36. 36.

    Penque BA, Su L, Wang J, Ji W, Bale A, Luh F, Fulbright RK, Sarmast U, Sega AG, Konstantino M, Spencer-Manzon M, Pierce R, Yen Y, Lakhani SA. A homozygous variant in RRM2B is associated with severe metabolic acidosis and early neonatal death. Eur J Med Genet. 2019;62:103574.

    PubMed  Article  Google Scholar 

  37. 37.

    Morel F, Rauch C, Petit E, Piton A, Theret N, Coles B, Guillouzo A. Gene and protein characterization of the human glutathione S-transferase kappa and evidence for a peroxisomal localization. J Biol Chem. 2004;279(16):16246–53.

    CAS  PubMed  Article  Google Scholar 

  38. 38.

    Thomson RE, Bigley AL, Foster JR, Jowsey IR, Elcombe CR, Orton TC. Tissue-specific expression and subcellular distribution of murine glutathione S-transferase class kappa. J Histochem Cytochem. 2004;52:653–62.

    CAS  PubMed  Article  Google Scholar 

  39. 39.

    Jin X, Xu Z, Cao J, Shao P, Zhou M, Qin Z, Liu Y, Yu F, Zhou X, Ji W, Cai W, Ma Y, Wang C, Shan N, Yang N, Chen X, Li Y. Proteomics analysis of human placenta reveals glutathione metabolism dysfunction as the underlying pathogenesis for preeclampsia. Biochim Biophys Acta, Proteins Proteomics. 2017;1865(9):1207–14.

    CAS  PubMed  Article  Google Scholar 

  40. 40.

    Jiang XL, Pang YW, Zhao SJ, Hao HS, Zhao XM, Du WH, Wang YC, Zhu HB. Thioredoxin-interacting protein regulates glucose metabolism and improves the intracellular redox state in bovine oocytes during in vitro maturation. Am J Physiol Endocrinol. 2020;318(3):E405–16.

    Article  Google Scholar 

  41. 41.

    Rubin CI, Atweh GF. The role of stathmin in the regulation of the cell cycle. J Cell Biochem. 2004;93:242–50.

    CAS  PubMed  Article  Google Scholar 

  42. 42.

    Dou YD, Zhao H, Huang T, Zhao SG, Liu XM, Yu XC, Ma ZX, Zhang YC, Liu T, Gao X, Li L, Lu G, Chan WY, Gao F, Liu HB, Chen ZJ. STMN1 promotes progesterone production via StAR up-regulation in mouse granulosa cells. Sci Rep. 2016;6:26691–2670.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  43. 43.

    Suzuki S, Iwamoto M, Hashimoto M, Suzuki M, Nakai M, Fuchimoto D, Sembon S, Eguchi-Ogawa T, Uenishi H, Onishi A. Generation and characterization of RAG2 knockout pigs as animal model for severe combined immunodeficiency. Vet Immunol Immunopathol. 2016;178:37–49.

    CAS  PubMed  Article  Google Scholar 

  44. 44.

    Funke L, Dakoji S, Bredt DS. Membrane-associated guanylate kinases regulate 324 adhesion and plasticity at cell junctions. Annu Rev Biochem. 2005;74:219–45.

    CAS  PubMed  Article  Google Scholar 

  45. 45.

    Zhang ZB, Tang JS, Di R, Liu QY, Wang XY, Gan SQ, Zhang XS, Zhang JL, Cuhu MX, Hu WP. Integrated hypothalamic transcriptome profiling reveals the reproductive roles of mRNAs and miRNAs in sheep. Front Genet. 2019;10:1296.

    CAS  PubMed  Article  Google Scholar 

  46. 46.

    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.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  47. 47.

    Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, Liu Y, Chen R, Zhao Y. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013;41(17):e166.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  48. 48.

    Kang YJ, Yang DC, Kong L, Hou M, Meng YQ, Wei L, Gao G. CPC2: a fast and accurate coding potential calculator based on sequence intrinsic features. Nucleic Acids Res. 2017;45(W1):W12–W6.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  49. 49.

    Bateman A, Birney E, Cerruti L, Durbin R, Etwiller L, Eddy SR, Griffiths JS, Howe KL, Marshall M, Sonnhammer EL. The pfam protein families database. Nucleic Acids Res. 2002;30(1):276–80.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  50. 50.

    Liu X, Liu KQ, Shan BS, Wei SJ, Li DF, Han HY, Wei W, Chen J, Liu HL, Zhang LF. A genome-wide landscape of mRNAs, lncRNAs and circRNAs during subcutaneous adipogenesis in pigs. J Anim Sci Biotechnol. 2018;9:76.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  51. 51.

    Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550–70.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  52. 52.

    Alessandro F, Irene B. Long non-coding RNAs: new players in cell differentiation and development. Nat Rev Genet. 2014;15(1):7–21.

    Google Scholar 

  53. 53.

    Chen C, Tan H, Bi J, Li Z, Rong T, Lin Y, Sun L, Li X, Shen J. Identification of competing endogenous RNA regulatory networks in vitamin a deficiency-induced congenital scoliosis by transcriptome sequencing analysis. Cell Physiol Biochem. 2018;48:2134–46.

    CAS  PubMed  Article  Google Scholar 

  54. 54.

    Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–504.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

Download references

Acknowledgements

We are grateful to Ran Di, Qiuyue Liu and Caihong Wei for their suggestions on experimental design.

Funding

This work was supported by National Natural Science Foundation of China (31772580), the Earmarked Fund for China Agriculture Research System (CARS-38), Agricultural Science and Technology Innovation Program of China (ASTIP-IAS13), Central Public-interest Scientific Institution Basal Research Fund (Y2017JC24), China Agricultural Scientific Research Outstanding Talents and Their Innovative Teams Program, China High-level Talents Special Support Plan Scientific and Technological Innovation Leading Talents Program (W02020274), Tianjin Agricultural Science and Technology Achievements Transformation and Popularization Program (201704020). Data analysis, interpretation and manuscript preparation were funded by the National Natural Science Foundation of China (31772580) and the Earmarked Fund for China Agriculture Research System (CARS-38). The funding bodies had no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Author information

Affiliations

Authors

Contributions

This study was designed by MXC and CYL, who performed data analysis and prepared figures, tables. CYL wrote the manuscript. MXC, ZJZ and CHR contributed to revision of the manuscript. XYH contributed to field experiment. All authors read and approved the final manuscript for publication.

Corresponding author

Correspondence to Mingxing Chu.

Ethics declarations

Ethics approval and consent to participate

Experimental animals were authorized by the Science Research Department (in charge of animal welfare issues) of the Institute of Animal Science, Chinese Academy of Agricultural Sciences (IAS-CAAS). The study complies with current laws of the country in which the experiments were performed.

Consent for publication

Not applicable.

Competing interests

All authors declare no conflicts of interest.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Figure S1.

Distribution of DE lncRNAs on chromosomes in MM_FP vs MM_LP.

Additional file 2: Figure S2.

Distribution of DE lncRNAs on chromosomes in MM_FP vs ww_FP.

Additional file 3: Figure S3.

Distribution of DE lncRNAs on chromosomes in MM_LP vs ww_LP.

Additional file 4: Figure S4.

Distribution of DE lncRNAs on chromosomes in ww_FP vs ww_LP.

Additional file 5: Figure S5.

Distribution of DE mRNAs on chromosomes in MM_FP vs MM_LP.

Additional file 6: Figure S6.

Distribution of DE mRNAs on chromosomes in MM_FP vs ww_FP.

Additional file 7: Figure S7.

Distribution of DE mRNAs on chromosomes in MM_LP vs ww_LP.

Additional file 8: Figure S8.

Distribution of DE mRNAs on chromosomes in ww_FP vs ww_LP.

Additional file 9: Figure S9.

Density distribution of candidate transcripts.

Additional file 10: Table S1.

Overview of DE mRNAs closely related to reproductive signal pathways.

Additional file 11: Supplementary material 1A.

Total set of known lncRNAs were identified from four groups. Supplementary material 1B. Total set of novel lncRNAs were identified from four groups.

Additional file 12: Supplementary material 2.

Total set of mRNAs were identified from four groups.

Additional file 13: Supplementary material 3A.

Total set of lncRNAs were up- and down- regulated in MM_FP vs MM_LP. Supplementary material 3B. Total set of lncRNAs were up- and down- regulated in MM_FP vs ww_FP. Supplementary material 3C. Total set of lncRNAs were up- and down- regulated in MM_LP vs ww_LP. Supplementary material 3D. Total set of lncRNAs were up- and down- regulated in ww_FP vs ww_LP.

Additional file 14: Supplementary material 4A.

Total set of mRNAs were up- and down- regulated in MM_FP vs MM_LP. Supplementary material 4B. Total set of mRNAs were up- and down- regulated in MM_FP vs ww_FP. Supplementary material 4C. Total set of mRNAs were up- and down- regulated in MM_LP vs ww_LP. Supplementary material 4D. Total set of mRNAs were up- and down- regulated in ww_FP vs ww_LP.

Additional file 15: Supplementary material 5A.

GO enrichment of differentially expressed lncRNA targets in MM_FP vs MM_LP. Supplementary material 5B. GO enrichment of differentially expressed lncRNA targets in MM_FP vs ww_FP. Supplementary material 5C. GO enrichment of differentially expressed lncRNA targets in MM_LP vs ww_LP. Supplementary material 5D. GO enrichment of differentially expressed lncRNA targets in ww_FP vs ww_LP.

Additional file 16: Supplementary material 6A.

GO enrichment of differentially expressed mRNAs in MM_FP vs MM_LP. Supplementary material 6B. GO enrichment of differentially expressed mRNAs in MM_FP vs ww_FP. Supplementary material 6C. GO enrichment of differentially expressed mRNAs in MM_LP vs ww_LP. Supplementary material 6D. GO enrichment of differentially expressed mRNAs in ww_FP vs ww_LP.

Additional file 17: Supplementary material 7A.

Total set of the top 20 KEGG enrichment pathways for differentially expressed lncRNA targets in MM_FP vs MM_LP. Supplementary material 7B. Total set of the top 20 KEGG enrichment pathways for differentially expressed lncRNA targets in MM_FP vs ww_FP. Supplementary material 7C. Total set of the top 20 KEGG enrichment pathways for differentially expressed lncRNA targets in MM_LP vs ww_LP. Supplementary material 7D. Total set of the top 20 KEGG enrichment pathways for differentially expressed lncRNA targets in ww_FP vs ww_LP.

Additional file 18: Supplementary material 8A.

Total set of the top 20 KEGG enrichment pathways for differentially expressed mRNAs in MM_FP vs MM_LP. Supplementary material 8B. Total set of the top 20 KEGG enrichment pathways for differentially expressed mRNAs in MM_FP vs ww_FP. Supplementary material 8C. Total set of the top 20 KEGG enrichment pathways for differentially expressed mRNAs in MM_LP vs ww_LP. Supplementary material 8D. Total set of the top 20 KEGG enrichment pathways for differentially expressed mRNAs in ww_FP vs ww_LP.

Additional file 19: Supplementary material 9A.

Co-expression details of DE lncRNA-mRNA after lncRNA targets coincided with DE mRNAs in MM_FP vs MM_LP. Supplementary material 9B. Co-expression details of DE lncRNA-mRNA after lncRNA targets coincided with DE mRNAs in MM_FP vs ww_FP. Supplementary material 9C. Co-expression details of DE lncRNA-mRNA after lncRNA targets coincided with DE mRNAs in MM_LP vs ww_LP. Supplementary material 9D. Co-expression details of DE lncRNA-mRNA after lncRNA targets coincided with DE mRNAs in ww_FP vs ww_LP.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Li, C., He, X., Zhang, Z. et al. Pineal gland transcriptomic profiling reveals the differential regulation of lncRNA and mRNA related to prolificacy in STH sheep with two FecB genotypes. BMC Genom Data 22, 9 (2021). https://doi.org/10.1186/s12863-020-00957-w

Download citation

Keywords

  • LncRNAs
  • RNA-Seq
  • Pineal gland
  • Prolificacy
  • Sheep
\