Overview of the sequencing data
We constructed cDNA libraries of miRNAs, mRNAs, and lncRNAs using the RNA from the ovary and testis. After filtering out low-quality transcripts, 5′ and 3′ adapters, and reads < 18 nt, a total of 113.5 M clean reads were produced by Illumina technology for miRNAs. The 21 and 22 nt length tanscripts were the most abundant (Fig. 1a), and 60.4% high-quality reads were mapped to the turtle genome (Pelsin-1.0, NCBI). We obtained 153.25 Gb clean reads for mRNA and lncRNA sequencing and approximately 84.41% to 87.72% clean reads were mapped to the turtle genome after removing low-quality reads (Fig. 1c). The length distribution of lncRNAs and mRNAs were shown in Fig. 1c and 1d. After mapping the genome, approximately 84.41% ~ 87.72% of the reads were mapped to inergenic regions in the P. sinensis reference genome (Fig. 1d, Additional 1).
Identification of differentional expression mRNAs, miRNAs, and lncRNAs
According to the miRNA expression profiles, we detected 10 446 novel miRNAs. A total of 633 miRNAs were significantly differently expressed between the ovaries and testes (P < 0.05), including 138 up-regulated miRNAs and 495 down-regulated miRNAs (Fig. 2a, 2d, Additional file 2). These DEmiRNAs belong to 438 families (Additional file 3). Among these DEmiRNAs, we identified a set of miRNAs which had been reported that might regulate animal reproduction, including miR-133, miR-138, miR-145, miR-143, miR-378, etc.
We detected 20 414 mRNAs, and 11 319 mRNAs were differently expressed based on sex including 5206 up-regulation mRNAs and 6113 down-regulation mRNAs (Fig 2b). A total of 28 500 lncRNAs with 10495 lncRNAs were differently expressed were detected including 1716 up-regulation lncRNAs and 8779 down-regulation lncRNAs between ovaries and testes (Fig 2c). Among the differentional expression lncRNAs and miRNAs, 3 miRNAs and 917 lncRNAs exhibited testis-specific expression, and 186 miRNAs and 4491 lncRNAs showed ovary-specific expression. Prediction of the potential targets of lncRNAs in cis and trans was performed to investigate the function of lncRNAs (Additional file 5). After searching for protein-coding genes 100 kb upstream and downstream, 3 904 DElncRNAs were found to correspond to regulating protein-coding genes in cis. The target genes included Foxl2, Cyp19a1, Gper, Esr, Dazl, and Sox30, etc, which suggest that the reproductive process might be regulated by the action of these lncRNAs on the protein-coding genes. Conversely, we identified 2 160 lncRNAs showing trans action by LncTar, including a set of genes that might regulate reproduction.
Functional analysis of DEmiRNAs and DElncRNAs
To annotate the molecular function of the differentially expressed miRNAs, both the RNA hybrid and miRanda software were used to improve the prediction of miRNA targets, resulting in 8 088 targets genes including 2 814 differentially expressed genes potentially regulated by 633 DEmiRNAs. GO categories of miRNAs and lncRNAs were assigned to all target genes based on the following three ontologies: cellular component, molecular function, and biological process (Fig. 3a, Additional file 6, 7, 8). Functions of target genes in the cellular component category were mainly focused on cell part, cell, and membrane. Based on the molecular function, the most abundant targets genes were focused on binding, followed by catalytic activity. Regarding the biological process, the most abundant target genes were focused on the single organism process, followed by cellular process, and biological regulation.
KEGG pathway enrichment analysis revealed that the DEmiRNAs were involved in 167 signalling pathways (Additional file 9). The identified metabolism networks were related to neuroactive ligand-receptor interaction and the regulation of the actin cytoskeleton. The most abundant target genes of DEmiRNAs focused on glyoxylate and dicarboxylate metabolism. We detected at least 12 pathways involved in reproductive biology, including oocyte meiosis, TGF-β signalling, ovarian steroidogenesis, GnRH signalling, Wnt signalling, cAMP signalling, steroid biosynthesis, MAPK signalling, p53 signalling, RNA polymerase, metabolism of xenobiotics by cytochrome P450, and mTOR signalling.
KEGG pathway enrichment analysis showed that the DElncRNAs were involved in 221 signalling pathways in a cis-regulatory manner and 225 signalling pathways in a trans-regulatory manner (Additional 10, 11). The KEGG pathway enrichment analysis revealed that the DGE lncRNAs were involved in oocyte meiosis, steroid hormone biosynthesis, Wnt signalling pathway, GnRH signalling pathway, p53 signalling pathway, apoptosis, MAPK signalling pathway, AMPK signalling pathway, TGFβ signalling pathway, cAMP signalling pathway, RIG-I-like receptor signalling pathway, mTOR signalling pathway, and insulin signalling pathway.
Validation of differentially expressed miRNAs and lncRNAs
To validate the sequencing data of miRNAs and lncRNAs, ten DEmiRNAs and ten DElncRNAs were randomly selected to test the relative expression in ovaries and testes. The expression of eight miRNAs and seven lncRNAs in ovaries and testes were consistent with the results of RNA sequencing. Among the miRNAs, novel-miR-1361, novel-miR-2322, novel-miR-6721, novel-miR-10042, novel-miR-10231, novel-miR-10322, and novel-miR-10468 were downregulated in testis, while novel-miR-1236 were upregulated in testes (Fig. 3a). Among the lncRNAs, MSTRG.435295.1, MSTRG.88998.1, MSTRG.127189.1, MSTRG.100955.1 were upregulated in testes, while MSTRG.129036.2, MSTRG.281180.2, MSTRG.561412.1 were downregulated in testis (Fig. 3b). The expression of these miRNAs and lncRNAs patterns among different groups were well-matched with the RNA-Seq data, which could guarante the accuracy of subsequent functional analysis.
Construction of ceRNA network
To construct the ceRNA networks, we screened miRNAs that included miRNA response elements, which could bind with both lncRNAs and mRNAs. We constructed a series of ceRNA networks of the mRNAs, miRNAs, and lncRNAs related with the DGE by integrating the expression profiles and regulated relationships among the mRNAs, lncRNAs, and miRNAs from the high-throughput sequencing data (Additional file 12). These networks included 102 DEmiRNAs, 635 DEmRNAs, and 1 621 DElncRNAs. The DEmiRNAs included novel-miR-227, novel-miR-9914, novel-miR-6375, novel-miR-1222, novel-miR-6721, novel-miR-2026, novel-miR-6671, novel-miR-642, novel-miR-6319, and novel-miR-42, etc. These ceRNA networks included a set of mRNAs regulating reproduction (Fig. 4a, 4b, Additonal file 12). For instance, Dazl mRNA and MSTRG.71049.8 shared a common binding site of miRNA novel-miR-1222. We also identified Wt1, CREB3l2, Gata4, Wnt2, Nr5a1, Hsd17, Igf2r, H2afz, Lin52, Trim71, Zar1, and Jazf1 in the ceRNA network. These miRNAs and mRNAs participate in regulating the reproductive process, including meiosis and spermatogenesis.