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

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. Supplementary Information The online version contains supplementary material available at 10.1186/s12863-020-00957-w.

Page 3/32 STH sheep are 2.61, 286.5%, respectively [7]. There are three genotypes based on effects of FecB mutation in STH sheep, namely FecB BB (with two-copy FecB mutation), FecB B+ (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 speci city 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 signi cantly enriched in the oxytocin signaling pathway [14]. Then, Feng et al. (2018) identi ed 5 lncRNAs and 76 mRNAs in ovaries of Hu sheep with high and low proli cacy, 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 FecB BB (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% to 93% of the reads were successfully mapped to the Ovis aries reference genome (Table 1).
GO analysis of the biological function of DE lncRNAs and mRNAs GO annotation enrichment was used to describe functions of the DE lncRNAs and mRNAs involved in cellular components, molecular function and biological processes, as shown in Fig. 3. Between MM_F_P and MM_L_P, targeted genes for DE lncRNAs were most enriched, and the terms were related to regulation of trans-membrane transport, antigen processing and presentation, immune system process. DE mRNAs were most enriched, the meaningful terms were related to the regulation of C-terminal protein methylation, C-terminal protein amino acid modi cation, post-translation protein modi cation, cellular macromolecular complex assembly and cellular metabolic process (Fig. 3A, Supplementary material 5A, 6A).
Between MM_F_P and ww_F_P, targeted genes for DE lncRNAs were enriched, the terms were related to regulation of protein complex assembly and biogenesis, protein complex subunit organization, spindle assembly involved in mitosis process. DE mRNAs were most enriched, the meaningful terms were related to regulation of secondary metabolic and biosynthetic process, viral protein processing, single-organism biosynthetic process (Fig. 3B, Supplementary material 5B, 6B).
Between MM_L_P and ww_L_P, targeted genes for DE lncRNAs were enriched, the terms were mainly related to regulation of single organism signaling, signal transduction, cellular response to stimulus and cellular communication. DE mRNAs were enriched, the meaningful terms were related to regulation of RNA methylation, metabolic process, organic substance metabolic process (Fig. 3C, Supplementary material 5C, 6C).
Between ww_F_P and ww_L_P, targeted genes for DE lncRNAs were enriched, the terms were related to regulation of immune response, glycerolipid metabolic process, cellular response to abiotic stimulus. DE mRNAs were enriched, the terms were related to regulation of nucleosome and chromatin assembly, nucleosome organization process (Fig. 3D, Supplementary material 5D, 6D). 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 coexpression of DE lncRNAs and DE target mRNAs, after screening the overlaps between target mRNAs and DE mRNAs in each comparison group, which indicated regulation of lncRNA and mRNA in reproduction (|Pearson correlation| 0.95). Between MM_F_P and MM_L_P, a total of 5 DE lncRNAs and 9 DE mRNAs lncRNAs and 10 DE mRNAs were involved in the network, and it consists of 10 edges (Fig. 8C, Supplementary material 9C). Between ww_F_P and ww_L_P, a total of 30 DE lncRNAs and 12 DE mRNAs were involved in the network, and it consists of 47 edges (Fig. 9, Supplementary material 9D).

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 pro les 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 rst genome-wide analysis of pineal gland in sheep with different genotypes, and might provide valuable resource for searching functional lncRNAs associated with sheep proli cacy.
In present study, we screened 21,282 lncRNAs and 43,674 coding transcripts. 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 in NC_019459.2, NC_019460.2, NC_019458.2 were greater than those in 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 1,360 (764+596) DE mRNAs in pineal gland at follicular and luteal phases between high and low proli cacy 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 speci c β 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 FecB BB 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, 105604037 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 detoxi cation 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 FecB BB 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 identi ed 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 nal lambing.

Conclusion
In summary, the pineal gland transcriptomic study reveals differential regulation of lncRNAs and mRNAs related to proli cacy 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 pro les provide a novel resource for elucidating regulatory mechanism underlying STH sheep proli cacy.

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 signi cant difference in age, weight, height, body length, chest circumference and tube circumference were selected for this experiment. 12 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 slaughtered on the 50 th hour after CIDR removal and the pineal tissues were collected (follicular phase, MM_F_P and ww_F_P, respectively). The other 6 sheep were slaughtered on the 7 th day after CIDR removal and the pineal tissues were collected (luteal phase, MM_L_P and ww_L_P, respectively). Obtained samples were stored immediately at -80℃ 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 rRNAdepleted 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 rstly 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].
Identi cation of potential lncRNA candidates Potential lncRNA candidates were identi ed as following: selecting transcripts length >200 nt with exon number ≥2. Transcripts that overlap with the database annotation exon region were removed. After calculating the expression level of each transcript using Cuffquant, transcripts of FPKM ≥0.5 were selected. The intersection of transcripts with no coding potential in software analysis results was used as the lncRNA data set. These transcripts were predicted by CNCI v.2.0 [47] . LncRNAs and mRNAs with P-value <0.05 and a fold change >2.0 were considered as differentially expression between two groups of data.

Bioinformatics analysis
The function of DE lncRNAs was predicted by GO and KEGG analysis of their cis-and trans-target mRNAs, which distance less than 100 Kb in genomic position were assumed to be cis-target genes, and Pearson correlation coe cient with the lncRNA-RNA pairs being ≥0.95 for trans-target mRNAs [52].
Statistical enrichment of DE lncRNA targets and DE mRNAs in GO annotation and KEGG pathway were evaluated, P-value ≤0.05 de ned as the signi cant threshold, signi cance 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 signi cant.    i-Note: "i-" represents the identical information with previous one in the same column.       Co-expression of DE lncRNA-mRNA after lncRNA targets coincided with DE mRNAs in ww_F_P vs ww_L_P.

Additional Files
Tangerine and green represent upregulated and downregulated, respectively. Octagons and triangles represent lncRNAs and mRNAs, respectively.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.