Overview of small RNA libraries
In order to identify differentially expressed miRNAs during the process of testicular development of Duroc and Meishan boars, six small RNA (sRNA) libraries of testis tissues of 20-, 75-, 270-day-old Duroc and Meishan boars (D20, D75, D270, M20, M75, and M270) were constructed and sequenced by the Illumina HiSeqTM 2000 platform. In total, 13,335,120, 13,051,493, 13,352,724, 13,606,755, 12,695,970, and 13,721,075 raw reads were generated in D20, D75, D270, M20, M75, and M270, respectively. After removing the low-quality sequences and adaptors, and then discarding the sequences shorter than 18 nt, 13,133,806, 12,892,394, 13,167,240, 13,441,877, 12,561,392, and 13,522,792 clean reads were obtained and used for further analysis (Table 1). A total of 1,078,105 and 878,097 unique sRNA from Duroc and Meishan boar testes were mapped to the porcine reference genome (Sscrofa10.2), respectively (Table 2). Read length distribution analyses of the six sRNA libraries showed that the dominant length of sRNAs was 22 nt, accounting for at least 36.24% of the population (Fig. 1). More 22 nt sRNAs were found in D20, D75, and M20 than in D270, M75, and M270. While few sRNAs with the length of 18 to 19 nt and 25 to 30 nt were detected, and sRNAs with the length of 25 to 30 nt may mainly represent Piwi-interacting RNAs (piRNA). These results are similar to those of previous study of pig [22].
Differentially expressed miRNAs across Duroc and Meishan boars
The expression profiles of known miRNAs from the six samples were analyzed and 36-139 significantly differentially expressed miRNAs (DE miRNAs) (P ≤ 0.05, |log2Ratio|≥1) were filtered in each pairwise comparison (Table 3; Additional file 9: Table S6). For example, the comparisons between the two breeds indicated that the number of DE miRNAs in the 75- and 270- day time points were substantially larger than that in the 20- day time point. These results suggested that significant differences existed during testicular development between Duroc and Meishan boars at the age of 75 and 270 days. More significantly DE miRNAs were detected in Meishan boars at age of 20 to 75 days than those in Duroc boars, which suggested that the testis developed at a faster rate in Meishan boars than in Duroc boars from 20 to 75 days. These findings were consistent with those in previous study of mRNA expression data of matched samples [17].
Venn diagram showed 45 significantly DE miRNAs were filtered from the four pairwise comparisons (D20 vs D270, D75 vs D270, M20 vs M75, and M20 vs M270) (Fig. 2a). The four pairs represented the comparisons before and after puberty since D270, M75, and M270 have reached puberty. Figure 2b showed two main sample branches (D20, D75, and M20 versus D270, M75, and M270), which indicated that the expression pattern of M20 was similar to those of D20 and D75, and that the expression pattern of D270 was similar to those of M75 and M270. The results were consistent with those of previous study of expression pattern of mRNA [17], demonstrating large differences existed in the process of testicular development between Meishan and Duroc boars. It could be concluded that miRNAs were pivotal factors regulating sexual function development.
Integrated analysis between differentially expressed miRNAs and target mRNAs in Duroc and Meishan boars at different stages
The mRNA expression data of six samples from our previous study were used for a pairwise integrated analysis [17]. Through the Trinity de novo assembly method, 20,525 non-redundant genes were obtained from the six samples (Additional file 1: Table S1), 19,310 (94.08%) and 18,241(88.87%) genes were mapped against Kyoto Encyclopedia of Genes and Genomes (KEGG) and (Gene Ontology) GO databases, respectively (Additional file 2: Table S2).
In order to better understand the potential roles of the miRNAs during the process of testicular development, computational target prediction was performed using Targetscan and miRanda. Then we performed integrated analyses of differentially expressed miRNAs and target mRNAs at the expression levels. A large number of correlated miRNA-mRNA pairs were detected in each pairwise comparisons (Fig. 3). The number of miRNA/mRNA-negative pairs between Duroc and Meishan boars at 75-day time point was obviously higher than that at 20- and 270-day time points, and more negative pairs were detected from 20 to 75 days in Meishan boars than those in Duroc boars (Fig. 3a). The previous study demonstrated that Meishan boars attained puberty and their testes generated sperms prior to 75 days earlier than Duroc boars [17]. These findings indicate that miRNAs as negative gene expression regulators significantly control the expression of genes involved in regulating the process of testicular development. Meanwhile, the number of miRNA/mRNA-positive pairs was also compared between Duroc and Meishan boars (Fig. 3b) with a similar tendency found in Duroc and Meishan boars at different ages. These results reveal that the subset of miRNAs may function as enhancers activating the transcription of genes which play an important role during the process of testicular development.
A representative miRNA-mRNA regulatory network of biological pathways was shown in Fig. 4. Eight DE miRNAs (ssc-mir-423-5p, ssc-mir-34c, ssc-mir-107, ssc-mir-196b-5p, ssc-mir-92a, ssc-mir-320, ssc-mir-10a-5p, and ssc-mir-181b) were selected from 45 miRNAs deriving from Fig. 2 serving as functional miRNAs in testicular development of Meishan and Duroc boars according to their annotations and the potential relationship between miRNAs and spermatogenesis and gonad development. GO analysis and KEGG functional annotation of potential target genes of eight DE miRNAs were performed to detect the functional characteristics of miRNAs. A large number of target genes were assigned to the functional categories related to sexual reproduction, male gamete generation, spermatogenesis, sperm development as well as meiosis, indicating that the eight miRNAs were highly involved in spermatogenesis and testis development. Five important pathways including GnRH, Wnt, p53, mTOR, and MAPK signaling pathway related to the regulation of male sexual function were enriched by functional genes including phospholipase C beta 1 (PLCβ1) involved in GnRH and Wnt signaling pathway with PLCβ1 being the target of ssc-mir-423-5p and ssc-mir-34c, serine/threonine/tyrosine interacting protein (STYX) involved in MAPK signaling pathway with STYX being the target of ssc-mir-320, ssc-mir-10a-5p, ssc-mir-92a and ssc-mir-107; cyclin D2 (CCND2), phosphatase and tensin homolog (PTEN), and cyclin B1 (CCNB1) with the latter 3 genes involved in the p53 signaling pathway, and so on. The representative miRNA-mRNA regulatory networks contained biological pathways regulating male sexual function, which illustrated a complex relationship and interaction between the two biomolecular types.
Verification of DE miRNAs and target verification of ssc-mir-423-5p
To evaluate our DE miRNAs library, the expression profiles of 3 DE miRNAs (ssc-mir-181b, ssc-mir-423-5p, ssc-mir-196b-5p), four DEGs from the miRNA-mRNA interaction networks, and ssc-mir-4334-3p from the Venn diagram in Fig. 2a, all of which were highly related to boar sexual function and reproduction, were further analyzed by quantitative real-time PCR (qRT-PCR) with specific primers. As shown in Additional file 3: Figure S1, the results of RNA-seq data and qRT-PCR data were identical. CYLD did not show consistent expression between RNA-seq and qRT-PCR data from Duroc and Meishan boars at age of 20 and 270 days, which was probably caused by the sensitivity of the different methods. In general, the results of qRT-PCR validated the RNA-seq results and demonstrated the reliability of our data.
ssc-mir-423-5p was one of the differentially expressed miRNAs and was selected as a candidate miRNA for analyzing male sexual function. PLCβ1 was predicted to be a target of ssc-mir-423-5p (Fig. 5a). The dual-luciferase reporter assay system analyzed the interaction between ssc-mir-423-5p and PLCβ1 gene. The analysis results indicated that Luciferase activity was significantly suppressed after we co-transfected ssc-mir-423-5p mimic (Additional file 4: Table S3). However, luciferase activity was not significantly changed when we co-transfected ssc-mir-423-5p mimic and pmirGLO- PLCβ1-3´-UTR -mut into Swine Testis (ST) cells (Fig. 5b). Meanwhile, ssc-mir-423-5p inhibitor significantly promoted luciferase activity after we co-transfected ssc-mir-423-5p inhibitor and luciferase reporter vector (pmirGLO- PLCβ1-3´-UTR), and luciferase activity was unchanged after we co-transfected ssc-mir-423-5p inhibitor and pmirGLO-PLCβ1-3´-UTR-mut into ST cells (Fig. 5c). qRT-PCR and western blotting analyses revealed that PLCβ1 mRNA and protein expression levels were significantly reduced after ssc-mir-423-5p mimic was transfected into ST cells, whereas the inhibition of ssc-mir-423-5p increased the expression of PLCβ1 mRNA and protein in ST cells (Fig. 5d-g, Additional file 5: Figure S2, Additional file 6: Figure S3). These results suggested that PLCβ1 was the target gene of ssc-mir-423-5p and the regulation of ssc-mir-423-5p in the process of spermatogenesis may be mediated by PLCβ1.