Genome-wide identification and characterization of long non-coding RNAs in developmental skeletal muscle of fetal goat
- Siyuan Zhan†1,
- Yao Dong†1,
- Wei Zhao1,
- Jiazhong Guo1,
- Tao Zhong1,
- Linjie Wang1,
- Li Li1Email author and
- Hongping Zhang1Email author
© The Author(s). 2016
Received: 31 May 2016
Accepted: 10 August 2016
Published: 22 August 2016
Long non-coding RNAs (lncRNAs) have been studied extensively over the past few years. Large numbers of lncRNAs have been identified in mouse, rat, and human, and some of them have been shown to play important roles in muscle development and myogenesis. However, there are few reports on the characterization of lncRNAs covering all the development stages of skeletal muscle in livestock.
RNA libraries constructed from developing longissimus dorsi muscle of fetal (45, 60, and 105 days of gestation) and postnatal (3 days after birth) goat (Capra hircus) were sequenced. A total of 1,034,049,894 clean reads were generated. Among them, 3981 lncRNA transcripts corresponding to 2739 lncRNA genes were identified, including 3515 intergenic lncRNAs and 466 anti-sense lncRNAs. Notably, in pairwise comparisons between the libraries of skeletal muscle at the different development stages, a total of 577 transcripts were differentially expressed (P < 0.05) which were validated by qPCR using randomly selected six lncRNA genes. The identified goat lncRNAs shared some characteristics, such as fewer exons and shorter length, with the lncRNAs in other mammals. We also found 1153 lncRNAs genes were neighbored 1455 protein-coding genes (<10 kb upstream and downstream) and functionally enriched in transcriptional regulation and development-related processes, indicating they may be in cis-regulatory relationships. Additionally, Pearson’s correlation coefficients of co-expression levels suggested 1737 lncRNAs and 19,422 mRNAs were possibly in trans-regulatory relationships (r > 0.95 or r < −0.95). These co-expressed mRNAs were enriched in development-related biological processes such as muscle system processes, regulation of cell growth, muscle cell development, regulation of transcription, and embryonic morphogenesis.
This study provides a catalog of goat muscle-related lncRNAs, and will contribute to a fuller understanding of the molecular mechanism underpinning muscle development in mammals.
KeywordsMuscle development LncRNA Goat Transcriptome cis-acting trans-acting Differential expression
Genome-wide transcriptional studies have revealed that large regions of eukaryotic genomes are transcribed into a heterogeneous population of non-coding RNAs (ncRNAs). Generally, ncRNAs shorter than 200 nucleotides are usually described as small/short ncRNA, such as microRNAs (miRNAs), PIWI-interacting RNAs (piRNAs), small interfering RNAs (siRNAs), and classical ncRNAs such as ribosomal RNAs (rRNAs), transfer RNAs (tRNAs), and small nucleolar RNAs (snoRNAs), whereas ncRNAs longer than 200 nucleotides are described as long ncRNAs (lncRNAs). In the past few years, an increasing number of lncRNAs have been discovered in mammal, including Homo sapiens [1, 2], Mus musculus [3–7], Bos taurus [8, 9], Sus scrofa [10–13], and Ovis aries . And accordingly unveiled that lncRNAs play critical roles in biological processes like transcriptional regulation [15–17], epigenetic modification [18–20], development [21–23], cell differentiation [24–26], as well as in some diseases [27–29].
As important economic animals worldwide, domestic goats (Capra hircus) are raised mainly for meat production. Thus, unveiling the molecular mechanisms underneath skeletal muscle formation and development is of vital interest. Muscle development is a complex process that requires the concerted expression and interaction of multiple factors . Several recent studies have indicated that lncRNAs play crucial roles in myogenic differentiation and myogenesis [31–35]. Nevertheless, currently the majority strategy for exploring molecular mechanisms underlying skeletal muscle growth and development in mammals [36–38] is detecting the expression and functions of coding genes like the MRF (myogenic regulatory factor) [39, 40], MEF2 (myocyte enhancer factor-2) [41, 42] families, and the paired box proteins , though high-throughput sequencing technologies is also employed to profile expression of miRNA and mRNA expression in goat [44, 45]. Therefore, information about skeletal muscle development-related lncRNAs is still limited especially in goats.
Here, we report the systematic identification and characterization of lncRNAs in fetal and postnatal goat skeletal muscle using an Illumina HiSeq 2500 platform. A total of 3981 lncRNA transcripts were identified and 577 of these transcripts were significantly differentially expressed in pairwise comparisons between RNA libraries of skeletal muscle at the different development stages. To the best of our knowledge, no other report on muscle lncRNAs and their biological functions in goat is currently available. Our results will provide a useful resource for better understanding the regulatory functions of lncRNAs in goat and for annotating the goat genome, as well as contribute to better comprehending skeletal muscle development in mammals.
Overview of RNA sequencing (RNA-seq)
To identify lncRNAs expressed in goat skeletal muscle development, we constructed 11 cDNA libraries (E45-1, E45-2, E45-3, E60-1, E60-2, E105-1, E105-2, E105-3, B3-1, B3-2, B3-3) from goat longissimus dorsi muscle samples at four developmental stages: three gestation stages at 45, 60, and 105 days of gestation (E45, E60, and E105), and one postnatal stage (B3). Three biological replicates for E45, E105, and B3, and two biological replicates for E60 were used. The libraries were sequenced using an Illumina HiSeq 2500 platform and 125 bp paired-end reads were generated. A total of 1,052,994,828 raw reads were generated in all 11 libraries. After discarding adaptor sequences and low-quality reads, we obtained 1,034,049,894 clean reads. The percentage of clean reads in each library ranged from 97.90 to 98.53 % (for details of the sequencing results see Additional file 1). We mapped the clean reads to the goat reference genome sequence (CHIR_1.0, NCBI). Approximately 75.20–86.60 % of the clean reads in all the libraries were mapped to the goat reference genome (Additional file 1). The mapped sequences in each library were assembled and a total of 56,710 unique assembled transcripts were obtained.
Identification of lncRNAs in goat skeletal muscle
Identification of differentially expressed lncRNAs
Enrichment analysis of nearest neighbor genes of the lncRNAs
To investigate the possible functions of the lncRNAs, we predicted the potential targets of lncRNAs in cis-regulatory relationships. We searched for protein-coding genes 10-kb upstream and downstream of all the identified lncRNAs. We found 1153 lncRNAs that were transcribed close to (<10 kb) 1455 protein-coding neighbors (Additional file 6). Gene Ontology (GO)  analysis of the cis lncRNA targets was performed to explore their functions. We found 88 GO terms that were significantly enriched (P < 0.05), and 12 of these terms were associated with regulation of gene expression. For example, the top 10 enriched terms included nucleotide binding, regulation of RNA metabolic process, DNA-dependent regulation of transcription, transcription regulator activity, and transcription factor activity (Additional file 7). These results suggest that one of the principal roles of lncRNAs may be transcriptional regulation of gene expression. Interestingly, we also found genes, including RBP4, PLN, MYLK2, RARA, CACNB4, NR2F2, CDK5 and PITX1, that were annotated with muscle development-related GO term, striated muscle tissue development (GO:0014706). These results suggest that muscle development may be regulated by the action of lncRNAs on these neighboring protein-coding genes. Pathway analysis  showed that the 1547 candidate cis target genes were enriched in 252 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, several of which were related to muscle development such as insulin signaling pathway, MAPK signaling pathway, TGF-beta signaling pathway, and PPAR signaling pathway (Additional file 7).
Enrichment analysis of co-expressed genes of lncRNAs
We also predicted the potential targets of lncRNAs in trans-regulatory relationships using co-expression analysis. A total of 288,020 interaction relationships (285,161 positive and 2859 negative correlations) were detected between 1747 lncRNA transcripts and 19,846 protein-coding transcripts that corresponded to 7718 protein-coding genes in the goat reference genome (Additional file 8). Functional analysis showed that the co-expressed genes were enriched in 446 GO terms (253 under biological process, 91 GO under cellular component, and 102 under molecular function) that encompassed a variety of biological processes (Additional file 9). Importantly, some of the terms were muscle development-related terms, including muscle cell development (GO:0055001), striated muscle cell development (GO:0055002), muscle contraction (GO:0006936), and muscle system process (GO:0003012). In addition, the co-expressed genes were enriched in 285 KEGG pathways, several of which were related to muscle development, including TGF-beta signaling, MAPK signaling, and PPAR signaling pathways (Additional file 9). These findings indicate that lncRNAs also regulate trans target genes.
Validation of differentially expressed lncRNAs
The identification and characterization of goat lncRNAs, particularly in fetal skeletal muscle development, have been very limited compared with lncRNAs in human [2, 48] and other model organisms [3, 49]. In goat skeletal muscle, the main focus has been on genes and miRNAs rather than on lncRNAs [44, 45, 50, 51]. In the present study, we identified a total of 3981 multiple exon lncRNAs in fetal and postnatal goat skeletal muscle. To the best of our knowledge, this is the first report to systematically identify lncRNAs from RNA-seq data during goat skeletal muscle development.
Non-coding and protein-coding genes are distinguished by their potential coding capability. In this study, we developed a highly stringent filtering pipeline (Fig. 1) to minimize the selection of false positive lncRNAs, which aimed to remove transcripts with evidence of protein-coding potential. We identified 3981 putative lncRNAs with high confidence across four development stages of goat skeletal muscle. In agreement with similar studies on different organisms, the identified putative lncRNAs had fewer exon numbers, shorter transcript lengths, and lower expression levels than protein-coding genes [11, 48, 49]. The number of putative lncRNAs detected in this study was more than that reported in previous studies in cattle and goat [9, 52]. Six randomly selected differentially expressed lncRNA transcripts were validated using qPCR, and the results were consistent with the results from the RNA-seq data. Together, these results confirmed that the identified lncRNAs were of high quality.
LncRNAs are a group of endogenous RNAs that function as regulators of gene expression, and are involved in developmental and physiological processes [23, 53, 54]. We detected 577 putative lncRNAs that were differentially expressed in pairwise comparisons between goat skeletal muscle at the different development stages. These lncRNAs may have specific biological roles in early muscle development in fetal goat. Skeletal muscle development from the fetal to the adult stage consists of a series of exquisitely regulated and orchestrated changes in the expression of many genes . In recent years, the roles of some lncRNAs in muscle biology have been reported. For example, the long intergenic ncRNA muscle differentiation linc-MD1 was the first muscle-specific lncRNA to be identified . Linc-MD1 is required for appropriate muscle differentiation, at least in part because it regulates the levels of myocyte enhancer factor 2C (MEF2C) and mastermind-like protein 1 (MAML1) by sponging endogenous miR-133 and miR-135 in the cytoplasm of muscle cells . In addition, the substantial disintegration of linc-MD1 in primary myoblasts of patients with Duchenne muscular dystrophy suggests that it is likely involved in the pathogenesis of this muscle disorder . Another study revealed that the lncRNA, lncMyoD, regulates skeletal muscle differentiation by blocking IMP2-mediated mRNA translation . Therefore, the differentially expressed lncRNAs reported here can be considered as important novel regulators of goat skeletal muscle biology.
Unlike miRNAs or proteins, the functions of ncRNAs cannot currently be inferred from their sequence or structure; therefore, in this study, we predicted the potential function of the detected lncRNAs using cis and trans methods. The cis nature of a lncRNA refers to its ability to act on a neighboring gene on the same allele from which it is transcribed; thus, this type of lncRNA commonly forms a feedback loop that regulates itself and its neighboring genes.
In the cis prediction, we searched for coding genes 10-kb upstream and downstream of all the identified lncRNAs. GO and KEGG analyses of the neighboring protein-coding genes revealed that major enriched pathways were associated with transition metal ion binding, nucleotide binding, zinc ion binding, regulation of RNA metabolic process, regulation of transcription, and transcription regulator activity. These results indicate the possible role of lncRNAs in transcriptional regulation of gene expression. Interestingly, we found some of the cis target protein-coding genes were involved in skeletal muscle tissue development (e.g., MYLK2, NR2F2, CDK5 and PITX1) (Additional file 6), implying that the corresponding lncRNAs play regulatory roles in skeletal muscle development. Several recent studies also indicated that lncRNAs were involved in cis-regulatory activity in muscle development; for example, the lncRNA Dum (developmental pluripotency-associated 2 (Dppa2) upstream binding muscle lncRNA) was identified in skeletal myoblast cells . Dum promotes myoblast differentiation and damage-induced muscle regeneration by silencing its neighboring gene, Dppa2, in cis through recruiting Dnmt1, Dnmt3a, and Dnmt3b . Similarly, a ChIP-seq study of the Yin Yang 1 (YY1) transcription factor, an important component of the PcG complex that negatively regulates myogenesis, identified a number of lncRNAs regulated by YY1 (YY1-associated muscle lncRNAs or Yams) . Among the Yams, Yam-1 displayed a cis effect on the expression of neighboring genes, including one that encodes miR-715, which targets and represses Wnt7b in skeletal muscle.
Many lncRNAs can also function in trans mode to target gene loci distant from where the lncRNAs are transcribed . In the co-expression analysis, we detected 1747 lncRNA transcripts that were related to protein-coding genes based on the expression correlation coefficient (r > 0.95 or < −0.95). GO enrichment analysis found that the enriched GO terms referred mainly to development process, transcriptional regulation, and biosynthetic process. Furthermore, a cluster of lncRNAs in the co-expression analysis often targeted protein-coding genes that were expressed specifically in muscle and were involved in muscle development (e.g., TNNT1, TNNT3, MYH1, MYH2, MyoG, and MYL3). This is an interesting observation, which indicates the functional complexity of lncRNAs and is worth investigating further. Mousavi et al.  found two lncRNAs in two enhancer regions of the MyoD gene that they named DDRRNA and CERNA, where DRR indicates distal regulatory regions and CE indicates core enhancer. The study showed that CERNA facilitated the occupancy of RNA polymerase II in cis by increasing chromatin accessibility, stimulating the expression of MyoD, while DDRRNA functions in trans to promote the expression of myogenin, a key member of the myogenic transcription factor family. More recently, Mueller et al.  identified a lncRNA transcribed upstream of MyoD named MUNC (MyoD upstream non-coding RNA), and demonstrated that one of the spliced isoforms of MUNC was DRRRNA. Consistent with the results of Mousavi et al. , experimentally decreasing MUNC expression blocked myoblast differentiation, further highlighting the role of enhancer-associated lncRNAs during myogenesis . These results suggest that lncRNAs exhibit regulatory function through cis-acting or trans-acting mechanisms in skeletal muscle biology and diseases.
All the studies mentioned above have demonstrated that lncRNAs are an integral part of the regulatory network of muscle biology. The present study provides evidence for the role of lncRNAs in skeletal muscle development in goat, which is a starting point for understanding the regulatory mechanisms in which they are involved. The identification of the lncRNAs has greatly improved the annotation of the goat reference genome. Further, we believe that the putative lncRNAs may contribute to a better understanding of the biological basis of regulatory interactions amongst mRNA, miRNA, and lncRNA.
We elucidated skeletal muscle lncRNA profiles of fetal and postnatal goats by RNA-seq analysis and identified and characterized caprine lncRNAs that may be involved in skeletal muscle development in goat. This study provides a catalog of goat muscle lncRNAs that will help in understanding their regulatory roles in goat muscle development. In future studies, we plan to investigate the functions of some of these lncRNAs to provide basic information that will add to the understanding of the regulatory mechanisms associated with goat muscle development at the molecular level.
The methods used in this study were performed in accordance with the guidelines of Good Experimental Practices adopted by the Institute of Animal Science (Sichuan Agricultural University, Chengdu, China). All surgical procedures involving goats were performed according to the approved protocols of the Biological Studies Animal Care and Use Committee, Sichuan Province, China.
Animal and tissue preparation
Jianzhou big-eared goats were used in this study. Three pregnant ewes at each developmental stage were subjected to caesarean section to collect the fetuses at 45, 60, 105 days of gestation, and three female lambs were collected at the third day after birth. Longissimus dorsi muscle samples were collected at these four developmental stages: three gestation stages (E45, E60, and E105) and one postnatal stage (B3). Three biological replicates for E45, E105, and B3, and two biological replicates for E60 were collected. The eleven samples were immediately frozen in liquid nitrogen for RNA extraction.
RNA extraction, library construction, and sequencing
Total RNA was isolated from the 11 libraries using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. RNA degradation and contamination were monitored on 1 % agarose gels. RNA purity was checked using a NanoPhotometer® spectrophotometer (IMPLEN, Los Angeles, CA, USA). RNA concentration was measured using a Qubit® RNA Assay Kit in a Qubit® 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA). RNA integrity was assessed using a RNA Nano 6000 Assay Kit in a Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). Only samples that had RNA Integrity Number (RIN) scores > 8 were used for sequencing. A total of 3 μg RNA per sample was used as input material for RNA sample preparation.
First, rRNA was removed using an Epicentre Ribo-zero rRNA Removal Kit (Epicentre, Madison, WI, USA), and the rRNA-free residue was obtained by ethanol precipitation. Subsequently, high strand-specificity libraries were generated using the rRNA-depleted RNA and a NEBNext Ultra Directional RNA Library Prep Kit for Illumina (NEB, Ipswich, MA, USA) following the manufacturer’s recommendations. Briefly, the rRNA-depleted RNA was fragmented using divalent cations under elevated temperature in NEBNext. First-strand cDNA was synthesized using random hexamer primers and M-MuLV reverse transcriptase (RNase H−). Subsequently, second-strand cDNA synthesis was performed using second-strand synthesis reaction buffer, DNA polymerase I, and RNase H. Remaining overhangs were converted into blunt ends by exonuclease/polymerase activity. After adenylation of the 3’ ends of the DNA fragments, NEBNext adaptors with hairpin loop structures were ligated to the fragments to prepare them for hybridization. To select cDNA fragments that are 150–200 bp in length, the fragments in each of the library were purified with an AMPure XP system (Beckman Coulter, Brea, CA, USA). Then 3 μl USER Enzyme (NEB, Ipswich, MA, USA) was used with size-selected, adaptor-ligated cDNA at 37 °C for 15 min followed by 5 min at 95 °C before PCR. The qPCRs were performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index (X) Primer. The PCR products were purified (AMPure XP system) and library quality was assessed on an Agilent Bioanalyzer 2100 system. Clustering of the index-coded samples was performed on a cBot Cluster Generation System using a TruSeq PE Cluster Kit v3-cBot-HS (Illumina, San Diego, CA, USA) according to the manufacturer’s instructions. After cluster generation, the libraries were sequenced on an Illumina HiSeq 2500 platform and 125-bp long paired-end reads were generated.
Clean data were obtained by removing reads containing adapters, reads containing over 10 % of poly(N), and low-quality reads (>50 % of the bases had Phred quality scores ≤ 10) from the raw data. The Phred score (Q20) and GC content of the clean data were calculated. All the downstream analyses were based on the high quality clean data. Goat reference genome and gene model annotation files were downloaded from NCBI database (CHIR_1.0, NCBI) . Index of the reference genome was built using Bowtie v2.0.6 [61, 62] and paired-end clean reads were aligned to the reference genome using TopHat v2.0.14 . The mapped reads from each library were assembled with Cufflinks v2.2.1 . Cufflinks was run with ‘min-frags-per-transfrag = 0’ and‘–library-type fr-firststrand’, and other parameters were set as default.
Filtering pipeline to identify multiple exon lncRNAs
We filtered the assembled novel transcripts from the different libraries to obtain putative lncRNAs following the steps in the pipeline (Fig. 1) as follows. (1) Single exon transcripts and transcripts < 200-bp long were removed. (2) The remaining transcripts that overlapped (>1 bp) with goat protein-coding genes were removed. (3) Transcripts that were likely to be assembly artifacts or PCR run-on fragments according to the class code annotated by cuffcompare  were removed. Among the cuffcompare classes, only transcripts annotated as “i”, “u”, and “x” representing novel intronic, intergenic, and antisense transcripts respectively, were retained. (4) The Coding Potential Calculator (CPC) , Coding-Non-Coding-Index (CNCI) , and phylogenetic codon substitution frequency (PhyloCSF)  tools were used to assess the coding potential of the remaining transcripts, and transcripts with CPC score > 0, CNCI score > 0, and PhyloCSF Max_score > 100 were removed. (5) The remaining transcripts that contained a known protein-coding domain were removed. To accomplish this, each transcript sequence was translated in all six reading frames, then HMMER was used to exclude transcripts with a translated protein sequence that had a significant hit in the Pfam (PfamA and PfamB) database, release 28.0 . (6) The remaining transcripts that belonged to known classes of small RNAs (including snRNA, tRNAs, and miRNAs) were removed using Rfam release 12.0 .
Cuffdiff 2.1.1  was used to calculate FPKM scores for the lncRNAs and coding genes in each library. Differentially expressed lncRNAs between any two libraries were identified by edgeR (release 3.2) . We used a P-adjust < 0.01 and an absolute value of the |log2(fold change)| > 2 as the threshold to evaluate statistical significance of lncRNA expression differences. We analyzed the clusters obtained by systematic analysis for all differentially expressed lncRNAs in eleven libraries using the Heatmaps software package in R .
Target gene prediction and functional enrichment analysis
Cis-acting lncRNAs target neighboring genes [74, 75]. We searched for coding genes 10-kb upstream and downstream of all the identified lncRNAs and then predicted their functional roles as follows. The names of the neighbor genes were used to form a gene list that was input into DAVID software  for GO analysis. A KEGG enrichment analysis of the predicted target genes was performed with KOBAS software  using a hypergeometric test. GO terms and KEGG pathways with P < 0.05 were considered significantly enriched.
We used the expression levels of the identified lncRNAs and the known protein-coding genes from the four different developmental stages to analyze the co-expression of lncRNAs and protein-coding genes. We calculated Pearson’s correlation coefficients between the expression levels of 1747 lncRNAs and 19,846 protein-coding transcripts with custom scripts (r > 0.95 or r < −0.95). Then, we performed a functional enrichment analysis of the candidate lncRNA target genes using DAVID and KOBAS software. Significance was calculated using the Expression Analysis Systematic Explorer (EASE) test method (P-value < 0.05 was considered significant).
Validation of differentially expressed lncRNAs by qPCR
Primers for the lncRNAs and internal control genes (Additional file 10) were designed using Primer-BLAST (http://www.ncbi.nlm.nih.gov/tools/primer-blast/). The expression levels of the lncRNAs were normalized to ACTB, YWHAZ, and HPRT1. Total RNA was converted to cDNA using a PrimeScript™ RT Reagent Kit with gDNA Eraser (TAKARA, Dalian, China), with oligo(dT) and random hexamer primers included in the kit. The qPCR was performed according to the SYBR Premix Ex Taq™ II instructions (TAKARA). The reaction volume contained 10 μl SYBR Premix Ex Taq™ II, 0.8 μl of 10 μM forward and reverse primers, 1.6 μl template cDNA, and dH2O to a final volume of 20 μl. The reactions were performed on a CFX96 Real-Time PCR System (Bio-Rad, CA, USA) as follows: 95 °C for 2 min, followed by 39 cycles of 95 °C for 10 s, and 10 s at the Tm indicated in Additional file 10: Table S8. Melting curve analysis was performed from 65 to 95 °C with increments of 0.5 °C. Amplifications were performed in triplicate for each sample. Relative gene expression levels were calculated using the 2-ΔΔCt method , and data were expressed as least square mean ± standard error of the mean (SEM).
All data were analyzed using one-way analysis of variance (ANOVA) to test homogeneity of variances via Levene’s test followed by Student’s t-test analyses in SAS software version 9.0 (SAS, Cary, North Carolina, USA). The significance level for the statistical analysis was P < 0.05.
We thank all contributors of the present study. We also thank Mr. Cheng Zufu and his colleagues at Gene Denovo Ltd., Co (Guangzhou, China) for assistance in data processing.
This work was supported by Sichuan Province Science and Technology Support Program (2014NZ0077, 2015NZ0112 and 16ZC2849).
Availability of data and materials
The raw data we obtained from RNA-seq are available in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) (accession number: SRR3567144). The transcript annotation information for the lncRNAs have been provided in form of a gtf file (Additional file 11).
SYZ conceived the research, performed data analysis, and drafted the manuscript. YD conducted the animal experiments and sampling. WZ conducted the qPCR validation. LL and JZG performed part of the data analysis. TZ and LJW participated in the experimental design and surgical processes. HPZ provided the experimental environment and edited the manuscript. All authors have read and approved the final manuscript.
The authors declare that they have no competing interests.
Consent for publication
Ethics approval and consent to participate
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- He C, Hu H, Wilson KD, Wu H, Feng J, Xia S, et al. Systematic characterization of long noncoding RNAs reveals the contrasting coordination of Cis- and trans-molecular regulation in human fetal and adult hearts. Circ Cardiovasc Genet. 2016;9(2):110–8.View ArticlePubMedGoogle Scholar
- Iyer MK, Niknafs YS, Malik R, Singhal U, Sahu A, Hosono Y, et al. The landscape of long noncoding RNAs in the human transcriptome. Nat Genet. 2015;47(3):199–208.View ArticlePubMedPubMed CentralGoogle Scholar
- Lv J, Liu H, Yu S, Liu H, Cui W, Gao Y, et al. Identification of 4438 novel lincRNAs involved in mouse pre-implantation embryonic development. Mol Genet Genomics. 2015;290(2):685–97.View ArticlePubMedGoogle Scholar
- Goff LA, Groff AF, Sauvageau M, Trayes-Gibson Z, Sanchez-Gomez DB, Morse M, et al. Spatiotemporal expression and transcriptional perturbations by long noncoding RNAs in the mouse brain. Proc Natl Acad Sci U S A. 2015;112(22):6855–62.View ArticlePubMedPubMed CentralGoogle Scholar
- Zhu JG, Shen YH, Liu HL, Liu M, Shen YQ, Kong XQ, et al. Long noncoding RNAs expression profile of the developing mouse heart. J Cell Biochem. 2014;115(5):910–8.View ArticlePubMedGoogle Scholar
- Lv J, Huang Z, Liu H, Liu H, Cui W, Li B, et al. Identification and characterization of long intergenic non-coding RNAs related to mouse liver development. Mol Genet Genomics. 2014;289(6):1225–35.View ArticlePubMedGoogle Scholar
- Liang M, Li W, Tian H, Hu T, Wang L, Lin Y, et al. Sequential expression of long noncoding RNA as mRNA gene expression in specific stages of mouse spermatogenesis. Sci Rep. 2014;4:5966.PubMedPubMed CentralGoogle Scholar
- 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.View ArticlePubMedPubMed CentralGoogle Scholar
- Billerey C, Boussaha M, Esquerre D, Rebours E, Djari A, Meersseman C, et al. Identification of large intergenic non-coding RNAs in bovine muscle using next-generation transcriptomic sequencing. BMC Genomics. 2014;15:499.View ArticlePubMedPubMed CentralGoogle Scholar
- Wang Y, Xue S, Liu X, Liu H, Hu T, Qiu X, et al. Analyses of long non-coding RNA and mRNA profiling using RNA sequencing during the pre-implantation phases in pig endometrium. Sci Rep. 2016;6:20238.View ArticlePubMedPubMed CentralGoogle Scholar
- Ran M, Chen B, Li Z, Wu M, Liu X, He C, et al. Systematic identification of long noncoding RNAs in immature and mature porcine testes. Biol Reprod. 2016;94(4):77.View ArticlePubMedGoogle Scholar
- Zhao Y, Li J, Liu H, Xi Y, Xue M, Liu W, et al. Dynamic transcriptome profiles of skeletal muscle tissue across 11 developmental stages for both Tongcheng and Yorkshire pigs. BMC Genomics. 2015;16:377.View ArticlePubMedPubMed CentralGoogle Scholar
- Zhao W, Mu Y, Ma L, Wang C, Tang Z, Yang S, et al. Systematic identification and characterization of long intergenic non-coding RNAs in fetal porcine skeletal muscle development. Sci Rep. 2015;5:8957.View ArticlePubMedPubMed CentralGoogle Scholar
- Bakhtiarizadeh MR, Hosseinpour B, Arefnezhad B, Shamabadi N, Salami SA. In silico prediction of long intergenic non-coding RNAs in sheep. Genome. 2016;59(4):263–75.View ArticlePubMedGoogle Scholar
- Chen YA, Aravin AA. Non-Coding RNAs in transcriptional regulation. Curr Mol Biol Rep. 2015;1(1):10–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Vance KW, Ponting CP. Transcriptional regulatory functions of nuclear long noncoding RNAs. Trends Genet. 2014;30(8):348–55.View ArticlePubMedPubMed CentralGoogle Scholar
- Bonasio R, Shiekhattar R. Regulation of transcription by long noncoding RNAs. Annu Rev Genet. 2014;48:433–55.View ArticlePubMedPubMed CentralGoogle Scholar
- Morlando M, Ballarino M, Fatica A, Bozzoni I. The role of long noncoding RNAs in the epigenetic control of gene expression. ChemMedChem. 2014;9(3):505–10.View ArticlePubMedGoogle Scholar
- Cao J. The functional role of long non-coding RNAs and epigenetics. Biol Proced Online. 2014;16:11.View ArticlePubMedPubMed CentralGoogle Scholar
- Mercer TR, Mattick JS. Structure and function of long noncoding RNAs in epigenetic regulation. Nat Struct Mol Biol. 2013;20(3):300–7.View ArticlePubMedGoogle Scholar
- Grote P, Herrmann BG. Long noncoding RNAs in organogenesis: making the difference. Trends Genet. 2015;31(6):329–35.View ArticlePubMedGoogle Scholar
- Mathieu EL, Belhocine M, Dao LT, Puthier D, Spicuglia S. Functions of lncRNA in development and diseases. Med Sci (Paris). 2014;30(8–9):790–6.View ArticleGoogle Scholar
- Dey BK, Mueller AC, Dutta A. Long non-coding RNAs as emerging regulators of differentiation, development, and disease. Transcription. 2014;5(4):e944014.View ArticlePubMedPubMed CentralGoogle Scholar
- 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.View ArticlePubMedPubMed CentralGoogle Scholar
- Fatica A, Bozzoni I. Long non-coding RNAs: new players in cell differentiation and development. Nat Rev Genet. 2014;15(1):7–21.View ArticlePubMedGoogle Scholar
- Dey BK, Pfeifer K, Dutta A. The H19 long noncoding RNA gives rise to microRNAs miR-675-3p and miR-675-5p to promote skeletal muscle differentiation and regeneration. Genes Dev. 2014;28(5):491–501.View ArticlePubMedPubMed CentralGoogle Scholar
- Tang JY, Lee JC, Chang YT, Hou MF, Huang HW, Liaw CC, et al. Long noncoding RNAs-related diseases, cancers, and drugs. Sci World J. 2013;2013:1–7.Google Scholar
- Ling H, Vincent K, Pichler M, Fodde R, Berindan-Neagoe I, Slack FJ, et al. Junk DNA and the long non-coding RNA twist in cancer genetics. Oncogene. 2015;34(39):5003–11.View ArticlePubMedPubMed CentralGoogle Scholar
- Liu XH, Sun M, Nie FQ, Ge YB, Zhang EB, Yin DD, et al. LncRNA HOTAIR functions as a competing endogenous RNA to regulate HER2 expression by sponging miR-331-3p in gastric cancer. Mol Cancer. 2014;13(1):92.View ArticlePubMedPubMed CentralGoogle Scholar
- Buckingham M. Myogenic progenitor cells and skeletal myogenesis in vertebrates. Curr Opin Genet Dev. 2006;16(5):525–32.View ArticlePubMedGoogle Scholar
- 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.View ArticlePubMedPubMed CentralGoogle Scholar
- 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.View ArticlePubMedGoogle Scholar
- Han X, Yang F, Cao H, Liang Z. Malat1 regulates serum response factor through miR-133 as a competing endogenous RNA in myogenesis. FASEB J. 2015;29(7):3054–64.View ArticlePubMedGoogle Scholar
- Ballarino M, Cazzella V, D’Andrea D, Grassi L, Bisceglie L, Cipriano A, et al. Novel long noncoding RNAs (lncRNAs) in myogenesis: a miR-31 overlapping lncRNA transcript controls myoblast differentiation. Mol Cell Biol. 2015;35(4):728–36.PubMedGoogle Scholar
- Legnini I, Morlando M, Mangiavacchi A, Fatica A, Bozzoni I. A feedforward regulatory loop between HuR and the long noncoding RNA linc-MD1 controls early phases of myogenesis. Mol Cell. 2014;53(3):506–14.View ArticlePubMedPubMed CentralGoogle Scholar
- Braun T, Gautel M. Transcriptional mechanisms regulating skeletal muscle differentiation, growth and homeostasis. Nat Rev Mol Cell Biol. 2011;12(6):349–61.View ArticlePubMedGoogle Scholar
- Buckingham M, Rigby PW. Gene regulatory networks and transcriptional mechanisms that control myogenesis. Dev Cell. 2014;28(3):225–38.View ArticlePubMedGoogle Scholar
- Eng D, Ma H-Y, Gross MK, Kioussi C. Gene networks during skeletal myogenesis. ISRN Dev Biol. 2013;2013:1–8.View ArticleGoogle Scholar
- Moncaut N, Rigby PW, Carvajal JJ. Dial M(RF) for myogenesis. FEBS J. 2013;280(17):3980–90.View ArticlePubMedGoogle Scholar
- Berkes CA, Tapscott SJ. MyoD and the transcriptional control of myogenesis. Semin Cell Dev Biol. 2005;16(4–5):585–95.View ArticlePubMedGoogle Scholar
- Naya FJ, Olson E. MEF2: a transcriptional target for signaling pathways controlling skeletal muscle growth and differentiation. Curr Opin Cell Biol. 1999;11(6):683–8.View ArticlePubMedGoogle Scholar
- Estrella NL, Desjardins CA, Nocco SE, Clark AL, Maksimenko Y, Naya FJ. MEF2 transcription factors regulate distinct gene programs in mammalian skeletal muscle differentiation. J Biol Chem. 2015;290(2):1256–68.View ArticlePubMedGoogle Scholar
- Buckingham M, Relaix F. The role of Pax genes in the development of tissues and organs: Pax3 and Pax7 regulate muscle progenitor cell functions. Annu Rev Cell Dev Biol. 2007;23:645–73.View ArticlePubMedGoogle Scholar
- Wang Y, Zhang C, Fang X, Zhao Y, Chen X, Sun J, et al. Identification and profiling of microRNAs and their target genes from developing caprine skeletal muscle. PLoS One. 2014;9(5):e96857.View ArticlePubMedPubMed CentralGoogle Scholar
- Wang YH, Zhang CL, Plath M, Fang XT, Lan XY, Zhou Y, et al. Global transcriptional profiling of longissimus thoracis muscle tissue in fetal and juvenile domestic goat using RNA sequencing. Anim Genet. 2015;46(6):655–65.View ArticlePubMedGoogle Scholar
- Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57.View ArticleGoogle Scholar
- Xie C, Mao X, Huang J, Ding Y, Wu J, Dong S, 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–322.View ArticlePubMedPubMed CentralGoogle Scholar
- 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.View ArticlePubMedPubMed CentralGoogle Scholar
- 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(3):577–91.View ArticlePubMedPubMed CentralGoogle Scholar
- Xu TS, Zhang XH, Gu LH, Zhou HL, Rong G, Sun WP. Identification and characterization of genes related to the development of skeletal muscle in the Hainan black goat. Biosci Biotechnol Biochem. 2012;76(2):238–44.View ArticlePubMedGoogle Scholar
- 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.View ArticlePubMedGoogle Scholar
- Ren H, Wang G, Chen L, Jiang J, Liu L, Li N, et al. Genome-wide analysis of long non-coding RNAs at early stage of skin pigmentation in goats (Capra hircus). BMC Genomics. 2016;17(1):67.View ArticlePubMedPubMed CentralGoogle Scholar
- Knoll M, Lodish HF, Sun L. Long non-coding RNAs as regulators of the endocrine system. Nat Rev Endocrinol. 2015;11(3):151–60.View ArticlePubMedPubMed CentralGoogle Scholar
- Ounzain S, Micheletti R, Beckmann T, Schroen B, Alexanian M, Pezzuto I, et al. Genome-wide profiling of the cardiac transcriptome after myocardial infarction identifies novel heart-specific long non-coding RNAs. Eur Heart J. 2015;36(6):353–368a.View ArticlePubMedGoogle Scholar
- Bismuth K, Relaix F. Genetic regulation of skeletal muscle development. Exp Cell Res. 2010;316(18):3081–6.View ArticlePubMedGoogle Scholar
- Gong C, Li Z, Ramanujan K, Clay I, Zhang Y, Lemire-Brachat S, et al. A long non-coding RNA, LncMyoD, regulates skeletal muscle differentiation by blocking IMP2-mediated mRNA translation. Dev Cell. 2015;34(2):181–91.View ArticlePubMedGoogle Scholar
- Lu L, Sun K, Chen X, Zhao Y, Wang L, Zhou L, et al. Genome-wide survey by ChIP-seq reveals YY1 regulation of lincRNAs in skeletal myogenesis. EMBO J. 2013;32(19):2575–88.View ArticlePubMedPubMed CentralGoogle Scholar
- Nie M, Deng ZL, Liu J, Wang DZ. Noncoding RNAs, emerging regulators of skeletal muscle development and diseases. Biomed Res Int. 2015;2015:1–17.Google Scholar
- Mousavi K, Zare H, Dell’Orso S, Grontved L, Gutierrez-Cruz G, Derfoul A, et al. eRNAs promote transcription by establishing chromatin accessibility at defined genomic loci. Mol Cell. 2013;51(5):606–17.View ArticlePubMedPubMed CentralGoogle Scholar
- Dong Y, Xie M, Jiang Y, Xiao N, Du X, Zhang W, et al. Sequencing and automated whole-genome optical mapping of the genome of a domestic goat (Capra hircus). Nat Biotechnol. 2013;31(2):135–41.View ArticlePubMedGoogle Scholar
- Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9(4):357–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10(3):R25.View ArticlePubMedPubMed CentralGoogle Scholar
- Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):R36.View ArticlePubMedPubMed CentralGoogle Scholar
- Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28(5):511–5.View ArticlePubMedPubMed CentralGoogle Scholar
- Roberts A, Pimentel H, Trapnell C, Pachter L. Identification of novel transcripts in annotated genomes using RNA-Seq. Bioinformatics. 2011;27(17):2325–9.View ArticlePubMedGoogle Scholar
- 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–349.View ArticlePubMedPubMed CentralGoogle Scholar
- 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.View ArticlePubMedPubMed CentralGoogle Scholar
- Lin MF, Jungreis I, Kellis M. PhyloCSF: a comparative genomics method to distinguish protein coding and non-coding regions. Bioinformatics. 2011;27(13):i275–282.View ArticlePubMedPubMed CentralGoogle Scholar
- Finn RD, Coggill P, Eberhardt RY, Eddy SR, Mistry J, Mitchell AL, et al. The Pfam protein families database: towards a more sustainable future. Nucleic Acids Res. 2016;44(D1):D279–285.View ArticlePubMedGoogle Scholar
- Nawrocki EP, Burge SW, Bateman A, Daub J, Eberhardt RY, Eddy SR, et al. Rfam 12.0: updates to the RNA families database. Nucleic Acids Res. 2015;43(Database issue):D130–137.View ArticlePubMedGoogle Scholar
- 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.View ArticlePubMedPubMed CentralGoogle Scholar
- Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.View ArticlePubMedGoogle Scholar
- Ploner A. Heatplus: heatmaps with row and/or column covariates and colored clusters. R package version 2.18.0. 2015.Google Scholar
- Ponjavic J, Oliver PL, Lunter G, Ponting CP. Genomic and transcriptional co-localization of protein-coding and long non-coding RNA pairs in the developing brain. PLoS Genet. 2009;5(8):e1000617.View ArticlePubMedPubMed CentralGoogle Scholar
- Orom UA, Derrien T, Beringer M, Gumireddy K, Gardini A, Bussotti G, et al. Long noncoding RNAs with enhancer-like function in human cells. Cell. 2010;143(1):46–58.View ArticlePubMedPubMed CentralGoogle Scholar
- 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–8.View ArticlePubMedGoogle Scholar