Determination of intramuscular fat content
Results for the IMF content of sheep are shown in Table 1. The IMF content of AFWS at 2, 4, 6 and 12 months increased gradually with growth in the two muscle tissues. In the longissimus dorsi muscle (LDM), the IMF content at Mth-4 was significantly higher than that at Mth-2 (P < 0.01) and was significantly lower than that at Mth-6 and Mth-12 (P < 0.01).No significant differences were detected between Mth-6 andMth-12. The same pattern occurred in the biceps femoris muscle (BFM). IMF content in the LDM was significantly higher than that in the BFM in the same month (P < 0.01). Thus, Mth-2 and Mth-12 were selected to conduct screening of differentially expressed mRNAs and lncRNAs.
Profiles of lncRNAs and mRNAs in sheep muscle
A total of 6 RNA expression profiles were generated in this experiment. The results are shown in Table 2. The average raw reading was 13.62 G. 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 filtered clean reads were compared with the reference genome with TopHat. The comparison rate in all six samples was greater than 88%, indicating that the experiment was free of contamination and that the experiment was suitable. The alignment data were assembled with String Tie and used for additional screening and analysis of lncRNAs.
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 S2A). We used circos (http://www.circos.ca) software to perform genomic mapping on the lncRNAs obtained by screening. 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 intronic lncRNAs, 1,710 intron lncRNAs and 4,046 intergenic lncRNAs (Fig. 1b, Additional file 2: Table S2B). The sequence information of all identified lncRNAs is shown in Additional file 3 (Table S3).
The structural characteristics and the expression levels of lncRNAs and mRNAs were different. In our study, the average length of lncRNAs was 868 nt, which was shorter than the average length of mRNAs (2,131 nt) (Fig. 1c). As illustrated in Fig. 1d, lncRNAs consisted of 1.7 exons on average, while mRNAs had 9.9 exons on average; 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 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 606 differentially expressed protein-coding genes (DEGs) were identified in muscle tissue (|log2(foldchange)| ≥ 1, P < 0.05). Of these DEGs, 154 mRNAs were up-regulated and 452 mRNAs were down-regulated (Fig. 2a, 2c). A summary of DEGs is provided in Additional file 4 (Table S4). Meanwhile, we identified 408 lncRNAs that were differentially expressed, of which 254 lncRNAs were up-regulated and 154 lncRNAs were down-regulated (Fig. 2b, 2d). The list of differentially expressed lncRNAs is provided in Additional file 5 (Table S5). To understand the overall distribution of differentially expressed genes, we created a heat map of differentially expressed mRNAs and lncRNAs (Fig. 2e, f).
We randomly selected 5 differentially expressed lncRNAs and 5 mRNAs whose expression levels were determined by qRT-PCR (Fig. 3a). The results of RNA sequencing are shown in Fig. 3b. By comparing the two sets of results above, the patterns of regulation of the two sets of genes detected by the two methods were consistent, indicating that the RNA sequencing data were accurate.
Enrichment analysis of differentially expressed mRNAs
GO functional enrichment analysis of DEGs revealed that these genes participated in a total of 419 significantly enriched functional classifications (P < 0.05), 282 of which were related to biological processes, 41 related to cellular components and 96 related to molecular functions (Additional file 6: Table S6A). The results of the biological process analysis revealed that these mRNAs were primarily involved in the regulation of processes related to metabolism (phospholipid metabolic processes and cholesterol metabolic processes), the immune response and lipid transport (reverse cholesterol transport and phospholipid efflux), as well as processes related to growth and development of fat and muscle, including skeletal system morphogenesis, cell maturation, fatty acid biosynthetic processes and lipid storage. The top 20 GO terms are displayed in Fig. 4a.
In addition, results of the KEGG pathway analysis showed that these DEGs were involved in 232 biological pathways (Additional file 6: Table S6B), 18 pathways of which were significantly enriched (Additional file 7: Table S7), including cholesterol metabolism (ko04979), arachidonic acid metabolism (ko00590) and glycine, serine and threonine metabolism (ko00260), which were related to fat metabolism. Systemic lupus erythematosus (ko05322) had the highest level of significance with 23 DEGs. Although multiple pathways were not significantly enriched, many fat-related DEGs were involved in these biological pathways. For example, SCD was involved in three pathways: the PPAR signaling pathway (ko03320), the AMPK signaling pathway (ko04152) and the biosynthesis of unsaturated fatty acids (ko01040). Moreover, the top 20 signaling pathways are shown in Fig. 4b. The results indicated that these pathways may have significantly contributed to the deposition of intramuscular fat.
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 183 differentially expressed lncRNAs regulated 218 differentially expressed mRNAs, 5 lncRNAs of which acted on 4 mRNAs through cis-regulation and 183 lncRNAs that acted on 218 mRNAs through trans-regulation (Additional file 8: Table S8). GO analysis of targets of lncRNAs revealed that these genes participated in a total of 1,840 GO terms, 546 of which were significantly enriched (P < 0.05) (Additional file 9: Table S9A). The results of biological process analysis showed that there were many genes involved in fat synthesis, transport, regulation and other processes. This finding is consistent with how SCD is known to be involved in unsaturated fatty acid biosynthetic processes and the involvement of BNIP3 in brown fat cell differentiation.
The results of KEGG pathway enrichment analysis of targeted genes showed that a total of 214 pathways were annotated according to the analysis of cis-regulation and trans-regulation of differentially expressed lncRNAs (Additional file 9: Table S9B). Of these pathways, 28 were significantly enriched (P < 0.05) (Additional file 10: Table S10). Most of these pathways were related to material metabolism and the FoxO signaling pathway. Biosynthesis of unsaturated fatty acids, phosphonate and phosphinate metabolism and ether lipid metabolism is known to be related to lipid deposition. In addition, although some pathways, such as the PPAR, Wnt and AMPK signaling pathways, were not significantly enriched, they still played an important role in lipid deposition.
Based on the comprehensive results of the KEGG enrichment analysis of target genes, we summarized DEGs related to lipid deposition (Table 3). Several DEGs in these critical pathways, such as SCD, ACAA2, FZD4, FADS2 and PLA2G4E, could be cis- or trans-regulated by several differentially expressed lncRNAs (Table 4). Based on GO terms and results of pathway enrichment analysis, we summarized lipid-related lncRNAs and their potential target genes (Additional file 11: Table S11). In sum, we suggest that these differentially expressed lncRNAs are likely involved in the regulation of the development and deposition of fats. However, the regulatory mechanisms underlying these lncRNAs require further study. Based on the patterns of enrichment of important pathways, we constructed a lncRNA-mRNA co-expression network (Fig. 5).