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).
Table 1
The six library simples information.
Sample name
|
Raw
reads
|
Clean
reads
|
Clean
bases
|
Mapped Reads
|
Mapping Ratio
|
Error
rate(%)
|
Q20(%)
|
Q30(%)
|
GC
content(%)
|
MFM2DY4_1
|
111656928
|
107238580
|
16.09G
|
80514662
|
75.08%
|
0.02
|
96.59
|
91.73
|
47.14
|
MFM2DY4_2
|
121457850
|
116315790
|
17.45G
|
88556504
|
76.13%
|
0.02
|
96.78
|
92.03
|
52.18
|
MFM2DY4_3
|
111275282
|
106746570
|
16.01G
|
82037293
|
76.85%
|
0.02
|
96.75
|
91.96
|
50.09
|
DFM2DY4_1
|
118575510
|
111359682
|
16.7G
|
92421276
|
82.99%
|
0.02
|
96.51
|
91.79
|
56.43
|
DFM2DY4_2
|
107494708
|
101760644
|
15.26G
|
82733818
|
81.3%
|
0.02
|
96.7
|
92.14
|
49.84
|
DFM2DY4_3
|
142313048
|
138324398
|
20.75G
|
115239405
|
83.31%
|
0.02
|
96.52
|
91.55
|
53.96
|
Identification and characterisation of lncRNA
All transcripts were screened in five 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 confirmed the accuracy of selection.
DE Analysis of lncRNAs and mRNAs
We use cufflinks to normalize transcript expression to FPKM value, and performed a differential expression transcript analysis between the Meishan and Duroc medium follicle samples. A total of 201 differentially expressed lncRNAs (DELs) [p value<0.05 and |log10(foldchange)| > 1.3] between the two breeds were detected, in which 117 DELs were downregulated, 84 were upregulated in Meishan follicles (Figure 2 A). Moreover, 675 DE protein-coding genes were detected, and 265 DEGs were upregulated and 410 were downregulated in Meishan follicles (Figure 2 B). DE lncRNAs and mRNAs were widely distributed on chromosomes, with larger numbers on chromosomes 1, 6 and 13 (Figures 2 A and B). 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).
Table 2
Differentially expressed lncRNAs (DELs) and their target genes (PTGs).
Up DELs
|
Numbers of PTGs
|
Down DELs
|
Number of PTGs
|
PTGs
|
Up
Regulated
|
Down
Regulated
|
PTGs
|
Up
Regulated
|
Down
Regulated
|
ALDBSSCT0000002624
|
4
|
4
|
0
|
LNC_000179
|
77
|
0
|
77
|
LNC_000811
|
15
|
13
|
2
|
LNC_000715
|
69
|
68
|
1
|
ALDBSSCT0000004325
|
21
|
18
|
3
|
LNC_001511
|
29
|
1
|
28
|
LNC_001167
|
48
|
48
|
0
|
LNC_000846
|
8
|
2
|
6
|
LNC_000252
|
36
|
35
|
1
|
LNC_001333
|
5
|
1
|
4
|
ALDBSSCT0000009309
|
26
|
26
|
0
|
LNC_000802
|
2
|
0
|
2
|
LNC_001490
|
26
|
10
|
16
|
ALDBSSCT0000010765
|
29
|
2
|
27
|
LNC_000718
|
2
|
2
|
0
|
ALDBSSCT0000010764
|
17
|
0
|
17
|
LNC_000460
|
20
|
19
|
1
|
LNC_000116
|
49
|
7
|
42
|
LNC_001307
|
17
|
10
|
7
|
ALDBSSCT0000001721
|
11
|
0
|
11
|
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, significantly. 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 significantly 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 significantly (p<0.05). They were all downregulated in Meishan compared with Duroc sows. LNC_000644 and ALDBSSCT0000006057 were correlated with CYP3A7 significantly (p<0.05) and were unregulated in Meishan. ALDBSSCT0000000051 was correlated with its target gene HSD17B8 significantly (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.
Screening of potential reproduction-related lncRNAs in PI3K-AKT signalling pathway in pig ovarian follicle
The PI3K-AKT signalling pathway is very important in porcine follicular development. In order to explore the regulatory roles of the DE-lncRNAs involved in ovarian follicle growth and the relationship between these lncRNAs and PI3K-AKT signalling pathway, we analysed the expression of lncRNAs and their PTGs in PI3K-AKT signalling pathway. We found 12 protein coding genes, namely, THBS1, ITGA3, ITGA6, ITGB1, ITGB4, ITGB7, PIK3C2B, AKT2, CREB3L3, IKBKG and TP53 in the PI3K-AKT signalling pathway that were regulated by 19 DE-lncRNAs. Some DE-lncRNAs and their PTGs showed significant positive correlation (p<0.05), as follows: extracellular matrix (ECM) THBS1 and lncRNA-ALDBSSCT0000004244; ITGA3 and the novel lncRNAs LNC_000857 and LNC_001052; ITGA6 and ALDBSSCT0000010443; ITGB4 and LNC_000536; ITGB7 and ALDBSSCT0000004232; PIK3CA and ENSSSCT00000036444; AKT2 and LNC_000751, LNC_000819 and LNC_000058; CREB3L3 and ALDBSSCT0000008900 and ALDBSSCT000001184; IKBKG and ALDBSSCT0000004325 and LNC_001172; and TP53 and ALDBSSCT0000002268. Some DE-lncRNAs and their PTGS shown significant negative correlation (p<0.05), as follows: ITGA6 and ALDBSSCT0000009356; ITGB1 and LNC_001556; PIK3C2B and ALDBSSCT0000000805; CREB3L3 and ALDBSSCT0000011847; and TP53 and LNC_000076 (Figure 6).