Determination of IMF content
Results for the IMF content of sheep are shown in Table 1. The IMF content of the longissimus dorsi muscle (LDM) at 2, 4, 6 and 12 months was 2.202 ± 0.006, 4.566 ± 0.178, 10.685 ± 0.690 and 11.163 ± 0.878, respectively. We found that the IMF content of LDM 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). The IMF content of LDM in Mth-12 was also significantly higher than that observed in Mth-2 (P < 0.01). No significant differences were detected between Mth-6 and Mth-12. The same pattern was observed for 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 for RNA sequencing (RNA-seq)
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. 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 HISAT2. The comparison rate in all six samples was greater than 88%, 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 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 bidirectional lncRNAs, 1,710 intronic lncRNAs and 4,046 intergenic lncRNAs (Fig. 1b, Additional file 2: Table S2B). 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,131 nt) (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 606 DEMs were identified in muscle tissue (log2 (fold change) ≥ 1 or log2 (fold change) ≤ -1 and P <0.05). Of these differentially expressed genes (DEGs), 154 were up-regulated and 452 were down-regulated (Fig. 2a, 2c). A summary of DEGs is provided in Additional file 4 (Table S4). 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 DELs is provided in Additional file 5 (Table S5). To illustrate the overall distribution of DEGs, we created clustering maps of DEMs and 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 that IGF2 was highly expressed at Mth-12. IGF2 has been reported to be a candidate gene, as it can play an important role in fat deposition in pig [29]. Some DEGs were also associated with lipid deposition in this study. As the target gene of miR-132-3p, UCP2 can regulate the differentiation of sheep precursor fat cells [30]. RNA-seq on the broiler pectoralis major muscles has revealed that ADIPOQ and CIDEC were key genes involved in lipid deposition [31]. SOCS2 can also act as a regulator of adipocyte size [32]. Transcriptome analyses examining IMF content in the LDM in heavy Iberian Pigs identified FOSB as a candidate gene and other regulatory factors [33]. We thus plan to focus on studying these genes in future analyses. Among the DELs, we found that MSTRG.8215.2, MSTRG.21380.1, MSTRG.21942.1, MSTRG.21719.1, MSTRG.8227.1, MSTRG.792.1 and MSTRG.19396.1 were highly expressed at Mth-12. We speculated that these novel lncRNAs might promote lipid deposition. Meanwhile, MSTRG.11343.4, MSTRG.13921.1, MSTRG.19788.2, MSTRG.2469.2, MSTRG.8912.2, MSTRG.2792.1 and MSTRG.21775.2 were highly expressed at Mth-2. We speculated that these novel lncRNAs might inhibit 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 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 top 20 GO terms are shown in Fig. 3a. The most significantly enriched GO terms were: toll-like receptor 3 signaling pathway (GO: 0034138), synapse pruning (GO: 0098883), reverse cholesterol transport (GO: 0043691), protein hetero dimerization activity (GO: 0046982), Nucleosome (GO: 0000786), extracellular region (GO: 0005576), DNA binding (GO: 0003677) and chromosome (GO: 0005694).
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), all of which were related to fat metabolism. Moreover, the top 20 signaling pathways are shown in Fig. 3b. Systemic lupus erythematosus (ko05322) showed the highest level of significance with 23 DEGs. 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 183 DELs regulated 218 DEMs, five lncRNAs of which acted on four 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). In GO annotation, these DEGs primarily played a role in biological processes. For example, positive regulation of phospholipid translocation (GO: 0061092), regulation of phospholipid catabolic process (GO: 0060696), Wnt signaling pathway, calcium modulating pathway (GO: 0007223), regulation of intracellular cholesterol transport (GO: 0032383), phospholipase C-activating dopamine receptor signaling pathway (GO: 0060158) and regulation of phospholipid biosynthetic process (GO: 0071071). The KEGG pathway enrichment analysis of target genes revealed a total of 214 annotated pathways (Additional file 9: Table S9B). Of these pathways, 28 were significantly enriched (P < 0.05) (Additional file 10: Table S10). Among them, eight pathways were related to lipid deposition and metabolism, including alpha-Linolenic acid metabolism (ko00592), FoxO signaling pathway (ko04068), Biosynthesis of unsaturated fatty acids (ko01040), Cell adhesion molecules (ko04514), Phosphonate and phosphinate metabolism (ko00440), Ether lipid metabolism (ko00565), Tight junction (ko04530) and Arachidonic acid metabolism (ko00590). Although some pathways were not significantly enriched, such as PPAR, Wnt, AMPK and mTOR signaling pathways, these pathways still played an important role in lipid deposition [31, 34-36]. Overall, the KEGG enrichment analysis of target genes revealed 16 critical pathways with 17 target genes (Table 3). We selected the DEMs and DELs with Pearson correlation coefficients ≥0.8 or ≤-0.8 (Additional file 11: Table S11). Among these genes, six were associated with lipid deposition, including SCD, ACAA2, FADS2, PLA2G4E, FZD4 and ULK1 [37-41], and were used to construct the lncRNA-mRNA co-expression network (Fig. 4). We speculate that MSTRG.792.1, MSTRG.8227, MSTRG.10679.1, MSTRG.21942.1, MSTRG.21380.1 and MSTRG.9270.1 might significantly contribute to the deposition of IMF. 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.