Identiﬁcation and Characterization of lncRNAs Expression Proﬁle Related to Goat Skeletal Muscle at Different Development Stages

Simple Summary: The skeletal muscle growth and development affects the production of goat meat. Here, we performed transcriptome sequencing of longissimus dorsi muscle of goats at two development stages (namely, 1-month-old and 9-month-old) by RNA sequencing, assessing the lncRNA expression proﬁle during skeletal muscle development. This study identiﬁed regulatory lncRNAs and pathways related to the development of skeletal muscle, providing a reference for future study on the molecular mechanism that regulates the skeletal growth and development. Abstract: LncRNAs are essential for regulating skeletal muscle. However, the expression proﬁle and function of lncRNAs in goat muscle remains unclear. Here, an average of ~14.58 Gb high-quality reads were obtained from longissimus dorsi tissues of 1-month-old (n = 3) and 9-month-old (n = 3) Wu’an black goats using RNA sequencing. Of a total of 3441 lncRNAs, 1281 were lincRNAs, 805 were antisense lncRNAs, and 1355 were sense_overlapping lncRNAs. These lncRNAs shared some properties with goats, such as fewer exons, shorter transcript, and open reading frames (ORFs) length. Among them, 36 differentially expressed lncRNAs (DE lncRNA) were identiﬁed, and then 10 random lncRNAs were validated by RT-qPCR. Furthermore, 30 DE lncRNAs were neighboring 71 mRNAs and several genes were functionally enriched in muscle development-related pathways, such as APC , IFRD1 , NKX2-5 , and others. Additionally, 36 DE lncRNAs and 2684 mRNAs were included in co-expression interactions. A lncRNA-miRNA-mRNA network containing 4 lncRNAs, 3 miRNAs, and 8 mRNAs was ﬁnally constructed, of which XR_001296113.2 might regulate PDLIM7 expression by interaction with chi-miR-1296 to affect skeletal muscle development. This study revealed the expression proﬁle of goat lncRNAs for further investigative studies and provides a fuller understanding of skeletal muscle development.


Introduction
As an important economically farm animal, the goat is raised for the utilization of meat, cashmere, and milk. With people's living standards having improved, the demand for goat meat has gradually increased. As a result, the 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 transcription have provided a series of valuable candidate genes that regulate muscle growth in goats [3,4]. Thus, uncovering the genetic mechanisms underneath muscle growth could help us to improve the meat

Classification and Differential Expression analysis
The potential lncRNAs were classified based on their association with the annotated protein-coding genes. In this study, the association between the potential lncRNAs and the annotated protein-coding genes was compared with cuffmerge [28]. Then, the classification and characteristics of lncRNAs were analyzed. StringTie was selected to analyze the expression levels of lncRNAs by calculating fragments per kilobase million (FPKM). The corrected p-value < 0.05 and an absolute value of the |log2FoldChange| ≥ 2 were set as the threshold to assess statistically significant differences of lncRNA expression.

Co-location and Co-expression analysis, and Functional Annotation Analysis
LncRNAs exert cis-regulatory effects on their co-localized genes [13]. However, lncRNAs regulate the expression of genes located on other chromosomes through a transacting mechanism [29]. To assess the potential function of DE lncRNAs, we predicted the possible target genes that were co-located and co-expressed with lncRNAs. In this study, coding genes located 100 kb/100 kb of DE lncRNA were classified as putative cis-acting targets. The trans-acting correlations between lncRNAs and genes (co-expression) were calculated with the Pearson's correlation coefficients method (|r| > 0.95 and p < 0.05), as previously reported [22,30]. The cis-and trans-acting putative target genes were annotated by Gene Ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) to explore the potential functions 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
To validate the reliability of the RNA-seq data, 10 DE lncRNAs were randomly selected to verify their expression patterns in skeletal muscle of goats at different development stages. Primers for the 10 lncRNAs and endogenous reference gene (Table 1) were designed by Primer 5, and the goat GAPDH was selected as the endogenous reference. Briefly, total RNA was extracted from the longissimus dorsi muscle samples. For RT-qPCR analysis, 1 μg of the total RNA was reverse transcribed using PrimeScript™ RT reagent kit (Takara, Beijing, China) according to the manufacturer's protocol. Then, the qPCR was performed by a RocheLight Cycler ® 480 Ⅱ system with SYBR Green qPCR Mix kit as previous study [31]. The expression levels were analyzed using comparative threshold (2 −ΔΔCt ). All experiment were performed in triplicate.

Classification and Differential Expression Analysis
The potential lncRNAs were classified based on their association with the annotated protein-coding genes. In this study, the association between the potential lncRNAs and the annotated protein-coding genes was compared with cuffmerge [28]. Then, the classification and characteristics of lncRNAs were analyzed. StringTie was selected to analyze the expression levels of lncRNAs by calculating fragments per kilobase million (FPKM). The corrected p-value < 0.05 and an absolute value of the |log2FoldChange| ≥ 2 were set as the threshold to assess statistically significant differences of lncRNA expression.

Co-Location and Co-Expression Analysis, and Functional Annotation Analysis
LncRNAs exert cis-regulatory effects on their co-localized genes [13]. However, lncR-NAs regulate the expression of genes located on other chromosomes through a trans-acting mechanism [29]. To assess the potential function of DE lncRNAs, we predicted the possible target genes that were co-located and co-expressed with lncRNAs. In this study, coding genes located 100 kb/100 kb of DE lncRNA were classified as putative cis-acting targets. The trans-acting correlations between lncRNAs and genes (co-expression) were calculated with the Pearson's correlation coefficients method (|r| > 0.95 and p < 0.05), as previously reported [22,30]. The cis-and trans-acting putative target genes were annotated by Gene Ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) to explore the potential functions 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
To validate the reliability of the RNA-seq data, 10 DE lncRNAs were randomly selected to verify their expression patterns in skeletal muscle of goats at different development stages. Primers for the 10 lncRNAs and endogenous reference gene (Table 1) were designed by Primer 5, and the goat GAPDH was selected as the endogenous reference. Briefly, total RNA was extracted from the longissimus dorsi muscle samples. For RT-qPCR analysis, 1 µg of the total RNA was reverse transcribed using PrimeScript™ RT reagent kit (Takara, Beijing, China) according to the manufacturer's protocol. Then, the qPCR was performed by a RocheLight Cycler ® 480 II system with SYBR Green qPCR Mix kit as previous study [31]. The expression levels were analyzed using comparative threshold (2 −∆∆Ct ). All experiment were performed in triplicate. To further predict the functional of lncRNA in skeletal muscle development, the ceRNA (lncRNA-miRNA-mRNA) networks were constructed based on the theory that lncRNA directly associates with miRNA and subsequently affects the activity of its mRNA [32]. Firstly, the correlation between lncRNA and miRNA or miRNA and mRNA was analyzed by correlation coefficient. All pairs with COR > 0.85 and adjusted p < 0.05 were selected as potential lncRNA-miRNA or miRNA-mRNA pairs. Then, the ceRNA networks were constructed with the DE lncRNAs, miRNAs, and mRNAs. Finally, the networks were visualized by Cytoscape software.

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

RNA-Seq Data Filtering, Mapping, and Transcript Assembly
To identify the function of lncRNAs in skeletal muscle growth, two cDNA libraries were constructed using longissimus dorsi samples of goats at two development stages: 1-month-old and 9-month-old. A total of 89.57 Gb of raw data were generated. After filtering out low-quality and adaptor sequences, an average of~14.58 Gb high-quality clean reads remained. The average GC percentage was 53.41%, with the quality scores of Q20 and Q30 above 96% and 91%, respectively ( Table 2). Approximately 72.47-85.25% of the high-quality reads were mapped to the goat reference genome (Table S1). The mapped sequences in the library were assembled and a total of 65,247 transcripts were obtained.

Identification and Confirmation of LncRNAs in Goat Longissimus Dorsi Tissue
We developed a highly stringent filtering pipeline to discard transcripts that did not display characteristics of lncRNAs. The assembled transcripts from the two libraries were filtered to obtain candidate lncRNAs. A total of 3441 lncRNAs were screened (Table S2), including 1281 (37.2%) long intergenic ncRNAs (lincRNAs), 805 (23.4%) antisense lncRNAs, and 1355 (39.4%) sense_overlapping lncRNAs ( Figure 2a). There was no sense_intronic lncR-NAs in this study. These 3441 putative lncRNA were encoded by 2675 genes. There were 1.3 transcripts on average per locus (Table S2). The lncRNA transcripts were widespread in chromosomes, including 29 autosomes and the X chromosome ( Figure S1), which reflected the function diversity of lncRNAs.

Identification and confirmation of lncRNAs in goat longissimus dorsi tissue
We developed a highly stringent filtering pipeline to discard transcripts that did not display characteristics of lncRNAs. The assembled transcripts from the two libraries were filtered to obtain candidate lncRNAs. A total of 3441 lncRNAs were screened (Table S2), including 1281 (37.2%) long intergenic ncRNAs (lincRNAs), 805 (23.4%) antisense lncRNAs, and 1355 (39.4%) sense_overlapping lncRNAs (Figure 2a). There was no sense_intronic lncRNAs in this study. These 3441 putative lncRNA were encoded by 2675 genes. There were 1.3 transcripts on average per locus (Table S2). The lncRNA transcripts were widespread in chromosomes, including 29 autosomes and the X chromosome (Figure S1), which reflected the function diversity of lncRNAs.
The Illumina RNA-seq also identified 66 novel mRNAs (Table S3) (Table S2). Importantly, the open reading frame (ORF) length of lncRNA was shorter compared with mRNA ( Figure 2d). The result showed that the coding potential of lncRNAs was lower.  The Illumina RNA-seq also identified 66 novel mRNAs (Table S3). As shown in Figure 2b,c, the transcript length and exon number of lncRNAs were both lower than that of the mRNA. The average length of lncRNA was 2001 bp with an average of 3.2 exons. The principal lncRNA transcripts with 2 exons accounted for 63.6% of the 3441 lncRNAs (Table S2). Importantly, the open reading frame (ORF) length of lncRNA was shorter compared with mRNA ( Figure 2d). The result showed that the coding potential of lncRNAs was lower.

Functional enrichment analysis
To evaluate the potential function of DE lncRNAs, we predicted the cis-and transpotential targets of DE lncRNAs. In the co-location analysis, 30 lncRNAs (1 annotated lncRNAs and 29 novel lncRNAs) were transcribed close to 71 protein-coding genes (Table  S5). GO enrichment analysis showed that 210 GO terms were significantly enriched (p < 0.05), including biological process (BP), cellular component (CC), and molecular function (MF) ( Table S6). Only the top 30 GO terms were shown in Figure 4a. According to GO enrichment analysis, 12 unique genes were enriched in muscle development-related terms, such as myoblast fate determination, skeletal system development, negative regulation of G1/S transition of mitotic cell cycle, and negative regulation of cell cycle G1/S phase transition. These genes, including IFRD1, CSRNP1, DYM, FLI1, MEPE, MUSTN1, TNFRSF11B, WDR48, APC, CRADD, MYO5B, and MYO16, may regulate skeletal muscle development (Table S6). The co-location interactions between lncRNAs and potential target genes related to muscle development were visualized by cis-regulatory network (Figure 5A). Moreover, 57 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were enriched through the pathway analysis (Table S7). The top 20 significantly enriched KEGG analyses were shown in Figure 4b. The potential target genes of DE lncRNAs associated with muscle development were involved in Glycerolipid metabolism, signaling pathways regulating pluripotency of stem cells, fatty acid metabolism, and biosynthesis of amino acids (Figure 4b).
Furthermore, we further predicted the potential targets of lncRNA in trans regulation based on Pearson′s correlation coefficients (|r| > 0.95). A total of 8438 interaction relationships were detected between 36 lncRNAs and 2684 mRNAs in goat reference genome (Table S8). The top 200 interaction relationships, which included 3 lncRNAs (novel lncRNAs) and 103 putative target genes, were selected for the subsequent functional cluster analysis. The putative trans-targets of DE lncRNAs were significantly enriched in 230 GO terms (p < 0.05) (Table S9), of which the top 30 are shown in Figure 4c. GO analysis indicated that NKX2-5, POU4F1, NPPA, IGFBP1, and TLE6 were enriched in muscle development-

Functional Enrichment Analysis
To evaluate the potential function of DE lncRNAs, we predicted the cis-and transpotential targets of DE lncRNAs. In the co-location analysis, 30 lncRNAs (1 annotated lncR-NAs and 29 novel lncRNAs) were transcribed close to 71 protein-coding genes (Table S5). GO enrichment analysis showed that 210 GO terms were significantly enriched (p < 0.05), including biological process (BP), cellular component (CC), and molecular function (MF) ( Table S6). Only the top 30 GO terms were shown in Figure 4a. According to GO enrichment analysis, 12 unique genes were enriched in muscle development-related terms, such as myoblast fate determination, skeletal system development, negative regulation of G1/S transition of mitotic cell cycle, and negative regulation of cell cycle G1/S phase transition. These genes, including IFRD1, CSRNP1, DYM, FLI1, MEPE, MUSTN1, TNFRSF11B, WDR48, APC, CRADD, MYO5B, and MYO16, may regulate skeletal muscle development (Table S6). The co-location interactions between lncRNAs and potential target genes related to muscle development were visualized by cis-regulatory network (Figure 5a). Moreover, 57 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were enriched through the pathway analysis (Table S7). The top 20 significantly enriched KEGG analyses were shown in Figure 4b. The potential target genes of DE lncRNAs associated with muscle development were involved in Glycerolipid metabolism, signaling pathways regulating pluripotency of stem cells, fatty acid metabolism, and biosynthesis of amino acids (Figure 4b).
scriptional regulation of POU4F1 and NKX2-5 or others. Finally, KEGG analysis indicated the putative trans-targets of DE lncRNAs were involved in inositol phosphate metabolism, phosphatidylinositol signaling system, Notch signaling pathway, and HIF-1 signaling pathway (Table S10, Figure 4d). Overall, lncRNAs and their putative target genes showed great potential in the regulation of skeletal muscle growth and development.  scriptional regulation of POU4F1 and NKX2-5 or others. Finally, KEGG analysis indicated the putative trans-targets of DE lncRNAs were involved in inositol phosphate metabolism, phosphatidylinositol signaling system, Notch signaling pathway, and HIF-1 signaling pathway (Table S10, Figure 4d). Overall, lncRNAs and their putative target genes showed great potential in the regulation of skeletal muscle growth and development.  Furthermore, we further predicted the potential targets of lncRNA in trans regulation based on Pearson's correlation coefficients (|r| > 0.95). A total of 8438 interaction relationships were detected between 36 lncRNAs and 2684 mRNAs in goat reference genome (Table S8). The top 200 interaction relationships, which included 3 lncRNAs (novel lncRNAs) and 103 putative target genes, were selected for the subsequent functional cluster analysis. The putative trans-targets of DE lncRNAs were significantly enriched in 230 GO terms (p < 0.05) (Table S9), of which the top 30 are shown in Figure 4c. GO analysis indicated that NKX2-5, POU4F1, NPPA, IGFBP1, and TLE6 were enriched in muscle development-related terms of canonical Wnt signaling pathway, muscle tissue morphogenesis, muscle organ morphogenesis, myotube differentiation, and so on. The co-expression interactions between lncRNAs and the putative target genes related to muscle development were visualized by trans-regulatory network (Figure 5b). The analysis suggested lncRNA TCONS_00085732 and TCONS_00005111 may affect muscle development through transcriptional regulation of POU4F1 and NKX2-5 or others. Finally, KEGG analysis indicated the putative trans-targets of DE lncRNAs were involved in inositol phosphate metabolism, phosphatidylinositol signaling system, Notch signaling pathway, and HIF-1 signaling pathway (Table S10, Figure 4d). Overall, lncRNAs and their putative target genes showed great potential in the regulation of skeletal muscle growth and development.

Verification of DE LncRNAs Expression Profile with RT-qPCR
Ten DE lncRNAs were randomly selected to explore their expression profile in skeletal muscle of goats at two development stages (1-month-old, 9-month-old) using RT-qPCR. Expression analysis revealed these 10 lncRNAs exhibited differential expression profile during skeletal growth in goats ( Figure 6). The expression of five lncRNAs, including TCONS_00134767, TCONS_00126170, TCONS_00124841, TCONS_00121766, and TCONS_00026838, increased during the process of skeletal muscle growth. In contrast, the expression of five other lncRNAs, including TCONS_00036756, TCONS_00066056, TCONS_00062751, XR_001296113.2, and TCONS_00150002, decreased. Importantly, the expression trends of the 10 significant differentially expressed lncRNAs were basically consistent with the RNA-seq results. The above analysis revealed that the pipeline we developed was reasonable to identify putative lncRNAs, and that the RNA-seq results were reliable.

Verification of DE lncRNAs expression profile with RT-qPCR
Ten DE lncRNAs were randomly selected to explore their expression profile in skeletal muscle of goats at two development stages (1-month-old, 9-month-old) using RT-qPCR. Expression analysis revealed these 10 lncRNAs exhibited differential expression profile during skeletal growth in goats ( Figure 6). The expression of five lncRNAs, including TCONS_00134767, TCONS_00126170, TCONS_00124841, TCONS_00121766, and TCONS_00026838, increased during the process of skeletal muscle growth. In contrast, the expression of five other lncRNAs, including TCONS_00036756, TCONS_00066056, TCONS_00062751, XR_001296113.2, and TCONS_00150002, decreased. Importantly, the expression trends of the 10 significant differentially expressed lncRNAs were basically consistent with the RNA-seq results. The above analysis revealed that the pipeline we developed was reasonable to identify putative lncRNAs, and that the RNA-seq results were reliable.

Construction and bioinformatics analysis of the ceRNA network
To elucidate the potential interaction of above lncRNAs, the lncRNA-miRNA-mRNA ceRNA networks were constructed by Cytoscape software (Figure S2). Considering the functional diversity of lncRNAs, we selected lncRNAs with muscle functions, which were based on the reported roles of mRNA in previous studies to construct the following

Construction and Bioinformatics Analysis of the ceRNA Network
To elucidate the potential interaction of above lncRNAs, the lncRNA-miRNA-mRNA ceRNA networks were constructed by Cytoscape software (Figure S2). Considering the functional diversity of lncRNAs, we selected lncRNAs with muscle functions, which were based on the reported roles of mRNA in previous studies to construct the following ceRNA networks diagram (Figure 7). The ceRNA networks comprised 3 miRNAs, 8 mRNAs, and 4 lncRNAs transcripts, which interacted with at least one miRNA.

Discussion
The growth and development of skeletal muscle is one of the most important factors affecting the meat production of livestock [36]. It is well known that skeletal muscle development involves a series of exquisitely regulated and orchestrated changes in the expression of many genes [37]. Importantly, the muscle development in goats involves two main stages-the embryonic and the postnatal stages. During the embryonic stage, muscle development is completed, and the number of muscle fibers generally does not change after birth. At postnatal stage, the muscle growth is mainly triggered by muscle fiber hypertrophy (the diameter and length of the myoblast) and an increase in intermuscular fat [38]. In general, kid goat (from newborn to 90-day-old) was related to muscle fibers' hypertrophy and regulation of myoblast proliferation [39]. In addition, goats grow very quickly from 0 to 7 months, but slowly after 18 months [40]. In this study, two growth stages (1-month-old and 9-month-old) were selected for researching the molecular mechanism of muscle growth and development. The lncRNAs expression profile of goat longissimus dorsi muscle at the two stages were detected by RNA-seq and bioinformatics analysis. We investigated the transcript structure and expression patterns of lncRNAs in goat skeletal muscle tissue, and then explored the potential functions of cis-and trans-potential target genes of lncRNAs. It also showed that miRNA forms 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_001296113.2 could regulate PDLIM7 by competing with chi-miR-1296. Meanwhile, lncRNA XR_001917947.1, XR_001917946.1, and XR_001917948.1 functioned as ceRNAs by regulating chi-miR-30b-3p, which affects smad3 expression.

Discussion
The growth and development of skeletal muscle is one of the most important factors affecting the meat production of livestock [36]. It is well known that skeletal muscle development involves a series of exquisitely regulated and orchestrated changes in the expression of many genes [37]. Importantly, the muscle development in goats involves two main stages-the embryonic and the postnatal stages. During the embryonic stage, muscle development is completed, and the number of muscle fibers generally does not change after birth. At postnatal stage, the muscle growth is mainly triggered by muscle fiber hypertrophy (the diameter and length of the myoblast) and an increase in intermuscular fat [38]. In general, kid goat (from newborn to 90-day-old) was related to muscle fibers' hypertrophy and regulation of myoblast proliferation [39]. In addition, goats grow very quickly from 0 to 7 months, but slowly after 18 months [40]. In this study, two growth stages (1-month-old and 9-month-old) were selected for researching the molecular mechanism of muscle growth and development. The lncRNAs expression profile of goat longissimus dorsi muscle at the two stages were detected by RNA-seq and bioinformatics analysis. We investigated the transcript structure and expression patterns of lncRNAs in goat skeletal muscle tissue, and then explored the potential functions of cis-and trans-potential target genes of lncRNAs.
Few studies have reported on the crucial function of lncRNAs expression profile in Wu'an black goat, especially in skeletal muscle growth and development. In this study, the cDNA libraries were generated with Illumina NovaSeq. An average of~14.58Gb highquality clean reads were obtained after the filtering process. Subsequently, the high-quality clean reads were mapped to the goat reference genome and assembled with StringTie, and a total of 65,247 library transcripts were obtained. LncRNAs can be single-or multiexon, making it difficult to distinguish putative lncRNAs from the plentiful sing-exon, lowly expressed and unreliable sequenced fragments [41,42]. To minimize the selection of false positive lncRNAs, we set up a relatively stringent filtering pipeline to obtain true lncRNAs with high confidence. Only multi-exon lncRNAs were selected for further exploration, as was done in other studies [23,43]. As a result, 3,441 putative lncRNAs with high confidence were identified. In addition, most lncRNAs were longer than 2000 bp in length and contained 2 exons, which was in agreement with previous studies in goats and sheep [21,23]. Moreover, our current data demonstrated that lncRNAs had fewer exons and shorter transcript length and ORFs length than mRNA, which was consistent with other studies [22,23]. These similarities suggested that the putative lncRNAs verified in our study were reliable. This was essential to expand our understanding of lncRNAs via association with multiple structural features.
In recent years, numerous studies have reported on the biological functions of lncR-NAs. For example, lnc-SEMT regulates IGF2 expression via competing with miR-125b to facilitate skeletal muscle growth and development in sheep [17]. Besides, lncR125b promotes the differentiation of goat skeletal muscle satellite cells by competing with miR-125b [16]. In addition, a series of lncRNAs, such as lncIRS1, lnc-RAM, and lncMD1, have been identified to have an essential effect on skeletal myogenesis [44][45][46]. LncRNAs are non-coding transcripts that act as regulators of gene expression and are involved in skeletal muscle development [47]. In the present study, we detected 36 DE lncRNAs in skeletal muscle of goats at two growth stages. These lncRNAs showed significantly different changes between 1-month-old and 9-month-old goats. They may have certain biological functions during skeletal muscle growth. Therefore, the DE lncRNAs verified in this study could be regarded as vital putative regulators of muscle biology. Furthermore, the expression trends of 10 randomly DE lncRNA (either up or down regulated) identified by RT-qPCR were consistent with the RNA-seq results. Together, these data provided a valuable resource for exploring the role of lncRNAs in postnatal muscle growth and provided new insights into the dynamic gene regulation of muscle biology in Wu'an black goat.
Unlike protein-coding genes, the biofunctions of lncRNAs cannot be directly speculated from their sequence or structure. Therefore, we attempted to uncover the function of lncRNAs based on their potential cis-and tans-acting target genes in our study [23,39,48]. Additionally, GO and KEGG analyses were carried out to further explore the function of lncRNAs. It has been known that some lncRNAs are thought to work in cis on neighboring genes, and other lncRNAs work in trans to regulate distantly located genes. The cis-regulating target genes play an important role in estimating the biological functions of lncRNAs. We searched for the potential cis-target genes located within 100 kb upstream and downstream of the identified lncRNAs, which was comparable with previous studies [19,30]. Of these potential cis-target genes, some were all enriched in the GO term of skeletal muscle system development, including APC, IFRD1, CSRNP1, TNFRSF11B, WDR48, and so on. The result implied that the corresponding lncRNAs had a vital role in regulating skeletal muscle development. Emerging studies have documented that lncRNAs were associated with cis-regulation 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 [40]. Co-location correlation analysis combining GO analysis indicated that TCONS_00062751 could affect muscle development by targeting APC. The lncRNA TCONS_00062751 was selected because of its cis-target APC enriching in the GO terms of negative regulation of G1/S transition of mitotic cell cycle and negative regulation of cell cycle G1/S phase transition. Moreover, it has been reported that APC is required for muscle stem cell proliferation and skeletal muscle tissue repair [49]. Furthermore, the potential cis-target genes of TCONS_00026388 and TCONS_00026389 mainly enriched myoblast fate determination, in which IFRD1 was significantly enriched. IFRD1 played a vital role in myoblast differentiation by regulating the expression of MyoD and NF-kappaB [50]. These results suggested that lncRNAs acted in cis on neighboring protein-coding genes to regulate muscle development in goats.
Still, some lncRNAs have been found to function in trans-acting to target gene loci far from the transcriptional location of lncRNAs [51]. For instance, LncRNA MUNC, encode 5 kb upstream of the MyoD transcription start site, facilitates the biofunction of MyoD in muscle biology [52]. In this study, the co-expression analysis between DE lncRNA and mRNA on different chromosomes were conducted based on the Pearson correlation coefficient. We found that several mRNAs were regulated by the trans-action of lncRNAs clustered in GO terms and KEGG pathways associated with muscle development. For example, NKX2-5 was enriched in terms of negative regulation of canonical Wnt signaling pathway, muscle tissue morphogenesis, and muscle organ morphogenesis and myotube differentiation. NKX2-5 has been reported to regulate the differentiation of skeletal myoblasts in vitro [39]. Based on the trans-regulatory interactions, TCONS_00085732 and TCONS_00005111 may affect muscle development by targeting NKX2-5. It was worth noting that POU4F1 was also the potential trans-action target of TCONS_00085732 and TCONS_00005111. However, we found no reports on the regulation of muscle development by POU4F1. The lncRNA TCONS_00085732 and TCONS_00005111 could regulate muscle development on account of its putative trans-target POU4F1, mainly enriching in the GO terms of muscle tissue morphogenesis and muscle organ morphogenesis.
Beyond that, another aspect worth considering is the low number of genes in "functionally enriched" categories (some terms have only one gene). Taking the co-location analysis as an example, the potential cis-targets of only the differentially expressed lncR-NAs (not all screened 3441 LncRNAs) were predicted in the present study. As shown in the manuscript, 71 putative target genes were enriched in 210 GO terms. Therefore, the GO terms with assigned 1 gene were reasonable. Similar results have been found in several other studies [30,53,54]. The above information validated that lncRNAs may be involved in goat skeletal muscle biology through cis-or trans-regulation. Overall, these genes associated with muscle development were derived from target predictions. Further validation is necessary to explore the function of lncRNAs in skeletal muscle development in goats.
Currently, the function roles and mechanisms of ceRNA are widely reported [55,56], in which lncRNAs act as molecular sponges that regulate the expression of target mRNAs upon binding miRNAs. In this study, the lncRNA-miRNA-mRNA ceRNA networks associated with muscle development in goats were constructed by sequencing and bioinformatic. The functional lncRNAs serve as key regulators in muscle biology [44][45][46]. Therefore, it was important to investigate the function of lncRNA in skeletal muscle growth. The ceRNA networks analysis showed that lncRNA XR_001296113.2 acted as a sponge of chi-miR-1296 to regulate the expression of PDLIM7. Similarly, XR_001917947.1 functioned as a ceRNA to regulate smad3 expression by sponging chi-miR-30b-3p. The potential function of lncRNAs can be inferred from co-located or co-expressed protein-coding genes. Obviously, PDLIM7 and smad3 attracted our attention in this study. Previous studies have reported that PDLIM7 is a member of the PDZ-LIM proteins family, which is known to regulate skeletal muscle development [57,58]. Smad3 plays an important role in skeletal muscle regeneration, as lack of smad3 signaling leads to impaired skeletal muscle regeneration [59]. Taken together, the constructed ceRNA networks may function in goat muscle. Nevertheless, the biological functions of the lncRNA-miRNA-mRNA interactions described in this study need to be validated by further studies. The discussions mentioned above all suggested that lncRNAs were an important component of the regulatory network in skeletal muscle growth and development. The bioinformatics analysis showed that these lncRNAs exhibited significant changes during skeletal muscle growth in goats, implying that they may have certain functions in myofiber growth.

Conclusions
This study used RNA-seq to systematically characterize the expression profile of lncR-NAs in skeletal muscle of goats at two development stages. The DE lncRNAs associated with muscle growth in goats were validated. The verified lncRNAs in our study shared many common features in their structure. The functional annotation of DE lncRNAs was performed using GO and KEGG pathway enrichment analysis. In addition, the visualization of lncRNA-associated ceRNA networks provide a more comprehensive understanding of candidate lncRNAs that may regulate skeletal muscle growth in goats. However, the lncRNA-associated ceRNA network involved in skeletal muscle development was inferred from the protein-coding genes, and its role needs to be validated by further studies.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/ani12192683/s1, Figure S1: Chromosome distribution of lncRNAs identified in goat skeletal muscle.; Figure S2: The lncRNA-miRNA-mRNA ceRNA networks diagram.; Table S1: The information of RNA-seq data.; Table S2: List of lncRNAs identified in goat skeletal muscle libraries.; Table S3: List of mRNAs identified in goat skeletal muscle libraries.; Table S4. List of DE lncRNAs from two muscle developmental stages.; Table S5: Protein-coding genes detected 100 kb upstream and downstream of DE lncRNAs.; Table S6: GO enrichment analysis of proteincoding genes co-located with lncRNAs.; Table S7: KEGG analysis of protein-coding genes co-located with lncRNAs.; Table S8: Co-expression analysis between protein-coding genes and lncRNAs.; Table S9: GO enrichment analysis of protein-coding genes co-expressed with lncRNAs.; Table S10: KEGG analysis of protein-coding genes co-expressed with lncRNAs.

Informed Consent Statement: Not applicable.
Data Availability Statement: The datasets supporting the conclusions of this article are included within the article and its supplementary files. The raw data obtained from RNA-seq can be found in the National Center for Biotechnology Information (NCBI) Sequence Read Archive repository, https://www.ncbi.nlm.nih.gov/sra/PRJNA749569 (accessed on 24 August 2022).