All animal experiments were performed according to protocols approved by the Institutional Animal Care and Use Committee at the South China Agricultural University (Guangzhou, People's Republic of China). All experiments were conducted in compliance with the regulations and guidelines established by that committee.

Fertilized eggs of a native Chinese yellow meat-type chicken known as Xinghua chickens were used in the this study. The Xinghua chicken is an indigenous chicken breed from Xinghua town, Fengkai County, Guangdong province, China is characterized by a relatively slow growth rate. The fertilized eggs were incubated in a humidified incubator at 37.5°C and 65% humidity. Skeletal leg muscles were resected at embryonic day 11 (E11), embryonic day 16 (E16), and 1 day after hatching (D1). All tissues were immediately frozen in liquid nitrogen and stored at −80°C until RNA extraction. The sex of the chicken was determined by PCR amplification of the CHD1 gene using sex-specific primers (Li Z. L. et al., 2012). Chickens with two bands of 600 bp and 450 bp were born as females, whereas those with one band of 600 bp were born as males. The sex-specific primer sequences are presented in Table 1.

Total RNA for RNA sequencing (RNA-seq) was isolated from six chicken leg muscle samples (two each from E11, E16, and D1), using Trizol Reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol. RNA quantity and quality were evaluated on an Agilent 2100 Bioanalyzer (Agilent Technologies, Waldbronn, Germany), and RNA integrity was further examined using agarose gel electrophoresis. Ribosomal RNA (rRNA) was removed from the total RNA using the Epicentre Ribo-Zero™ rRNA Removal Kit (Epicentre, Madison, Wisconsin, USA) following the manufacturer's instructions. High-throughput RNA-seq was performed on the Illumina Hiseq 2500 platform (Illumina, San Diego, CA, USA). The Illumina sequencing raw reads were cleaned by removing empty reads, adapter sequences, reads with over 10% N sequences, and low quality reads, in which the number of bases with a quality value Q ≤ 10 was >50%. In addition, rRNA reads were identified, by blasting against the rRNA database (http://www.arb-silva.de/) using Simple Object Access Protocol (SOAP) software, and removed from the dataset. The filtered reads were mapped to the chicken reference genome (ftp://ftp.ensembl.org/pub/release-73/fasta/gallus_gallus/dna/) using the Tophat2 software. The mapped reads were assembled and transcripts were constructed using the Cufflinks 2.0.2 software. The Ensembl, NCBI RefGene, and UCSC databases were chosen as annotation references for lncRNA analyses. The Coding Potential Calculator (CPC), Coding-Non-Coding Index (CNCI) and Phylogenetic codon substitution frequency (PhyloCSF) were applied to estimate coding ability of potential lncRNAs. An RNA length ≥ 200 nt, CPC score ≤ 0, CNCI score ≤ 0, and PhyloCSF score ≤ −20 were used to evaluate the coding potential of transcripts. The RefSeq and Ensembl databases were used for gene analysis. The expression levels of the transcripts were expressed as FPKM (fragments per kilobase of transcript per million mapped reads) values using the Cuffnorm program. Differentially expressed lncRNAs and mRNAs for each comparison (E11 vs. E16, E11 vs. D1, and E16 vs. D1) were identified using Cuffdiff (the differential expression analysis tool packaged with the Cufflinks software), with a q < 0.05 and |(fold change)| ≥ 2 as the cut-off points. The sequencing data obtained from RNA-seq were released to the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database under the accession number GSE91060.

Differentially expressed lncRNAs were selected for cis- and trans-target gene prediction, and these findings were integrated with differentially expressed mRNA data to improve the accuracy of target prediction. To classify lncRNA cis-target genes, we used the UCSC Genome Bioinformatics tool to identify differentially expressed genes (DEGs) located ~300 kb upstream and downstream of differentially expressed lncRNAs, and determined the potential that the lncRNA was cis-acting. To classify lncRNA trans-target genes, the BLAST (e < 1E − 5) software was used to assess the impact of lncRNA binding on complete mRNA molecules. The RNAplex program (−e −20) was then used to identify possible trans-target genes of the lncRNAs. Differentially expressed lncRNAs and their corresponding differentially expressed cis- and trans-target genes were used to construct lncRNA-gene interaction networks using the Cytoscape program.

Gene interaction network and biological process enrichment analyses for the cis- and trans-target genes of differentially expressed lncRNAs were performed using IPA (www.ingenuity.com). Ingenuity Pathway Analysis is a web-based software from Ingenuity Systems® (Qiagen, Redwood City, CA, USA) that enables investigators to analyze data derived from high-throughput sequencing. IPA provides a comprehensive database of known networks and pathways that are continuously being updated based on published work on gene functions and interactions. Each gene identifier was imported and mapped to its corresponding gene object in the Ingenuity Pathways Knowledge Base (IPKB) and superimposed onto a global molecular network, based on the information stored in the IPKB. IPA uses computational algorithms to identify local networks that are particularly enriched for the genes of interest. A network score that is used to rank networks is also provided. This numerical score, which is based on the hyper-geometric distribution, is calculated using a right-tailed Fisher's exact test and is represented as the negative log of the P-value. Functional enrichment analysis was also performed on these networks to understand the significance of the biological functions and/or disease phenotypes of the genes in the network. The predicted lncRNA target genes were entered into the Molecule Annotation System (MAS) 3.0 (http://bioinfo.capitalbio.com/mas3), which is based on the Kyoto Encyclopedia of Genes and Genomes database (KEGG) (CapitalBio, Beijing, China). The MAS uses GO data as an advantage to identify the molecular functions represented in the gene profile. The MAS 3.0 software was used to perform GO and KEGG analyses.

Single-strand cDNA was synthesized from 1 μg of total RNA using the PrimeScript™ RT reagent kit with gDNA Eraser (Takara, Osaka, Japan). The qPCR reactions were performed on a Bio-Rad CFX96 Real-Time Detection system to detect the RNA expression level of each gene, using the iTaq™ Universal SYBR® Green Supermix kit (Bio-Rad Laboratories Inc., Hercules, CA, USA), for which the chicken β-actin gene was used as an internal control. Data analysis was carried out with the 2^−ΔΔCt method as described previously (Livak and Schmittgen, 2001). Real-time qPCR primers were designed using the Premier Primer 5.0 software (Premier Biosoft International, Palo Alto, CA, USA). The qPCR primer sequences are presented in Table 1.

Quantitative expressions of the data are presented as mean ± S.E.M of three independent experiments. In each experiment, the assay was performed in triplicate. Statistical significance of the data was tested by performing one-sample or paired t-tests when required. When applicable, the specific types of tests and the P < 0.05 are indicated in the figure legends.

The Illumina Hiseq 2500 platform was used to perform RNA-seq for the six cDNA libraries, and 125 bp paired-end reads were generated. The number of raw reads generated from each library was exceeded 100 million reads (Table 2). After filtering the low quality reads, most of the clean reads still comprised more than 90% of the raw data, with the exception of the E16-2 group (85.13%) and D1-1 group (89.47%). More than 73.70% of the clean reads were perfectly mapped to the chicken reference genome. The unique mapped reads ranged from 95.62 to 96.54% of the total mapped reads (Table 2).

After a stringent filtering process to discard transcripts that did not have the characteristics of lncRNAs, RNA-seq yielded 1995 lncRNA transcripts. The length of the lncRNAs ranged from 204 to 32984 bp, with 30.78 and 25.31% of lncRNAs having a length of 204–1000 and 1000–2000 bp, respectively (Figure 1A). The lncRNA genes had an average length of 2941 bp and 2.9 exons (Figures 1A,B). The 1995 lncRNA transcripts were distributed among chicken chromosomes 1–28, 33, MT, Z, and W. Quite notably, 316 lncRNAs were located in the chromosome 1 (Figure 1C).

To investigate the key lncRNAs and mRNAs involved in chicken skeletal muscle development, RNA-seq (with rRNA removed) was performed to detect the differentially expressed lncRNAs (DE-lncRNAs) and DEGs at three stages of leg muscle development in the chicken (E11, E16, and D1). DE-lncRNAs and DEGs were identified by Cuffdiff analysis (q < 0.05; |(fold change)| ≥ 2). Among three comparisons, E11 vs. E16, E11 vs. D1, and E16 vs. D1, 129, 132, and 45 DE-lncRNAs were detected, respectively. Among these DE-lncRNAs, 26, 24, and 6 lncRNAs were uniquely expressed in one of the two samples in the E11 vs. E16 (Table S1), E11 vs. D1 (Table S2), and E16 vs. D1 (Table S3) comparisons, respectively. Of these DE-lncRNAs, six, including lnc00060568, lnc00057952, lnc00019955, lnc00064868, lnc00002354, and lnc00029488, were identified as differentially expressed among all three comparisons, whereas 61, 54, and 7 lncRNAs were found exclusively in the E11 vs. E16, E11 vs. D1, and E16 vs. D1 comparisons, respectively (Figure 2A).

Comparisons of gene expression levels among the three developmental stages revealed that there were 1798 DEGs (q < 0.05; |(fold change)| ≥ 2) between E11 and E16 (E11 vs. E16), 3072 DEGs between E11 and D1 (E11 vs. D1), and 1211 DEGs between E16 and D1 (E16 vs. D1). Among these DEGs, 27, 50, and 5 genes were uniquely expressed in one of the two samples in the E11 vs. E16 (Table S4), E11 vs. D1 (Table S5), and E16 vs. D1 (Table S6) comparisons, respectively. Of these DEGs, 386 genes were detected in all three comparisons, whereas 281, 1019, and 110 genes were exclusively differentially expressed in E11 vs. E16, E11 vs. D1, and E16 vs. D1 comparisons, respectively (Figure 2B).

To address the question of how lncRNAs function in concert with their target genes (mRNAs) to regulate chicken muscle development, and to identify the key molecular players in this process, cis- and trans-targets of the differentially expressed lncRNAs were predicted, and the possible regulatory networks of interactions between lncRNAs and their target genes (mRNAs) were constructed. Previous studies have confirmed that lncRNAs regulate the expression of neighboring protein-coding genes through cis-acting mechanisms (Bu et al., 2012; Han et al., 2012). In addition, lncRNAs can regulate the expression of genes that are located on other chromosomes; such regulation is called trans-acting (Han et al., 2012). In the present study, bioinformatics analysis was used to predict potential cis- and trans-target genes of DE-lncRNAs, and further construct a lncRNA-gene interaction network between the DE-lncRNAs and their corresponding differentially expressed cis- and trans-target genes. Construction of the lncRNA-gene network consisted of three steps (Figure 3). Firstly, the genes within 300 kb upstream and downstream of the lncRNA were identified using the UCSC Genome Browser. Among these genes, those that were coexpressed with the lncRNA were classified as potential lncRNA cis-target genes. Secondly, BLAST and RNAplex software were used to identify trans-target genes. To narrow the field and improve the accuracy of target prediction, only those putative trans-target genes that were coexpressed with the lncRNAs were classified as trans-lncRNA target genes. Thirdly, Cytoscape was used to construct a lncRNA-gene interaction network based on the DE-lncRNAs and their corresponding differentially expressed cis- and trans-target genes. For the E11 vs. E16 comparison, the lncRNA-gene interaction network comprised 211 network nodes and 148 lncRNA-gene connections among 83 lncRNAs and 128 protein-coding genes (Figure 4, Table S7). For the E11 vs. D1 comparison, the lncRNA-gene interaction network comprised 305 network nodes and 242 lncRNA-gene connections among 102 lncRNAs and 203 protein-coding genes (Figure 5, Table S8). For the E16 vs. D1 comparison, the lncRNA-gene interaction network comprised 46 network nodes and 28 lncRNA-gene connections among 21 lncRNAs and 25 protein-coding genes (Figure 6, Table S9).

We then analyzed the expression correlation between the network lncRNAs and their corresponding target genes. In the network constructed from DE-lncRNAs and target genes identified from the E11 vs. E16 comparison, 92 lncRNA-gene connections were positively correlated, whereas another 56 connections were negatively correlated (Figure 4, Table S7). For the E11 vs. D1 comparison, 138 lncRNA-gene connections were positively correlated and 104 connections were negatively correlated (Figure 5, Table S8). For the E16 vs. D1 comparison, 16 lncRNA-gene connections were positively correlated and 12 connections were negatively correlated (Figure 6, Table S9). Directionality analysis showed that for all three contrast networks, the number of positive correlations between lncRNA-gene pairs was higher than the number of negative correlations.

In total, 128, 203, and 25 DEGs that were also targets of DE-lncRNAs belonged to the E11 vs. E16, E11 vs. D1, and E16 vs. D1 comparison lncRNA-gene networks, respectively. To identify the functions of the DE-lncRNAs target genes and their potential network connections, we used IPA to determine which gene networks might be affected by these DE-lncRNAs target genes. For the E11 vs. E16 comparison, we identified a total of ten networks, and the top four were related to cancer, organismal injury and abnormalities, cell cycle, developmental disorders, hereditary disorders, metabolic disease, reproductive system disease, cellular development, cellular growth and proliferation, and connective tissue development and function (Figures 7A–D, Table 1). In total, 23 (Figure 7A), 22 (Figure 7B), 20 (Figure 7C), and 15 (Figure 7D) DEGs belonged to these top four gene-gene networks. For the E11 vs. D1 comparison, 16 networks were identified and the top two were related to cell morphology, organ morphology, reproductive system development and function, embryonic development, organismal development, and tissue development (Figures 7E,F, Table 1). A total of 22 (Figure 7E), and 19 (Figure 7F) DEGs belonged to these top two gene-gene networks. In addition, for the E16 vs. D1 comparison, we identified a total of four networks, and the top two were related to nervous system development and function, cellular development, cellular growth and proliferation, developmental disorders, hereditary disorders, and ophthalmic disease (Figures 7G,H, Table 1). A total of 11 (Figure 7G), and 10 (Figure 7H) DEGs belonged to these top two gene-gene networks.

We then conducted GO and KEGG pathway analyses to gain insight into the functions of DE-lncRNA target genes. GO analysis was done to uncover enriched (q < 0.01) biological process terms for each contrast. A total of 107 enriched GO terms, including 43 for E11 vs. E16 (Figure S1A, Table S10); 51 for E11 vs. D1 (Figure S1B, Table S11); and 13 for E16 vs. D1 (Figure S1C, Table S12), were identified. Among these terms, 71 were unique GO terms that each appeared only once among all three comparisons. Many of the GO terms that were enriched in more than one comparison were related to muscle maintenance, lysine transport, ornithine transport, arginine transport, and peptide biosynthesis. A total of seven DE-lncRNA target genes involved in muscle maintenance, embryonic skeletal joint morphogenesis, post-embryonic development and organ growth were identified in the E11 vs. E16, and E11 vs. D1 comparisons, including well-known genes that affect muscle growth, such as dystrophin (DMD) and delta like non-canonical notch ligand 1 (DLK1).

KEGG pathway analysis of DE-lncRNAs target genes was done to identify pathways that were enriched in DE-lncRNA target genes. For the E11 vs. E16, E11 vs. D1, and E16 vs. D1 comparisons, 18 (Table S13), 17 (Table S14), and one (Table S15) enriched pathways (q < 0.001) were detected, respectively. For the E11 vs. E16 comparison, 11 DE-lncRNAs target genes played a role in the most significantly enriched pathways, including fatty acid metabolism, valine, leucine and isoleucine degradation, fatty acid elongation in mitochondria, and benzoate degradation via CoA ligation. In contrast, for the E11 vs. D1 comparison, the mitogen-activated protein kinase (MAPK) signaling pathway, mismatch repair, beta-Alanine metabolism, and chondroitin sulfate biosynthesis were identified as the most significantly enriched pathways, with 13 DE-lncRNAs target genes being associated with these pathways.

To confirm the accuracy and reproducibility of the DE-lncRNA and DEG expression levels obtained from RNA-seq, four well-known DE-lncRNA target genes that affect muscle development, including the TEA domain transcription factor 4 (TEAD4) (Jacquemin et al., 1996), DMD (Ghahramani et al., 2010), fibroblast growth factor 13 (FGF13) (Itoh and Ornitz, 2008), and DLK1 (Andersen et al., 2013), and their corresponding lncRNA regulators, lnc00003323, lnc00005738, lnc00023378, and lnc00073545, were selected for real-time qRT-PCR validation. The qRT-PCR primers for these lncRNAs and genes are shown in Table 1. All four lncRNAs and their target genes had similar expression patterns in comparison to the RNA-seq data (Figure 8), indicating the reliability of our RNA-seq data.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.