Characterization of the ovary transcriptome and identification of mRNA and lncRNAs
To identify key differences associated with reproductive efficiency in sows with extreme phenotypes, we constructed 16 cDNA libraries from ovarian tissues during the follicular and luteal phases of the estrous cycle. A total of 1,509,728,794 raw reads were generated from 16 porcine ovary samples. After removing adapter sequences and low-quality sequences, 1,500,044,340 clean reads were retained and used for further analysis. In each sample, the percentage of clean reads ranged from 99.20% to 99.52% (Additional file 1: Table S1-1). In addition, most clean reads were aligned to the reference genome (Sscrofa10.2), accounting for 81.49% to 84.48% (Table S1-2).A total of 1,122,470 transcripts were assembled by Cuffmerge and Scripture [31]. According to the characteristics of lncRNAs, we used four tools (CPC, CNCI, CPAT and PFAM) to discard potential coding transcripts. In the end, 24,447 lncRNA transcripts were identified (Additional file 2: Table S2). These included 6,392 anti-sense lncRNAs (26.15%) and 18,055 intergenic lncRNAs (73.85%) (TableS2). In addition, 27,370 mRNAs were identified by mapping Illumina RNA-seqreads (Additional file 3: Table S3).
microRNA sequencing and identification
Small RNA sequencing data from ovarian tissue from 16 sows was generated on Illumina HiSeq and provided a total of 433,792,963 raw reads. After filtering out the low-quality sequences, including adaptor dimmers and less than 18 nt, 363,273,750 clean reads were ultimately achieved. The percentage of clean reads ranged from 74.75% to 88.31% in each small RNA library (Additional file 4: Table S4-1). The length distribution of clean reads showed that most of the reads were 20–24 nt, and 22 nt was the most abundant length identified. Such reads accounted for 36.30% of the total sequences (Additional file 4: Table S4-2). In total, 216 known miRNAs (Additional file 5: Table S5-1) and 1,724 novel miRNAs (Additional file 5: Table S5-2) were identified, while 198 known miRNAs were expressed in all four groups (Fig. 1).
Genomic features and expression patterns of lncRNAs
Overall, 24,447 lncRNAs and 27,370 mRNAs were detected in the ovaries of all 16 individual sows. In order to examine the differences in features between lncRNAs and mRNAs in ovarian tissues, their lengths were compared. The average length of lncRNAs was 2,955 bp, which was longer than that of the mRNAs (Fig. 2A). We also observed that the number of exons of lncRNAs was lower than that of the mRNAs, which tend to contain 2.3 exons (Fig. 2B). The ORFs of the lncRNAs were shorter than those of the mRNAs (Fig. 2C). Lastly, their expression levels were also compared (Fig. 3); In general, lncRNAs had lower expression levels.
lncRNAs can act in cis or trans to positively or negatively regulate gene expression; however, cis-acting lncRNAs are restricted to the chromosome from which they are transcribed [32]. Several studies also demonstrated that some lncRNAs have a high correlation with expression of neighboring gene [32, 33]. To further explore the relationship between lncRNAs and their neighboring coding genes in ovarian tissues, we searched for neighboring protein-coding genes (<10k) of all the identified lncRNAs and analyzed gene pairs formed by lncRNAs and their neighboring genes. We identified 4,044 protein-coding genes: coding gene pairs (873 in divergent) and 1,664 lncRNA: coding gene pairs (195 in divergent) (Fig. 4). We observed that the expression pattern of lncRNAs with their neighboring gene pairs (average Pearson correlation: 0.20) was similar to the coding gene pairs (average Pearson correlation: 0.28) and exhibited a significantly higher correlation than random coding gene pairs (average Pearson correlation: 0.041, P <0.01) (Fig. 5A). We observed that there was a relatively low correlation between divergent lncRNAs: coding gene pairs (average correlation: 0.19) than divergent coding gene pairs (average Pearson correlation: 0.30, P <0.05), and a higher correlation compared with random coding gene pairs (average Pearson correlation: 0.013, P <0.01) (Fig. 5B). This result indicated that the correlation between lncRNAs and their neighboring gene was higher than random coding gene pairs.
Identification of differentially expressed mRNAs, lncRNAs and miRNAs between the high and low fertility groups
From the expression profiles, differentially expressed mRNAs, lncRNAs and miRNAs in the ovaries of Large White pigs were obtained by comparing LH vs. LL and FH vs. FL (Table 1). A total of 956 (345 up-regulated and 611 down-regulated) lncRNA transcripts were differentially expressed in LH vs. LL (P <0.05), while 415 (247 up-regulated and 168 down-regulated) were differentially expressed in FH vs. FL (P <0.05) (Additional file 6: Table S6-1 and 2). We also identified 457 mRNA transcripts that were differentially expressed between the LH and LL groups (Table 1) (Additional file 6: Table S6-3). Among these transcripts, 334 were annotated as known genes. In the FH vs. FL comparison, we found that 475 mRNAs were differentially expressed, while 316 mRNAs were annotated (Additional file 6: Table S6-4). Analyses of the small RNA sequencing data showed that 29 and 11 known miRNAs were differentially expressed when comparing LH vs. LL and FH vs. FL, respectively (Additional file 6: Table S6-5 and 6).
Function enrichment analysis of the lncRNAs
To investigate the function of the differentially expressed lncRNAs in each comparison, the potential targets of lncRNAs were predicted in this study. GO analysis revealed that there were 18 and 15 significant GO terms (corrected P <0.05) in LH vs. LL and FH vs. FL, respectively (Additional file 7: Table S7-1 and 2). We noticed that three significant GO terms were common in all four groups: catalytic activity, single-organism metabolic process and vitamin D metabolic process. The KEGG analysis indicated that a total of 23 and 14 significant pathways were found in LH vs. LL and FH vs. FL, respectively (Additional file 7: Table S7-3 and 4). Importantly, it was observed that “ovarian steroidogenesis” and “lysosome” were the specific enrichment pathways in LH vs. LL, while steroid biosynthesis was the common pathway in the four comparison groups.
Target prediction of miRNAs and construction of miRNA–mRNA networks
To understand the biological functions of differentially expressed miRNAs on fertility, we predicted the potential target genes of these miRNAs in two comparisons. We found that there were 13,458 putative target sites for 122 miRNAs in LH vs. LL and 4,466 target sites for 46 miRNAs in FH vs. FL (Additional file 8: Table S8-1 and 2). Furthermore, GO and pathway enrichment analyses were performed. GO analysis of the target genes revealed that there were 410 and 236 significant GO terms (corrected P <0.05) in LH vs. LL and FH vs. FL comparisons (Additional file 9: Table S9-1 and 2). KEGG pathway analysis revealed that a total of 97 and 31 significant pathways (Hypergeometric Distribution Hypothesis Test, P <0.05) were identified in LH vs. LL and FH vs. FL comparisons, respectively (Additional file 9: Table S9-3 and 4). Among these KEGG pathways, multiple pathways were closely involved in the reproductive process, such as the Insulin signaling pathway, MAPK signaling pathway, Estrogen signaling pathway, GnRH signaling pathway, PI3K-Akt signaling pathway, Ras signaling pathway, Cytokine-cytokine receptor interaction, Jak-STAT signaling pathway and Lysosome pathway in LH vs. LL, and the Notch signaling pathway, TGF-beta signaling pathway and Steroid biosynthesis in FH vs. FL. It is worth noting that the Wnt signaling pathway, Insulin secretion and Adherens junction were common in LH vs. LL and FH vs. FL.
Additionally, we aimed to illustrate negative interactions between differentially expressed miRNAs and mRNAs in the porcine ovary that might lead to differences in fertility; thus, regulatory networks of miRNA–mRNA pairs were constructed (Fig. 6). Of these negative interactions, three miRNAs (ssc-miR-1307, ssc-miR-1343 and ssc-miR-671-5p) targeted multiple mRNAs, but several miRNAs targeted only one mRNA. For example, up-regulated ssc-miR-1307 targets 12 genes, but down-regulated ssc-miR-361-3p targets only one gene (Fig. 6A). Moreover, down-regulated progestin and adipoQ receptor 7 (PAQR7) (FH vs. FL) was regulated by two differentially expressed miRNAs including ssc-miR-885-5p and ssc-miR-671-5p (Fig. 6B).
Construction of lncRNA-miRNA-mRNA networks
To explore the role and relation of lncRNAs and miRNAs mediation in pig fertility, differentially expressed lncRNAs were selected by miRanda analysis [34]. The lncRNA–miRNA negative pairs between differently expressed lncRNAs and miRNAs were selected to construct the co-expression network. In the LH vs. LL comparison, we found that the key miRNAs interacted with 19 lncRNAs (Fig. 7A). In FH vs. FL group, the key miRNAs interacted with 7 lncRNAs (Fig. 7B). A total of 19 and 7 lncRNA–miRNA pairs were identified in LH vs. LL and FH vs. FL, respectively. It is worth noting that most lncRNAs were targeted by the same miRNA. Among these key miRNAs, ssc-miR-1343 and ssc-miR-671-5p had more interactions than others. Ssc-miR-1343 is the key miRNAs targeted with nine key lncRNAs (TCONS_00009287, TCONS_00196796, TCONS_00309415, TCONS_00309419, TCONS_00372560, TCONS_00521720, TCONS_00521721, TCONS_00703255, and TCONS_00814106) through MREs, and ssc-miR-671-5p targeted with six key lncRNAs (TCONS_00019076, TCONS_00229497, TCONS_00429823, TCONS_00651713, TCONS_00702922, and TCONS_00817482), which may be key regulators related to fertility.
Based on the above data, we integrated the lncRNA–miRNA interactions and miRNA–mRNA interactions to establish lncRNA–miRNA–mRNA networks and then visualized using the Cytoscape software (Fig. 8). The network of LH vs. LL was composed of 44 nodes and 40 edges, and the nodes included 4 miRNAs, 14 lncRNAs and 27 mRNAs, which could be the important nodes involved in the ceRNA network during the luteal phase of the estrous cycle (Fig. 8A). In this network, some of them have been reported to be reproduction-associated molecules such as NUMBL, ILF3, GRIK4, SLC9A1, TGFBR1and LOXL4. We noticed that nine lncRNAs were interrelated with ssc-miR-1343 and may act as ceRNA to inhibit target miRNAs and mediated, related hub genes translation such as NUMBL, ILF3, TGFBR1, TMEM8B, PRR14, TSHZ2, and CAMKV. In addition, we found that TCONS_00309450 and TCONS_00429684 may serve as ceRNA to mediate GRIK4 by sponging ssc-miR-1249. In the FH vs. FL group, there were 23 nodes and 22 edges, consisting of 1 miRNA, 6 lncRNAs and 16 mRNAs (Fig. 8B). These six lncRNAs may serve as ceRNA to mediate the corresponding gene transcripts by sponging ssc-miR-671-5p. We also found that several genes, such as GRIN2D and FZD5, were mainly involved in the “cAMP signaling pathway”, “Calcium signaling pathway”, “Wnt signaling pathway” and “mTOR signaling pathway”, implying that they might be acting as reproduction related genes. Thus, we hypothesize that these lncRNAs, miRNAs and mRNAs play critical roles in fertility regulation.
RT-qPCR Verification
In the lncRNA–miRNA–mRNAs interaction networks, we selected 14 and 5 key nodes to validate expression levels in LH vs. LL and FH vs. FL, respectively, using RT-qPCR (Fig. 9). The expression of each miRNA was significantly higher in the LL groups compared to the LH groups. In contrast, the expression of each lncRNA or mRNA was significantly lower in the LL groups than in the LH groups (Fig. 9A-D). We also found that the expression of ssc-miR-671-5p was higher in the FL groups compared to the FH groups (Fig. 9E). These results suggested that the post-transcriptional regulatory functions of miRNAs negatively correlated with their targets and that these differentially expressed miRNAs and lncRNAs may contribute to the fertility differences in sows with extreme phenotypes during the F and L phases of the estrous cycle. In brief, the results demonstrated that the expression patterns of 19 differentially expressed genes were consistent between the RT-qPCR data and the RNA-Seq data, implying that the accuracy of RNA-Seq data was reliable (Fig. 9).