Identification and co-expression analysis of long noncoding RNAs and mRNAs involved in the deposition of intramuscular fat in Aohan fine-wool sheep

Abstract

Background: Intramuscular fat (IMF) content has become one of the most important indicators for measuring meat quality, and levels of IMF are affected by various genes. Long non-coding RNAs (lncRNAs) are widely expressed non-coding RNAs that play an important regulatory role in a variety of biological processes; however, research on the lncRNAs involved in sheep IMF deposition is still in its infancy. Aohan fine-wool sheep (AFWS), one of China's most important meat-hair, dual-purpose sheep breed, provides a great model for studying the role of lncRNAs in the regulation of IMF deposition. We identified lncRNAs by RNA sequencing in longissimus dorsi muscle (LDM) samples of sheep at two ages: 2 months (Mth-2) and 12 months (Mth-12). 

Results: We identified a total of 26,247 genes and 6,935 novel lncRNAs in LDM samples of sheep. Among these, 606 mRNAs and 408 lncRNAs were differentially expressed. We then compared the structural characteristics of lncRNAs and mRNAs. We obtained target genes of differentially expressed lncRNAs (DELs) and performed enrichment analyses using Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG). We found that target mRNAs were primarily enriched in lipid metabolism, lipid transport, regulation of primary metabolic processes and developmental pathways. Based on the results of important KEGG pathways, we obtained six candidate lncRNAs that potentially regulate lipid deposition and constructed an lncRNA-mRNA co-expression network that included MSTRG.792.1- SCD , MSTRG.8227.1- ACAA2 , MSTRG.10679.1- FADS2 , MSTRG.21942.1- PLA2G4E , MSTRG.21380.1- FZD4 and MSTRG.9270.1- ULK1 . We speculated that these candidate lncRNAs might play a role by regulating the expression of target genes. We randomly selected five mRNAs and five lncRNAs to verify the accuracy of the sequencing data by qRT-PCR.

Conclusions: Our study provided a list of the lncRNAs and mRNAs related to intramuscular lipid deposition in sheep and laid a foundation for future research on regulatory mechanisms.

Background

High-quality lamb meat is becoming increasingly popular as living standards improve and dietary patterns change. Currently, evaluations of the meat quality of livestock have revealed that the content of intramuscular fat (IMF) is lower in carcass fats, yet IMF has a critically important influence on the edibility and flavor of muscle meat [1]. Indeed, the quantity of IMF has become one of the most critical parameters of meat quality indicators, as it is considered to be positively related to meat quality and texture [2, 3]. When a certain amount of fat is deposited between the muscle bundles and muscle fibers, the marbled section of the meat has a high score, and the meat is fresh and juicy, which is often considered ideal [4,5]. The selective deposition of fat can improve production efficiency and play a key role in improving meat quality. This practice is also a major focus and challenge of modern livestock breeding [6]. Therefore, ensuring the appropriate deposition of IMF in lean meat can enhance the future quality of sheep meat. 

Studies have shown that intramuscular lipid deposition is affected by multiple genes and signaling pathways, such as the FAS, FAM134B, and HSL genes and the Wnt and AMPK signaling pathways [7-9]. Recently, long non-coding RNAs (lncRNAs) have received increased attention for their wide-ranging functions. LncRNAs refer to a class of non-coding RNAs longer than 200 nt in length [10]. Most lncRNAs have significant temporal and spatial expression specificity [11, 12] and have low sequence conservation among species [13-15]. LncRNAs can be divided into five types based on their positions relative to neighboring protein-coding genes: intronic lncRNAs, bidirectional lncRNAs, sense lncRNAs, intergenic lncRNAs and antisense lncRNAs [16].

LncRNAs can regulate various life activities of the body, including epigenetic regulation, transcriptional regulation and post-transcriptional regulation [17-19]. The most common regulation methods of lncRNAs include cis-regulation of the transcription of neighboring protein-coding genes and the trans-regulation of non-adjacent genes. In addition, lncRNAs can interact with miRNAs to affect the post-transcriptional translation of related mRNAs [20-22]. Studies have shown that lncRNAs can play direct or indirect roles in the process of lipid accumulation [23]. SRA (steroid receptor RNA activator) is one of the earliest discovered lncRNAs and plays an important role in fat metabolism. SRA can bind to peroxisome proliferator-activated receptor gamma (PPARγ) and enhance PPARγ activity, thereby promoting the differentiation of pre-adipocytes [24]. A study of the expression levels of lncRNAs in the IMF of Jinhua and Landrace pigs revealed a total of 119 differentially expressed lncRNAs (DELs), six of which were involved in fat deposition and lipid metabolism-related pathways [25]. Furthermore, an analysis of transcriptome data from IMF in Inner Mongolia goats revealed that 1,472 lncRNAs were involved in adipocyte growth regulation and morphological changes of adipocytes [26]. Another study has shown that lncRNAs can play a key regulatory role in fat deposition in sheep tails [27]. Overall, these findings demonstrate that lncRNAs can regulate lipid deposition through a variety of regulatory mechanisms. However, few studies have assessed the roles of lncRNAs in intramuscular lipid deposition in sheep. 

Aohan fine-wool sheep (AFWS) is an important meat-hair, dual-purpose sheep breed in China that grows rapidly early in development. The elimination of male lambs for fat lamb production can increase both hair and meat gains as well as improve the overall benefits provided by fine wool sheep [28].Exploring the developmental characteristics of IMF deposition and selecting candidate genes for AFWS provide references for future studies and applications in sheep breeding, improve the quality of mutton and accelerate the breeding process. The goal of our study was to systematically identify the profiles of differentially expressed mRNAs (DEMs) and DELs during intramuscular lipid deposition in sheep through high-throughput sequencing. We hoped that by studying the relationship between lncRNAs and lipid deposition, our findings would shed light on the mechanisms underlying selective muscle lipid deposition in sheep.

Results

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.  

Discussion

IMF content increased gradually with growth, as significant differences were detected between Mth-2 and Mth-12. These findings were consistent with a previous study showing that the IMF content of sheep increased from 0 to 6 months but remained stable thereafter until 12 months of age [42]. Furthermore, these findings are consistent with the characteristics of muscle growth and the development of experimental sheep. The sheep switched to a fattening phase after weaning at 2 months. The weight of sheep increased rapidly between the ages of 4 to 6 months, after which weight gain stabilized. IMF is an important feature contributing to meat quality. Therefore, we selected the LDM samples at Mth-2 (less lipid deposition) and Mth-12 (more lipid deposition) for RNA-seq to provide a robust test of gene expression differences.

Overall, we identified a total of 26,247 genes and 6,935 predicted novel lncRNAs in LDM samples of sheep by RNA-seq. Among these, 606 mRNAs (154 up-regulated and 452 down-regulated) and 408 lncRNAs (254 up-regulated and 154 down-regulated) were differentially expressed. To further characterize the mechanisms underlying DEGs, we performed GO and KEGG analysis of DEMs. In GO annotation, these DEGs primarily played a role in biological processes. These processes were closely related to fat formation and deposition, such as reverse cholesterol transport (GO: 0043691), positive regulation of cholesterol efflux (GO: 0010875) and positive regulation of cholesterol esterification (GO: 0010873) [43-45]. In addition, we found that many genes were enriched in biological processes, such as signal transmission, organ development, biosynthesis and cell proliferation, and these are also important processes in muscle development [46-48]. KEGG pathway analysis revealed that the DEMs were significantly enriched in the immune system, inflammatory response and infectious diseases pathways, demonstrating that signal transmission between adipocytes and immune cells can greatly affect the function of adipose tissue [49]. This result was consistent with the fact that inflammatory cell infiltration has been documented to commonly occur in adipose tissue and stimulate the activation of the immune defense system [50]. In addition, we also found that many pathways related to lipid metabolism (Cholesterol metabolism and Arachidonic acid metabolism) were significantly enriched for many DEMs, such as ADIPOQ, CIDEC, SOCS2, SCD and ACAA2, that have been reported to participate in lipid metabolism [51-55]. Many DEMs, such as C1QC, AGPAT4, SNCA and CAV3, were found to be related to lipid deposition for the first time [56-58]. Based on the GO and KEGG analysis, we obtained DEM expression profiles that affected the IMF deposition of sheep. However, further research is required to identify the underlying mechanisms.

In our study, 183 DELs participated in the regulation of mRNAs, five lncRNAs of which acted on four mRNAs through cis-regulation and 183 lncRNAs on 218 mRNAs through trans-regulation. The functionality of lncRNA is reflected through the study of their target genes [59]. We performed GO and KEGG enrichment analysis on these target genes. We focused on the GO terms related to lipid deposition. These included terms under biological processes, such as lipoprotein metabolic process (GO: 0042157), long-chain fatty acid biosynthetic process (GO: 0042759) and unsaturated fatty acid biosynthetic process (GO: 0006636), as well as molecular functions, such as Wnt-protein binding (GO: 0017147), phosphatidic acid binding (GO: 0070300) and phosphodiesterase I activity (GO:0004528). We found that the identified lncRNAs might be related to these GO terms; however, further research is required to identify the precise mechanisms.

The KEGG enrichment analysis revealed that target genes were significantly enriched in lipid metabolism pathways, such as alpha-linolenic acid metabolism (ko00592), phosphonate and phosphinate metabolism (ko00440) and ether lipid metabolism (ko00565). Cell adhesion molecules, tight junction and the FoxO signaling pathway also have important regulatory effects on lipid deposition [60-62]. Although some signaling pathways, such as the “Wnt signaling pathway,” “MAPK signaling pathway,” “AMPK signaling pathway” and “PPAR signaling pathway” were not significantly enriched in our study, they are critically important in the process of lipid deposition [63-66]. The above analysis and KEGG pathways led to the identification of six novel lncRNAs: MSTRG.21380.1, MSTRG.21942.1, MSTRG.792.1, MSTRG.8227.1, MSTRG.10679.1 and MSTRG.9270.1. These lncRNAs might play important roles in lipid deposition and deserve further study.

To facilitate future studies of the mechanisms underlying lncRNAs, we constructed the lncRNA-mRNA co-expression network based on important KEGG pathways. The network contained 127 lncRNAs and eight mRNAs. We obtained six pairs of lncRNAs-mRNAs, including MSTRG.792.1-SCD, MSTRG.8227.1-ACAA2, MSTRG.10679.1-FADS2, MSTRG.21942.1-PLA2G4E, MSTRG.21380.1-FZD4 and MSTRG.9270.1-ULK1.Target genes of these lncRNAs have been reported to be involved in lipid deposition. For example, a study examining gene expression differences in metabolism and function between intramuscular and subcutaneous adipocytes in cattle found that SCD was highly expressed in adipocytes and closely associated with fat formation [53]. ACAA2 can play an important role in fatty acid metabolism by promoting the differentiation of sheep precursor adipocytes into adipocytes [38, 67]. The inhibition of FADS2 in essential fatty acid deficiency induces hepatic lipid accumulation via impairment of very low-density lipoprotein secretion [68]. PLA2G4E was significantly enriched in alpha-linolenic acid metabolism pathways (linolenic acid is closely related to the flavor of lamb) and is associated with lipid metabolism or cholesterol metabolism [40, 69]. FZD4 is highly expressed during fat production [41], and ULK1 participates in lipid metabolism [62, 63]. These results provide information for future studies examining how lncRNAs regulate IMF deposition in sheep. The specific regulatory mechanisms require further study and testing.

Conclusions

Our study systematically identified mRNA and lncRNA expression profiles during intramuscular lipid deposition in sheep. We obtained a total of 606 DEMs and 408 DELs and identified some important lncRNAs related to lipid deposition through GO and KEGG enrichment analysis. In addition, co-expression network analysis of lncRNAs and mRNAs involving 127 lncRNAs and eight mRNAs was conducted based on significant KEGG pathways. Six pairs of lncRNA-mRNA, including MSTRG.792.1-SCD, MSTRG.8227.1-ACAA2, MSTRG.10679.1-FADS2, MSTRG.21942.1-PLA2G4E, MSTRG.21380.1-FZD4 and MSTRG.9270.1-ULK1, were selected for further research. Our study provided a list of the lncRNAs and mRNAs related to intramuscular lipid deposition and laid a foundation for future research on the regulatory mechanisms of lncRNA on sheep muscle lipid deposition.

Methods

Sample preparation

All experimental sheep came from the AFWS Stud Farm (Chifeng, Inner Mongolia, China). All sheep were fed under the same feeding and management conditions. A total of 12 healthy AFWS rams (3 individuals for each stage) at 2, 4, 6 and 12 months of age were killed for the sample collection. AFWS rams were obtained from 12 ewes of a similar age and weight that were in estrus simultaneously and were artificially inseminated from the same ram. The 12 healthy AFWS rams were placed in a closed chamber and anesthetized with sodium pentobarbital at a dose of 25 mg/kg by intravenous injection. Rams were anesthetized and eventually sacrificed in the enclosed chamber by having it filled with 20% carbon dioxide every minute until the gas concentration had reached 80%. Experimental animal handling procedures were performed following published protocols [72, 73]. Samples of the LDM and BFM were collected, placed in RNAase-free Eppendorf tubes and stored immediately in liquid nitrogen. Likewise, the LDM and BFM samples (150 g for each muscle) were quickly collected and stored at -20 °C for the determination of IMF content using Soxhlet petroleum-ether extraction. 

Determination of IMF content

After removing the white intermuscular fat from the muscle samples, the samples were minced thoroughly with a meat grinder, loaded into glassware and dried at 105° C until completely dry. Samples were then weighed after crushing (marked as z), wrapped with quantitative filter paper and baked at 105 ℃ until samples were dry and their weight did not change. Samples were then weighed in paper bags after drying (marked as x). The dried paper bag was then placed in the Soxhlet extraction bottle, and the ether reflux device was used to reflux the sample at 65 ℃ until the drops are transparent. The paper bag was then placed in a fume hood to fully volatilize the ether, followed by drying at 105 ℃ until the weight did not change(marked as y). The measurement was repeated three times for each sample. The following formula was used to calculate the IMF content: IMF content (%) = (x–y) / z × 100%.

RNA extraction and quality assessment

Total RNA of the longissimus dorsi muscles was extracted using Trizol reagent (Invitrogen, Carlsbad, CA, USA) per the manufacturer's instructions. RNA purity was measured at an OD 260/280 with a NanoDrop ND-2000 instrument (Thermo Fisher Scientific, MA, USA). RNA integrity (RIN) was evaluated by 1% agarose gel electrophoresis and Aglian 2100 (Agilent, Santa Clara, CA, USA). RNA samples with OD₂₆₀/OD₂₈₀ ratio greater than 1.8 and RIN value greater than 7.5 were selected for sequencing.

Library preparation and sequencing

First-strand complementary DNA (cDNA) was synthesized using random hexamer primers and M-MuLV reverse transcriptase (RNase H-) [74], with rRNA-depleted RNA used as a template. Second-strand cDNA was then synthesized with dNTPs, DNA polymerase I and RNase H. Next, T4 DNA polymerase and Klenow DNA polymerase were used to repair and modify the ends to add an A base and ligate the sequencing adapter. The cDNA products were then purified using AMPure XP beads (Beckman Coulter, Brea, CA, USA). Finally, uracil DNA glycosylase (NEB, Ipswich, MA, USA) was used to degrade the U-containing chain to remove second-strand cDNA. The purified first-strand cDNA was enriched by PCR to obtain a cDNA library. The quality of the libraries was assessed using an Agilent 2100 Bioanalyzer, and sequencing was performed using paired-end sequencing (2*150 bp) with the Illumina HiSeq 4000 platform (LC Sciences, Houston, TX, USA).

Mapping of reads and transcriptome assembly

Cutadapt was used to remove the reads that were contaminated by adapters, low-quality bases and undetermined bases. The clean reads were mapped to Ovis aries Ensembl release 96 (ftp://ftp.ensembl.org/pub/release-96/fasta/ovis_aries/dna/) using HISAT2 [75]. Mapped reads of each sample were assembled by StringTie [76]. The R package “edge” was used for difference statistics and visual drawing. Gffcompare was used to combine all transcripts from samples to reconstruct a comprehensive transcriptome. StringTie was used to determine the expression level for all transcripts by calculating FPKM (FPKM = [total exon fragments/mapped reads (millions)×exon length(kb)]). The R package edgeR was used to select the differentially expressed transcripts that satisfied the condition of log₂ (fold change) ≥ 1 or log₂ (fold change) ≤ -1 and P < 0.05.

Identification of lncRNAs

Known transcripts and transcripts less than 200 bp in length were removed from the data set. CPC (Coding Potential Calculator) and CNCI (Coding-Non-Coding Index) were then used to screen lncRNAs [77, 78]. When the CPC software score was less than 0.5 and the CNCI software score was less than 0, a transcript was considered a novel lncRNA.

Enrichment analysis of differentially expressed mRNAs

We used the Gene Ontology database (http://www.geneontology.org) and the Kyoto Encyclopedia of Genes and Genomes (http://www.kegg.jp/kegg) to annotate DEGs. The genes were mapped to GO terms and KEGG pathways based on annotation information and then the hypergeometric test was performed. The clustering map was drawn by the R package. GO terms and KEGG pathways were defined as significantly enriched when P < 0.05.

Prediction of lncRNA target genes

Based on the cis-and trans-regulation mechanisms of lncRNA, we identified the protein-coding genes (100-kb upstream and downstream) located on the same chromosome as the lncRNA that was a target for cis-regulation. RIsearch was used to predict the free energy of lncRNA-mRNA gene combinations on different chromosomes; combinations of lncRNA and mRNA with free energies below -11 kcal/mol were identified as trans target genes of lncRNA [79]. The results of the cis and trans-regulation were used to calculate Pearson correlations between lncRNA and mRNA expression. Cytoscape was used to plot the co-expression network.

Verification of sequencing data

We randomly selected five lncRNAs and five mRNAs to validate their expression using SYBR Green PCR Master Mix (Takara, Dalian, China). Primer 5 was used to design primers for the candidate genes. The sequences of the primers used are listed in Additional file 12 (Table S12). A 20-μL PCR mixture consisted of 10 μL SYBR® Premix Ex Taq II (2×), 0.5 μL forward primer (10 μM/L), 0.5 μL reverse primer (10 μM/L), 1 μL cDNA and 8 μL ddH₂O. The PCR parameters were as follows: 95 °C for 30 s; 40 cycles of 95 °C for 5 s; 60 °C for 30 s; 72 °C for 30 s; and 72 °C for 5 min. Three replicates were conducted for each sample. The 2^-ΔΔCt method was used to quantify relative expression levels [80].

Statistical analysis

Data on IMF content were expressed as means± standard deviation. One-way analyses of variance in SPSS 17.0 were used to analyze experimental results. Independent sample t-tests were used to compare the IMF content of muscles at the same age. All the data from the qRT-PCR were obtained using atleast three independent replicates. Differences were deemed statistically significant if p-values are smaller than 0.05.

Abbreviations

AFWS: Aohan fine wool sheep; IMF: Intramuscular fat; LncRNAs: long noncoding RNAs; Mth-2, Mth-4, Mth-6, Mth-12: Rams at 2, 4, 6, 12 months; LDM: Longissimus dorsi muscle; BFM: Biceps femoris muscle; RNA-seq: RNA sequencing; DEGs: Differentially expressed genes; DEMs: Differentially expressed mRNAs; DELs: Differentially expressed lncRNAs ; ORF:open reading frame; GO: Gene ontology; KEGG: Kyoto encyclopedia of genes and genomes; qRT-PCR: Quantitative real-time PCR; CPC: Coding Potential Calculator; CNCI: Coding-Non-Coding Index; RIN :RNA integrity

Declarations

Ethics approval and consent to participate

All the experimental operation have conforming to the Guidelines for Experimental Animals of the Ministry of Science and Technology (Beijing, China) and were approved by the Experimental Animal Ethics Committee of Qingdao Agricultural University. The managements of laboratory animal has in keeping with 《Laboratory Animal-Requirements of Environment and Housing Facilities》(GB 14925-2001). The written informed consentto participate was obtained from the AFWS Stud Farm in Inner Mongolia Autonomous Region. All efforts were made to minimize suffering.

Consent for publication

Not applicable.

Availability of data and materials

Additional data can be found in supplementary files. The RNA-Seq data is being uploaded and can be uploaded to the database before publication.

Competing interests

The authors declare that they have no competing interests.

Funding

The work was supported by the Earmarked Fund for Modern China Wool & Cashmere Technology Research System (CARS-39); the Project of Shandong Province Agricultural Variety Program (2019LZGC012); the National Natural Science Foundation of China (31402047）, and the Project of Shandong Province Higher Educational Science and Technology Program(J18KA136). The funding body had no role in the design of the study, collection, analysis, interpretation of data, decision to publish or writing the manuscript.

Authors’ contributions

NL and JNH designed this study; FHH, RRZ, LRL, LLL, and QL participated in sample collection; FHH, RRZ, LLL and QL performed qRT-PCR validation; FHH, JL, NL, LRL and JNH analyzed the RNA-Seq data; FHH wrote the manuscript with contribution from JL, NL, LRL and JNH. All authors reviewed and approved the final manuscript.

Acknowledgments

The authors thank the staff of the Laboratory of Animal Genetics, Breeding and Reproduction, Qingdao Agricultural University for providing the research laboratories.

Authors' information

¹College of Animal Science and Technology, Qingdao Agricultural University, Qingdao 266109, China. 

²Qufu Animal Husbandry and Veterinary Technical Service Center, Qufu 273100, China.

³China Animal Health and Epidemiology Center, Qingdao 266032, China.

References

1. Anderson F, Pannier L, Pethick DW, Gardner GE. Intramuscular fat in lamb muscle and the impact of selection for improved carcass lean meat yield. Animal, 2015;9(6): 1081-1090.
2. Dodson MV, Allen RE, Du M, Bergen WG, Velleman SG, Poulos SP, Fernyhough-Culver M, Wheeler MB, Duckett SK, Young MRI, Voy BH, Jiang Z, Hausman GJ. INVITED REVIEW: Evolution of meat animal growth research during the past 50 years: Adipose and muscle stem cells. Journal of Animal Science, 2015;93(2): 457-481.
3. Mohrhauser DA, Taylor AR, Underwood KR, Pritchard RH, Wertz-Lutz AE, Blair AD. The influence of maternal energy status during midgestation on beef offspring carcass characteristics and meat quality. Journal of Animal Science, 2015;93(2): 786-793.
4. Urban T, Mikolasova R, Kuciel J, et al.A study of associations of the HFABP genotypes with fat and meat production of pigs [J].Journal of Applied Genetics, 2002;43(4):505-510.
5. Wang YH, Byrne KA, Reverter A, et al. Transcriptional profiling of skeletal muscle tissue from two breeds of cattle[J].Mammalian Genome, 2005;16(3):201-210.
6. Wood JD, Enser M, Fisher AV, et al. Fat deposition, fatty acid composition and meat quality: A review. Meat Sci. 2008;78(4):343–358.
7. Luan ZJ, Liu KD, He JN, Cheng M, Qu XX, Liu N. The FAM134B, PPARγ, HSL and FAS gene expression patterns and their association with intramuscular fat contentin in sheep [J]. Acta Veterinaria et Zootechnica Sinica. 2016;47(12):2379-2389.
8. Okamura M, Kudo H, Wakabayashi K, et al. COUP-TFII acts downstream of Wnt/beta-catenin signal to silence PPARgamma gene expression and repress adipogenesis. Proc Natl Acad Sci U S A. 2009;106(14):5819-5824.
9. Xu Z, Yu S, Hsu CH, Eguchi J, Rosen ED. The orphan nuclear receptor chicken ovalbumin upstream promoter-transcription factor II is a critical regulator of adipogenesis. Proc Natl Acad Sci U S A. 2008;105(7):2421-2426.
10. Wang J, Hua L, Chen J, et al. Identification and characterization of long non-coding RNAs in subcutaneous adipose tissue from castrated and intact full-sib pair Huainan male pigs. BMC Genomics. 2017;18(1):542.
11. Cabili MN, Trapnell C, Goff L, et al. Integrative annotation of human large intergenicnoncoding RNAs reveals global properties and specific subclasses[J]. Genes Dev, 2011;25,1915-1927.
12. Zhu J, Fu H, Wu Y, et al. Function of lncRNAs and approaches to lncRNA-protein interactions [J]. Science China Life Sciences, 2013;56,876-885.
13. Huang W, Long N, Khatib H. Genome-wide identification and initialCharacterization of bovinelong noncoding RNAs from EST data [J]. Anim Genet, 2012;43(6):674-682.
14. Koufariotis LT, Chen YP, Chamberlain A, Vander Jagt C, Hayes BJ. A catalogue of novel bovine long noncoding RNA across 18 tissues. PLoS One. 2015;10(10):e0141225.
15. Weikard R, Hadlich F, Kuehn C. Identification of noveltranscripts and noncodingRNAsinbovine skin by deep next generation sequencing [J]. BMC Genomics, 2013;14:789.
16. Ponting CP, Oliver PL, Reik W. Evolutionand functions of long noncoding RNAs [J]. Cell, 2009;136(4):629-641.
17. Zhang C, Peng G. Non-coding RNAs: an emerging player in DNA damageresponse. Mutat Res Rev Mutat Res. 2015;763:202-211.
18. Li M, Xu M, Wen S, Bai B, Chen R, Liu L. One novel long noncoding RNA lnc10 in Drosophila. J Genet Genomics. 2014;41(2):79–82.
19. Wapinski O, Chang HY. Long noncoding RNAs and human disease. Trends Cell Biol. 2011;21(6):354–361.
20. Ørom UA, Derrien T, Beringer M, et al. Long noncoding RNAs with enhancer-like function in human cells. Cell. 2010;143(1):46–58.
21. Pruneski JA, Hainer SJ, Petrov KO, Martens JA. The Paf1 complex represses SER3 transcription in Saccharomyces cerevisiae by facilitating intergenic transcription-dependent nucleosome occupancy of the SER3 promoter. Eukaryot Cell. 2011;10(10):1283–1294.
22. Li M, Sun X, Cai H, et al. Long non-coding RNA ADNCR suppresses adipogenic differentiation by targeting miR-204. Biochim Biophys Acta. 2016;1859(7):871–882.
23. Chen Z. Progress and prospects of long noncoding RNAs in lipid homeostasis. Mol Metab. 2015;5(3):164–170.
24. Xu B, Gerin I, Miao H, et al. Multiple roles for the non-coding RNA SRA in regulation of adipogenesis and insulin sensitivity. PLOS One. 2010;5(12):e14199.
25. Miao Z, Wang S, Zhang J, et al. Identification and comparison of long non-conding RNA in Jinhua and Landrace pigs. Biochem Biophys Res Commun. 2018;506(3):765–771.
26. Zheng ZQ, Du C, Fu SY, et al. Identification and characterization analysis of lncRNA from RNA-seq data of intramuscular adipocytes in goats[J].Acta Veterinaria et Zoote chnicaSinica, 2014;45(12):1924-1931.
27. Ma L, Zhang M, Jin Y, et al. Comparative Transcriptome Profiling of mRNA and lncRNA Related to Tail Adipose Tissues of Sheep. Front Genet. 2018;9:365.
28. Liu K, Ge S, Luo H, Yue D, Yan L. Effects of dietary vitamin E on muscle vitaminE and fatty acid content in Aohan fine-wool sheep. J Anim Sci Biotechnol. 2013;4(1):21.
29. Liu GL, Jiang SW, Xiong YZ, Zheng R, Qu YC. Yi Chuan Xue Bao. 2003;30(12):1107‐1112.
30. Shi Tao. MiR-132-3p targeting UCP2 regulates the differentiation of sheep precursor fat cells [D]. Shanxi :Shanxi Agricultural University, 2018.
31. Liu L, Liu X, Cui H, Liu R, Zhao G, Wen J. Transcriptional insights into key genesand pathways controlling muscle lipid metabolism in broiler chickens. BMC Genomics. 2019;20(1):863.
32. Val CH, de Oliveira MC, Lacerda DR, et al. SOCS2 modulates adipose tissue inflammation and expansion in mice. J Nutr Biochem. 2020;76:108304.
33. Muñoz M, García-Casco JM, Caraballo C, et al. Identification of Candidate Genes andRegulatory Factors Underlying Intramuscular Fat Content Through Longissimus Dorsi Transcriptome Analyses in Heavy Iberian Pigs. Front Genet. 2018;9:608.
34. Du M, Yin J, Zhu MJ. Cellular signaling pathways regulating the initial stage of adipogenesis and marbling of skeletal muscle. Meat Sci. 2010;86(1):103‐109.
35. Kwak HJ, Jeong MY, Um JY, Park J. β -Lapachone Regulates Obesity through Modulating Thermogenesis in Brown Adipose Tissue and Adipocytes: Role of AMPK Signaling Pathway. Am J Chin Med. 2019;47(4):803‐822.
36. Cai R, Tang G, Zhang Q, et al. A Novel lnc-RNA, Named lnc-ORA, Is Identified byRNA-Seq Analysis, and Its Knockdown Inhibits Adipogenesis by Regulating the PI3K/AKT/mTOR Signaling Pathway. Cells. 2019;8(5):477.
37. Stachowiak M, Flisikowski K. Analysis of allele-specific expression of seven candidategenes involved in lipid metabolism in pig skeletal muscle and fat tissues reveals allelic imbalance of ACACA, LEP, SCD, and TNF. J Appl Genet. 2019;60(1):97-101.
38. Zhang Y, Wang Y, Wang X, et al. Acetyl-coenzyme A acyltransferase 2 promote the differentiation of sheep precursor adipocytes into adipocytes [published online ahead of print, 2018 Nov 28]. J Cell Biochem. 2018;10.1002/jcb.28080.
39. Kothapalli KSD, Park HG, Guo X, et al. A novel FADS2 isoform identified in human milk fat globule suppresses FADS2 mediated Δ6-desaturation of omega-3 fatty acids. Prostaglandins Leukot Essent Fatty Acids. 2018;138:52-59.
40. Song YB, An YR, Kim SJ, et al. Lipid metabolic effect of Korean red ginseng extract in mice fed on a high-fat diet. J Sci Food Agric. 2012;92(2):388‐396.
41. Tian XW. The role of Wnt signaling molecules in fat metabolism and its protein properties [D]. Jinan: Shandong Normal University, 2006.
42. Hao CL, Li QF, Qiao Y, et al. Correlation between H-FABP and PPARγ gene expression levels and IMF content in muscle tissue of Hu sheep [J]. Chinese Agricultural Science, 2008; 41 (11): 3776-3783.
43. Zhang M, Li F, Sun JW, et al. LncRNA IMFNCR Promotes Intramuscular Adipocyte Differentiation by Sponging miR-128-3p and miR-27b-3p. Front Genet. 2019;10:42.
44. Pauli A, Valen E, Lin MF, Garber M, Vastenhouw NL, Levin JZ, et al. Systematic identification of long noncoding RNAs expressed during zebrafish embryogenesis. Genome Res. 2012;22:577–91.
45. Lee HJ, Jang M, Kim H, et al. Comparative Transcriptome Analysis of Adipose Tissues Reveals that ECM-Receptor Interaction Is Involved in the Depot-Specific Adipogenesis in Cattle. PLoS One. 2013;8(6):e66267.
46. Qiao Y. Study on developmental changes of gene expression related to IMF depositionin different parts of Hu Lambs [D]. Nanjing: Nanjing Agricultural University, 2007.
47. Boone C, Mourot J, Gregoire F, et al.The adipose conversion process:Regulation by extracellular and intracellular factors[J]. Reproduction Nutrition Development. 2000;40(4):325-358.
48. Da Costa AS, Pires VM, Fontes CM, Mestre Prates JA. Expression of genes controlling fat deposition in two genetically diverse beef cattle breeds fed high or low silage diets. BMC Vet Res. 2013;9:118.
49. Engin AB. Adipocyte-Macrophage Cross-Talk in Obesity. Adv Exp Med Biol. 2017;960:327–343.
50. Morimoto K, Tanaka T, Nitta Y, Ohnishi K, Kawashima H, Nakatani T. NEDD9 crucially regulates TGF-β-triggered epithelial-mesenchymal transition and cell invasion in prostate cancer cells: involvement in cancer progressiveness. Prostate. 2014;74(8):901–910.
51. Zhao NQ, Li XY, Wang L, et al. Palmitate induces fat accumulation by activating C/EBPβ-mediated G0S2 expression in HepG2 cells. WorldJ Gastroenterol.2017;23(43):7705-7715.
52. Huang JZ, Huang LM, Zeng QJ, et al. Distribution and quantitative analysis of CIDEa and CIDEc in broiler chickens: accounting for differential fat deposition between strains. Br Poult Sci. 2018;59(2):173-179.
53. Hudson NJ, Reverter A, Griffiths WJ, et al. Gene expression identifies metabolic and functional differences between intramuscular and subcutaneous adipocytes in cattle. BMC Genomics. 2020;21(1):77.
54. Li B, Qiao L, An L, et al. Transcriptome analysis of adipose tissues from two fat-tailed sheep breeds reveals key genes involved in fat deposition. BMC Genomics. 2018;19(1):338.
55. Peng M, Li S, He Q, Zhao J, Li L, Ma H. Proteomics reveals changes in hepatic proteins during chicken embryonic development: an alternative model to study human obesity. BMC Genomics. 2018;19(1):29.
56. Dias A, Santos D, Coelho T, et al. C1QA and C1QC modify age-at-onset in familial amyloid polyneuropathy patients. Ann Clin Transl Neurol. 2019;6(4):748‐754.
57. Minakaki G, Menges S, Kittel A, et al. Autophagy inhibition promotes SNCA/alpha-synuclein release and transfer via extracellular vesicles with a hybrid autophagosome-exosome-like phenotype. Autophagy. 2018;14(1):98‐119.
58. Feng XJ, Ma LX, Jiao C, et al. Nerve injury elevates functional Cav3.2 channels in superficial spinal dorsal horn. Mol Pain. 2019;15:1744806919836569.
59. Sheik Mohamed J, Gaughwin PM, Lim B, Robson P, Lipovich L. Conserved long noncoding RNAs transcriptionally regulated by Oct4 and Nanog modulate pluripotency in mouse embryonic stem cells. RNA. 2010;16(2):324–337.
60. Bosanská L, Michalský D, Lacinová Z, et al. The influence of obesity anddifferent fatdepots on adipose tissue gene expression and protein levels of celladhesion molecules. Physiol Res. 2010;59(1):79–88.
61. Hussain M , Umair Ijaz M , Ahmad MI , et al. Meat proteins in a high-fat diet havea substantial impact on intestinal barriers through mucus layer and tight junction protein suppression in C57BL/6J mice. Food Funct. 2019;10(10):6903‐6914.
62. Brown KM, Sharma S, Baker E, Hawkins W, van der Merwe M, Puppa MJ. Delta-6-desaturase (FADS2) inhibition and omega-3 fatty acids in skeletal muscle protein turnover. Biochem Biophys Rep. 2019;18:100622.
63. Ross SE, Hemati N, Longo KA, et al. Inhibition of adipogenesis by Wnt signaling. Science. 2000;289(5481):950–953.
64. Wu W, Zhang J, Zhao C, Sun Y, Pang W, Yang G. CTRP6 Regulates Porcine Adipocyte Proliferation and Differentiation by the AdipoR1/MAPK Signaling Pathway. J Agric Food Chem. 2017;65(27):5512–5522.
65. Yao C, Pang D, Lu C, et al. Data Mining and Validation of AMPK Pathway as a Novel Candidate Role Affecting Intramuscular Fat Content in Pigs. Animals (Basel). 2019;9(4):137.
66. Li B, Qiao L, An L, et al. Transcriptome analysis of adipose tissues from two fat-tailed sheep breeds reveals key genes involved in fat deposition. BMC Genomics. 2018;19(1):338.
67. Symeou S, Tzamaloukas O, Banos G, Miltiadou D. ACAA2and FASNpolymorphisms affect the fatty acid profile of Chios sheep milk. J Dairy Res. 2020;87(1):23‐26.
68. Hayashi Y, Shimamura A, Ishikawa T, Fujiwara Y, Ichi I. FADS2 inhibition in essential fatty acid deficiency induces hepatic lipid accumulation via impairment of very low-density lipoprotein (VLDL) secretion. Biochem Biophys Res Commun. 2018;496(2):549–555.
69. Zhang QE, Ao CJ. Factors affecting fatty acid composition of mutton [J]. Animal Husbandry and Veterinary Medicine. 2007;39(3):19-21.
70. Kim S, Kim N, Park S, et al. Tanycytic TSPO inhibition induces lipophagy to regulate lipid metabolism and improve energy balance. Autophagy. 2019;1‐21.
71. Wang X, Wu R, Liu Y, et al. m6A mRNA methylation controls autophagy and adipogenesis by targeting Atg5 and Atg7 . Autophagy. 2019;1‐15.
72. Wen S, Yan MR, Cao C, Yang SB, Zhang SB. Analysis of application of different doses of pentobarbital sodium to anesthesia of experimental dogs. J Acta Academiae Medicine Zunyi. 2009;32(05):465–466.
73. Raj AB, Leach MC, Morton DB. Carbon dioxide for euthanasia of laboratory animals. Comp Med. 2004;54(5):470–471.
74. Adiconis X, Borges-Rivera D, Satija R, et al. Comparative analysis of RNA sequencing methods for degraded or low-input samples. Nat Methods. 2013;10(7):623–629.
75. Kim T, Seo HD, Hennighausen L, Lee D, Kang K. Octopus-toolkit: a workflow to automate mining of public epigenomic and transcriptomic next-generation sequencing data. Nucleic Acids Res. 2018;46(9):e53.
76. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–295.
77. Kong L, Zhang Y, Ye ZQ, et al. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007;35(Web Server issue):W345–W349.
78. Sun L, Luo H, Bu D, et al. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013;41(17):e166.
79. Xie C, Mao X, Huang J, et al. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 2011;39(Web Server issue):W316–W322.
80. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods. 2001;25(4):402–408.

Tables

Table 1  IMF content of sheep (%)

IMF content in different parts of muscles among sheep with the same age. ** indicates that means were highly significantly different (P < 0.01); * indicates significant differences (P < 0.05); different lowercase letters indicate that means differ significantly (P < 0.05) between the same muscles groups of sheep of different ages; different capital letters indicate that means were highly significantly different (P < 0.01) between the same muscles groups of sheep of different ages.

Table 2 Statistical data derived from RNA Sequencing

Table 3 Critical mRNAs based on the KEGG pathways related to lipid deposition

Additional Files

Additionnal file 1: Table S1. Summary of protein-coding genes identified in the libraries. (XLSX 7373 kb)

Additionnal file 2: Table S2. (A). Summary of lncRNAs identified in the libraries.

(B). Type statistics of lncRNAs identified in the libraries. Class code: “x” represents the antisense lncRNAs, “o” represents the sense lncRNAs, “j” represents the bidirectional lncRNAs, “i” represents the intronic lncRNAs, and “u” represents the intergenic lncRNAs. (XLSX 5253 kb)

Additional file 3: Table S3. Sequence information of all expressed lncRNAs found in the study. (FA 157696 kb)

Additional file 4: Table S4. The summary of differentially expressed protein-coding genes.  (XLSX 212 kb)

Additional file 5: Table S5. The summary of differentially expressed lncRNAs. (XLSX 77 kb)

Additional file 6: Table S6. (A). GO enrichment analysis of the differentially expressed mRNAs. S gene number: the number of significant differentially expressed mRNAs which match to a GO term; TS gene number: the number of significant differentially expressed mRNAs which have GO annotations; B gene number: the number of detected mRNAs which match to a GO term; TB gene number: the number of all detected mRNAs which have GO annotations. (B). KEGG enrichment analysis (P < 0.05) of the differentially expressed mRNAs.  (XLSX 77 kb)

Additional file 7: Table S7. Significantly enriched KEGG pathways (P < 0.05) of the differentially expressed mRNAs. (XLSX 10 kb) S gene number: the number of significant differentially expressed mRNAs that match to a pathway; TS gene number: the number of significant differentially expressed mRNAs that have pathway annotations; B gene number: the number of detected mRNAs that match to a pathway term; TB gene number: the number of all detected mRNAs that have pathway annotations.

Additional file 8: Table S8. The differentially expressed target mRNAs of the differentially expressed lncRNA in either trans- or cis- regulatory roles. (XLSX 5919 kb)

Additional file 9: Table S9. (A). GO enrichment analysis of differentially expressed target mRNAs of differentially expressed lncRNAs in the study. (B) .KEGG enrichment analysis of differentially expressed target mRNAs of differentially expressed lncRNAs in the study. (XLSX 138 kb)

Additional file 10: Table S10. The significantly enriched KEGG pathway analysis (P < 0.05) of target genes of differentially expressed lncRNAs. (XLSX 11 kb)

Additional file 11: Table S11. Critical mRNAs and their lncRNAs that were related to the lipid deposition. (XLSX 102 kb)

Additional file 12: Table S12. Primers used in the qRT-PCR analysis. (XLSX 99 kb)