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

DOI: https://doi.org/10.21203/rs.3.rs-18941/v1

Abstract

Background: Intramuscular fat (IMF) content has become one of the most important indicators for measuring meat quality, and its level is 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 lncRNAs involved in sheep intramuscular fat deposition is still in its infancy. Aohan fine-wool sheep (AFWS), China's representative meat-hair, dual-purpose sheep breed, provides a great model for studying the role of lncRNAs in the regulation of intramuscular fat deposition. We identified lncRNAs by RNA sequencing in sheep longissimus dorsi muscle(LDM) samples at two ages: 2 months (Mth-2) and 12 months (Mth-12).

Results: We identified a total of 26,247 genes and 6,935 predicted 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 targeted genes of differentially expressed lncRNAs and performed an enrichment analysis using Gene Ontology(GO) and the Kyoto Encyclopedia of Genes and Genomes(KEGG). We found these targeted mRNAs were primarily enriched in lipid metabolism, lipid transport, regulation of primary metabolic processes and developmental pathways, such as alpha-linolenic acid metabolism, biosynthesis of unsaturated fatty acids, phosphonate and phosphinate metabolism and cell proliferation. Based on the results of this enrichment analysis, we obtained candidate lncRNAs that potentially regulate lipid deposition and constructed a lncRNA-mRNA co-expression network. We speculated that these lncRNAs have important regulatory roles in intramuscular fat deposition. We randomly selected five mRNAs and five lncRNAs to verify the accuracy of sequencing results by qRT-PCR.

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

Background

With improvements in living standards and changes in patterns of consumption, high-quality lamb meat is becoming increasingly popular. 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 in 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 difficulty 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 FAS, FAM134B, HSL genes and Wnt, AMPK signaling pathways [79]. Recently, long non-coding RNAs (lncRNAs) have received increased attention because of their wide range of functions. LncRNAs refer to a class of non-coding RNAs longer than 200 nt in length that do not encode proteins [10]. Most lncRNAs have significant temporal and spatial expression specificity [11, 12] and have low sequence conservation among species [1315]. According to their relative positions with neighboring protein-coding genes, they can be divided into five types: sense lncRNAs, antisense lncRNAs, intronic lncRNAs, intergenic lncRNAs and bidirectional lncRNAs [16].

LncRNAs can regulate various life activities of the body in a variety of ways, including epigenetic regulation, transcriptional regulation and post-transcriptional regulation [1719]. 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, affecting post-transcriptional translation of related mRNAs [2022]. Studies have shown that lncRNAs can play a direct or indirect role 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 preadipocytes [24]. A study of the expression levels of lncRNAs in the intramuscular fat of Jinhua and Landrace pigs revealed a total of 119 differentially expressed lncRNAs, six of which were involved in fat deposition and lipid metabolism-related pathways [25]. Furthermore, an analysis of transcriptome data from intramuscular fat 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 a representative meat-hair, dual-purpose sheep breed in China, whose lambs grow quickly early in their development. The elimination of male lambs for fat lamb production can increase both hair and meat gains and increase the overall benefits of fine wool sheep [28]. Exploring the developmental characteristics of IMF deposition and selecting candidate genes for AFWS has many potential benefits, such as providing a reference for future studies and applications in sheep breeding, improving the quality of mutton and accelerating the breeding process of sheep. The goal of our study was to systematically identify the profiles of differentially expressed mRNAs and lncRNAs 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 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).  

Discussion

Muscle tissue is the most important tissue for IMF deposition [29], so we initially studied IMF deposition in two types of muscle tissue: the LDM and the BFM. The sheep used in the experiment switched to the 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. Therefore, sheep at 2, 4, 6 and 12 months of age were selected for experimental study of IMF deposition to enhance the specificity of the results. Studies have shown that the deposition rate of IMF had clear differences; specifically, the IMF content of the LDM has been documented to be higher than that of the psoas and BFM from 0 to 6 months [30]. The same conclusion was made in our study. We found that the lipid deposition of AFWS initially increased and then stabilized. This finding was consistent with a previous study showing that the IMF content of sheep from 0 to 6 months increased but thereafter remained stable until 12 months of age [31]. Therefore, to explore the regulatory mechanisms underlying intramuscular fat genes, we selected the LDM samples at Mth-2 (less lipid deposition) and Mth-12 (more lipid deposition) for RNA sequencing to provide a robust test of gene expression differences.

Changes in the regulation and expression of genes must accompany changes in the deposition of lipids as sheep move from the breeding stage to the maturity stage [10, 27, 28]. Many genes related to lipid deposition were found to be differentially expressed inour study. In previous studies, RNA sequencing on the broiler pectoralis major muscles revealed that ADIPOQ and CIDEC were key genes in lipid deposition [32]. SOCS2 can also act as a regulator of adipocyte size [33]. Another 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 [34]. The aforementioned findings indicate that these genes might also play an important role in lipid deposition in sheep.

GO and KEGG enrichment analysis of significantly different genes revealed that many of the significantly enriched GO terms and pathways were related to fat metabolism and development, such as reverse cholesterol transport, phospholipid efflux and positive regulation of the fatty acid metabolic process. These processes were all closely related to the formation and deposition of fat [3537]. There is no doubt that these genes have been the focus of our attention. In addition, we found that many genes were enriched in biological processes, such as signal transmission, organ development, biosynthesis and cell proliferation, and these were also important processes in muscle development [38]. The pathway results showed that signal transmission can occur between cells or between cells, especially the "session" between adipocytes and immune cells, which seemed to greatly affect the function of adipose tissue [39]. Many genes that were significantly enriched in KEGG pathways were associated with the immune system, inflammatory response and infectious diseases. This observation is consistent with the fact that inflammatory cell infiltration has been documented to commonly occur in adipose tissue and thereby stimulate the activation of the immune defense system [40]. Given that the experimental sample included the LDM, IMF deposition was likely to have been accompanied by skeletal muscle growth and development, such as locomotory behavior, positive regulation of cell motility and actin-myosin filament sliding. Both adipocytes and myoblasts are regulated by genes during the cell proliferation process [4143]. For example, genes related to corpus callosum development and cell maturation deserve increased attention. In sum, analysis based on GO and KEGG showed that the process of lipid deposition was accompanied by the continuous regulation of metabolic, immune and developmental functions at the cellular and organ levels.

Overall, 408 differentially expressed lncRNAs (254 up-regulated and 154 down-regulated lncRNAs) were obtained by RNA sequencing. In recent years, RNA sequencing technology has been widely used to explore new genes, providing a basis for more efficient mechanistic research [44]. A previous study showed that ribosome-removed high-throughput sequencing technology can be used to identify differentially expressed lncRNAs before and after bovine adipocyte differentiation. Specifically, 16 lncRNAs were found to be differentially expressed before and after adipocyte differentiation, and lncRNA-ADNCR inhibited bovine adipocyte differentiation through competitive binding to miR-204 [22]. Therefore, the differentially expressed lncRNAs reported in our study might also potentially be important regulators of intramuscular lipid deposition.

An increasing number of studies have shown that lncRNA can regulate the expression levels of genes in a variety of ways [45]. In our study, 183 differentially expressed lncRNAs can participate in the regulation of mRNAs, 5 lncRNAs of which acted on 4 mRNAs through cis-regulation and 183 lncRNAs on 218 mRNAs through trans-regulation. In cis-regulation, the target gene HSPA6 was reported to be differentially expressed in gestational diabetes mellitus adipose tissue samples [46] and affects the proliferation and differentiation of various cells in the body [47, 48]. Therefore, HSPA6 may be involved in the regulation of sheep adipocytes. We also identified many target genes related to lipid deposition, such as ACAA2 [49], SCD [50], FADS2 [51] and PLA2G4E [52]. We further performed GO and KEGG enrichment analysis on these target genes. The results showed that target genes were significantly enriched in lipid metabolism pathways, such as alpha-linolenic acid metabolism, phosphonate and phosphinate metabolism and ether lipid metabolism. The target gene PLA2G4E was significantly enriched in alpha-linolenic acid metabolism pathways. Previous studies have characterized the correlation between fatty acids and the flavor of mutton in more than 20 species of sheep. Generally, the flavor of mutton is significantly related to the presence of oleic acid (18: 1) and linolenic acid (18: 3) (r = 0.33, r = 0.33) [53], indicating that the genes in this pathway are important for determining the flavor of lamb during fat formation.

In addition, target genes were significantly enriched in pathways that are known to be involved in lipid synthesis and metabolism, such as cell adhesion molecules, tight junction and the FoxO signaling pathway. The cell adhesion molecules pathway mediates cell-to-cell contact and cell-to-extracellular matrix interactions [54]. To date, the mRNAs enriched in this pathway have not been reported for adipocytes. For example, research on CACNA1D has been mostly conducted on mutants that lead to central nervous system disease [55]. However, cell adhesion has a great effect on the cell proliferation and differentiation process. Whether these candidate genes play an important role in lipid deposition remains to be studied [54]. In addition, although various pathways, such as the “Wnt signaling pathway,” “MAPK signaling pathway,” “AMPK signaling pathway” and “PPAR signaling pathway” were not significantly enriched, they were related to fat development and deposition [5659]. Therefore, the differentially expressed lncRNAs involved in the aforementioned pathways may play an important role in the process of lipid deposition and deserve further study.

Conclusions

In sum, our study systematically identified mRNA and lncRNA expression profiles during intramuscular lipid deposition in sheep via RNA sequencing. We identified some differentially expressed mRNAs and lncRNAs related to lipid deposition through GO and KEGG enrichment analysis. Co-expression analysis of lncRNAs and mRNAs suggests that some lncRNAs are associated with lipid metabolism, such as those involved in Wnt signaling pathways, PPAR signaling pathways, biosynthesis of unsaturated fatty acids and glycerophospholipid metabolism. Our results not only provide insight into the identification of genes that play a role in IMF deposition but also lay a foundation for future research on the regulatory mechanisms underlying 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 selected for study and they were obtained from 12 ewes of the similar age and weight whose estrus were synchronized, and artificially insemination from a 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. under rams anesthesia, the enclosed chamber was filled with 20 % carbon dioxide every minute. When the gas concentration had reached 80 %, the rams died Experimental animal handling procedures were performed following published protocols [60, 61]. 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 determination of IMF content using Soxhlet petroleum-ether extraction. 

Library preparation and sequencing

First-strand complementary DNA was synthesized using random hexamer primers and M-MuLV reverse transcriptase (RNase H-)[62], while rRNA-depleted RNA was used as a template. Second-strand complementary DNA 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 cDNA library was obtained by enriching purified first-strand cDNA. The quality of the libraries was identified using an Agilent 2100 Bioanalyzer, and they were sequenced using paired-end sequencing (2*150 bp) with Illumina HiSeq 4000 technology (LC Sciences, Houston, TX, USA).

Mapping of reads and transcriptome assembly

Raw data were generally provided in FASTQ format. First, the clean data were obtained using Cutadapt (v 1.10) to remove unqualified data. Reads were removed if they (1) were contaminated by sequencing adapters; (2) contained poly-N; or (3) were low quality with Q-values less than 10. The clean reads were mapped to ovis_aries ensembl release 96(ftp://ftp.ensembl.org/pub/release-96/fasta/ovis_aries/dna/) using TopHat [63]. The number of reads following comparisons and their distribution and density on chromosomes were counted according to the genome annotation file. The mapped reads of each sample were assembled and quantified by StringTie [64]. The R package “edge” was used for difference statistics and visual drawing.

Identification of lncRNAs

Known transcripts with protein-coding capabilities and transcripts less than 200 bp in length were removed from the data set. CPC (Coding Potential Calculator) [65] and CNCI (Coding-Non-Coding Index) [66] were then used to screen lncRNAs. When the CPC software score was less than 0.5 and the CNCI software score was less than 0, a transcript was considered a lncRNA that could not encode proteins.

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. We performed GO analyses using GO seqplatform, while KEGG pathways were analyzed using KOBAS software [67]. The software TM4 (http://www. tm4.org/mev.html) was used to screen for DEGs and to perform cluster analysis. GO terms and KEGG pathways were defined as significantly enriched when P < 0.05.

Prediction of lncRNA target genes

Based on cis-and trans-regulation mechanisms of lncRNA, we identified the protein-coding gene (100 kb upstream and downstream) located on the same chromosome as the lncRNA that was a target for cis-regulation. Genes that were located on a different chromosome from lncRNAs and with free energy used to form secondary structures less than -11 kcal/mol were identified as trans targets for lncRNAs [68]. GO terms and KEGG pathways were considered significantly enriched when P < 0.05.

Verification of sequencing results

To verify the accuracy of the sequencing results, we randomly selected 5 lncRNAs and 5 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® PremixExTaq II (2×), 0.5 μL forward primer (10 μM/L), 0.5 μL reverse primer (10 μM/L), 1 μL cDNA and 8 μL dd H2O. 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 set for each sample to be tested. The 2-ΔΔCt method was used to calculate and analyze results [69].

Statistical analysis

Data on IMF content were expressed as means. SPSS 17.0 was used to analyze the experimental results with one-way analysis of variance. Independent sample t-tests were used to compare the IMF content of muscles at the same age. The results of RNA sequencing were analyzed using SPSS 17.0 and R software (2016.09.29). Co-expression networks of lncRNAs-mRNAs were constructed using Cytoscape software (version 3.4.0; The Cytoscape Consortium, San Diego, CA, USA).

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;   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

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

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

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

3China 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.
  1. Wapinski O, Chang HY. Long noncoding RNAs and human disease. Trends Cell Biol. 2011;21(6):354–361.
  2. Ørom UA, Derrien T, Beringer M, et al. Long noncoding RNAs with enhancer-like function in human cells. Cell. 2010;143(1):46–58.
  3. 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.
  4. 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.
  5. Chen Z. Progress and prospects of long noncoding RNAs in lipid homeostasis. Mol Metab. 2015;5(3):164–170.
  6. 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.
  7. 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.
  8. 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.
  9. 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.
  10. 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.
  11. 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.
  12. Qiao Y. Study on developmental changes of gene expression related to IMF depositionin different parts of Hu Lambs [D]. Nanjing: Nanjing Agricultural University, 2007.
  13. 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.
  14. 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.
  15. Val CH, de Oliveira MC, Lacerda DR, et al. SOCS2 modulates adipose tissue inflammation and expansion in mice. J Nutr Biochem. 2020;76:108304.
  16. 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.
  17. 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.
  18. 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.
  19. 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.
  20. Wang GQ, Wang Y, Xiong Y, et al. Sirt1 AS lncRNA interacts with its mRNA to inhibit muscle formation by attenuating function of miR-34a. Sci Rep. 2016;6:21865.
  21. Engin AB. Adipocyte-Macrophage Cross-Talk in Obesity. Adv Exp Med Biol. 2017;960:327–343.
  22. 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.
  23. Cui HX, Zheng MQ, Liu RR, Zhao GP, Chen JL, Wen J. Liver dominant expression of fatty acid synthase (FAS) gene in two chicken breeds during intramuscular-fat development.  Mol Biol Rep. 2012;39(4):3479-3484.
  24. 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.
  25. Ren H, Li L, Su H, et al. Histological and transcriptome-wide level characteristics of fetal myofiber hyperplasia during the second half of gestation in Texel and Ujumqin sheep. BMC Genomics. 2011;12:411.
  26. Stark R, Grzelak M, Hadfield J. RNA sequencing: the teenage years. Nat Rev Genet. 2019;20(11):631–656.
  27. 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.
  28. Deng X, Yang Y, Sun H, Qi W, Duan Y, Qian Y. Analysis of whole genome-wide methylation and gene expression profiles in visceral omental adipose tissue of pregnancies with gestational diabetes mellitus. J Chin Med Assoc. 2018;81(7):623–630.
  29. Shin SS, Song JH, Hwang B, et al. HSPA6 augments garlic extract-induced inhibition of proliferation, migration, and invasion of bladder cancer EJ cells; Implication for cell cycle dysregulation, signaling pathway alteration, and transcription factor-associated MMP-9 regulation. PLoS One. 2017;12(2):e0171860.
  30. McCullagh KJ, Cooney R, O'Brien T. Endothelial nitric oxide synthase induces heat shock protein HSPA6 (HSP70B') in human arterial smooth muscle cells. Nitric Oxide. 2016;52:41–48.
  31. Zhang Y, Wang Y, Wang X, et al. Acetyl-coenzyme A acyltransferase 2 promote the differentiation of sheep precursor adipocytes into adipocytes. J Cell Biochem. 2018;10.1002/jcb.28080.
  32. Ackerman D, Tumanov S, Qiu B, et al. Triglycerides Promote Lipid Homeostasis during Hypoxic Stress by Balancing Fatty Acid Saturation. Cell Rep. 2018;24(10):2596–2605.e5
  33. 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.
  34. 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.
  35. Zhang QE, Ao CJ. Factors affecting fatty acid composition of mutton [J]. Animal Husbandry and Veterinary Medicine 2007;39(3):19-21.
  36. 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.
  37. Pinggera A, Striessnig J. Cav 1.3 (CACNA1D) L-type Ca2+ channel dysfunction in CNS disorders. J Physiol. 2016;594(20):5839–5849.
  38. Ross SE, Hemati N, Longo KA, et al. Inhibition of adipogenesis by Wnt signaling. Science. 2000;289(5481):950–953.
  39. 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.
  40. 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.
  41. 59. 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.
  42. 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.
  43. Raj AB, Leach MC, Morton DB. Carbon dioxide for euthanasia of laboratory animals. Comp Med. 2004;54(5):470–471.
  44. 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.
  45. 63. Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009;25(9):1105–1111.
  46. 64. 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.
  47. 65. 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.
  48. 66. 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.
  49. 67. Wenzel A, Akbasli E, Gorodkin J. RIsearch: fast RNA-RNA interaction search using a simplified nearest-neighbor energy model. Bioinformatics. 2012;28(21):2738–2746.
  50. 68. 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.
  51. 69. 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 (%)

Age

Longissimus dorsi IMF(%)

Biceps femoris  IMF(%)

Mth-2

2.202±0.006Aa**

2.012±0.058Aa

Mth-4

4.566±0.178Bb**

3.390±0.149Bb

Mth-6

10.685±0.690Cc**

7.925±0.378Cc

Mth-12

11.163±0.878Cc*

8.867±0.188Cc

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



Table 2  Statistical data derived from RNA Sequencing

 

Mth-2-1

Mth-2-2

Mth-2-3

Mth-12-1

Mth-12-2

Mth-12-3

Raw reads

90397074

89528496

88813860

90619672

95989952

89577658

(13.59G)

(14.40G)

(13.44G)

(13.56G)

(13.43G)

(13.32G)

Valid reads

85098620

84675502

81473556

85252798

91240846

85201776

(12.79G)

(13.69G)

(12.78G)

(12.76G)

(12.70G)

(12.22G)

Valid Ratio%

94.14

94.58

91.74

94.08

95.05

95.11

Q20%

99.98

99.98

99.97

99.98

99.98

99.98

Q30%

98.02

97.98

97.83

98.03

98.15

98.17

GC content%

49.5

51

53

49

49

49

Mapped reads

76713303

74525753

69966338

77129925

83264100

77731307

(90.15%)

(88.01%)

(85.88%)

(90.47%)

(91.26%)

(91.23%)

Expressed genes

24537

24229

24295

24605

24563

24077

Unique lncRNAs

6495

6492

6507

6532

6541

6425

 



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

Critical pathways

critical mRNAs 

FoxO signaling pathway

HOMER2,AGAP2,FBXO32,BNIP3,ARAF

Biosynthesis of unsaturated fatty acids

FADS2,SCD

Phosphonate and phosphinate metabolism

CHPT1

Cell adhesion molecules (CAMs)

CDH5,SIGLEC1

Ether lipid metabolism

CHPT1,PLA2G4E

PPAR signaling pathway

FADS2,SCD

Fatty acid elongation

ACAA2

Glycerophospholipid metabolism

CHPT1,PLA2G4E

mTOR signaling pathway

ULK1,FZD4

Fatty acid degradation

ACAA2

NOD-like receptor signaling pathway

MFN2

Wnt signaling pathway

FZD4

Cell cycle

CDKN2C,ANAPC13

AMPK signaling pathway

ULK1,SCD

Insulin signaling pathway

ARAF,PPP1CB

Apoptosis

TUBA4A,CAPN1

Jak-STAT signaling pathway

FHL1,IL10RA

cAMP signaling pathway

DRD1,PPP1CB

MAPK signaling pathway

RPS6KA5,ARAF,PLA2G4E

Phospholipase D signaling pathway

PLA2G4E

Cytokine-cytokine receptor interaction

IL10RA

alpha-Linolenic acid metabolism

PLA2G4E,FADS2

Tight junction

CTTN,TUBA4A

Arachidonic acid metabolism

PLA2G4E

Linoleic acid metabolism

PLA2G4E

Retinol metabolism

CYP1A1

Inositol phosphate metabolism

TPI1

Phosphatidylinositol signaling system

PPIP5K1

 

 

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” represent the antisense lncRNAs, “o” represent the sense lncRNAs,“j” represent the intronic lncRNAs, “i” represent the intron lncRNAs, and “u” represent 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 KEGG pathway analysis (P < 0.05) of targeted genes of differentially expressed lncRNAs. (XLSX 11 kb)

Additional file 11: Table S11. LncRNAs and its potential target genes that are involved in some key pathways and GO terms related to the lipid deposition. (XLSX 550 kb)

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