Integrated transcriptome analysis reveals roles of long non- coding RNAs (lncRNAs) in caprine skeletal muscle mass and meat quality

doi:10.21203/rs.3.rs-1960027/v1

Download PDF

Research Article

Integrated transcriptome analysis reveals roles of long non- coding RNAs (lncRNAs) in caprine skeletal muscle mass and meat quality

https://doi.org/10.21203/rs.3.rs-1960027/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background

Long non-coding RNAs (lncRNAs) play important roles in growth and development of skeletal muscle. However, there was limited information in goats. In this study, expression profiles of lncRNAs in Longissimus dorsi muscle from Liaoning cashmere (LC) goats and Ziwuling black (ZB) goats with divergent meat yield and meat quality were compared using RNA-sequencing. Based on our previous microRNAs (miRNAs) and mRNAs profiles obtained from the same tissues, the target genes and binding miRNAs of differentially expressed lncRNAs were obtained. Subsequently, lncRNA-mRNA interaction networks and a ceRNA network of lncRNA-miRNA-mRNA were constructed.

Results

A total of 136 differentially expressed lncRNAs were identified between the two breeds. 15 cis target genes and 143 trans target genes were found for differentially expressed lncRNAs, and they were enriched in muscle contraction, muscle system process, muscle cell differentiation, and p53 signaling pathway. A total of 69 lncRNA-trans target gene pairs were constructed, with close relationship with muscle development, intramuscular fat deposition and meat tenderness. A total of 16 lncRNA-miRNA-mRNA ceRNA pairs were identified, of which some reportedly associated with skeletal muscle development and fat deposition were found.

Conclusion

The study identified some crucial lncRNAs related to muscle development, intramuscular fat deposition and meat tenderness, which will provide an improved understanding of the roles of lncRNAs in caprine meat yield and meat quality.

long non-coding RNA

meat yield

meat quality

ceRNA

goat

Domestic goats (Capra hircus) are one of the most essential farm animals economically. Goat meat is increasingly favored by consumers for higher protein, lower fat as well as unique flavor and palatability when compared to meat produced from other domestic animals [1]. It is well known that as other animals, caprine meat yield and quality are regulated by both mRNAs and non-coding RNAs including long non-coding RNAs (lncRNAs), so the identification of the RNAs that regulate skeletal muscle offers an opportunity to improve caprine meat production performance.

The lncRNAs are a class of novel non-coding RNAs with greater than 200 nucleotides in length. They are primarily transcribed by RNA polymerase II and widely distributed in mammalian cells [2]. Although the vast majority of lncRNAs expressed at lower levels compared to protein-coding genes [3], they were found to play crucial roles in regulation of mRNA expression at transcriptional, post-transcriptional or epigenetic levels. For example, some lncRNAs transcribed from genome sequence containing cis-regulatory DNA elements can regulate expression of their neighboring genes in cis by affecting their transcription [4]. The lncRNAs can also trans-regulate the expression of the target genes by interacting with trans-acting factors [4]. In addition, some lncRNAs act as competing endogenous RNAs (ceRNAs) to sequester microRNAs (miRNAs), eventually leading to increase the expression of the target gene by miRNAs at post-transcriptional level [5]. In this context, lncRNAs are considered to widely involve in embryonic development, organ morphogenesis, as well as cell cycles, including differentiation, proliferation, and apoptosis [6–7].

It has been confirmed that lncRNAs played important roles in the growth and development of skeletal muscle. For example, lncR-125b accelerated differentiation of caprine muscle satellite cells by increasing expression of insulin growth factor 2 (IGF2) derived from a ceRNA for miR-125b [5]. LncRNA-Six1 promoted chicken muscle growth by cis-regulating expression of the target gene SIX homeobox 1 (Six1) [8]. Additionally, lncMAAT was found to regulate muscle atrophy by negatively regulating the expression of miR-29b [9]. Meanwhile, the expression profiles of lncRNAs in skeletal muscle tissues were also investigated and lncRNAs were found to be differentially expressed in samples with different genetic backgrounds. However, these studies in livestock have mainly been focused on sheep [10–12], cattle [13–15], chicken [16–18], and pigs [19–21].

Despite there being some studies reporting lncRNAs expression profiles of skeletal muscle tissues in goats, these studies were all performed with muscle samples collected from different developmental periods. For example, a total of 547 lncRNAs were differentially expressed in skeletal muscles of Anhui white goats among five fetal stages and two kid stages, and they were involved in signaling pathways closely associated with muscle development, including structure formation, p53 signaling pathway, and MAPK signaling pathway [22]. Additionally, 577 and 648 differentially expressed lncRNAs were found in skeletal muscle tissues between embryonic and postnatal stages of Jianzhou big-eared goats and Dazu black goats, respectively [23–24]. However, there have been no reports on comparison of lncRNAs expression profiles in skeletal muscle tissues between different goat breeds.

In this study, Liaoning cashmere (LC) goats and Ziwuling black (ZB) goats were selected for the investigation. The two breeds are all dual-purpose breed for meat and cashmere fiber, and have significant difference in meat production performance. Briefly, LC goats have higher carcass weight, net meat weight, muscle fiber size, and intramuscular fat content, but have poorer meat tenderness compared to ZB goats (P < 0.05) [25]. RNA-Seq analysis of the same Longissimus dorsi muscle tissues as those used in the study revealed that differentially expressed genes, miRNAs and circular RNAs (circRNAs) were responsible for these phenotypic differences [26–28]. However, the biological mechanism by lncRNAs regulate meat production performance difference is still unclear. Accordingly, the expression profiles of lncRNAs in Longissimus dorsi muscle tissues between LC and ZB goats were compared using RNA sequencing (RNA-Seq), and differentially expressed lncRNAs were then screened. We also analyzed functional enrichment of the target genes, and constructed lncRNA-mRNA interaction networks and a ceRNA network of lncRNA-miRNA-mRNA, with the aim of uncovering the possible function of lncRNAs in muscle growth and development of goats.

Identification and characterization of lncRNAs in caprine skeletal muscle

For the ten Longissimus dorsi muscle tissue samples (five LC goats and five ZB goats), their clean reads and mapped results to the Caprine Genome Assembly ARS1 have been described in our previous study [26]. In this study, a total of 2,302 lncRNAs were identified, including 1,767 known caprine lncRNAs and 535 novel lncRNAs (Supplementary File 1). Of all the lncRNAs identified, 1,945 lncRNAs were co-expressed in the two goat breeds, while 178 and 179 lncRNAs were specifically expressed in LC and ZB goats, respectively (Fig. 1A). Of the five lncRNA types classified according to their location relative to protein-coding genes, intergenic lncRNA was the most common with a proportion of 54.4%, followed by antisense lncRNAs (15.3%) and bidirectional lncRNAs (15.3%). Sense lncRNAs (3.4%) and intronic lncRNAs (2.3%) were the least. In addition, 9.3% lncRNAs were defined as other type of lncRNAs that were not classified into the five types (Fig. 1B).

In our previous study, the mRNA expression profiles of the same Longissimus dorsi muscle tissues as those used in the study has been reported [26]. In this study, characteristics in Longissimus dorsi muscle between lncRNAs and mRNAs identified were accordingly compared. The average transcript length of known lncRNAs identified in this study was 1,668 bp, which was less than novel lncRNAs and mRNAs with an average length of 5,065 and 3,674 bp, respectively (Fig. 1C). The average expression level of lncRNAs was far less than mRNA transcripts (Fig. 1D). The average exon number of known and novel lncRNAs was 3.8 and 2.6, respectively, which was much less than average 12.7 exons for mRNAs (Fig. 1E). Additionally, the average length of open reading frame (ORF) of known and novel lncRNAs was shorter than mRNAs (Fig. 1F). This was consistent with the results that lncRNAs had much lower average coding potential score than mRNAs (Fig. 1G).

Screening and validation of differentially expressed lncRNAs

A total of 136 differentially expressed lncRNAs were identified in comparing Longissimus dorsi muscle tissue from LC goats with Longissimus dorsi muscle tissue from ZB goats. These included 62 up-regulated lncRNAs and 74 down-regulated lncRNAs in the Longissimus dorsi muscle of LC goats (Supplementary File 2).

The results from RT-qPCR showed that the relative expression levels of the 12 differentially expressed lncRNAs selected were in accordance with those obtained from RNA-Seq (Fig. 2). This confirmed the reliability and repeatability of the RNA-Seq results.

Prediction and function annotation of the target genes of differentially expressed lncRNAs

To predict the target genes of differentially expressed lncRNAs, a total of 391 up-regulated genes and 222 down-regulated genes in the same Longissimus dorsi muscle as those used in the study, were re-identified in LC goats when comparing to ZB goats (Supplementary File 3). In this study, 15 cis target genes were predicted for 18 differentially expressed lncRNAs, resulting in 22 lncRNA-cis target gene pairs (Supplementary File 4). Additionally, 143 trans target genes were predicted for 58 differentially expressed lncRNAs, resulting in 203 lncRNA-trans target gene pairs (Supplementary File 4). However, there were no target genes predicted for the remaining 60 differentially expressed lncRNAs.

The target genes in cis and trans were significantly enriched in 155 biological process (GO-BP) terms, 19 cellular components (GO-CC) terms, and 19 molecular function (GO-MF) terms (Supplementary File 5). It was noteworthy that some significant GO terms directly related to the growth and development of skeletal muscle were found (Fig. 3A), including muscle contraction (P = 0.002), regulation of muscle contraction (P = 0.004), muscle system process (P = 0.014), regulation of muscle system process (P = 0.022), muscle cell differentiation (P = 0.040), and striated muscle cell differentiation (P = 0.041). In addition, connective tissue development process related to meat tenderness was also significantly enriched by several target genes (P = 0.010) (Fig. 3A).

KEGG pathway analysis showed that the target genes were significantly enriched in 19 signaling pathways (Supplementary File 5) and the top fifteen significant pathways with the lowest P value are shown in Fig. 3B. Systemic lupus erythematosus was the most enriched pathway (P = 2.685E-05), followed by alcoholism (P = 2.623E-04) and allograft rejection (P = 0.007). The p53 signaling pathway (P = 0.023) was noteworthy as it played dual roles in the growth and development of muscle and adipose tissues (Fig. 3B).

Construction of lncRNA-mRNA interaction networks

Given unknown roles for 15 cis target genes in meat yield and quality, their interaction network with lncRNAs was not constructed. Of 203 lncRNA-trans target gene pairs described above, based on the roles of the target genes in meat yield and quality, a total of 69 lncRNA-trans target gene pairs were further selected to construct lncRNA-mRNA interaction networks, including 38 pairs related to muscle development (Fig. 4A), 18 pairs associated with intramuscular fat deposition (Fig. 4B) and 13 pairs related to meat tenderness (Fig. 4C). For example, the target genes CDK1, GREM1, MEGF10, and DNER were involved in muscle cell differentiation process (Supplementary File 5). The target genes RORC, KLK7 and RYR3 were relate to adipogenesis [29–31], while the target genes GREM1 and LOXL2 were involved in connective tissue development process affecting meat tenderness [32].

Construction of a lncRNA-miRNA-mRNA ceRNA network

In addition to the roles in regulating expression of the target genes, lncRNAs can also function as miRNA sponges to increase expression level of the target genes by miRNAs. In this study, a total of 16 lncRNA-miRNA-mRNA ceRNA pairs were identified (P < 0.05) (Fig. 5; Supplementary File 6). It was noteworthy that several ceRNA networks related to the growth and development of skeletal muscle were found including XR_001917125.1-miR-34-3p-NCAM2, XR_001917125.1-miR-460-5p-NCAM2 and XR_001918832.1-miR-200c-CRHBP. A ceRNA network related to fat deposition was also found including XR_001918832.1-miR-200b-CRHBP.

Meat production performance is a kind of complex economic traits that include multiple important traits influencing meat yield and meat quality. For example, carcass weight and net meat weight directly reflect meat yield. Meat tenderness is regarded as important palatability trait that significantly affects consumer acceptance for meat, while intramuscular fat influences sensory quality of meat, including flavor, juiciness and tenderness [33]. An increasing number of evidences suggest that meat yield and quality traits can be directly regulated by some non-coding RNAs. Given that lncRNAs account for 80% of the total number of non-coding RNAs, the role of lncRNAs in skeletal muscle are worthy of further investigation.

In this study, a total of 2,302 lncRNAs were identified in Longissimus dorsi muscle of goats, and this was less than those reported in muscle tissue of Jianzhou big-eared goats with an average of 2,739 lncRNAs identified [23]. This may partly reflect breed-specific expression pattens of lncRNAs. Additionally, the number of lncRNAs identified in this study was more than those for a study in Longissimus dorsi muscle of sheep [12], but less than those reported in cattle by Yan et al. [14]. The differences may be related to species-specific expression of lncRNAs. The characteristics of lncRNAs identified in this study are in concordance with previous studies in skeletal muscle tissue. For example, intergenic lncRNA was found to be the most common in skeletal muscle of Dazu black goats [24], rabbit [34] and cattle [13]. Our observation that known lncRNAs had shorter transcript length than mRNAs, is consistent with findings in Anhui white goats [22], cattle [13] and donkey [35]. As expected, lncRNAs found in this study had much lower expression levels and coding potential score, shorter ORF length and fewer exon number when comparing to mRNA transcripts. The phenomenon has also been reported in skeletal muscle of Anhui white goats [22] and Jianzhou big-eared goats [23].

The regulation of the target genes in expression is one of the main functions of lncRNA. In lncRNA-mRNA interaction networks constructed in the study (Fig. 4A, B), up-regulated MSTRG.262.1, XR_001917125.1 and XR_001917605.1 in LC goats may regulate both myogenesis and adipogenesis as their some trans target genes GREM1, MEGF10, and DNER were involved in muscle cell differentiation and striated muscle cell differentiation process (P < 0.05), while other trans target genes CDK1 and RRM2 were enriched in p53 signaling pathway (P < 0.05) that played dual roles in muscle cell differentiation [36] and adipogenesis [37]. In addition, MSTRG.262.1 would target MYOZ2 and ANKRD2, while XR_001917125.1 would target MYOZ2. The genes MYOZ2 and ANKRD2 were found to be up-regulated in the same Longissimus dorsi muscle tissue of LC goats as those used in the study (Supplementary File 3) and were also reported to positively regulate muscle cell differentiation [38–39]. It was therefore inferred that the three up-regulated lncRNAs in LC goats may be responsible for the meat production and intramuscular fat content differences between LC goats and ZB goats.

It was noteworthy that MSTRG.262.1 and XR_001917125.1 also appeared in the interaction network related to meat tenderness. The lncRNAs were co-expressed with myosin light chain 6B (MYL6B) and collagen alpha-1(III) chain (LOC102176755) (Fig. 4C). MYL6B was a crucial muscle structure component and its nucleotide sequence variations were found to affect beef tenderness [40]. LOC102176755 encodes a type of collagen that has a strong negative effect on meat tenderness [41]. The GO analysis result further supported the effect of XR_001917125.1 on meat tenderness as its other two up-regulated trans target genes (GREM1 and LOXL2) in Longissimus dorsi muscle of LC goats (Supplementary File 3) were involved in connective tissue development process (P < 0.05) (Supplementary File 5). Connective tissue is mainly composed of collagen, which its content was highly positively correlated with the shear force value of beef (r = 0.95) [41]. In this context, the two lncRNAs may also regulate meat tenderness difference between the two breeds.

Two up-regulated lncRNAs (XR_001917386.1, XR_001918614.1) and one down-regulated lncRNA (XR_001297059.2) in LC goats caught our attention as their trans target genes (SMPX and RYR3) were enriched in muscle contraction (P < 0.01) and muscle system process (P < 0.05) (Supplementary File 5). Of the two target genes, up-regulated SMPX (Supplementary File 3) would be trans-regulated by up-regulated XR_001917386.1 and XR_001918614.1 in LC goats. SMPX was also known as small muscle protein, in that it was a positive regulator of muscle fiber size [42]. We therefore speculated that the up-regulation of SMPX in skeletal muscle of LC goats may be partly caused by trans-regulation of XR_001917386.1 and XR_001918614.1, eventually resulting in a higher muscle fiber size in LC goats. The lncRNA XR_001297059.2 was down-regulated in LC goats with higher intramuscular fat content (Supplementary File 2). This was not surprising as its trans target gene RYR3 played a negative role in adipogenesis [43] and was also down-regulated in LC goats (Supplementary File 3).

In interaction network related to muscle development, NR5A2 was common trans target gene of seven differentially expressed lncRNAs (Fig. 4A). NR5A2 has been reported to regulate muscle morphogenesis and glucose metabolism in muscle cells [44–45]. Similarly, RORC would be collectively trans-regulated by six differentially expressed lncRNAs in interaction network related to intramuscular fat deposition (Fig. 4B) and meat tenderness (Fig. 4C). RORC was enriched in connective tissue development process related to meat tenderness (P < 0.05) (Supplementary File 5). The variations in RORC were also associated with intramuscular fat and marbling score of cattle [29–30]. These suggest that these multiple lncRNAs may play roles in skeletal muscle development and intramuscular fat deposition in goats by collectively trans-regulating NR5A2 or RORC.

Other differentially expressed lncRNAs involved in the lncRNA-mRNA interaction networks may also regulate the phenotype differences of the two goat breeds in meat yield and meat quality. For example, up-regulated MSTRG.12645.1 in Longissimus dorsi muscle of LC goats compared to ZB goats would be co-expressed with KLK7, which was also identified as an up-regulated gene in LC goats (Supplementary File 3). KLK7 led to an increase of body fat mass in mice [31]. The up-regulated XR_001918823.1 in LC goats would target SLC7A8 that was up-regulated in Longissimus dorsi muscle of LC goats compared to ZB goats (Supplementary File 3). SLC7A8 was reported to promote muscle growth by weaking proteolysis [46]. However, 15 cis target genes predicted for differentially expressed lncRNAs obtained in this study were not found in the three lncRNA-mRNA interaction networks. This suggest that these differentially expressed lncRNAs screened in this study play roles in skeletal muscle development and meat quality mainly through a trans-regulatory mechanism.

The lncRNA can also act as a sponge for miRNA to relieve the repression of the target mRNA by miRNA, with an accompanying increase in expression level of mRNAs in mammalian cells [5, 8]. In the study, up-regulated XR_001918832.1 in LC goats was predicted to be a sponge of miR-200b and miR-200c to regulate the expression of CRHBP (Fig. 5). The miR-200b was reported to suppress proliferation of C2C12 myoblast [47] and differentiation of ovine preadipocytes [48], while miR-200c inhibited differentiation of C2C12 myoblast [49]. It was therefore speculated that up-regulated XR_001918832.1 in LC goats may sponge-absorb miR-200b/miR-200c to positively regulate caprine muscle development and intramuscular fat. Similarly, up-regulated XR_001917125.1 in LC goats would increase the expression of NCAM2 by relieving the repression effect of miR-34-3p and miR-460-5p (Fig. 5). Given that NCAM2 promoted myogenesis of mice [50], up-regulation of XR_001917125.1 would result in a higher skeletal muscle mass in LC goats compared to ZB goats. These suggest that ceRNA mechanisms may partly explain the phenotype differences in skeletal muscle mass and intramuscular fat between LC goats and ZB goats. However, the lncRNA-miRNA-mRNA ceRNA networks predicted in the study need to be further corroborated.

In summary, we identified 136 differentially expressed lncRNAs in muscle tissues between LC goats and ZB goats. The target genes of the differentially expressed lncRNAs were related to phenotype differences in meat production, intramuscular fat deposition and meat tenderness between the two caprine breeds. Several lncRNAs may responsible for the differences via ceRNA mechanism, including XR_001918832.1-miR-200b/miR-200c-CRHBP and XR_001917125.1-miR-34-3p/miR-460-5p-NCAM2.

Ethics statement

All animal procedures in this study were approved by Animal Experiment Ethics Committee of Gansu Agricultural University with an approval number of GSAU-ETH-AST-2021-028.

Experimental animals and sample collection

Five healthy nine-month-old LC goats and five healthy nine-month-old ZB goats raised under the same environmental conditions and nutrition program, were selected from the Yongfeng Goat Breeding Company (Huan County, Gansu Province, China). All the experimental animals were slaughtered by arterial bleeding, while under general anesthesia induced intravenously with Propofol (B.Braun, Melsungen, Germany; 5 mg kg^− 1 body weight). Their meat yield and quality were measured after slaughter and are presented in Supplementary File 7. The Longissimus dorsi muscle was collected from the region between 12th and 13th ribs in the left half carcass from each goat. These tissue samples were immediately stored in liquid nitrogen for RNA extraction.

RNA-Sequencing

Total RNA from the ten muscle samples was isolated using Trizol reagent kit (Invitrogen, Carlsbad, CA, USA). Its concentration and RNA Integrity Number (RIN) were measured using a Nanodrop 2000 (Thermo Scientific, MA, USA) and Agilent 2100 Bioanalyzer (Agilent, CA, USA), respectively. Only samples with an RIN value > 7 were used for subsequent cDNA library construction. The Ribosomal RNA (rRNA) was removed from total RNA using a Ribo-Zero Gold rRNA Removal Kit (Illumina, CA, United States), and the remaining Longissimus dorsi muscle RNA was used to generate complementary DNA (cDNA) libraries using a NEBNext Ultra RNA Library Prep Kit (New England Biolabs, MA, United States). The cDNA libraries obtained were paired-end sequenced using an Illumina HiSeqTM 4000 sequencer (Illumina, CA, United States) at Gene Denovo Biotechnology Co., Ltd (Guangzhou, China).

Identification and characterization of lncRNAs

The clean reads were obtained by removing low-quality reads with quality scores < Q20, reads with > 10% unknown nucleotides, and adaptor reads from raw reads using Fastp v0.18.0. The clean reads were aligned against Caprine Genome Assembly ARS1 using HISAT2 v2.1.0 and the known lncRNAs were then identified according to the alignment results. Subsequently, the mapped results from HISAT2 v2.1.0 were assembled to transcripts using Stringtie v1.3.4 and the novel lncRNAs were then identified. Briefly, the transcripts with smaller than 200 nucleotides in length or those with a single exon were discarded. The coding potential of the remaining transcripts were predicted by CNCI v2.0 [51] and CPC v0.9-r2 [52] and the predicted results from the two kinds of software were intersected. These transcripts with coding potential score < 0 were considered to be novel lncRNAs. The identified caprine lncRNAs were characterized by analyzing types, transcript length, exon number, and coding potential score to mRNAs identified from the same muscle samples as those used in the study [26], using in-house Perl scripts.

Screening of differentially expressed lncRNAs

The expression level of each annotated lncRNA was normalized by calculating Fragments Per Kilobase of transcript per Million mapped reads (FPKM) using StringTie v1.3.4. The differential expression analysis of lncRNAs between the two goat breeds was performed by DESeq v2.0. The lncRNAs with Fold change > 2 and P-value < 0.05 were defined as differentially expressed lncRNAs.

Prediction and function annotation of the target genes of differentially expressed lncRNAs

The lncRNAs have been reported to affect the expression of their neighboring genes in cis, while they also influence co-expressed genes in trans. In our previous study, differentially expressed genes were screened from the same caprine Longissimus dorsi muscle tissues as those used in the study, using the thresholds of |fold change| > 2.0 and false discovery rate (FDR) value < 0.05 [26]. In order to obtain more target genes of lncRNAs, we adjusted the thresholds to |fold change| > 2.0 and P value < 0.05 to re-identify differentially expressed genes. The differentially expressed genes screened above were used to search cis and trans target genes. For screening of cis target genes, the coding genes located within 100 kb upstream and downstream of differentially expressed lncRNAs were searched. The search results were compared with differentially expressed genes obtained. The overlapped genes were chosen as cis target genes of differentially expressed lncRNAs. Subsequently, co-expression analysis was performed in expression levels between differentially expressed lncRNAs identified in this study and differentially expressed genes obtained. Briefly, Pearson's coefficients in expression levels between the differentially expressed lncRNAs and the differentially expressed genes were calculated, only the differentially expressed genes with |r| > 0.95 and P-value < 0.05 were chosen as trans target genes of differentially expressed lncRNAs.

To further investigate the function of differentially expressed lncRNAs, their cis and trans target genes were used to perform Gene Ontology (GO) and Kyoto Encyclopedia of genes and genomes (KEGG) functional enrichment analysis using the Gene Ontology database (http://www.geneontology.org/) and KEGG database (http://www.kegg.jp/kegg/) [53], respectively. The significant GO terms and KEGG pathways (P < 0.05) were screened based on hypergeometric test.

Construction of lncRNA-mRNA interaction networks

To better understand how these differentially expressed lncRNAs regulate the phenotype differences in meat yield and meat quality between LC and ZB goats, their target genes related to skeletal muscle development, intramuscular fat deposition, and meat tenderness were selected based on GO and KEGG enrichment analysis results as well as their functions reported in literatures. The lncRNA-mRNA interaction networks were finally constructed using Cytoscape v3.5.1.

A ceRNA network of lncRNA-miRNA-mRNA

Based on our previous Illumina HiSeq miRNA data (SRR16760528-SRR16760537) [28] and mRNA data (SRR13008213-SRR13008222) [26] obtained from the same caprine Longissimus dorsi muscle tissues as those used in the study using small RNA sequencing and RNA-Seq, respectively, a lncRNA-associated ceRNA network was constructed according to ceRNA hypothesis as follows [54]: (1) The target genes and binding lncRNAs of differentially expressed miRNAs between LC goats and ZB goats found from small RNA sequencing were firstly predicted using miReap v0.2, Miranda v3.3a, and TargetScan v7.0, and the predicted results were then overlapped with differentially expressed genes and differentially expressed lncRNAs obtained from RNA-Seq. These overlapped lncRNAs and target genes were preliminary used to construct the ceRNA network. (2) Correlation in expression between the lncRNA and the miRNA or between the miRNA and the mRNA was evaluated using the Spearman's rank correlation coefficient (SCC). Pairs with SCC < -0.7 were selected as negatively co-expressed lncRNA-miRNA or miRNA-mRNA pairs. (3) For lncRNA-miRNA and miRNA-mRNA pairs screened above, correlation in expression between the lncRNA and the mRNA was evaluated using the Pearson correlation coefficient (PCC), only those pairs with PCC > 0.9 were selected as co-expressed lncRNA-mRNA pairs. (4) Hypergeometric cumulative distribution function test was used to test whether a specific common miRNA sponge between a lncRNA and a mRNA was significant. Only pairs with P < 0.05 were selected as candidate ceRNA pairs. The lncRNA-miRNA-mRNA network was therefore visualized using Cytoscape v3.5.1.

Validation of differentially expressed lncRNAs

To validate the accuracy of RNA-Seq results, 12 differentially expressed lncRNAs appeared in the lncRNA-mRNA interaction networks or the ceRNA network were selected for reverse transcription-quantitative PCR (RT-qPCR) analysis. These included seven up-rgulated lncRNAs (XR_001917386.1, XR_001917605.1, MSTRG.262.1, MSTRG.12645.1, XR_001918614.1, XR_001918921.1, and XR_001917125.1) and five down-regulated lncRNAs (XR_001917948.1, XR_001919911.1, XR_001919918.1, XR_001919910.1, and MSTRG.12980.1) in Longissimus dorsi muscle of LC goats compared to ZB goats. The RNA samples extracted originally for RNA-Seq were used to synthesize cDNA using SuperScript™ II reverse transcriptase (Invitrogen, CA, USA). The RT-qPCR reaction was conducted in triplicate using 2 × ChamQ SYBR qPCR Master Mix (Vazyme, Nanjing, China). The expression levels of the lncRNAs were normalized by caprine GAPDH [24] and β-actin [23] and then calculated using 2^−ΔΔCt method. The primer information involved was listed in Table 1.

Table 1

PCR primers used for RT-qPCR
RNA	Forward (5'→3')	Reverse (5'→3')
XR_001917386.1	TATGGCCAGAGTTGGGAAAG	GGTCATGAGAGGGAGCTGAG
XR_001917948.1	CGCTCCTCTGTCTTCCTCAT	ACATTCACTCAGGGCTCCAA
XR_001917605.1	GGACCTTAGTGTGAGTGCCT	CCAGTGTCTTCATGCTGCTC
XR_001919911.1	AGATCTGTGGTTGCCTGTGA	ACACCCCTTCTTCTCTGCTC
XR_001919918.1	ACTTGTCTGGAATTCTGAGTTGC	TCATGAACACAAAGACCACAGA
XR_001919910.1	GAGCAGAGAAGAAGGGGTGT	TCGCCCTGTTAGTCTCCATC
MSTRG.262.1	TCTCTCTTCATCCCCTTGCC	TGCTCCTGACATCGTTTCCT
MSTRG.12645.1	TGGCAGAGTGTTCAGTCCAT	CTGGCGGGTTACAGTTCATG
MSTRG.12980.1	TGCTTTGGTCGTGTTGTCTG	CTCCCCTGTGGCTGTAAAGA
XR_001918614.1	TGATGTCAGGAGCAGTGGAG	CGTCTGCTTACACGTTTCCA
XR_001918921.1	TGGCAGCATCCTTTTCTTCT	GGGTGCTTCTCCTGCTAGTG
XR_001917125.1	CGCTGACTCTTCCACAATGA	CTTCCAGGTCAAAGGCTGAG
GAPDH	ACACTGAGGACCAGGTTGTG	GACAAAGTGGTCGTTGAGGG
β-actin	AGCCTTCCTTCCTGGGCATGGA	GGACAGCACCGTGTTGGCGTAA

LncRNA: long non-coding RNA

LC: Liaoning cashmere

ZB: Ziwuling black

MiRNAs: MicroRNAs

IGF2: Increasing expression of insulin growth factor 2

Six1: SIX homeobox 1

RNA-Seq: RNA sequencing

FPKM: Fragments Per Kilobase of transcript per Million mapped reads

FDR: False discovery rate

GO: Gene Ontology

KEGG: Kyoto Encyclopedia of genes and genomes

SCC: Spearman's rank correlation coefficient

PCC: Pearson correlation coefficient

RT-qPCR: Reverse transcription-quantitative PCR

MYL6B: Myosin light chain 6B

Acknowledgements

Not applicable.

Author contributions

JS and JW conceived and designed the experiments. JS, YL, JH, XL and SL performed the experiments. ZH, ML and ZZ analyzed the data. YZ, SY, LW and YG contributed reagents, materials and tools and collected the samples. JS and JW wrote the manuscript and revised the manuscript. All authors read and approved the final manuscript.

Funding

This research was funded by the fund for the Fuxi Young Talent Fund of Gansu Agricultural University (Gaufx-02Y02), the Science and Technology project of Lanzhou city (2021-1-162), the Central Guidance on Local Science and Technology Development Fund of Gansu Province, Lanzhou City Overseas Expertise Introduction Base for Molecular Breeding of Mutton Sheep, and the Projects of Gansu Agricultural University (GSAU-ZL-2015-033).

Availability of data and materials

The datasets generated in this study have been deposited at the Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra) under accession number PRJNA675157.

Ethics approval and consent to participate

Ethical approval by the Ethics Committee of Gansu Agricultural University, was obtained (Ethic approval file No. GSAU-ETH-AST-2021-028). All experimental procedures and sample collection methods were performed in accordance with approved guidelines and regulations to ensure animal welfare. In this study, written informed consent was obtained from Yongfeng Goat Breeding Company to use these animals. Meanwhile, the study is in accordance with ARRIVE guidelines.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Teixeira A, Silva S, Rodrigues S. Advances in sheep and goat meat products research. Adv Food Nutr Res. 2019, 87: 305-370.
Agliano F, Rathinam VA, Medvedev AE, Vanaja SK, Vella AT. Long noncoding RNAs in Host-Pathogen interactions. Trends immunol. 2019, 40: 492-510.
Yao RW, Wang Y, Chen LL. Cellular functions of long noncoding RNAs. Nat Cell Biol. 2019, 21: 542.
Kopp F, Mendell JT, Functional classification and experimental dissection of long noncoding RNAs. Cell. 2018, 172: 393-407.
Zhan S, Qin C, Li D, Zhao W, Nie L, Cao J, Guo J, Zhong T, Wang L, Li L, et al. A novel long noncoding RNA, lncR-125b, promotes the differentiation of goat skeletal muscle satellite cells by sponging miR-125b. Front Genet. 2019, 10: 1171.
Cheng P, Lu P, Guan J, Zhou Y, Zou L, Yi X, Cheng H. LncRNA KCNQ1OT1 controls cell proliferation, differentiation and apoptosis by sponging miR-326 to regulate c-Myc expression in acute myeloid leukemia. Neoplasma. 2020, 67: 238-248.
Zhu S, Chen CY, Hao Y. LncRNA KCNQ1OT1 acts as miR-216b-5p sponge to promote colorectal cancer progression via up-regulating ZNF146. J Mol Histol. 2021, 52: 479-490.
Cai B, Li Z, Ma M, Wang Z, Han P, Abdalla BA, Nie Q, Zhang X. LncRNA-Six1 encodes a micropeptide to activate Six1 in cis and is involved in cell proliferation and muscle growth. Front Physiol. 2017, 8: 230.
Li J, Yang T, Tang H, Sha Z, Chen R, Chen L, Yu Y, Rowe GC, Das S, Xiao J. Inhibition of lncRNA MAAT controls multiple types of muscle atrophy by cis- and trans-regulatory actions. Mol Ther. 2021, 29: 1102-1119.
Ren C, Deng M, Fan, Y, Yang H, Zhang G, Feng X, Li F, Dan W, Wang, F, Zhang Y. Genome-Wide analysis reveals extensive changes in lncrnas during skeletal muscle development in Hu sheep. Genes. 2017, 8: 191.
Li CY. Li X, Liu Z, Ni W, Zhang X, Hazi W, Ma Q, Zhang Y, Cao Y, Qi J, et al. Identification and characterization of long non-coding RNA in prenatal and postnatal skeletal muscle of sheep. Genomics. 2019, 111: 133-141.
Li Q, Liu R, Zhao H, Di R, Lu Z, Liu E, Wang Y, Chu M, Wei C. Identification and characterization of long noncoding rnas in ovine skeletal muscle. Animals. 2018, 8: 127.
Li H, Huang K, Wang P, Feng T, Shi D, Cui K, Luo C, Shafique L, Qian Q, Ruan J. Comparison of long non-coding RNA expression profiles of cattle and buffalo differing in muscle characteristics. Front Genet. 2020a, 11: 98.
Yan XM, Zhang Z, Liu JB, Li N, Yang GW, Luo D, Zhang Y, Yuan B, Jiang H, Zhang JB. Genome-wide identification and analysis of long noncoding RNAs in longissimus muscle tissue from Kazakh cattle and Xinjiang brown cattle. Anim Biosci. 2021, 34: 1739-1748.
Huang C, Ge F, Ma X, Dai R, Dingkao R, Zhaxi Z, Burenchao G, Bao P, Wu X, Guo X. Comprehensive analysis of mRNA, lncRNA, circRNA, and miRNA expression profiles and their ceRNA networks in the Longissimus dorsi muscle of cat-tle-yak and yak. Front Genet. 2021, 12: 772557.
Li T, Wan S, Wu R, Zhou X, Zhu D, Zhang Y. Identification of long non-protein coding RNAs in chicken skeletal muscle using next generation sequencing. Genomics. 2012, 99: 292-298.
Li Z, Ouyang H, Zheng M, Cai B, Han P, Abdalla BA, Nie Q, Zhang X Integrated analysis of long non-coding RNAs (lncRNAs) and mRNA expression profiles reveals the potential role of lncRNAs in skeletal muscle development of the chicken. Front Physiol. 2017, 7: 687.
Liu J, Zhou Y, Hu X, Yang J, Lei Q, Liu W, Han H, Li F, Cao D. Transcriptome analysis reveals the profile of long non-coding RNAs during chicken muscle development. Front Physiol. 2021, 12: 660370.
Wang J, Ren Q, Hua L, Chen J, Zhang J, Bai H, Li H, Xu B, Shi Z, Cao H. Comprehensive analysis of differentially expressed mRNA, lncRNA and circRNA and their ceRNA networks in the Longissimus dorsi muscle of two different pig breeds. Int J Mol Sci. 2019, 20: 1107.
Li R, Li B, Jiang A, Cao Y, Hou L, Zhang Z, Zhang X, Liu H, Kim KH, Wu W. Exploring the lncRNAs related to skeletal muscle fiber types and meat quality traits in pigs. Genes. 2020b, 11: 883.
Hou X, Wang L, Zhao F, Liu X, Gao H, Shi L, Yan H, Wang L, Zhang L. Genome-wide expression profiling of mRNAs, lncRNAs and circRNAs in skeletal muscle of two different pig breeds. Animals. 2021,11: 3169.
Ling Y, Zheng Q, Sui M, Zhu L, Xu L, Zhang Y, Liu Y, Fang F, Chu M, Ma Y. Comprehensive analysis of lncRNA reveals the temporal-specific module of goat skeletal muscle development. Int J Mol Sci. 2019, 20: 3950.
Zhan S, Dong Y, Zhao W, Guo J, Zhong T, Wang L, Li L, Zhang H. Genome-wide identification and characterization of long non-coding RNAs in developmental skeletal muscle of fetal goat. BMC Genomics, 2016, 17: 666.
Huang CN, Liu CL, Zeng SQ, Liu CB, Si WJ, Yuan Y, Ren LX, He YM, Zhang WY, Zhang HY. Identification of differentially expressed long non-coding RNAs and messenger RNAs involved with muscle development in Dazu black goats through RNA sequencing. Anim Biotechnol. 2022, 5: 1-9.
Wang JQ, Shen JY, Liu X, Li SB, Luo YZ, Zhao ML, Hao ZY, Ke N, Song YZ, Qiao LQ. Comparative analysis of meat production traits, meat quality, and muscle nutrient and fatty acid contents between Ziwuling black goats and Liaoning cashmere goats. Acta Pratacult Sin. 2021, 30: 166-177.
Shen J, Hao Z, Wang J, Hu J, Liu X, Li S, Ke N, Song Y, Lu Y, Hu L. Comparative transcriptome profile analysis of Longissimus dorsi muscle tissues from two goat breeds with different meat production performance using RNA-Seq. Front Genet. 2021, 11: 619399.
Shen J, Zhen H, Li L, Zhang Y, Wang J, Hu J, Liu X, Li S, Hao Z, Li M Identification and characterization of circular RNAs in Longissimus dorsi muscle tissue from two goat breeds using RNA-Seq. Mol Genet Genomics. 2022a, 297: 817-831.
Shen J, Hao Z, Luo Y, Zhen H, Liu Y, Wang J, Hu J, Liu X, Li S, Zhao Z Deep small RNA sequencing reveals important miRNAs related to muscle development and intramuscular fat deposition in Longissimus dorsi muscle from different goat breeds. Front Vet Sci. 2022b, 9: 911166.
Austin S, Medvedev A, Yan ZH, Adachi H, Hirose T, Jetten AM. Induction of the nuclear orphan receptor RORgamma during adipocyte differentiation of D1 and 3T3-L1 cells. Cell Growth & Differentiation. 1998, 9: 267-276.
Barendse W, Bunch RJ, Harrison BE. The effect of variation at the retinoic acid receptor-related orphan receptor C gene on intramuscular fat percent and marbling score in Australian cattle. J Anim Sci. 2010, 88: 47-51.
Kunath A, Weiner J, Krause K, Rehders M, Pejkovska A, Gericke M, Biniossek ML, Dommel S, Kern M, Ribas-Latre. A Role of Kallikrein 7 in Body Weight and Fat Mass Regulation. Biomedicines. 2021, 9: 131.
Gur S, Epon M, Epin S. Influence of growth rate in two growth periods on intramuscular connective tissue and palatability traits of beef. Czech J Anim Sci. 2003, 48: 113-119.
Hocquette JF, Gondret F, Baéza E, Médale F, Jurie C, Pethick DW. Intramuscular fat content in meat-producing animals: development, genetic and nutritional control, and identification of putative markers. Animal. 2010, 4: 303-319.
Kuang L, Lei M, Li C, Zhang X, Ren Y, Zheng J, Guo Z, Zhang C, Yang C, Mei X. Identification of long non-coding RNAs related to skeletal muscle development in two rabbit breeds with different growth rate. Int J Mol Sci. 2018, 19: 2046.
Shi T, Hu W, Hou H, Zhao Z, Shang M, Zhang L. Identification and comparative analysis of long non-coding RNA in the skeletal muscle of two Dezhou donkey strains. Genes. 2020, 11: 508.
Park JW, Lee JH, Han JS, Shin SP, Park TS. Muscle differentiation induced by p53 signaling pathway-related genes in myostatin-knockout quail myoblasts. Mol Biol Rep. 47: 9531-9540.
Krstic J, Reinisch I, Schupp M, Schulz TJ, Prokesch A. P53 functions in adipose tissue metabolism and homeostasis. Int J Mol Sci. 2018, 19: 2622.
Bean C, Fa Cc Hinello N, Faulkner G, Lanfranchi G. The effects of ANKRD2 alteration indicate its involvement in cell cycle regulation during muscle differentiation. BBA-Mol Cell Res. 2008, 1783: 1023-1035.
Wang YN, Yang WC, Li PW, Wang HB, Zhang YY, Zan LS. Myocyte enhancer factor 2A promotes proliferation and its inhibition attenuates myogenic differentiation via myozenin 2 in bovine skeletal muscle myoblast. Plos One. 2018, 13: e0196255.
Rosa AF, Moncau CT, Poleti MD, Fonseca LD, Balieiro JCC, Silva SLE, Eler JP. Proteome changes of beef in Nellore cattle with different genotypes for tenderness. Meat Sci. 2017, 138: 1.
Dransfield E, Martin JF, Bauchart D, Abouelkaram S, Lepetit J, Culioli J, Jurie C, Picard B. Meat quality and com-position of three muscles from French cull cows and young bulls. 2003, 76: 387-399.
Kemp TJ, Sadusky TJ, Simon M, Brown R, Eastwood M, Sassoon DA, Coulton GR Identification of a novel stretch-responsive skeletal muscle gene (Smpx). Genomics. 2001, 72: 260-271.
Tsai SH, Chang EY, Chang YC, Hee SW, Tsai YC, Chang TJ, Chuang LM Knockdown of RYR3 enhances adi-ponectin expression through an ATF3-dependent pathway. Endocrinology. 2013, 154: 1117-1129.
Bolado-Carrancio A, Riancho JA, Sainz J, Rodríguez-Rey JC Activation of nuclear receptor NR5A2 increases Glut4 expression and glucose metabolism in muscle cells. Biochem Bioph Res Co. 2014, 446: 614-619.
Sheela SG, Lee WC, Lin WW, Chung BC Zebrafish ftz-f1a (nuclear receptor 5a2) functions in skeletal muscle organization. Dev Biol. 2005, 286: 377-390.
Liu YX, Ma XM, Xiong L, Wu XY, Liang CN, Bao PJ, Yu QL, Yan P. Effects of intensive fattening with total mixed rations on carcass characteristics, meat quality, and meat chemical composition of yak and mechanism based on serum and transcriptomic profiles. Front Vet Sci. 2020, 7: 599418.
Yao CX, Wei QX, Zhang YY, Wang WP, Xue LX, Yang F, Zhang SF, Xiong CJ, Li WY, Wei ZR. MiR-200b targets GATA-4 during cell growth and differentiation. RNA Biol. 2013, 10: 465-480.
Jin X, Wang J, Hu J, Liu X, Li S, Lu Y, Zhen H, Li M, Zhao Z, Luo Y. MicroRNA-200b regulates the proliferation and differentiation of ovine preadipocytes by targeting p27 and KLF9. Animals. 2021, 11: 2417.
D'Agostino M, Torcinaro A, Madaro L, Marchetti L, Sileno S, Beji S, Salis C, Proietti D, Imeneo G, C Capogrossi M, De Santa F. Role of miR-200c in myogenic differentiation impairment via p66Shc: implication in skeletal muscle regenera-tion of dystrophic mdx mice. Oxid Med Cell Longev. 2018: 4814696.
Ronzoni FL, Giarratana N, Crippa S, Quattrocelli M, Cassano M, Ceccarelli G, Benedetti L, Van Herck J, Cusella De Angelis MG, Vitale M. Guide cells support muscle regeneration and affect Neuro-Muscular junction organization. Int J Mol Sci. 2021, 22: 1939.
Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, Liu Y, Chen R, Zhao Y. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013, 41: e166.
Kong L, Zhang Y, Ye ZQ, Liu XQ, Zhao SQ, Wei L, Gao G CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007, 35: W345-349.
Minoru K, Michihiro A, Susumu G, Masahiro H, Mika H, Masumi I, Toshiaki K, Shuichi K, Shujiro O, Toshiaki T. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008, 36: 480-484.
Quan J, Kang Y, Luo Z, Zhao G, Liu Z. Integrated analysis of the responses of a circRNA-miRNA-mRNA ceRNA network to heat stress in rainbow trout (Oncorhynchus mykiss) liver. BMC Genomics. 2021, 22: 48.

No competing interests reported.

Download PDF

Version 1

posted

You are reading this latest preprint version

Integrated transcriptome analysis reveals roles of long non- coding RNAs (lncRNAs) in caprine skeletal muscle mass and meat quality

Status:

Version 1

Abstract

Background

Results

Conclusion

Figures

Background

Results

Identification and characterization of lncRNAs in caprine skeletal muscle

Screening and validation of differentially expressed lncRNAs

Prediction and function annotation of the target genes of differentially expressed lncRNAs

Construction of lncRNA-mRNA interaction networks

Construction of a lncRNA-miRNA-mRNA ceRNA network

Discussion

Conclusion

Methods

Ethics statement

Experimental animals and sample collection

RNA-Sequencing

Identification and characterization of lncRNAs

Screening of differentially expressed lncRNAs

Prediction and function annotation of the target genes of differentially expressed lncRNAs

Construction of lncRNA-mRNA interaction networks

A ceRNA network of lncRNA-miRNA-mRNA

Validation of differentially expressed lncRNAs

Abbreviations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1