Identification and characterization of lncRNAs expression profile related to goat skeletal muscle at different development stages

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

Background

Long noncoding RNAs (lncRNAs) act vital roles in muscle growth and development, however, the reports of their expression characteristics and effects in goat skeletal muscle growth are few.

Result

Here, the RNA sequencing reads of goat longissimus doris tissues collected from kid (1 month old) and youth (9 months old) were mapped to the reference genome. As a result, we obtained an average of ~ 14.58 Gb high-quality clean reads. A total of 3,441 lncRNAs were identified, including 1,281 lincRNAs, 805 antisense lncRNAs and 1,355 sense_overlapping lncRNAs. These lncRNAs shared some properties in goat, such as shorter transcript length, fewer exons and shorter open reading frame (ORF). Among them, 36 differentially expressed lncRNAs (DE lncRNA) were identified, of which 10 lncRNAs were randomly selected to validate by RT-qPCR. Moreover, 79 co-location interactions were found between 26 lncRNAs and 71 genes by co-location analysis and functional annotation, and these genes were functionally involved in skeletal muscle related biological processes, such as skeletal system development, negative regulation of G1/S transition of mitotic cell cycle, negative regulation of cell cycle G1/S phase transition, myoblast fate determination and myosin complex. Additionally, trans regulation analysis suggested 36 lncRNAs and 2,684 mRNAs were included in co-expression networks. These co-expressed mRNAs were significantly enriched in muscle development related GO terms. The key lncRNA-mRNA interactions could reveal the function of lncRNA in goat and may be valuable in regulating skeletal muscle growth and development. In the last, the ceRNA (lncRNA-miRNA-mRNA) networks containing 4 lncRNAs, 3 miRNAs, and 8 mRNAs were constructed. Importantly, muscle development related mRNAs were included in the ceRNA networks, such as IFRD1, NKX2-5 and other genes.

Conclusion

This study revealed the lncRNAs expression profile in skeletal muscle growth and development and can facilitate further study on their function in goat muscle biology.

Goat

skeletal muscle

lncRNAs

RNA sequencing

regulation networks

As an important economically farm animals, the goats are raised for the utilization of meat, cashmere, and goat milk. With people’s living standards improve, the demand for goat meat is gradually increasing. As a result, lower meat production has hindered the development of the goat industry. Postnatal muscle growth is positively correlated with muscle fiber diameter, with larger muscle fiber diameters resulting in faster muscle growth rates [1, 2]. Recently, the studies of genome-wide transcriptional have provided a series of valuable candidate genes for regulating muscle growth rate in goat [3, 4]. Thus, unveiling the genetic mechanisms underneath muscle growth rates help us improve the meat production. The Wu’an black goats, a unique black goat breed in the Hebei province of China, have attracted increasing attention as the better meat performance. Good meat quality with rich in proteins, low content of fat and cholesterol as well as nice meat flavor are its outstanding characteristics. The Wu’an black goats with good meat quality showed the high breeding value [5, 6]. LncRNAs are transcription RNAs longer than 200 bp in length, with a complex structure and none protein-coding ability [7]. These molecules are widely distributed in different species, such as animals, plants, yeast and even viruses [8–11]. Moreover, lncRNAs have time-, space-, and tissue-specific expression, and show effects on multiple physiological processes, including cell proliferation, differentiation, and apoptosis [12]. However, the function of long noncoding RNAs (lncRNAs) in revealing skeletal muscle growth and development is still unclear and understanding the molecular mechanism of their action is of crucial importance.

Plenty of lncRNAs were verified to act crucial roles in skeletal muscle development of multiple species, such as goat, sheep, and pig [13–15]. Importantly, it has been confirmed that lncRNAs have the regulatory function in the growth and development of sheep muscle, which is closely related to goat. For example, to explore the functional lncRNAs during Hu sheep muscles growth, the genome-wide analysis was used to sequence the skeletal muscles at three key development stages (fetus, lamb, and adult). In total, there 6,924 lncRNAs were generated, of which the differentially expressed lncRNAs (DE lncRNA) are involved in the essential bio-function and processes, including skeletal muscle development [16]. Moreover, the lncRNA temporal expression patterns were identified in sheep longissimus doris across the gestational to postnatal stages, and the functional lncRNAs regulating myogenesis differentiation were described under deep sequencing data [17]. Similarly, the expression profiles of lncRNAs were verified in ovine (Texel and Ujumqin) gastrocnemius muscle at fetal (days 85 and 120 of gestation), newborn and adult stage [18]. In recent years, there has been increasing evidences indicate the crucial functions of lncRNAs in goats skeletal muscle development [9, 19]. For example, transcriptional sequencing analysis provided 577 DE lncRNAs at fetal (45, 60, and 105 days of gestation) and postnatal (3 days after birth) stage of Jianzhou big-eared goats, and which may play vital role in goat skeletal muscle development [20]. Although the genome-wide research of goat lncRNAs related to skeletal muscle development have been carried out. The information about functional lncRNAs related to skeletal muscle development is still limited. We hope to elucidate the functional mechanisms of DE lncRNAs in goat skeletal muscle tissue from a new perspective.

Here, we report the expression pattern and the potential role of the DE lncRNAs from the longissimus dori of Wu’an black goat at two developmental stages (kid: 1 month; youth: 9 months) using RNA sequencing. Moreover, the biological functions of DE lncRNAs molecules were annotated through GO analysis and KEGG pathway enrichment analysis. Finally, to clarify the molecular mechanism underpinning skeletal muscle development, we constructed the lncRNA-miRNA-mRNA ceRNA networks. The top priority of the study will facilitate a better understanding of transcriptomic changes in goat skeletal muscle development and will provide a reference for future studies on the molecular mechanism that improve the skeletal growth rates of goat.

RNA-Seq read filtering, mapping, and transcript assembly

To identified lncRNAs function in goat skeletal muscle growth, we built two cDNA libraries from longissimus dori samples at two grwoth stages: 1-month-old and 9-month-old. Three biological replicates were used. A total of 89.57Gb of raw data were generated. Low-quality and adaptor sequences were filtered out, and an average of ~14.58Gb high-quality RNA-seq data remained. For further analysis, we found the GC content averaged ranged from 53.71% to 64%. Moreover, the quality scores of Q20 and Q30 were above 96% and 91%, respectively (Table 2). The high-quality data were first mapped to the goat reference genome (CHIR_1.0, NCBI). Approximately 72.47%-85.25% of the high-quality reads were mapped to the goat reference genome (Additional file 1). The mapped sequences in the libraries were assembled and a total of 65,247 library transcripts were acquired

Identification and confirmation of lncRNAs in goat longissimus doris

We set a highly stringent filtering pipeline to discard transcripts that did not have all the characteristics of lncRNAs. The assembled transcripts from the two libraries were filtered to get candidate lncRNAs. A total of 3,441 lncRNAs were screened (Additional file 2), which included 1,281(37.2%) long intergenic ncRNAs (lincRNAs), 805 (23.4%) antisense lncRNAs, and 1,355 (39.4%) sense_overlapping lncRNAs (Fig. 2a). There was no sense_intronic lncRNAs. These 3,441 putative lncRNA were encoded by 2,675 genes. That is 1.3 transcripts on average per locus (Additional file 2). The lncRNA transcripts were widespread in chromosomes including 29 autosomes and the X chromosome (Additional file 3), which reflected the function diversity of lncRNAs.

The Illumina RNA-seq also produced 66 mRNAs. As shown in Fig. 2b and 2c, the lncRNAs transcripts length and exon number both were lower than that of the mRNA. The average length of lncRNA was 2,001bp with an average of 3.2 exons. Importantly, the principal lncRNA transcripts with 2 exons accounted for 63.6% of the 3,441 lncRNAs (Additional file 2). Moreover, compared to the mRNA, the average length of lncRNA open reading frame (ORF) was shorter (Fig. 2d). The result showed that lncRNA has lower coding potential.

Differential expression and functional enrichment analysis

To dissect the crucial lncRNAs involved in goat skeletal muscle growth, we explored the DE lncRNAs (P < 0.05, |log2FoldChange | > 2) for 1-month-old vs. 9-month-old group. The lncRNA expression profiles were assessed with FPKM. In this study, 36 lncRNAs were differentially expressed in goat skeletal muscle at these two developmental stages (Additional file 4). Compared with 9-month-old stage, 8 DE lncRNAs were upregulated and 28 were downregulated lncRNAs in 1-month-old stage (P < 0.05) (Fig. 3).

To evaluate the potential function of lncRNAs, we predicted the possible target gene associated with lncRNA in location (co-location) and expression (co-expression) relationships. In the co-location analysis, 30 lncRNAs (1 annotated lncRNAs and, 29 novel lncRNAs) were transcribed close to 71 protein coding genes (Additional file 5). Besides, the GO analysis of these 71 targets gene showed that 210 GO terms were significantly enriched (P < 0.05), including biological process (BP), cellular component (CC) and molecular function (MF) (Additional file 6). Notably, the significant enriched GO terms that were identified include skeletal system development, negative regulation of G1/S transition of mitotic cell cycle, negative regulation of cell cycle G1/S phase transition, myoblast fate determination and myosin complex (Additional file 6, Fig. 4a). Only the top 30 GO terms are shown in Additional file 6 and Fig. 4a. Additionally, 57 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were enriched through the pathway analysis. The top 20 significantly enriched KEGG analyses were shown in Additional file 6 and Fig. 4b. The DE lncRNAs target genes that were related to muscle development including Glycerolipid metabolism, signaling pathways regulating pluripotency of stem cells, fatty acid metabolism and biosynthesis of amino acids (Fig. 4b).

Furthermore, in the co-expression analysis, the top 200 target genes were chosen for next functional cluster analysis. The results showed that 3 lncRNAs (novel lncRNAs) corresponded to 103 target genes (Additional file 7). In total, 230 GO terms were obviously enriched (P < 0.05), of which the top 30 were shown in Additional file 8 and Fig. 4c. The skeletal muscle growth- and development-related GO terms, called regulation of transcription elongation from RNA polymerase I promoter, negative regulation of canonical Wnt signaling pathway, muscle tissue morphogenesis, muscle organ morphogenesis, cGMP biosynthetic process and negative regulation of myotube differentiation, were found (Additional file 8). Pathway analysis indicated that these target genes were involved in inositol phosphate metabolism, phosphatidylinositol signaling system, Notch signaling pathway, and HIF-1 signaling pathway (Additional file 8, Fig. 4d). Overall, DE lncRNAs and their target genes showed great potential in the regulation of skeletal muscle growth and development.

lncRNA-mRNA interaction network associated with goat skeletal muscle development

To investigate how lncRNAs interact with their cis- or trans- acting target genes, we constructed the lncRNA-mRNA regulatory networks associated with skeletal muscle development (Additional file 6， file 8). A total of 23 lncRNA gene pairs were obtained. The network provided candidate lncRNAs with muscle growth and development. Consequently, cis-regulatory network visualization revealed 13 interactions between 12 lncRNA and 12 gene that was related to muscle development (Fig. 5a), while trans-regulatory networks showed 10 interaction between2 lncRNAs and 5 genes (Fig. 5b).

Verification the expression profiles of DE lncRNAs with RT-qPCR

Ten DE lncRNAs were randomly selected to explore their expression profile at two experimental time (1 month, 9 months) using RT-qPCR. The data showed that the10 lncRNAs exhibited differential expression patterns during the process of goat skeletal growth and development (Fig. 6). Importantly, the expression trend of 10 lncRNA between RNA-seq and RT-qPCR were similar. The results reveal that the pipeline we set was reasonable in identifying putative lncRNAs, suggesting that the RNA-seq were credible.

The ceRNA networks construction and bioinformatics analysis

To further elucidate the potential role of lncRNA in goat skeletal muscle growth and development, we constructed ceRNA networks of DE mRNA, miRNA and lncRNAs (P < 0.05, |log2FoldChange| > 2) by Cytoscape software (Additional file 9), and found some ceRNA networks associated with muscle development based on the given muscle development-related mRNAs and DE lncRNA (Fig. 7). In the networks, there were 3 miRNAs,8 muscle development related mRNAs, 4 lncRNAs transcript, which had at least one predicted target miRNA. Moreover, miRNAs form the center of the network with lncRNA as the bait, and mRNA as the target, suggesting that lncRNA acts as a sponge of miRNA to regulate gene expression. For instance, lncRNA XR_001917947.1 act as a sponge for miRNAs, such as chi-miR-30b-3p, may to regulated mRNA smad3. This ceRNA network may provide valuable information for goat skeletal muscle myogenesis.

Skeletal muscle growth and development mainly affects the yield of meat in livestock [28]. We all known that skeletal muscle development contains a series of exquisitely regulated and orchestrated changes in the expression of many genes [29]. Importantly, the muscle development in goat involves two major stages. During the embryonic stage, muscle development is completed, and the number of muscle fibers generally does not change after birth. Postnatal muscle growth is mainly triggered by muscle fiber hypertrophy (the diameter and length of the muscle cells) and increased intermuscular fat [30]. Generally, kid goat (the newborn to 90-day-old stage) was related to muscle fibers hypertrophy and regulation of myoblast proliferation [31]. In addition, goats grow very quickly from 0 to 7 months, but slowly after 18 months [32]. In this study, the choice of development stages (1-month-old and 9-month-old) was suited for researching the molecular mechanism of muscle growth and development. We detected goat longissimus dori lncRNAs profile across two development stages by RNA-sequencing and bioinformatics analysis. Briefly, six healthy Wu’an black goats in the period of kid (1-month-old) and youth (9-months-old) were selected randomly. We investigated the transcript structure and expression patterns of goat skeletal muscle lncRNAs and further explored the function of lncRNA by the target genes predicted by cis- and trans-regulation.

Few studies have reported the crucial function of lncRNAs expression profiles in Wu’an black goat, especially in skeletal muscle growth and development. In this study, RNA-seq cDNA libraries were generated with Illumina NovaSeq. We acquired an average of ~ 14.58Gb high-quality clean reads after filtering process. Subsequently, sequence reads that mapped to the goat reference genome were assembled with String Tie. Thus, a total of 65,247 library transcripts were acquired. LncRNAs can be either single- or multi-exon, while it is difficult to distinguish putative lncRNAs from the plentiful sing-exon, lowly expressed and unreliable fragments assembled from RNA-seq [33, 34]. To minimize the selection of false positive lncRNAs, we set a relatively stringent filtering pipeline to verify true lncRNAs with high confidence. Only multi-exon lncRNAs were selected for further exploration, as was done in other studies [20, 35]. As a result, 3,441 putative lncRNAs with high confidence were identified. Moreover, most lncRNAs were longer than 2,000bp and contained 2 exons, which was agree with previous studies in goat and sheep [18, 20]. Consistent with other studies, our current data demonstrate that lncRNAs have shorter transcript length, fewer exons, and shorter open reading frames (ORFs) than mRNA [19, 20]. It is vital to expanding our understanding of lncRNAs via association with multiple structural features. These similarities verified that the putative lncRNAs verified in our study were credible.

In recent years, numerous studies have reported the biological functions of lncRNAs. For instance, lnc-SEMT regulates IGF2 expression via competing with miR-125b to facilitate sheep skeletal muscle growth and development [14]. Moreover, lncR125b facilitates the differentiation of goat skeletal muscle satellite cells via competing with miR-125b [13]. Additionally, a series of lncRNAs, such as lncIRS1, lnc-RAM and lncMD1, have been identified to have an essential effect on skeletal myogenesis [36–38]. LncRNAs are non-coding transcripts that act as regulators of gene expression and involve in skeletal muscle biological process [39]. In the present study, we detected 36 DE lncRNAs at two development stages of goat skeletal muscles. They showed different significant changes between 1 month and 9 months old goat. These lncRNAs may have certain biological functions during skeletal muscle growth and development. Consequently, the DE lncRNAs showed in this study can be regarded as vital putative regulators of muscle biology. Furthermore, the expression trend of 10 randomly DE lncRNA (either up or down regulated) identified by RT-qPCR was consistent with the RNA-seq data. Together, the information validated that the identified lncRNAs were of high quality. Importantly, the data provided a valuable resource for exploring the role of lncRNA in postnatal muscle growth and new insight into the dynamic gene regulation of muscle biology in Wu’an black goat.

Unlike protein-coding genes, the biofunctions of lncRNAs cannot directly be speculate from their sequence or structure. Therefore, we attempted to uncover the function of lncRNAs with intersection of their cis- and tans- acting target genes in our study. [20, 31, 40]. Cis-acting refer to the co-location interactions between lncRNAs and mRNA. We searched the genes 100 kb upstream and downstream of the identified lncRNAs comparable with previous studies [16, 41]. The cis-regulating target genes act vital role in estimating the biological function of lncRNAs. Emerging studies have document that lncRNAs were related to cis-regulatory in skeletal muscle biology. For instance, the lncRNA Dum can facilitate myoblast differentiation and damage-induced muscle regeneration via silencing its cis-acting target gene, Dppa2 [32]. Thus, further research was need to validate the biological function of lncRNAs on muscle development.

Still, some lncRNAs were found to function in trans-acting to target gene loci distant from where the lncRNAs are transcribed [42]. For instance, LncRNA MUNC, encode 5 kb upstream of the MyoD transcription start site, facilitates the biofunction of MyoD in muscle biology [43]. According to co-expression analysis, we chose the DE lncRNA and mRNA that are on different chromosome. As a result, a total of 10 pairs of lncRNAs and mRNA were generated to have trans-acting interactions in 1 month vs. 9 months group. The information above validates that lncRNA may involve in goat skeletal muscle biology via cis- or trans- acting regulation.

For functional exploration of lncRNAs, GO and KEGG analyses was carried out based on the cis- and trans-acting target gene. IFRD1 plays vital role in myoblast differentiation via modulate MyoD and NF-kappaB expression [44]. In the co-location analysis, IFRD1, also the cis-target gene predicted by TCONS_00026389 and TCONS_00026388, is enriched for the term myoblast fate determination and myoblast fate commitment. Furthermore, some mRNA modulated by lncRNA trans-action are also clustered in terms and pathways associated with muscle development. For example, NKX2-5 is enriched in the terms of negative regulation of canonical Wnt signaling pathway, muscle tissue morphogenesis and muscle organ morphogenesis. It has been reported that NKX2-5 could regulate differentiation of skeletal myoblasts in vitro [45]. These muscle development-related genes are all from target prediction. Overall, these results validate that lncRNA may involve in the goat skeletal muscle development.

Currently, the ceRNA mechanism is widely reported [46, 47], in which lncRNAs function as molecular sponges that compete with mRNA for target miRNAs. In our study, a lncRNA-associated ceRNA network related to goat muscle development were constructed. We firstly identified 4 lncRNAs, 3 miRNAs, and DE mRNA between 1 month and 9 months stage. According to target prediction and expression correlation analysis, we obtained 15 ceRNAs of 4 lncRNs, 3 miRNAs and 8 mRNA that were related to muscle development. Then, a ceRNA network with these 15 lncRNA-miRNA-mRNA ceRNAs was constructed. Additionally, the lncRNA-trans-regulation networks indicate that the lncRNAs interact with 1 or 2 miRNAs with their predicted target genes related to muscle development. Collectively, the result suggests the potential molecular regulation of lncRNAs in goat skeletal muscle growth and development. Nevertheless, the biological functions of the lncRNA-miRNA-mRNA interaction described in this study need further studies to verify. All the studies referred above emphasize that lncRNAs are a vital part of the regulatory network of skeletal muscle growth and development. The bioinformatics results suggested that the lncRNAs showed significant changes during goat skeletal muscle growth, which means they may have certain functional during myofiber growth.

The study indicates the systematic description of lncRNAs in goat skeletal muscle growth at two development stages using RNA-seq. The DE lncRNAs associated with goat muscle growth were validated. The verified lncRNAs in our study shared lots of features in goat. The function of DE lncRNAs were annotated via GO and KEGG pathway enrichment analysis. Furthermore, the visual lncRNA-associated ceRNA network provide a fuller understanding of candidate lncRNAs that could contribute to goat skeletal muscle growth rates. However, the ceRNA network and the ceRNA interaction is results of prediction based on gene function and their roles still needs further study to validate.

Animal Preparation and Sample Collection

The animals used in this study were the Wu’an black goats from Yutian Black Goat farm (Wu’an, Hebei, China). A total of six health female Wu’an black goats were selected covering two groups at different growth stages: 1-month-old and 9-month- old. The animal samples (longissimus dorsi muscle) were harvested and immediately frozen in liquid nitrogen and stored at -80 ºC for RNA sequencing analysis.

RNA Isolation, lncRNA-seq Library Construction, and Sequencing

Total RNA was isolated from the 6 muscle tissue samples with TRIzol reagent (Invitrogen, Carlsbad, CA, USA). RNA quantity and purity were detected with the NanoDrop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA). The Ribo-zeroTM Gole Kit (Illumina, San Diego, CA, USA) were selected to remove Ribosomal RNAs from the total RNA according to the manufacturer’s instructions. RNA-seq complementary DNA (cDNA) libraries were generated and then sequenced on an Illumina NovaSeq platform at Novogene (Beijing, China).

Raw data in fastq format were treated by inhouse script as previous. After filtering process, high quality clean reads were generated. Here we filtered out reads containing adapter, poly-N and low-quality reads from the raw data. Meanwhile, the Phred score of Q20 and Q30 as well as GC content of the high-quality clean data were calculated. Based on the high-quality clean data, all the downstream analyses were completed. We downloaded the goat reference genome and gene model annotation files from the genome website. HISAT2 was chosen to align the high-quality clean reads to the reference genome. Based on it, the mapped transcripts were assembled using String Tie following a reference-based approach. Finally, both known genes and new transcripts from the results of HISAT2 alignment were constructed by Cuffmerge program and annotated by cuffcompare.

Prediction of Multiple Exon lncRNAs

The assembled new transcripts were filtered from the lncRNA-seq libraries to acquire the putative lncRNAs. The identification criteria of credible lncRNA candidates was in line with the protocol shown in Fig. 1. We filtered out the assembled transcripts based on the steps as follows: (1) Low-confidence single exon transcripts or transcripts shorter than 200bp long were discarded. (2) The remaining transcripts that overlapped with the exon region of the database annotated by cuffcompare were removed. (3) The transcripts that contained protein-coding potential also were filtered. The software programs of Coding Potential Calculator [21] and Coding-Non-Coding Index [22] as well as Protein Families Database [23] were used to predicted the coding potential of transcripts. Only the transcripts that have no protein coding domain belonged to the credible lncRNAs.

Classification and Differential Expression analysis

Comparison of potential lncRNA with the known transcripts based on their location in the genome using cuffmerge [24].Then, the lnRNAs classification and characteristic were analyzed. In line with the comparison results, lncRNA can be separated into three categories: lincRNAs, antisense lncRNAs, sense_overlapping lncRNAs. String Tie was selected to examine lncRNAs expression levels by computing fragments per kilo bases per million reads (FPKMs). Differentially expressed lncRNA (DE lncRNA) between the two libraries were verified by edgeR. The corrected P-value < 0.05 and an absolute value of the |log2FoldChange| > 2 were set as the threshold to assess statistically significance differences of lncRNA expression.

Co-location and Co-expression analysis, and Functional Annotation Analysis

The co-location and co-expression correlations of DE lncRNA and DE mRNA were processed based on the FPKMs to investigate the function of lncRNAs. LncRNAs exert cis-acting effects on their co-localized genes [25]. In this research, DE mRNA located 100 kb/100 kb of DE lncRNA were identified as cis-acting target gene. The trans-acting correlations between lncRNAs and genes (co-expression) was calculated with the Pearson’s correlation coefficient (|r| > 0.95 and P < 0.05).

The cis- and trans-acting target genes were annotated by Gene Ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) to explore the main function of lncRNAs. GO terms and KEGG pathways with a corrected P-value < 0.05 were defined as significantly enriched.

Verification of lncRNA expression pattern with RT-qPCR

We randomly selected 10 DE lncRNAs, which taken part in the process of skeletal muscle growth and development. Primers for the 10 lncRNAs and endogenous reference gene (Table 1) were designed by Primer 5, and the goat GAPDH was chosen as the endogenous reference. For RT-qPCR analysis, the cDNA as the template was converted from total RNA by PrimeScript™ RT reagent kit (Takara, Beijing, China). Then the qPCR was performed by a RocheLight Cycler^®480 Ⅱ system with SYBR Green qPCR Mix kit as our previous study [26]. The expression level was analyzed with the method of comparative threshold (2^-ΔΔCt).

Construction of the ceRNA (lncRNA-miRNA-mRNA) network

To further describe the functional role of lncRNA in skeletal muscle development, the ceRNA (lncRNA-miRNA-mRNA) network was build relied on the theory that lncRNA directly associated with miRNA and affect the activity of its mRNAs [27]. Firstly, the correlation between lncRNA and miRNA or miRNA and mRNA were analyzed by Pearson correlation coefficient. All pairs with |r| > 0.95 and adjusted P < 0.05 were chosen as potential lncRNA-miRNA or miRNA-mRNA pairs. Then, DE lncRNAs, miRNAs, and mRNAs were chosen to construct the ceRNA network. Finally, the network was visualized by Cytoscape software.

Statistical Analysis

RT-qPCR data and graphs were generated by GraphPad Prism 6.0 (San Diego, CA, USA). The results are presented as means ± SEMs. The unpaired two-tailed t-test was performed to determine the statistical difference. All experiments were performed a minimum of three replicates. Minimal standard of statistical significance was set at P < 0.05 or P < 0.01.

lncRNA

long chain non-coding RNA

ORF

open reading frame

BP

biological process

CC

cellular component

MF

molecular function

KEGG

Kyoto Encyclopedia of Genes and Genomes.

Ethics approval and consent to participate

All animal experiments in this study were carried out in adherence to the Regulations for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China, 2004). The experimental protocol was approved by the respective Animal Research and Ethics Committees of Hebei University of Engineering (AEE16015). The study was carried out in accordance with ARRIVE guidelines.

Consent for publication

Not applicable.

Availability of data and materials

The datasets supporting the conclusions of this article are included within the article and its additional files. The raw data we obtained from RNA-seq can be found in the National Center for Biotechnology Information (NCBI) Sequence Read Archive with the accession number PRJNA749569.

Competing interests

The authors declare that they have no competing interests.

Funding

This research was funded by the National Natural Science Foundation of China Youth Program (32102509)， Natural Science Foundation of Hebei Province of China for Youths (C2019402261 and C2019402293) and Natural Science Foundation of Hebei Province of China for Youths (C2019402293).

Author’ contributions

YFL designed the study. YHG supervised the work. HYH, WTL, JNL, YKF, HZ and JQY participated in performed experiments and processed the sample as well as analyzed the data. HYH and YFL prepared the first draft of the manuscript. All authors reviewed, edited, and supported the final version of the manuscript for submission.

Acknowledgements

We thank all contributors of the present study. We also thank Mr Guo Yanrong and his colleagues at Yutian Black Goat farm (Wuan, Hebei, China) for providing animals and assistance in animals processing.

Zeng C, Shi H, Kirkpatrick LT, Ricome A, Park S, Scheffler JM, et al. Driving an oxidative phenotype protects myh4 null mice from myofiber loss during postnatal growth. Front Physiol. 2021; 12:785151.
Kim DH, Choi YM, Lee J, Shin S, Kim S, Suh Y, et al. Differential expression of mstn isoforms in muscle between broiler and layer chickens. Animals (Basel). 2022; 12(5).
Yang BG, Yuan Y, Zhou DK, Ma YH, Mahrous KF, Wang SZ, et al. Genome-wide selection signal analysis of australian boer goat reveals artificial selection imprinting on candidate genes related to muscle development. Anim Genet. 2021; 52(4):550–555.
Shen J, Hao Z, Wang J, Hu J, Liu X, Li S, et al. Comparative transcriptome profile analysis of longissimus dorsi muscle tissues from two goat breeds with different meat production performance using rna-seq. Front Genet. 2020; 11:619399.
Li L, Zhang J, Zhu JQ, Gu S, Sun Q, Zhou GM, et al. Genetic diversity of nine populations of the black goat (capra hircus) in sichuan, pr china. Zoolog Sci. 2006; 23(3):229–34.
Liu K, Jin Y, Zhang FY, Zhang YZ, Quan XQ, Han QD, et al. The complete mitochondrial DNA sequence of yimeng black goat (capra hircus) and its potential application in mutton discrimination. Mitochondrial DNA B Resour. 2020; 5(3):2303–2305.
Prensner JR and Chinnaiyan AM The emergence of lncrnas in cancer biology. Cancer Discov. 2011; 1(5):391–407.
Shi X, Sun M, Liu H, Yao Y, and Song Y Long non-coding rnas: A new frontier in the study of human diseases. Cancer Lett. 2013; 339(2):159–66.
Sendler E, Johnson GD, Mao S, Goodrich RJ, Diamond MP, Hauser R, et al. Stability, delivery and functions of human sperm rnas at fertilization. Nucleic Acids Res. 2013; 41(7):4104–17.
Ma L, Bajic VB, and Zhang Z On the classification of long non-coding rnas. RNA Biol. 2013; 10(6):925–33.
Wang Z and Li X The role of noncoding rna in hepatocellular carcinoma. Gland Surg. 2013; 2(1):25–9.
Ming GF, Wu K, Hu K, Chen Y, and Xiao J Nampt regulates senescence, proliferation, and migration of endothelial progenitor cells through the sirt1 as lncrna/mir-22/sirt1 pathway. Biochem Biophys Res Commun. 2016; 478(3):1382–8.
Zhan S, Qin C, Li D, Zhao W, Nie L, Cao J, 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.
Wei C, Wu M, Wang C, Liu R, Zhao H, Yang L, et al. Long noncoding rna lnc-semt modulates igf2 expression by sponging mir-125b to promote sheep muscle development and growth. Cell Physiol Biochem. 2018; 49(2):447–462.
Ma X, Fu D, Chu M, Ding X, Wu X, Guo X, et al. Genome-wide analysis reveals changes in polled yak long non-coding rnas in skeletal muscle development. Front Genet. 2020; 11:365.
Ren C, Deng M, Fan Y, Yang H, Zhang G, Feng X, et al. Genome-wide analysis reveals extensive changes in lncrnas during skeletal muscle development in hu sheep. Genes (Basel). 2017; 8(8).
Yuan C, Zhang K, Yue Y, Guo T, Liu J, Niu C, et al. Analysis of dynamic and widespread lncrna and mirna expression in fetal sheep skeletal muscle. PeerJ. 2020; 8:e9957.
Li Q, Liu R, Zhao H, Di R, Lu Z, Liu E, et al. Identification and characterization of long noncoding rnas in ovine skeletal muscle. Animals (Basel). 2018; 8(7).
Ling Y, Zheng Q, Sui M, Zhu L, Xu L, Zhang Y, et al. Comprehensive analysis of lncrna reveals the temporal-specific module of goat skeletal muscle development. Int J Mol Sci. 2019; 20(16).
Zhan S, Dong Y, Zhao W, Guo J, Zhong T, Wang L, et al. Genome-wide identification and characterization of long non-coding rnas in developmental skeletal muscle of fetal goat. BMC Genomics. 2016; 17:666.
Kong L, Zhang Y, Ye ZQ, Liu XQ, Zhao SQ, Wei L, 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-9.
Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, et al. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013; 41(17):e166.
Sonnhammer EL, Eddy SR, Birney E, Bateman A, and Durbin R Pfam: Multiple sequence alignments and hmm-profiles of protein domains. Nucleic Acids Res. 1998; 26(1):320–2.
Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, et al. Differential gene and transcript expression analysis of rna-seq experiments with tophat and cufflinks. Nat Protoc. 2012; 7(3):562–78.
Bonasio R and Shiekhattar R Regulation of transcription by long noncoding rnas. Annu Rev Genet. 2014; 48:433–55.
Liu Y, Zhou Z, He X, Tao L, Jiang Y, Lan R, et al. Integrated analyses of mirna-mrna expression profiles of ovaries reveal the crucial interaction networks that regulate the prolificacy of goats in the follicular phase. BMC Genomics. 2021; 22(1):812.
Thomson DW and Dinger ME Endogenous microrna sponges: Evidence and controversy. Nat Rev Genet. 2016; 17(5):272–83.
Berri C, Le Bihan-Duval E, Debut M, Sante-Lhoutellier V, Baeza E, Gigaud V, et al. Consequence of muscle hypertrophy on characteristics of pectoralis major muscle and breast meat quality of broiler chickens. J Anim Sci. 2007; 85(8):2005–11.
Tripathi AK, Patel AK, Shah RK, Patel AB, Shah TM, Bhatt VD, et al. Transcriptomic dissection of myogenic differentiation signature in caprine by rna-seq. Mech Dev. 2014; 132:79–92.
Braun T and Gautel M Transcriptional mechanisms regulating skeletal muscle differentiation, growth and homeostasis. Nat Rev Mol Cell Biol. 2011; 12(6):349–61.
Zhou J, Xiong Q, Chen H, Yang C, and Fan Y Identification of the spinal expression profile of non-coding rnas involved in neuropathic pain following spared nerve injury by sequence analysis. Front Mol Neurosci. 2017; 10:91.
Wang L, Zhao Y, Bao X, Zhu X, Kwok YK, Sun K, et al. Lncrna dum interacts with dnmts to regulate dppa2 expression during myogenic differentiation and muscle regeneration. Cell Res. 2015; 25(3):335–50.
Cabili MN, Trapnell C, Goff L, Koziol M, Tazon-Vega B, Regev A, et al. Integrative annotation of human large intergenic noncoding rnas reveals global properties and specific subclasses. Genes Dev. 2011; 25(18):1915–27.
Derrien T, Johnson R, Bussotti G, Tanzer A, Djebali S, Tilgner H, et al. The gencode v7 catalog of human long noncoding rnas: Analysis of their gene structure, evolution, and expression. Genome Res. 2012; 22(9):1775–89.
Yue B, Li H, Liu M, Wu J, Li M, Lei C, et al. Characterization of lncrna-mirna-mrna network to reveal potential functional cernas in bovine skeletal muscle. Front Genet. 2019; 10:91.
Li Z, Cai B, Abdalla BA, Zhu X, Zheng M, Han P, et al. Lncirs1 controls muscle atrophy via sponging mir-15 family to activate igf1-pi3k/akt pathway. J Cachexia Sarcopenia Muscle. 2019; 10(2):391–410.
Yu X, Zhang Y, Li T, Ma Z, Jia H, Chen Q, et al. Long non-coding rna linc-ram enhances myogenic differentiation by interacting with myod. Nat Commun. 2017; 8:14016.
Cesana M, Cacchiarelli D, Legnini I, Santini T, Sthandier O, Chinappi M, et al. A long noncoding rna controls muscle differentiation by functioning as a competing endogenous rna. Cell. 2011; 147(2):358–69.
Knoll M, Lodish HF, and Sun L Long non-coding rnas as regulators of the endocrine system. Nat Rev Endocrinol. 2015; 11(3):151–60.
Zhang S, Qin C, Cao G, Xin W, Feng C, and Zhang W Systematic analysis of long noncoding rnas in the senescence-accelerated mouse prone 8 brain using rna sequencing. Mol Ther Nucleic Acids. 2016; 5(8):e343.
Liu Y, Qi B, Xie J, Wu X, Ling Y, Cao X, et al. Filtered reproductive long non-coding rnas by genome-wide analyses of goat ovary at different estrus periods. BMC Genomics. 2018; 19(1):866.
Nie M, Deng ZL, Liu J, and Wang DZ Noncoding rnas, emerging regulators of skeletal muscle development and diseases. Biomed Res Int. 2015; 2015:676575.
Mueller AC, Cichewicz MA, Dey BK, Layer R, Reon BJ, Gagan JR, et al. Munc, a long noncoding rna that facilitates the function of myod in skeletal myogenesis. Mol Cell Biol. 2015; 35(3):498–513.
Micheli L, Leonardi L, Conti F, Maresca G, Colazingari S, Mattei E, et al. Pc4/tis7/ifrd1 stimulates skeletal muscle regeneration and is involved in myoblast differentiation as a regulator of myod and nf-kappab. J Biol Chem. 2011; 286(7):5691–707.
Riazi AM, Lee H, Hsu C, and Van Arsdell G Csx/nkx2.5 modulates differentiation of skeletal myoblasts and promotes differentiation into neuronal cells in vitro. J Biol Chem. 2005; 280(11):10716–20.
Wang H, Zhong J, Zhang C, Chai Z, Cao H, Wang J, et al. The whole-transcriptome landscape of muscle and adipose tissues reveals the cerna regulation network related to intramuscular fat deposition in yak. BMC Genomics. 2020; 21(1):347.
Zhao R, Li FQ, Tian LL, Shang DS, Guo Y, Zhang JR, et al. Comprehensive analysis of the whole coding and non-coding rna transcriptome expression profiles and construction of the circrna-lncrna co-regulated cerna network in laryngeal squamous cell carcinoma. Funct Integr Genomics. 2019; 19(1):109–121.

Table 1 Primer information of RT-qPCRd in the study

Oligo Name

Sequence

TCONS_00036756

F: 5’ AGGCTGCAATCCACGCTAA 3’

R: 5’ TCCAACTCTGTGTGACCCCAT 3’

TCONS_00066056

F:5’ AGAAAATGAATCCCTGGAGTGTG 3’

R:5’ AACGCTGACCACCATGATGAC 3’

TCONS_00062751

F：5' CTGGGAGAATACAAAGGGGG 3’

R: 5’ GGATCTACGGGCCTTTTGTCT 3’

XR_001296113.2

F:5’ GCCGCCGTGAAGACTATTG 3’

R:5’ CCATGAAGCCAGGGTACAAAC 3’

TCONS_00150002

F: 5’ CCCCGAATGTAAGCAATGAG 3’

R:5’ AGGAGACCTACCGCTACCTGAG 3’

TCONS_00134767

F: 5’ CCCAACAAAGTGCCCAGAC 3’

R:5’ GGAGAAGACGGCGTTATGC 3’

TCONS_00126170

F: 5’ GCTAGTCCCAGACAGCATTCAT 3’

R: 5’ GGTGTTGTTCTCGCCTGGAA 3’

TCONS_00124841

F: 5’ CCCTTACCACAGGCACCACT 3’

R: 5’ CAGGTGAGAAGGTGTGTTCTGG 3’

TCONS_00121766

F: 5’ TGTCCCCAACCTCGGTATCT 3’

R: 5’ GGTCAAACCTCTGAGCCTCG 3’

TCONS_00026838

F: 5’CTTCTCCTTGCTTGGCACCT 3’

R: 5’ CAGGTGCCAAGCAAGGAGA 3’

GAPDH

F: 5’ CACGGCACAGTCAAGGCAG 3’'

R: 5’ AGATGATGACCCTCTTGGCG 3’

Table 2 Data summary of RNA-seq in goat muscle

Sample name	Raw reads	Clean reads	Raw bases(G)	Clean bases(G)		Error rate (%)	Q20 (%)	Q30 (%)	GC Content (%)
mon1_1	98721676	97036650	14.81		14.56	0.03	97.47	93	52.99
mon1_2	104194360	101986738	15.63		15.3	0.03	97.29	92.64	55.08
mon1_3	99110136	97525086	14.87		14.63	0.03	96.95	91.87	53.06
mon9_1	98408274	96645866	14.76		14.5	0.03	97.4	92.79	51.64
mon9_2	90534786	88937344	13.58		13.34	0.03	97.48	92.97	53.35
mon9_3	106156574	101082860	15.92		15.16	0.03	96.92	91.93	54.31

No competing interests reported.

Additionalfile1.xls
Additional file 1：Table S1. The information of RNA-seq data
Additionalfile2.xls
Additional file 2：Table S2. List of lncRNAs identified in the goat skeletal muscle libraries.
Additionalfile3.pdf
Additional file 3：Fig. S1. Chromosome distribution of lncRNAs identified in goat skeletal muscle.
Additionalfile4.xls
Additional file 4：Table S3. List of differentially expressed lncRNAs from two skeletal muscle developmental stages.
Additionalfile5.xls
Additional file 5：Table S4. Protein-coding genes detected 100 Kb upstream and downstream of the differentially expressed lncRNAs.
Additionalfile6.xls
Additional file 6：Table S5. The top 30 GO terms of functional enrichment analysis
Additionalfile7.xlsx
Additional file 7：Table S6. Co-expression analysis between protein-coding genes and lncRNAs
Additionalfile8.xlsx
Additional file 8：Table S7. Functional enrichment analysis of protein-coding genes co-expressed with lncRNAs.
Additionalfile9.pdf
Additional file 9：Fig.S2. The lncRNA-miRNA-mRNA ceRNA networks.

Identification and characterization of lncRNAs expression profile related to goat skeletal muscle at different development stages

Status:

Version 1

Abstract

Background

Result

Conclusion

Figures

Background

Results

Discussions

Conclusions

Methods

Abbreviations

Declarations

References

Tables

Additional Declarations

Supplementary Files

Status:

Version 1