LncRNAs and their regulatory networks in breast muscle tissue of Chinese Gushi chicken during postnatal late development


 Background: Chicken skeletal muscle is an important economic product. Late stages of chicken development are the main period that affects meat production. LncRNAs play important roles in controlling the epigenetic process of growth and development. However, studies on the role of lncRNAs in the late stages of chicken breast muscle development are still lacking. In this study, in order to study the expression characteristics of lncRNAs during chicken muscle development, 12 cDNA libraries were constructed from Gushi chicken breast muscle samples from 6-, 14-, 22-, and 30-week-old chickens.Results: Therefore, a total of 1,252 new lncRNAs and 1,376 annotated lncRNAs were identified. Furthermore, 53, 61, 50, 153, 117, and 78 DE-lncRNAs were found in the W14 vs. W6, W22 vs. W14, W22 vs. W6, W30 vs. W6, W30 vs. W14, and W30 vs. W22 comparison groups, respectively. After GO enrichment analysis of the DE-lncRNAs, several muscle development-related GO terms were found in the W22 vs. W14 comparison group. Furthermore, it was found that the MAPK signaling pathway was one of the most frequently enriched pathways in the different comparison groups. In addition, 12 common target DE-miRNAs of DE-lncRNAs were found in different comparison groups, some of which were muscle-specific miRNAs, such as gga-miR-206, gga-miR-1a-3p, and miR-133a-3p. Interestingly, the precursors of four newly identified miRNAs were found to be homologous to lncRNAs. Furthermore, we found some ceRNA networks associated with muscle development-related GO terms. For example, the ceRNA networks containing the DYNLL2 gene with 12 lncRNAs that targeted 2 miRNAs. In addition, we also constructed PPI networks, such as IGF-I-EGF and FZD6-WNT11, etc.Conclusions: This study for the first time revealed the dynamic changes of lncRNA expression in Gushi chicken breast muscle at different periods, and reveal the MAPK signaling pathway play a vital role in muscle development. Furthermore, MEF2C and its target lncRNA may be involved in muscle regulation through the MAPK signaling pathway. This research provided valuable resources for elucidating post-transcriptional regulatory mechanisms to promote the development of chicken breast muscles after hatching.

exons, and fewer ORF numbers than mRNAs ( Figure 1D-F) and lower average transcript abundance ( Figure 1G).

Characteristics of differentially expressed lncRNAs
To gain insight into the key lncRNAs involved in chicken breast muscle development, we  Figure S1). The clustering results showed that the intragroup repeats of each group were clustered together, indicating that the intragroup differences were smaller than the intergroup differences, which proved that the data were reliable. In addition, we found that 22 weeks clustered close to 14 weeks, followed by 6 weeks, and the farthest distance was from 30 weeks. LncRNAs are usually expressed in a tissue-and developmental stage-specific manner. During muscle development, lncRNA expression levels showed dynamic global changes. Moreover, after identifying DE-lncRNAs, we analyzed the chromosome distribution information of the DE-lncRNAs and found that DE-lncRNAs were distributed in almost all chromosomes but not in chromosomes 22 and 27, and the largest number was found in chromosome 1 ( Figure S2). In addition, we selected several lncRNAs for data validation. The lncRNA expression level was determined and showed a similar pattern compared to the RNA-seq data (Figure 3), thus indicating that the RNA-seq data were authentic.
To investigate the possible functions of the DE-lncRNAs in breast muscle between the different developmental stages, we conducted Gene Ontology (GO, http://www.geneontology.org/) enrichment analysis to uncover the enriched biological process terms associated with DE-lncRNAs target DEGs for each comparison group. Only the top 20 GO terms for the W14 vs. W6, W22 vs. W14, and W30 vs. W22 comparison groups are shown in Figure 4 and Figure S3. The cis-targets of all lncRNAs were predicted with a 100 kb upstream and downstream range. The GO enrichment analysis of the cistargets of lncRNAs showed that only one growth-and development-related GO term called positive regulation of embryonic development, was found in the W22 vs. W14 comparison group ( Figure 4A). In addition, we predicted the regulation in trans between lncRNAs and genes by the Pearson correlation coefficient r > 0.95. In the GO analysis of the transtargets of lncRNAs, we only found several muscle development-related GO terms in the To further understand how DE-lncRNA-targeted DEGs play roles in regulating chicken muscle development, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/) pathway analysis for each comparison. In the KEGG pathway analysis of the cis-targets of lncRNAs, the phagosome pathway was identified as the most significantly enriched pathway for the W14 vs. W6 comparison group ( Figure S4A). Furthermore, for the W22 vs. W14 comparison group, the endocytosis pathway was identified as the most significantly enriched pathway ( Figure S4B). Additionally, the focal adhesion and cytokine-cytokine receptor interaction pathways were the top two pathways for the W30 vs. W22 comparison group ( Figure S4C). Moreover, the KEGG pathway analysis of the trans-targets of lncRNAs showed that the propanoate metabolism pathway and fatty acid metabolism pathway were the top two pathways for the W14 vs. W6 comparison group ( Figure S4D). In the W22 vs. W14 comparison group, the two top pathways were arginine and proline metabolism and the MAPK signaling pathway ( Figure 5A). In addition, for the W30 vs. W22 comparison group, there were two top pathway terms, the MAPK signaling pathway and the regulation of actin cytoskeleton pathway ( Figure 5B). We found that the MAPK signaling pathway was one of the most frequently enriched pathways in both the W22 vs. W14 and W30 vs. W22 comparison groups.

Interactions between lncRNAs and mRNAs during breast muscle development
To explore how lncRNAs interact with their target genes to regulate chicken muscle development and to identify key molecular players in the process, we first predicted the cis-and trans-targets of DE-lncRNAs and then constructed the regulatory networks between DE-lncRNAs and their target genes. A total of 13,460 cis-regulatory interaction relationships were detected between 2,309 lncRNAs and 7,783 mRNAs (Table S2). In addition, 13,343 trans-regulatory interaction relationships were detected between 733 lncRNAs and 2,190 mRNAs (Table S3) we found a total of 10 interaction relationships between 3 genes and 10 lncRNAs ( Figure   6A). In addition, in the networks of trans-target DEGs of DE-lncRNAs of the muscle development-related GO terms we found 7 interaction networks between 3 genes and 7 lncRNAs ( Figure 6B). Furthermore, we also generated the lncRNA-mRNA networks of the frequently enriched MAPK signaling pathway, which had a total of 25 interaction relationships between 11 genes and 25 cis-regulating lncRNAs and 27 interaction relationships between 8 genes and 17 trans-regulating lncRNAs ( Figure 6C). Interestingly, we found that the networks contained MEF2C and its targeting lncRNAs (ALDBGALT0000008862, ALDBGALT0000008865, LNC_001247, etc.) were not only in the muscle development-related GO terms but also in the MAPK signaling pathway.

LncRNA-miRNA interactions
In our previous study [17], we found 388 known miRNAs and 31 novel miRNAs. To explore the interactions between lncRNAs and miRNAs, we predicted the target relationship between lncRNAs and miRNAs in different comparison groups (Table S4). Twelve common target DE-miRNAs of DE-lncRNAs were found in different comparison groups. It is important that some of them were muscle-specific miRNAs, such as gga-miR-206, gga-miR-1a-3p, and miR-133a-3p. Only the lncRNAs (FPKM>1) of these miRNA targets are shown in Figure 7. Then, we predicted the pre-miRNAs with homology to lncRNAs. Unexpectedly, the precursors of four newly identified miRNAs were found to be homologous to lncRNAs, and the precursors were temporarily named gga-miR-N1, gga-miR-N2, gga-miR-N3 and gga-miR-N4 ( Figure 8, Table S5), for example, the pre-miRNA of gga-miR-N1 exactly matches ALDBGALT0000008009 at its position from 309 to 374, and the pre-miRNA of gga-miR-N2 exactly matches lnc_000010 at its position from 1527 to 1587, and it also matches lnc_000011 from 1101 to 1161. These lncRNAs may form miRNA precursors through intracellular shearing and then could be processed to generate specific miRNAs that regulate the expression of target genes.
Among those, myosin heavy polypeptide 11 (MYH11) is related to muscle cell differentiation, which was involved in 8 ceRNA networks containing gga-miR-194 and 8 target genes, and only some of the genes associated with muscle development-related GO terms and the MAPK signaling pathway and their corresponding lncRNAs are shown in Fig. 6. Interestingly, we found that the networks containing MEF2C and its targeting lncRNAs were not only in the muscle development-related GO terms but also in the MAPK signaling pathway. Research has shown that MEF2 factors are involved in muscle differentiation, and MyoD and MEF2 family members interact to activate transcription and myogenesis [26]. Thus, we speculate that these lncRNAs, such as ALDBGALT0000008862, ALDBGALT0000008865, and LNC_001247, are involved in muscle differentiation.
LncRNAs have significant similarities to classic mRNAs during transcription, so miRNAs can not only target mRNAs but also regulate target lncRNAs and reduce their structural and functional stability [27]. It has been speculated that miRNA-192 and miRNA-204 directly suppress lncRNA HOTTIP and interrupt GLS1-mediated glutaminolysis in hepatocellular carcinoma [28]. Interestingly, miRNAs can also enhance lncRNA expression through a number of mechanisms [29]. miR-140 binding sites were identified in NEAT1, and the authors found that mature miR-140 could physically interact with NEAT1 in the nucleus, thereby promoting the expression of NEAT1 [30]. The predicted miRNA-lncRNA networks are shown in Fig. 7 and suggest that miRNAs may play a role in promoting or inhibiting lncRNA expression. Interestingly, several muscle-specific miRNAs, such as gga-miR-206, gga-miR-1a-3p, and miR-133a-3p have been identified [31]. Therefore, we speculate that the target lncRNAs of these miRNAs may also be involved in the muscle development process. In addition, lncRNAs play many complex functions in gene expression and signal transduction and can also be used as precursors of miRNAs. LncRNAs can produce specific miRNAs through intracellular RNA splicing to affect the posttranscriptional regulation of mRNAs. It has been found that the H19 RNA is a miRNA precursor and generates the exonic microRNA miR-675 [32]. In this study, we found 5 pairs of homologous lncRNAs and pre-miRNAs, including ALDBGAL0000008009 with gga-miR-N1 and lnc_000010 and lnc_000011 with gga-miR-N2. We believe that the miRNAs produced by these lncRNAs play a functional role.
There is increasing evidence that supports the ceRNA hypothesis that lncRNAs regulate target genes by competitively adsorbing miRNAs [33]. For instance, lncRNA-Unigene56159 can directly bind to miR-140-5p and act as a ceRNA of miR-140-5p to regulate the expression of its target gene Slug, thereby affecting HCC cell migration and invasion [34].
In thyroid cancer, lncRNA-Gas5 regulates gene of phosphate and tension homology deleted on chromsome ten (PTEN) expression through a ceRNA mechanism as a sponge for miR-222-3p [36]. To fully identify how lncRNA-associated ceRNA networks affect breast muscle development in different developmental stages, we predicted and successfully constructed a network of ceRNAs associated with DE-lncRNAs in different comparison groups. For example, in the W14 vs. W6 comparison group, there were ceRNA networks containing the DYNLL2 gene with 12 lncRNAs that targeted 2 miRNAs, gga-miR-148a-3p and gga-miR-130b-3p (Fig. 9A). DYNLL2 is involved in a variety of cellular processes and interacts with myosin 5a (myo5a) to participate in cargo transport [37]. Therefore, we speculate that the above 12 lncRNAs may regulate DYNLL2 by adsorbing miR-148a-3p and gga-miR-130b-3 to regulate muscle development. Furthermore, the ceRNA network analysis showed that ANKRD1 is involved in 13 ceRNA networks containing two miRNAs (miR-148a-3p and miR-10b-5p) and 12 lncRNAs in the W14 vs. W6 comparison group (Fig. 9A). Moreover, ANKRD1 was involved in 8 ceRNA networks containing gga-miR-92-3p with 8 lncRNAs in the W30 vs.
W22 comparison group (Fig. 9C). Studies have shown that ANKRD1 can induce gene expression in cultured skeletal muscle cells and trigger signaling via myogenic regulators (MRFs) during myogenesis [38]. Therefore, we think that the above 20 lncRNAs can also play roles in skeletal muscle cells to regulate the ANKRD1 gene by adsorbing miR-148a-3p, miR-10b-5p, and gga-miR-126-5p.
In addition, we also demonstrated the interactions between DE-lncRNA target genes in Gushi chicken breast development. The interaction between IGF-I and epidermal growth factor (EGF) was identified in the W14 vs. W6 comparison group. It has been reported that the expression of the IGF gene is enhanced during muscle hypertrophy and that locally produced IGF-I may play roles in skeletal muscle growth [39]. Heparin-binding epidermal growth factor (HB-EGF) promotes airway smooth muscle (ASM) cell proliferation through the MAPK pathway [40]. Thus, we surmise that IGF-I and EGF can interact to play a role in skeletal muscle cell proliferation. Moreover, the PPI networks based on the W22 vs. Therefore, the interactions between these genes may eventually affect muscle development, and the above results indicate that there are complex intergenic interactions in the development of chicken breast muscle.

Conclusions
In conclusion, we first described the lncRNA profile of Gushi chicken breast muscle development at 6, 14, 22, and 30 weeks. In the W22 vs. W14 comparison groups, some GO terms related to muscle development were found, further indicating that between 14 and 22 weeks, changes in the expression of some crucial lncRNAs and their target genes affected the growth and development of chicken breast muscles during these important stages. In addition, the MAPK signaling pathway has been found to play a vital role in muscle development. MEF2C and its target lncRNA, such as ALDBGALT0000008862, ALDBGALT0000008865, and LNC_001247 may be involved in muscle regulation through the MAPK signaling pathway. These results provide insight and valuable resources for further research on the molecular mechanism of postpartum skeletal muscle development.

Animals and samples preparation
The experimental animals were Gushi chickens which from the Animal Center of Henan Agricultural University in this study. A total of 300 one-day-old female Gushi chickens were raised in cages of the same environment with standard conditions for pure breeding conservation of Gushi chickens. In this study, three healthy chickens were randomly selected at 6 weeks, 14 weeks, 22 weeks, and 30 weeks of age, respectively. Therefore, twelve chickens were used in this study, our sample size is sufficient, and the remaining healthy chickens are still used for the pure breeding conservation of Gushi chickens.
These chickens were fed a two-stage feeding in which 18.5% crude protein and 12.35

Identification of lncRNAs
The combined transcript sets were screened for lncRNAs, and the lncRNA screening process was divided into the following five steps: the first step was to select the transcripts with exon number ≥ 2; the second step was to select the transcripts with transcript length > 200 bp; the third step was to use Cuffcompare software to filter out the transcripts that overlapped with exon regions annotated in the database, and the lncRNAs in the database that overlapped with the exon regions of splicing transcripts were annotated as lncRNAs in subsequent analyses; the fourth step was to calculate the expression of each transcript by Cuffquant, whereby transcripts with FPKM (fragments per kilobase of exons per million mapped fragments) ≥ 0.5 expression level were selected; and the fifth step was to screen the protein coding potential, in which the four databases CPC, PFAM, PhyloCSF and CNCI were used to remove potential protein coding transcripts.
In this study, the resulting transcripts with no protein coding potential in the software analyses resulted in the lncRNA dataset, and the transcripts that were predicted to have protein coding potential by at least one coding potential prediction software were set as TUCPs.

Differential expression analysis
Cuffdiff (v2.1.1) was used to calculate the FPKM of lncRNAs in each sample, and cuffdiff used a model based on a negative binomial distribution to provide statistical procedures for determining differential expression in digital transcript or gene expression data [43].
Based on Illumina sequencing data, FPKM values were used to assess the expression levels of lncRNAs in the libraries constructed from breast muscle. The FC for each lncRNA between two discretionary stages was calculated according to comparisons among four comparison groups (W6, W14, W22, and W30). Differential expression analysis among the four groups was performed using the DESeq R package (1.8.3), and lncRNAs and genes found by DESeq with q-value < 0.05 and |FC|≥1.7 were considered to be differentially expressed.

Cis-and trans-targeting analyses
Differentially expressed lncRNAs were selected for cis-and trans-target gene predictions and were integrated with differentially expressed gene data to improve the veracity of target prediction. In the present study, DEGs located ∼100 kb upstream and downstream of DE-lncRNAs were classified as cis-acting target genes. In addition, we predicted the regulation in trans between lncRNAs and genes by the Pearson correlation coefficient, r > 0.95. The relationship between lncRNAs and their target genes was demonstrated by