Transcriptome Analysis Reveals that Long Noncoding RNAs Contribute to Development Differences of Medium-sized Ovarian Follicle between Meishan and Duroc Sows

Background: Ovulation rate is an extremely important factor of litter size in sows. It differs greatly among pig breeds of different genetic backgrounds. Long non-coding RNAs (lncRNAs) can regulate follicle development, granulosa cell(GC) growth and hormone secretion, which in turn affects sow litter size. Results: In our research, we identied 3554 lncRNAs and 25491 mRNAs in M2 follicle from Meishan and Duroc pigs. lncRNAs sequence and open reading frame(ORF) length is shorter than mRNAs, and it has fewer exons, lower abundance and conserved than protein-coding RNAs. Furthermore, 201 lncRNAs were differentially expressed in breeds, and quantitative trait loci (QTL) analysis of differential expression (DE) lncRNAs were performed. 127 DE lncRNAs are located in 119 reproduction trait-related loci. In addition, the lncRNAs potential target genes (PTGs) in cis or trans were predicted. Gene ontology(GO) and KEGG pathway analysis revealed that some PTGs were include some follicular development and hormone secretion-related biological processes or pathways, such as regulation of progesterone biosynthetic process, oestrogen metabolic process and ovarian steroidogenesis and PI3K-Akt signalling pathway. Furthermore (cid:0) We also screened 19 differentially expressed lncRNAs of PI3K-Akt signalling pathway as candidates. Conclusions: This study provided a new signicance on the roles of lncRNAs in follicular growth and development and porcine reproduction. We analysed the RNA data of six libraries of M2 follicles from three Meishan and three Duroc pigs, and obtained 111,275,282 to 142,313,048 raw reads and 101,760,644 to 138,324,398 clean reads from per sample, respectively. The clean reads were used in the after discarding transcripts with adapters, the low quality reads or other possible contaminants. Tophat2 software was used to perform reference genomic alignment analysis on Clean reads. The results showed that the six sequencing libraries’ mapping ratio were all greater than 70%, indicating that the sequencing accuracy is high and can be used for subsequent analysis (Table 1).

Identi cation and characterisation of lncRNA All transcripts were screened in ve steps. Finally 3,554 lncRNAs were obtained, including 1,997 known lncRNAs and 1,557 novel lncRNAs ( Figures 1A and B). Among the lncRNAs, lncRNA and antisense-lncRNA accounted for 89.3% and 10.7%, respectively ( Figure 1C). In addition, we obtained 25,491 mRNAs. The lengths of lncRNA sequences were shorter than those of mRNAs. The lncRNA transcripts average length was 1417 bp, while the mRNAs was 2423 bp ( Figure 1D). lncRNAs had fewer exons. The exon number of lncRNAs ranged from 1 to 10, whereas that of mRNAs from 1 to 30 ( Figure 1E). The lncRNA have shorter ORF than mRNA. The ORF lengths of lncRNAs mostly ranged from 1 bp to 500 bp, whereas those of the mRNAs were mostly from 1 bp to 2000 bp ( Figure 1F). The expression abundance of lncRNA were lower than mRNA ( Figure 1G). The lncRNA were less conserved compared with the protein coding genes, as revealed through a phastCon analysis ( Figure 1H). The structural feature comparative analysis of selected lncRNAs con rmed the accuracy of selection.

DE Analysis of lncRNAs and mRNAs
We use cu inks to normalize transcript expression to FPKM value, and performed a differential expression transcript analysis between the To verify the accuracy of sequencing, 6 DE-lncRNAs (ENSSSCT00000018610, ALDBSSCT0000001721, ALDBSSCT0000000051, LNC_000116, ALDBSSCT0000011300 and ALDBSSCT0000006152) and 6 DE-mRNAs (COL3A1, LRP8, ENSSSCT00000009222, SEPP1, COL5A2 and CYP19A1) were selected randomly, and their relative expression levels in the DFM2DAY4 and MFM2DAY4 groups were detected using qPCR. The expression analysis of the six lncRNA and six mRNA are displayed in the Figure 3, which are consistent with the RNA-Seq analysis results both in lncRNA and mRNA ( Figure 3).

QTL Analysis of DE lncRNAs
To explore the preliminary DELs function, we mapped the differential expressions(DE) of lncRNAs transcri[pts to the QTL database, and performed QTL analysis. The result showed that a total of 1446 QTLs that 145 DELs were located in, the 127 DE lncRNAs are located in 119 reproduction trait-related loci ( Figure 4A). By studying the distribution of QTLs on the chromosome, 119 QTLs related to reproduction deposition were found to be distributed in chromosomes 1, 2, 3, 4, 7, 8, 13, 15 and 18 ( Figure 4B). The 127 lncRNAs associated with reproduction QTLs could affect corpus luteum number (53), litter size (18), androstenone (23), total number born alive (4), number of stillborn (3), FSH concentration (3) and number of viable embryos (2) ( Figure 4C).

Prediction of lncRNA potential target genes
In order to further explore the regulatory functions of DE lncRNAs, we predicted their PTGs of lncRNAs in 2 ways-cis and trans [13,14]. In our study, we found that lncRNA may regulate multiple coding genes. For PTGs regulated by lncRNAs in cis, we searched protein-coding transcripts located in 100 kb upstream or downstream of the DELs as its cis-regulatory target genes. We obtained a total of 320 co-localised target genes of 118 DELs. In trans way analysis, we obtained 2930 PTGs of 127 DELs. We only showed the PTGs s of 24 DELs. The PTGs numbers for each DEL were varied. For example, the maximum number among the lncRNA is LNC_000179, and it had 77 PTGs, the second is LNC_000715 with 69 target genes. Some lncRNAs, such as LNC_000718 and LNC_000802, only had 2 target genes ( Table 2). Function Enrichment Analysis for lncRNAs To study the regulatory role of DE lncRNAs in Meishan and Duroc M2 follicles, we predicted the function of DE target genes by using GO and KEGG to speculate the lncRNA functions. The GO and KEGG analysis results revealed that PTGs participated in 1063 biological processes and 111 pathways, signi cantly. Many biological processes were involved in follicular development and ovulation, such as the regulation of progesterone biosynthetic, oestrogen metabolism, negative regulation of cell proliferation, ITGA3-ITGB1-THBS1 complex, cellular response to transforming growth factor beta stimulus, meiotic cell cycle and steroid catabolic process (p < 0.05) ( Figure 5A;). The KEGG pathway analysis shown the PTGs were signi cantly involved in ovarian steroidogenesis, PI3K-Akt, MAPK, Wnt, BMP and TNF signalling pathway (p < 0.05) ( Figure 5B).
Some PTGs that participated in oestrogen metabolic process and ovarian steroidogenesis signalling pathway are highlighted, such as CYP1A1, CYP19A1 and HSD3B1. A gene that participated in oestrogen metabolic process and ovarian steroidogenesis signalling pathway was CYP3A7. HSD17B8 participated in oestrogen metabolic process. ALOX5, LHCGR, IGF1, GNAS, CYP2J2 and CYP17A1 participated in ovarian steroidogenesis signalling pathway. In addition, one protein-coding transcripts may be regulated by multiple lncRNAs, such as DE-lncRNA ALDBSSCT0000001929, ALDBSSCT0000006256 and ALDBSSCT0000002430, which were correlated with their target gene CYP1A1 signi cantly (p<0.05). They were all downregulated in Meishan compared with Duroc sows. LNC_000644 and ALDBSSCT0000006057 were correlated with CYP3A7 signi cantly (p<0.05) and were unregulated in Meishan. ALDBSSCT0000000051 was correlated with its target gene HSD17B8 signi cantly (p<0.05). ALDBSSCT0000000051 was downregulated, whereas HSD17B8 was unregulated ( Figure 5C). Therefore, we speculated that some DE-lncRNAs participate in the development of porcine follicles by positively or negatively regulating their target genes, which are related to hormone secretion and metabolism.

Discussion
The follicle is a crucial tissue in mammal reproduction, and its development determines ovulation rate. Follicle groth and development is a complex and precise process, and long no coding RNA play an important role in thisfollicle development process. Some research evidence indicated the important roles of lncRNAs in animal reproduction with the development of gene chip and sequencing technologies. A. D. Macaulay used confocal transmission electron microscopy and RNA-Seq to determine that the cumulus cells around bovine oocytes transported a large number of nutrients and substances, e.g., mRNA and lncRNA on the oocytes of adult cows [15]. Sun detected 92 differentially expressed transcripts between Erhualian and large white pigs ovulatory follicles by using microarray [16]. Liu identi ed 2076 lncRNAs and 25491 mRNAs in Duroc's ovaries on days 0, 2, and 4 of follicular growth and development [17]. Meishan and Drouc sows show a huge difference in ovulation rate and ovarian follicle development and are good animal models for uncovering the molecular mechanisms underlying these differences. Therefore, we used RNA-Seq to identify the lncRNAs expressed in the M2 follicles of Meishan and Drouc sows. In our research, the newly identi ed lncRNAs in pig ovarian follicle have many obvious characteristics, it has shorter and fewer exons, longer exon length and lower abundance, and they are less conserved than mRNAs. The characteristics of these lncRNAs were compliance with those observed in other studies [18,19]. This current research is a meaningful resource for further studies on the functional lncRNAs and mRNA in pig ovarian follicles.
We obtained 145 differentially expressed lncRNAs QTL analysis, and the result showed that 127 DE lncRNAs were located in known QTLs of corpus luteum number, litter size, androstenone, total number born alive, number of stillborn, FSH concentration and number of viable embryos. The previous studies showed that QTLs controlling ovulation rate are located in pig chromosome 8 [20,21]. QTLs of luteum and ovulation rate are in chromosome 3 [22]. The QTLs of ovulation rate are in chromosomes 4, 7, 13 and 15. These QTLs are closely related to pig reproduction [23,24]. In this study, we found that the QTLs enriched most of the differential lncRNAs, i.e., ovarian follicle lncRNAs, distributed in chromosomes 3, 7, 8, 13 and 15. The results indicated that the DE lncRNAs may be participate in the regulation of porcine reproductive traits and may be related to the follicular development in the two pig groups.
LncRNAs have a low expression abundance, and some lncRNAs are still unclearly. Thus, understanding the function of lncRNA is di cult. We obtained the target genes of DELs, and functional analysis found that these target genes participate in multiple biological processes and molecular signalling. CYP19A1 is the target gene of novel lncRNA-LNC_000757, it's expressed in granulosa cells(GCs) and is rate-limiting enzyme for oestrogen production. It can catalyse the conversion of androgens to oestrogens, thereby increasing the expression level of oestrogen, knockdown of CYP19A1 can result in the inability to ovulate and loss of corpus luteum in female mice [25]. CYP1A1, an important metabolic enzyme in the oestrogen metabolic process and ovarian steroidogenesis signalling pathway, can reduce the level of oestrogen and lead to the acceleration of follicular atresia [26,27]. DE-lncRNA ALDBSSCT0000001929, ALDBSSCT0000006256 and ALDBSSCT0000002430 targeted and were positively correlated (p<0.05) with CYP1A1. The expression level of CYP1A1 was lower in Meishan, whereas that of CYP19A1 was higher in Meishan compared with Duroc. These results may indicate that the level of oestrogen in the M2 follicles of Meishan sows was higher than that in Duroc. In addition, IGF1 and HSD3B1 can stimulate and promote granulosa cells proliferation and promote follicular development [28,29]. We also found HSD17B8, ALOX5, LHCGR, GNAS, CYP2J2 and CYP17A1 participate in ovarian steroidogenesis signalling pathway. These protein-coding genes may be candidate genes for porcine follicle development and are targeted by one or more lncRNAs ( Figure 5C). Therefore, we speculated that lncRNAs may participate in porcine follicular development by regulating their target genes.
Studies have con rmed that the PI3K-AKT signalling pathway plays an major roles in the ovarian function regulation and follicular development [30][31][32]. PI3K-AKT signalling pathway regulates the primordial follicle maintenance and activation of s and promotes apoptosis of ovarian granulosa cells in humans and mice [33,34]. KEGG analysis revealed that some target mRNAs of the DE lncRNAs were related to the PI3K-AKT signalling pathways. The THBS1-ITGA-ITGB complex members belong to the TGF-β family. TGF-β signalling can regulate the PI3K-Akt signalling pathway. The lncRNA-ALDBSSCT0000004 244 targets THBS1, which is an ECM. The expression of THBS1 is driven by oestradiol (E2), and its abnormal expression causes the apoptosis of granulosa cells and the acceleration of follicular atresia [35]. ITGA3 is the target gene of novel lncRNA LNC_000857 and LNC_001052. The expression of ITGA3 in M2 follicles of Meishan sows was higher than in those of Duroc. However, at present, no reports of ITGA3 related to animal reproduction exist. ALDBSSCT0000010443 and ALDBSSCT0000009356 target ITGA6, which may be participate in the regulation of cumulus expansion and oocyte maturation [36]. ITGB1, ITGB4 and ITGB7 are the target genes of LNC_001556, LNC_000536 and ALDBSSCT0000004232, respectively. ITGB1 is associated with the apoptosis process [37] and is regulated by progesterone and oestrogen [38]. The expressions of PI3K and AKT is the golden standard for the activation of PI3K-AKT signalling pathway. We found that ENSSSCT000 00036444 and ALDBSSCT0000000805 target PIK3CA and PIK3C2B, and LNC_000751, LNC_000819 and LNC_000058 target AKT2. PIK3CA and PIK3C2B are the catalytic subunits of PI3K, and the knockout of PIK3CA and PIK3C2B causes early embryo death in mice [39]. PI3K can react with a variety of growth factors, thereby phosphorylating PIP2 to PIP3, activating PDK1, and then indirectly or directly activating Akt [40]. AKT2 is a subtype of AKT, and a downstream signalling core molecule of the PI3K-AKT classic signalling pathway [41,42]. Activated AKT2 can phosphorylate its downstream signalling molecules and produces the cAMP that can activate CREB, TP53 and IKK protein for nuclear gene transcription and expression, thereby regulating cell proliferation [43,44]. ALDBSSCT0000011847 and ALDBSSCT0000008900 target CREB3L3. ALDBSSCT0000004325 and LNC_001172 target the IKBKG, TP53 and ALDBSSCT0000002268 and were both signi cantly positively correlated in pairs (p<0.05), LNC_000076 and ALDBSSCT0000002268 target TP53. CREB3L3 is a transcription factor of CREB, and its function is similar to that of CREB. Its abnormal expression can attenuate the up-regulation of Egr1 by GnRH receptor activation, which in turn affects the expression level of LH-β and the growth of granular cells [45]. IKK can induce phosphorylation of IkB (inhibitory protein of NF-kB), dissociate the NF-kB/IkB dimer, activate NF-Kb, initiate the inhibition of apoptosis pathway and maintain porcine follicular development [46]. TP53 is a critical factor for cell survival, and suppression of p53 in oocytes can promote follicular growth and development [47]. In this study, we screened 19 DE-lncRNA candidates and their 12 target gene in the PI3K-AKT signalling pathway. We suggest that lncRNA may participate in follicular growth and development by regulating their target genes, but the speci c functions still require systematic research.

Conclusions
We identi ed lncRNA and mRNA expressed in M2 follicles in Meishan and Duroc sows and found that some lncRNAs participate in the follicular development by regulating their target genes. These ndings in transcriptome provide a valuable resource for follicular development and reproduction-related transcripts. The interactions between DE-lncRNAs and their PTGs and enriched pathways provide clues for further research on the role of follicular growth and development in pig. However, the lncRNAs function and molecular mechanism in follicular development remain unclear. Nevertheless, our study provides new insights into lncRNAs associated with follicular development and reproduction in pig.

Animals and sample collection
Three multiparous Meishan cyclic sows were maintained at the Agricultural Experimental Station of Yangzhou University. Three multiparous Duroc cyclic sows were raised at the Animal Experimental Station of the Shihezi University. Observe the sows every day and determine natural oestrous cycle, day 0 was the rst oestrus day). Sows were injected with PGF2α according to the pigs' weight analogue (cloprostenol, Ningbo Second Hormone Factory) on the 14th day of the oestrous cycle to induce luteal regression and to synchronise the follicular growth phase. 4 days later, the sows were electric stunning and quick bleeding, then ovaries and medium follicles (follicles diameter: 5.0-6.9 mm, M2) were harvested , snap-frozen in LN2 and then stored at -80 ℃ refrigerator.

RNA-Seq preparation and sequencing analysis
Total RNAs of all the follicle samples were extracted using TRIzol reagent (Invitrogen). RNA degradation and contamination was monitored on 1% agarose gels. The total RNA concentration, integrity and purity were detected using qubit RNA assay kit in Qubit 2.0 urometer (Life Technologies), RNA Nano 6000 assay kit by a bioanalyser 2100 (Agilent Technologies) and NanoPhotometer spectrophotometer (IMPLEN).
The ribosomal RNAs of all samples were removed using the Ribo-zero rRNA Removal Kit (Epicentre). 6 strand-speci c RNA-Seq pools for the M2 follicles of the six sows were constructed according to the manufacturer's instructions. Purify the samples library fragments using the AMPure XP system (Beckman Coulter, USA) to remove the preferred cDNA fragments of 150-200 bp. The blunt end cDNA fragments were augmented with A base and ligated to the sequencing adapter. The nal products were puri ed (AMPure XP System), and the quality of library was assessed by the system Agilent Bioanalyser 2100. The samples pools were then analysed using one lane of 100-200nt pairedend HiSeq 4000 platform. Quality control (QC) of RNA-Seq reads was performed using Fast Q C.

Transcriptome assembly
We ltering reads with adapter and low-quality, poly-N reads from raw data through inhouse perl scripts, the clean reads were obtained. Clean data with high quality is the basis for all subsequent analyses. Reference genome les were downloaded from the Ensembl (Sus scrofa 10.2). Pigs reference genome annotation index was built using Bowtie v2.0.6 [48]. Appropriate parameters were set using Tophat2 v 2.0.9 [49]. The Scripture and Cu inks [50] were used to assemble and splice the aligned sequences, which can be as small as possible. The transcript set, Cu inks has speci c parameters for the chain-speci c library, and the directional information of the transcript chain were accurately provided.
lncRNA Identi cation The following steps to identify lncRNAs from the nonredundant transcriptome. (1) Transcripts has single exon or two exons were ltered out.
(2) Transcripts length with less than 200 bp were removed. (3) Any transcript with the FPKM (a fragment per kilobase of transcript per million mapped read) score lower than 0.5 in every pools was discarded. (4) The remaining transcripts were blasted in pig known annotation lncRNA database -ALDB database [51] using Cuffcompare. Only the transcript of lncRNAs whose splice sites were congruent between our results and those in ALDB were immediately brought in as a known lncRNA. (5) Transcripts of any known protein-coding were discarded, and transcripts that belonged to pre-miRNA, snRNA, rRNA, snoRNA and pseudogenes were removed. (6) The CNCI, CPC, and phyloCSF tool were used to calculate the transcripts that hascoding potential. Transcript with a CNCI score of <0 [52], CPC score of <0 [53], Pfam-scan E-value of <0.001 [54] and phyloCSF Max-score of ≤100 [55], as well as the intersection result of each tools were de ned as novel found lncRNA transcript.
Differentially Expressed lncRNAs and mRNA Analysis. qRT-PCR Veri cation Total RNAs were extracted using TRIzol (Invitrogen, CA, USA) and cDNA was synthesized using a RT-PCR kit (TaKaRa, Japan), qPCR reactions were performed using SYBR Green (TaKaRa Biotech, Dalian, China) according to the manufacturer's protocol. The reaction was conducted by combining 12.5 µl of 2× Real Master Mix (TaKaRa Biotechnology), 2 μL of cDNA, 1 μL each of the upstream and downstream primers, and 8.5 μL of RNase-free ddH2O water. Speci c primers were designed using the Primer Premier 5.0 program (Table 3) and con rmed with BLAST. The expression levels of gene were normalized to linear GAPDH levels using 2 (−ΔΔCt) method [56] and the statistical difference was analysed using SPSS17.0. The correlation between the results of RNA-Seq and qPCR was calculated using a correlation test.

QTL analysis of Differentially Expressed lncRNAs
To predict the functions of DELs, QTLs analysis was performed. The location information of DELs was obtained from the transcriptome le and the QTL data of pig were downloaded from Animal QTLdb. Next, BEDTools and the "intersectBed" command were used for QTL analysis [57].

Prediction of PTGs of lncRNAs
Two methods were used to predicte the PTGs of lncRNAs. We identi ed the PTGs regulated by lncRNAs in cis, which were de ned as proteincoding genes located at 100 Kb upstream and downstream of the lncRNA, by BEDTools 2.17.0 [58]. The trans regulation of a lncRNA and its PTG was identi ed by the expression level correlation analysis of each pair of lncRNA and PTG. According to the Pearson's correlation coe cients (|r| > 0.95), the PTGs were selected to construct a lncRNA-mRNA co-expression network.

GO and KEGG Pathway Enrichment Analysis
The prediction of the lncRNA and PTGs' function such as gene enrichment and pathway analysis were performed using gene ontology (GO) (http://www.geneontology.org/) [59] and Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg) [60]. To further explore the main biological functions of differentially expressed genes, KOBAS software was used to detect the statistical enrichment of lncRNA PTGs or DEGs in KEGG pathways. The enrichment ndings with a q-value of less than or equal to 0.05 was considered statistically signi cant.

Declarations
Ethics approval and consent to participate Consent to publish