Pineal gland transcriptomic profiling reveals the differential regulation of lncRNA and mRNA related to prolificacy in STH sheep with two FecB genotypes
BMC Genomic Data volume 22, Article number: 9 (2021)
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.
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.
All of these differential lncRNAs and mRNAs expression profiles in pineal gland provide a novel resource for elucidating regulatory mechanism underlying STH sheep prolificacy.
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 . 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  and GDF9  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 . To date, this mutation has been detected in diverse sheep species such as Booroola Merino sheep (Australia) , Small Tail Han (STH) and Hu sheep (China) . 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 . 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 . 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 . LncRNA is proposed to be the largest transcript class in mammalian transcriptome , 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 , signal transduction , immune regulation . 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 . Then, Feng et al. (2018) identified 5 lncRNAs and 76 mRNAs in ovaries of Hu sheep with high and low prolificacy, respectively . 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 . 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 . 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.
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).
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.
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).
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).
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).
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).
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).
Studies have found that lncRNA is involved in multiple reproductive functions such as spermatogenesis , placentation , signaling pathway of sex hormone response [20, 21] and gonadgenesis . It is known that the melatonin synthesized in pineal gland is closely related to the estrus cycle . 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 . 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 . 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 . 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 . 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 . 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 , 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 . Over-expression of STMN1 stimulates progesterone production by modulating the promoter activity of Star and Cyp11a1 in mouse granulosa cells . Besides, RAG2 is indispensable for generation of antigen receptor diversity in immune cells . 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 . 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.
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.
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 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 . 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) . 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 . The mapped reads of each sample were assembled by StringTie v. 1.3.1 .
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 , CPC v. 2.81  and PFAM v.1.3  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) . For experiments with three biological replicates, the R package DEseq2 was used to identify differentially expressed transcripts after a negative binomial distribution . 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.
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 . 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 .
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) .
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).
Small tailed Han sheep
Long noncoding RNA
The fragments per kilobase of transcript per million reads mapped
Controlled internal drug releasing device
Bartlewski PM, Baby TE, Giffin JL. Reproductive cycles in sheep. Anim Reprod Sci. 2011;124:259–68.
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.
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.
Davis GH. Fecundity genes in sheep. Anim Reprod Sci. 2004;82-83:247–53.
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.
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.
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.
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.
Wilusz JE, Sunwoo H, Spector DL. Long noncoding RNAs: functional surprises from the RNA world. Genes Dev. 2009;23(13):1494–504.
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.
Peng WX, Koirala P, Mo YY. LncRNA-mediated regulation of cell signaling in cancer. Oncogene. 2017;36(41):5661–7.
Mattick JS. Long noncoding RNAs in cell and developmental biology. Semin Cell Dev Biol. 2011;22(4):327.
ChenYG SAT, Chang HY. Gene regulation in the immune system by long noncoding RNAs. Nat Immunol. 2017;18(9):962–72.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Dalum JV, Melum VJ, Wood SH, Hazlerigg DG. Maternal photoperiodic programming: melatonin and seasonal synchronization before birth. Front Endocrinol (Lausanne). 2020;10:901.
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.
Nam JW, Bartel DP. Long noncoding RNAs in C. elegans. Genome Res. 2012;22:2529–40.
Nishiwaki-Ohkawa T, Yoshimura T. Molecular basis for regulating seasonal reproduction in vertebrates. J Endocrinol. 2016;229:117–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.
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.
Jin JM, Yang WX. Molecular regulation of hypothalamus-pituitary-gonads axis in males. Gene. 2014;551:15–25.
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.
Naor Z. Signaling by G-protein-coupled receptor (GPCR): studies on the GnRH receptor. Front Neuroendocrinol. 2009;30(1):10–29.
Torii S, Yamamoto T, Tsuchiya Y, Nishida E. ERK MAP kinase in G1 cell cycle progression and cancer. Cancer Sci. 2006;97:697–702.
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.
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.
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.
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.
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.
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.
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.
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.
Rubin CI, Atweh GF. The role of stathmin in the regulation of the cell cycle. J Cell Biochem. 2004;93:242–50.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Alessandro F, Irene B. Long non-coding RNAs: new players in cell differentiation and development. Nat Rev Genet. 2014;15(1):7–21.
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.
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.
We are grateful to Ran Di, Qiuyue Liu and Caihong Wei for their suggestions on experimental design.
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.
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
All authors declare no conflicts of interest.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Distribution of DE lncRNAs on chromosomes in MM_FP vs MM_LP.
Distribution of DE lncRNAs on chromosomes in MM_FP vs ww_FP.
Distribution of DE lncRNAs on chromosomes in MM_LP vs ww_LP.
Distribution of DE lncRNAs on chromosomes in ww_FP vs ww_LP.
Distribution of DE mRNAs on chromosomes in MM_FP vs MM_LP.
Distribution of DE mRNAs on chromosomes in MM_FP vs ww_FP.
Distribution of DE mRNAs on chromosomes in MM_LP vs ww_LP.
Distribution of DE mRNAs on chromosomes in ww_FP vs ww_LP.
Density distribution of candidate transcripts.
Overview of DE mRNAs closely related to reproductive signal pathways.
Total set of known lncRNAs were identified from four groups. Supplementary material 1B. Total set of novel lncRNAs were identified from four groups.
Total set of mRNAs were identified from four groups.
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.
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.
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.
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.
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.
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.
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.
About this article
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