RNA-Seq read filtering, mapping, and transcript assembly
To identified lncRNAs function in goat skeletal muscle growth, we built two cDNA libraries from longissimus dori samples at two grwoth stages: 1-month-old and 9-month-old. Three biological replicates were used. A total of 89.57Gb of raw data were generated. Low-quality and adaptor sequences were filtered out, and an average of ~14.58Gb high-quality RNA-seq data remained. For further analysis, we found the GC content averaged ranged from 53.71% to 64%. Moreover, the quality scores of Q20 and Q30 were above 96% and 91%, respectively (Table 2). The high-quality data were first mapped to the goat reference genome (CHIR_1.0, NCBI). Approximately 72.47%-85.25% of the high-quality reads were mapped to the goat reference genome (Additional file 1). The mapped sequences in the libraries were assembled and a total of 65,247 library transcripts were acquired
Identification and confirmation of lncRNAs in goat longissimus doris
We set a highly stringent filtering pipeline to discard transcripts that did not have all the characteristics of lncRNAs. The assembled transcripts from the two libraries were filtered to get candidate lncRNAs. A total of 3,441 lncRNAs were screened (Additional file 2), which included 1,281(37.2%) long intergenic ncRNAs (lincRNAs), 805 (23.4%) antisense lncRNAs, and 1,355 (39.4%) sense_overlapping lncRNAs (Fig. 2a). There was no sense_intronic lncRNAs. These 3,441 putative lncRNA were encoded by 2,675 genes. That is 1.3 transcripts on average per locus (Additional file 2). The lncRNA transcripts were widespread in chromosomes including 29 autosomes and the X chromosome (Additional file 3), which reflected the function diversity of lncRNAs.
The Illumina RNA-seq also produced 66 mRNAs. As shown in Fig. 2b and 2c, the lncRNAs transcripts length and exon number both were lower than that of the mRNA. The average length of lncRNA was 2,001bp with an average of 3.2 exons. Importantly, the principal lncRNA transcripts with 2 exons accounted for 63.6% of the 3,441 lncRNAs (Additional file 2). Moreover, compared to the mRNA, the average length of lncRNA open reading frame (ORF) was shorter (Fig. 2d). The result showed that lncRNA has lower coding potential.
Differential expression and functional enrichment analysis
To dissect the crucial lncRNAs involved in goat skeletal muscle growth, we explored the DE lncRNAs (P < 0.05, |log2FoldChange | > 2) for 1-month-old vs. 9-month-old group. The lncRNA expression profiles were assessed with FPKM. In this study, 36 lncRNAs were differentially expressed in goat skeletal muscle at these two developmental stages (Additional file 4). Compared with 9-month-old stage, 8 DE lncRNAs were upregulated and 28 were downregulated lncRNAs in 1-month-old stage (P < 0.05) (Fig. 3).
To evaluate the potential function of lncRNAs, we predicted the possible target gene associated with lncRNA in location (co-location) and expression (co-expression) relationships. In the co-location analysis, 30 lncRNAs (1 annotated lncRNAs and, 29 novel lncRNAs) were transcribed close to 71 protein coding genes (Additional file 5). Besides, the GO analysis of these 71 targets gene showed that 210 GO terms were significantly enriched (P < 0.05), including biological process (BP), cellular component (CC) and molecular function (MF) (Additional file 6). Notably, the significant enriched GO terms that were identified include skeletal system development, negative regulation of G1/S transition of mitotic cell cycle, negative regulation of cell cycle G1/S phase transition, myoblast fate determination and myosin complex (Additional file 6, Fig. 4a). Only the top 30 GO terms are shown in Additional file 6 and Fig. 4a. Additionally, 57 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were enriched through the pathway analysis. The top 20 significantly enriched KEGG analyses were shown in Additional file 6 and Fig. 4b. The DE lncRNAs target genes that were related to muscle development including Glycerolipid metabolism, signaling pathways regulating pluripotency of stem cells, fatty acid metabolism and biosynthesis of amino acids (Fig. 4b).
Furthermore, in the co-expression analysis, the top 200 target genes were chosen for next functional cluster analysis. The results showed that 3 lncRNAs (novel lncRNAs) corresponded to 103 target genes (Additional file 7). In total, 230 GO terms were obviously enriched (P < 0.05), of which the top 30 were shown in Additional file 8 and Fig. 4c. The skeletal muscle growth- and development-related GO terms, called regulation of transcription elongation from RNA polymerase I promoter, negative regulation of canonical Wnt signaling pathway, muscle tissue morphogenesis, muscle organ morphogenesis, cGMP biosynthetic process and negative regulation of myotube differentiation, were found (Additional file 8). Pathway analysis indicated that these target genes were involved in inositol phosphate metabolism, phosphatidylinositol signaling system, Notch signaling pathway, and HIF-1 signaling pathway (Additional file 8, Fig. 4d). Overall, DE lncRNAs and their target genes showed great potential in the regulation of skeletal muscle growth and development.
lncRNA-mRNA interaction network associated with goat skeletal muscle development
To investigate how lncRNAs interact with their cis- or trans- acting target genes, we constructed the lncRNA-mRNA regulatory networks associated with skeletal muscle development (Additional file 6， file 8). A total of 23 lncRNA gene pairs were obtained. The network provided candidate lncRNAs with muscle growth and development. Consequently, cis-regulatory network visualization revealed 13 interactions between 12 lncRNA and 12 gene that was related to muscle development (Fig. 5a), while trans-regulatory networks showed 10 interaction between2 lncRNAs and 5 genes (Fig. 5b).
Verification the expression profiles of DE lncRNAs with RT-qPCR
Ten DE lncRNAs were randomly selected to explore their expression profile at two experimental time (1 month, 9 months) using RT-qPCR. The data showed that the10 lncRNAs exhibited differential expression patterns during the process of goat skeletal growth and development (Fig. 6). Importantly, the expression trend of 10 lncRNA between RNA-seq and RT-qPCR were similar. The results reveal that the pipeline we set was reasonable in identifying putative lncRNAs, suggesting that the RNA-seq were credible.
The ceRNA networks construction and bioinformatics analysis
To further elucidate the potential role of lncRNA in goat skeletal muscle growth and development, we constructed ceRNA networks of DE mRNA, miRNA and lncRNAs (P < 0.05, |log2FoldChange| > 2) by Cytoscape software (Additional file 9), and found some ceRNA networks associated with muscle development based on the given muscle development-related mRNAs and DE lncRNA (Fig. 7). In the networks, there were 3 miRNAs,8 muscle development related mRNAs, 4 lncRNAs transcript, which had at least one predicted target miRNA. Moreover, miRNAs form the center of the network with lncRNA as the bait, and mRNA as the target, suggesting that lncRNA acts as a sponge of miRNA to regulate gene expression. For instance, lncRNA XR_001917947.1 act as a sponge for miRNAs, such as chi-miR-30b-3p, may to regulated mRNA smad3. This ceRNA network may provide valuable information for goat skeletal muscle myogenesis.