Determination of IMF content
Profiles of lncRNAs and mRNAs in sheep muscle
A total of six RNA expression profiles were generated in this experiment. The results are shown in Table 2. The average raw reading was 13.62 G (1G means 1*109 base). After preprocessing the raw data, the average value of the filtered data obtained in each library was 12.82 G. The data obtained from the six expression profiles were relatively average with Q20 ≥ 99% and G/C contents ranging from 49 to 53%, indicating that the quality of the filtered data was reliable. The comparison rate between the filtered clean reads and the reference genome was greater than 88% in all six samples, indicating that the experiment was free of contamination and that the experimental results were robust.
An average of 24,384 expressed genes was identified in the six libraries, and a summary of the protein-coding genes identified is provided in Additional file 1 (Table S1). An average of 6,499 unique lncRNAs was identified in the libraries. The information associated with all identified lncRNAs is shown in Additional file 2 (Table S2). We found that the number of reads was positively related to the length of the chromosome (Fig. 1a). Based on the locations of novel lncRNAs in the genome, we identified 525 antisense lncRNAs, 304 sense lncRNAs, 350 bidirectional lncRNAs, 1,710 intronic lncRNAs and 4,046 intergenic lncRNAs (Fig. 1b). The sequence information for all identified lncRNAs is shown in Additional file 3 (Table S3).
The structural characteristics and the expression levels of lncRNAs and mRNAs were different. The average length of lncRNAs was 868 nt, which was shorter than the average length of mRNAs (2,131nt) (Fig. 1c). LncRNAs consisted of 1.7 exons on average, while mRNAs had 9.9 exons on average (Fig. 1d); thus, lncRNAs had fewer exons than mRNAs. Meanwhile, lncRNAs had lower expression levels relative to mRNAs (Fig. 1e). Moreover, the length of the open reading frame (ORF) of lncRNAs tended to be shorter than that of mRNAs (Fig. 1f, 1g). Overall, lncRNAs were characterized by shorter lengths, fewer exons, lower expression levels and shorter ORF length distributions compared with mRNAs.
Identification of differentially expressed mRNAs and lncRNAs
A total of 199 DEMs were identified in muscle tissue (log2 (fold change) ≥ 1 or log2 (fold change) ≤ -1 and q < 0.05). Of these differentially expressed genes (DEGs), 70 were up-regulated and 129 were down-regulated (Fig. 2a, 2c). A summary of DEGs is provided in Additional file 4 (Table S4A). We identified 61 lncRNAs that were differentially expressed, of which 25 lncRNAs were up-regulated and 36 lncRNAs were down-regulated (Fig. 2b, 2d). Interestingly, 58 DELs were novel lncRNAs, we will focus on these novel lncRNAs in subsequent research. The list of DELs is provided in Additional file 4 (Table S4B). To illustrate the distribution of DEGs, we created clustering maps of top 100 DEMs and all the DELs (Fig. 2e, f). Red indicates that the gene had a higher expression level, and blue indicates that the gene had a lower level of expression.
Among the DEMs, we found many genes that have been shown to be related to lipid deposition. IGF2, CAPN6, UCP2 and SOCS2 were highly expressed at Mth-2. IGF2 and CAPN6 were reported to affect the deposition of intramuscular fat and play an important role in meat efficiency [29,30]. UCP2 and SOCS2 can regulate the proliferation and differentiation of preadipocytes [31,32]. FOSB, SCD and CMYA5 were highly expressed at Mth-12. They have all been found to be potential candidate gene affecting meat quality [33-35]. Among the DELs, we found that MSTRG.13347.2, MSTRG.16157.3, MSTRG.11343.1, MSTRG.11343.4, MSTRG.10678.1 were highly expressed at Mth-2. We speculated that these novel lncRNAs might inhibit lipid deposition. Meanwhile, MSTRG.3004.3, MSTRG.21053.3, MSTRG.14887.1, MSTRG.790.1, MSTRG.10518.3 were highly expressed at Mth-12. We speculated that these novel lncRNAs might promote lipid deposition. However, the regulatory mechanisms underlying these lncRNAs require further study.
Enrichment analysis of differentially expressed mRNAs
GO functional enrichment analysis of DEGs revealed that these genes participated in a total of 346 significantly enriched functional classifications (P < 0.05), 235 of which were related to biological processes, 30 related to cellular components and 81 related to molecular functions (Additional file 5: Table S5A). The top 25 of biological processes, top 15 of cellular components, top 10 of molecular functions were shown in Fig. 3a. The most significantly enriched GO terms were: DNA binding (GO:0003677), extracellular region (GO:0005576), calcium ion binding (GO:0005509), extracellular space (GO:0005615), oxidation-reduction process (GO:0055114), oxidoreductase activity (GO:0016491).
In addition, results of the KEGG pathway analysis showed that these DEGs were involved in 126 biological pathways (Additional file 5: Table S5B), 18 pathways of which were significantly enriched, including arachidonic acid metabolism (ko00590), linoleic acid metabolism(ko00591), steroid hormone biosynthesis(ko00140), and retinol metabolism(ko00830), all of which were related to lipid metabolism. Moreover, the top 20 signaling pathways are shown in Fig. 3b. The results indicated that these pathways may have significantly contributed to the deposition of IMF.
Comprehensive analysis of candidate lncRNAs and mRNAs
To understand the potential function of novel lncRNAs, we performed cis-regulation and trans-regulation analyses on candidate lncRNAs. A total of 61 DELs regulated 49 DEMs, all the lncRNAs that acted on 49 mRNAs through trans-regulation (Additional file 6: Table S6). GO analysis of targets of lncRNAs revealed that these genes participated in a total of 422 GO terms, 180 of which were significantly enriched (P<0.05) (Additional file 7: Table S7A). In GO annotation, these DEGs primarily played a role in biological processes. For example, glyceraldehyde-3-phosphate metabolic process (GO:0019682), positive regulation of phospholipid translocation (GO:0061092), glyceraldehyde-3-phosphate biosynthetic process (GO:0046166), cell-cell junction assembly (GO:0007043), muscle cell cellular homeostasis (GO:0046716) and actin filament polymerization (GO:0030041). The top 25 of biological processes, top 15 of cellular components, top 10 of molecular functions were shown in Figure. 3c. The KEGG pathway enrichment analysis of target genes revealed a total of 47 annotated pathways (Additional file 7: Table S7B), of these pathways, one was significantly enriched (P<0.05), namely Tight junction(ko04530). Although some pathways were not significantly enriched, such as Wnt signaling pathway(ko04310), AMPK signaling pathway(ko04152), mTOR signaling pathway(ko04150), Cell adhesion molecules (CAMs)(ko04514) and MAPK signaling pathway(ko04010) and Jak-STAT signaling pathway, these pathways have been reported in the literature to play an important role in lipid deposition [36-41]. Overall, the significantly enriched pathways and GO terms involve 49 target genes. Of these, seven target genes are associated with lipid deposition. There were lncRNAs whose pearson correlation coefficients ≥0.9 or ≤-0.9 associated with these seven target genes (Table 3). Of these, seven have high expression levels: MSTRG.16157.3, MSTRG.21053.3, MSTRG.19941.2, MSTRG.2469.2, MSTRG.4051.3, MSTRG.21381.1, MSTRG.12864.1. A network describing the connections between the source genes and lncRNAs (whose pearson correlation coefficients ≥0.8 or ≤-0.8) was constructed (Fig. 4). However, the regulatory mechanisms underlying these lncRNAs require further study.
Validation of lncRNA and mRNA expression by qRT-PCR
To validate the expression levels of DELs and DEMs, we randomly selected five DELs and five DEMs and detected their expression levels by qRT-PCR (Fig. 5a). The results of RNA-seq are shown in Fig. 5b. Comparison of the two sets of results above revealed consistent regulatory trends of genes detected by the two methods, indicating that the RNA-seq data were accurate.