Integrated transcriptome analysis reveals roles of long non-coding RNAs (lncRNAs) in caprine skeletal muscle mass and meat quality

Long non-coding RNAs (lncRNAs) play important roles in the growth and development of skeletal muscle. However, there is limited information on goats. In this study, expression profiles of lncRNAs in Longissimus dorsi muscle from Liaoning cashmere (LC) goats and Ziwuling black (ZB) goats with divergent meat yield and meat quality were compared using RNA-sequencing. Based on our previous microRNA (miRNA) and mRNA profiles obtained from the same tissues, the target genes and binding miRNAs of differentially expressed lncRNAs were obtained. Subsequently, lncRNA-mRNA interaction networks and a ceRNA network of lncRNA-miRNA-mRNA were constructed. A total of 136 differentially expressed lncRNAs were identified between the two breeds. Fifteen cis target genes and 143 trans target genes were found for differentially expressed lncRNAs, and they were enriched in muscle contraction, muscle system process, muscle cell differentiation, and p53 signaling pathway. A total of 69 lncRNA-trans target gene pairs were constructed, with close relationship with muscle development, intramuscular fat deposition, and meat tenderness. A total of 16 lncRNA-miRNA-mRNA ceRNA pairs were identified, of which some reportedly associated with skeletal muscle development and fat deposition were found. The study will provide an improved understanding of the roles of lncRNAs in caprine meat yield and meat quality.


Introduction
Domestic goats (Capra hircus) are one of the most essential farm animals economically. Goat meat is increasingly favored by consumers for higher protein, lower fat, and unique flavor and palatability when compared to meat produced from other domestic animals (Teixeira et al. 2019). It is well known that as in other animals, caprine meat yield and quality are regulated by both mRNAs and non-coding RNAs including long non-coding RNAs (lncRNAs), so the identification of the RNAs that regulate skeletal muscle offer an opportunity to improve caprine meat production performance.
The lncRNAs are a class of novel non-coding RNAs with greater than 200 nucleotides in length. They are primarily transcribed by RNA polymerase II and widely distributed in mammalian cells (Agliano et al. 2019). Although the vast majority of lncRNAs expressed at lower levels compared to protein-coding genes , they were found to play crucial roles in the regulation of mRNA expression at transcriptional, post-transcriptional, or epigenetic levels. For example, some lncRNAs transcribed from genome sequences containing cis-regulatory DNA elements can regulate the expression of their neighboring genes in cis by affecting their transcription (Kopp and Mendell 2018). The lncRNAs can also trans-regulate the expression of the target genes by interacting with trans-acting factors (Kopp and Mendell 2018). In addition, some lncRNAs act as competing endogenous RNAs (ceRNAs) to sequester microR-NAs (miRNAs), eventually leading to increased expression of the target gene by miRNAs at post-transcriptional level (Zhan et al. 2019). In this context, lncRNAs are considered to widely involve in embryonic development, organ morphogenesis, and cell cycles, including differentiation, proliferation, and apoptosis (Cheng et al. 2020;Zhu et al. 2021).
It has been confirmed that lncRNAs played important roles in the growth and development of skeletal muscle. For 63 Page 2 of 13 example, lncR-125b accelerated differentiation of caprine muscle satellite cells by increasing expression of insulin growth factor 2 (IGF2) derived from a ceRNA for miR-125b (Zhan et al. 2019). LncRNA-Six1 promoted chicken muscle growth by cis-regulating expression of the target gene SIX homeobox 1 (Six1) . Additionally, lncMAAT was found to regulate muscle atrophy by negatively regulating the expression of miR-29b . Meanwhile, the expression profiles of lncRNAs in skeletal muscle tissues were also investigated and lncRNAs were found to be differentially expressed in samples with different genetic backgrounds. However, these studies in livestock have mainly been focused on sheep (Ren et al. 2017;Li et al. 2019;Li et al. 2018), cattle (Li et al. 2020a;Yan et al. 2021;Huang et al. 2021), chicken (Li et al. 2012;Li et al. 2017;Liu et al. 2021a), and pigs Li et al. 2020b;Hou et al. 2021).
Despite there being some studies reporting lncRNA expression profiles of skeletal muscle tissues in goats, these studies were all performed with muscle samples collected from different developmental periods. For example, a total of 547 lncRNAs were differentially expressed in skeletal muscles of Anhui white goats among five fetal stages and two kid stages, and they were involved in signaling pathways closely associated with muscle development, including structure formation, p53 signaling pathway, and MAPK signaling pathway (Ling et al. 2019). Additionally, 577 and 648 differentially expressed lncRNAs were found in skeletal muscle tissues between embryonic and postnatal stages of Jianzhou big-eared goats and Dazu black goats, respectively (Zhan et al. 2016;Huang et al. 2022). However, there have been no reports on comparison of lncRNA expression profiles in skeletal muscle tissues between different goat breeds.
In this study, Liaoning cashmere (LC) goats and Ziwuling black (ZB) goats were selected for the investigation. The two breeds are all dual-purpose breeds for meat and cashmere fiber, and have significant differences in meat production performance. Briefly, LC goats have higher carcass weight, net meat weight, muscle fiber size, and intramuscular fat content, but have poorer meat tenderness compared to ZB goats (P < 0.05) . RNA-Seq analysis of the same Longissimus dorsi muscle tissues as those used in the study revealed that differentially expressed genes, miRNAs, and circular RNAs (circRNAs) were responsible for these phenotypic differences (Shen et al. , 2022a. However, the biological mechanism by which lncRNAs regulate meat production performance difference is still unclear. Accordingly, the expression profiles of lncRNAs in Longissimus dorsi muscle tissues between LC and ZB goats were compared using RNA sequencing (RNA-Seq), and differentially expressed lncRNAs were then screened. We also analyzed functional enrichment of the target genes, and constructed lncRNA-mRNA interaction networks and a ceRNA network of lncRNA-miRNA-mRNA, with the aim of uncovering the possible function of lncRNAs in muscle growth and development of goats.

Ethics statement
All animal procedures in this study were approved by the Animal Experiment Ethics Committee of Gansu Agricultural University with an approval number of GSAU-ETH-AST-2021-028.

Experimental animals and sample collection
Five healthy 9-month-old LC goats and five healthy 9-month-old ZB goats raised under the same environmental conditions and nutrition program were selected from the Yongfeng Goat Breeding Company (Huan County, Gansu Province, China). Their meat yield and quality were measured after slaughter and are presented in Supplementary File 1. The Longissimus dorsi muscle was collected from the region between 12th and 13th ribs in the left half carcass from each goat. These tissue samples were immediately stored in liquid nitrogen for RNA extraction.

RNA-sequencing
Total RNA from the ten muscle samples was isolated using a Trizol reagent kit (Invitrogen, Carlsbad, CA, USA). Its concentration and RNA Integrity Number (RIN) were measured using a Nanodrop 2000 (Thermo Scientific, MA, USA) and Agilent 2100 Bioanalyzer (Agilent, CA, USA), respectively. Only samples with an RIN value > 7 were used for subsequent cDNA library construction. The ribosomal RNA (rRNA) was removed from total RNA using a Ribo-Zero Gold rRNA Removal Kit (Illumina, CA, United States), and the remaining Longissimus dorsi muscle RNA was used to generate complementary DNA (cDNA) libraries using a NEBNext Ultra RNA Library Prep Kit (New England Biolabs, MA, USA). The cDNA libraries obtained were pairedend sequenced using an Illumina HiSeq TM 4000 sequencer (Illumina, CA, United States) at Gene Denovo Biotechnology Co., Ltd (Guangzhou, China).

Identification and characterization of lncRNAs
The clean reads were obtained by removing low-quality reads with quality scores < Q20, reads with > 10% unknown nucleotides, and adaptor reads from raw reads using Fastp v0.18.0 (Chen et al. 2018). The clean reads were aligned against Caprine Genome Assembly ARS1 using HISAT2 v2.1.0 (Kim et al. 2015), and the known lncRNAs were then identified according to the alignment results. Subsequently, the mapped results from HISAT2 v2.1.0 were assembled into transcripts using Stringtie v1.3.4 (Pertea et al. 2015) and the novel lncRNAs were then identified. Briefly, the transcripts with smaller than 200 nucleotides in length or those with a single exon were discarded. The coding potentials of the remaining transcripts were predicted by CNCI v2.0 (Sun et al. 2013) and CPC v0.9-r2 (Kong et al. 2007), and the predicted results from the two kinds of software were intersected. These transcripts with coding potential score < 0 were considered to be novel lncRNAs. The identified caprine lncRNAs were characterized by analyzing types, transcript length, exon number, and coding potential score to mRNAs identified from the same muscle samples as those used in the study , using in-house Perl scripts.

Screening of differentially expressed lncRNAs
The expression level of each annotated lncRNA was normalized by calculating Fragments Per Kilobase of transcript per Million mapped reads (FPKM) using StringTie v1.3.4. The computational formula of FPKM is: FPKM = 10 6 C NL∕10 3 . Given FPKM to be the expression of a specific lncRNA transcript, C to be the number of fragments mapped to the specific lncRNA transcript, N to be total number of fragments that mapped to reference genome, and L to be number of bases on the specific lncRNA transcript. The lncRNAs with FPKM > 0 were considered to be expressed in skeletal muscle. The differential expression analysis of lncRNAs between the two goat breeds was performed based on the negative binomial distribution model in DESeq v2.0 (Love et al. 2014). The lncRNAs with Fold change > 2 and P-value < 0.05 were defined as differentially expressed lncRNAs.

Prediction and functional annotation of the target genes of differentially expressed lncRNAs
The lncRNAs have been reported to affect the expression of their neighboring genes in cis, while they also influence coexpressed genes in trans. In our previous study, differentially expressed genes were screened from the same caprine Longissimus dorsi muscle tissues as those used in the study, using the thresholds of |fold change| > 2.0 and false discovery rate (FDR) value < 0.05 . In order to obtain more target genes of lncRNAs, we adjusted the thresholds to |fold change| > 2.0 and P value < 0.05 to re-identify differentially expressed genes. The differentially expressed genes screened above were used to search cis and trans target genes. For screening of cis target genes, the coding genes located within 100 kb upstream and downstream of differentially expressed lncRNAs were searched. The search results were compared with differentially expressed genes obtained. The overlapped genes were chosen as cis target genes of differentially expressed lncRNAs. Subsequently, co-expression analysis was performed in expression levels between differentially expressed lncRNAs identified in this study and differentially expressed genes obtained. Briefly, Pearson's coefficients in expression levels between the differentially expressed lncRNAs and the differentially expressed genes were calculated; only the differentially expressed genes with |r| > 0.95 and P-value < 0.05 were chosen as trans target genes of differentially expressed lncRNAs.
To further investigate the function of differentially expressed lncRNAs, their cis and trans target genes were used to perform Gene Ontology (GO) and Kyoto Encyclopedia of genes and genomes (KEGG) functional enrichment analysis using the Gene Ontology database (http:// www. geneo ntolo gy. org/) (Ashburner et al. 2000) and KEGG database (http:// www. kegg. jp/ kegg/) (Kanehisa et al. 2008), respectively. The significant GO terms and KEGG pathways (P < 0.05) were screened based on the hypergeometric test.

Construction of lncRNA-mRNA interaction networks
To better understand how these differentially expressed lncRNAs regulate the phenotype differences in meat yield and meat quality between LC and ZB goats, their target genes related to skeletal muscle development, intramuscular fat deposition, and meat tenderness were selected based on GO and KEGG enrichment analysis results as well as their functions reported in literatures. The lncRNA-mRNA interaction networks were finally constructed using Cytoscape v3.5.1.

A ceRNA network of lncRNA-miRNA-mRNA
Based on our previous Illumina HiSeq miRNA data (SRR16760528-SRR16760537) (Shen et al. 2022b) and mRNA data (SRR13008213-SRR13008222) ) obtained from the same caprine Longissimus dorsi muscle tissues as those used in the study using small RNA sequencing and RNA-Seq, respectively, a lncRNA-associated ceRNA network was constructed according to ceRNA hypothesis as follows (Quan et al. 2021): (1) The target genes and binding lncRNAs of differentially expressed miRNAs between LC goats and ZB goats found from small RNA sequencing were firstly predicted using miReap v0.2, Miranda v3.3a, and TargetScan v7.0, and the predicted results were then overlapped with differentially expressed genes and differentially expressed lncRNAs obtained from RNA-Seq. These overlapped lncRNAs and target genes were preliminary used to construct the ceRNA network.
(2) Correlation in expression between the lncRNA and the miRNA or between the miRNA and the mRNA was evaluated using the Spearman's rank correlation coefficient (SCC). Pairs with SCC < −0.7 were selected as negatively co-expressed lncRNA-miRNA or miRNA-mRNA pairs. (3) For lncRNA-miRNA and miRNA-mRNA pairs screened above, correlation in expression between the lncRNA and the mRNA was evaluated using the Pearson correlation coefficient (PCC); only those pairs with PCC > 0.9 were selected as coexpressed lncRNA-mRNA pairs. (4) Hypergeometric cumulative distribution function test was used to test whether a specific common miRNA sponge between a lncRNA and an mRNA was significant. Only pairs with P < 0.05 were selected as candidate ceRNA pairs. The lncRNA-miRNA-mRNA network was therefore visualized using Cytoscape v3.5.1.

Identification and characterization of lncRNAs in caprine skeletal muscle
For the ten Longissimus dorsi muscle tissue samples (five LC goats and five ZB goats), their clean reads and mapped results to the Caprine Genome Assembly ARS1 have been described in our previous study . The correlations between the five samples within Liaoning cashmere goats and Ziwuling black goats are above 0.92 and 0.93, respectively. In this study, a total of 2302 lncRNAs were identified, including 1767 known caprine lncRNAs and 535 novel lncRNAs (Supplementary File 2). Of all the lncRNAs identified, 1945 lncRNAs were co-expressed in the two goat breeds, while 178 and 179 lncRNAs were specifically expressed in LC and ZB goats, respectively (Fig. 1A). Of the five lncRNA types classified according to their location relative to proteincoding genes, intergenic lncRNA was the most common with a proportion of 54.4%, followed by antisense lncR-NAs (15.3%) and bidirectional lncRNAs (15.3%). Sense lncRNAs (3.4%) and intronic lncRNAs (2.3%) were the least. In addition, 9.3% lncRNAs were defined as another type of lncRNAs that were not classified into the five types (Fig. 1B). In our previous study, the mRNA expression profiles of the same Longissimus dorsi muscle tissues as those used in the study have been reported . In this study, characteristics in Longissimus dorsi muscle between lncR-NAs and mRNAs identified were accordingly compared. The average transcript length of known lncRNAs identified in this study was 1668 bp, which was less than novel lncRNAs and mRNAs with an average length of 5065 and 3674 bp, respectively (Fig. 1C). The average expression level of lncR-NAs was far less than mRNA transcripts (Fig. 1D). The average exon number of known and novel lncRNAs was 3.8 and 2.6, respectively, which was much less than average 12.7 exons for mRNAs (Fig. 1E). Additionally, the average length of open reading frame (ORF) of known and novel lncRNAs was shorter than mRNAs (Fig. 1F). This was consistent with the results that lncRNAs had much lower average coding potential score than mRNAs (Fig. 1G).

Screening and validation of differentially expressed lncRNAs
A total of 136 differentially expressed lncRNAs were identified in comparing Longissimus dorsi muscle tissue from LC goats with Longissimus dorsi muscle tissue from ZB goats. These included 62 upregulated lncRNAs and 74 downregulated lncRNAs in the Longissimus dorsi muscle of LC goats (Supplementary File 3).
The results from RT-qPCR showed that the relative expression levels of the 12 differentially expressed lncRNAs selected were in accordance with those obtained from RNA-Seq (Fig. 2). This confirmed the reliability and repeatability of the RNA-Seq results.

Prediction and functional annotation of the target genes of differentially expressed lncRNAs
To predict the target genes of differentially expressed lncR-NAs, a total of 391 upregulated genes and 222 downregulated genes in the same Longissimus dorsi muscle as those used in the study were re-identified in LC goats when compared to ZB goats (Supplementary File 4). In this study, 15 cis target genes were predicted for 18 differentially expressed lncRNAs, resulting in 22 lncRNA-cis target gene pairs (Supplementary File 5). Additionally, 143 trans target genes were predicted for 58 differentially expressed lncRNAs, resulting in 203 lncRNA-trans target gene pairs (Supplementary File 5). However, there were no target genes predicted for the remaining 60 differentially expressed lncRNAs.
The target genes in cis and trans were significantly enriched in 155 biological process (GO-BP) terms, 19 cellular components (GO-CC) terms, and 19 molecular function (GO-MF) terms (Supplementary File 6). It was noteworthy that some significant GO terms directly related to the growth and development of skeletal muscle were found (Fig. 3A), including muscle contraction (P = 0.002), regulation of muscle contraction (P = 0.004), muscle system process (P = 0.014), regulation of muscle system process (P = 0.022), muscle cell differentiation (P = 0.040), and striated muscle cell differentiation (P = 0.041). In addition, connective tissue development process related to meat tenderness was also significantly enriched by several target genes (P = 0.010) (Fig. 3A).
KEGG pathway analysis showed that the target genes were significantly enriched in 19 signaling pathways (Supplementary File 6) and the top fifteen significant pathways with the lowest P value are shown in Fig. 3B. Systemic lupus erythematosus was the most enriched pathway (P = 2.685E−05), followed by alcoholism (P = 2.623E−04) and allograft rejection (P = 0.007). The p53 signaling pathway (P = 0.023) was noteworthy as it played dual roles in the growth and development of muscle and adipose tissues (Fig. 3B).

Construction of lncRNA-mRNA interaction networks
Given unknown roles for 15 cis target genes in meat yield and quality, their interaction network with lncRNAs was not constructed. Of 203 lncRNA-trans target gene pairs described above, based on the roles of the target genes in meat yield and quality, a total of 69 lncRNA-trans target gene pairs were further selected to construct lncRNA-mRNA interaction Fig. 2 RT-qPCR validation of 12 differentially expressed lncRNAs identified using RNA-Seq. These included seven upregulated lncRNAs and five downregulated lncRNAs in Liaoning cashmere (LC) goats compared to Ziwuling black (ZB) goats networks, including 38 pairs related to muscle development (Fig. 4A), 18 pairs associated with intramuscular fat deposition (Fig. 4B), and 13 pairs related to meat tenderness (Fig. 4C). For example, the target genes CDK1, GREM1, MEGF10, and DNER were involved in muscle cell differentiation process (Supplementary File 6). The target genes RORC, KLK7, and RYR3 were related to adipogenesis (Tsai et al. 2013;Barendse et al. 2010;Kunath et al. 2021), while the target genes GREM1 and LOXL2 were involved in connective tissue development process affecting meat tenderness (Gur et al. 2003).

Construction of a lncRNA-miRNA-mRNA ceRNA network
In addition to the roles in regulating expression of the target genes, lncRNA can also acts as the sponge of miRNA to relieve the inhibitory effect of miRNA on its target mRNA, and finally increase the expression level of the target mRNA. In this study, a total of 16 lncRNA-miRNA-mRNA ceRNA pairs were identified (P < 0.05) (Fig. 5, Supplementary File 7). It was noteworthy that several ceRNA networks related to the growth and development of skeletal muscle were found including XR_001917125.1-miR-34-3p-NCAM2, XR_001917125.1-miR-460-5p-NCAM2, and XR_001918832.1-miR-200c-CRHBP. A ceRNA network related to fat deposition was also found including XR_001918832.1-miR-200b-CRHBP.

Discussion
Meat production performance is a kind of complex economic trait that include multiple important traits influencing meat yield and meat quality. For example, carcass weight and net meat weight directly reflect meat yield. Meat tenderness is regarded as an important palatability trait that significantly Fig. 3 GO and KEGG analysis of the target genes of differentially expressed lncRNAs identified between Liaoning cashmere (LC) goats and Ziwuling black (ZB) goats. A The important GO terms related to the growth and development of skeletal muscle and meat tenderness. B The top fifteen significant pathways enriched by the target genes of differentially expressed lncRNAs. The left side Y-axis represents the number of the target genes of differentially expressed lncRNAs, while the Y-axis on the right side shows the value of -Log10 (P-value) affects consumer acceptance of meat, while intramuscular fat influences sensory quality of meat, including flavor, juiciness and tenderness (Hocquette et al. 2010). An increasing number of evidences suggest that meat yield and quality traits can be directly regulated by some non-coding RNAs.
Given that lncRNAs account for 80% of the total number of non-coding RNAs, the role of lncRNAs in skeletal muscle is worthy of further investigation.
In this study, a total of 2302 lncRNAs were identified in Longissimus dorsi muscle of goats, and this was less than Fig. 4 LncRNA-mRNA interaction networks related to muscle development (A), intramuscular fat deposition (B), and meat tenderness (C). The red triangles and inverted triangles represent upregulated and downregulated lncRNAs in Longissimus dorsi muscle of Liaoning cashmere (LC) goats compared to Ziwuling black (ZB) goats, respectively. The green circles represent the target genes of the differentially expressed lncRNAs those reported in muscle tissue of Jianzhou big-eared goats with an average of 2739 lncRNAs identified (Zhan et al. 2016). This may partly reflect breed-specific expression patterns of lncRNAs. Additionally, the number of lncRNAs identified in this study was more than those for a study in Longissimus dorsi muscle of sheep ), but less than those reported in cattle by Yan et al. (2021). The differences may be related to species-specific expression of lncRNAs. The characteristics of lncRNAs identified in this study are in concordance with previous studies in skeletal muscle tissue. For example, intergenic lncRNA was found to be the most common in skeletal muscle of Dazu black goats (Huang et al. 2022), rabbit (Kuang et al. 2018) and cattle (Li et al. 2020a). Our observation that known lncR-NAs had shorter transcript length than mRNAs is consistent with findings in Anhui white goats (Ling et al. 2019), cattle (Li et al. 2020a) and donkey (Shi et al. 2020). As expected, lncRNAs found in this study had much lower expression levels and coding potential score, shorter ORF length, and fewer exon number when compared to mRNA transcripts. The phenomenon has also been reported in skeletal muscle of Anhui white goats (Ling et al. 2019) and Jianzhou bigeared goats (Zhan et al. 2016).
The regulation of the target genes in expression is one of the main functions of lncRNA. In lncRNA-mRNA interaction networks constructed in the study (Fig. 4A,  B), upregulated MSTRG.262.1, XR_001917125.1 and XR_001917605.1 in LC goats may regulate both myogenesis and adipogenesis as their some trans target genes GREM1, MEGF10, and DNER were involved in muscle cell differentiation and striated muscle cell differentiation process (P < 0.05), while other trans target genes CDK1 and RRM2 were enriched in p53 signaling pathway (P < 0.05) that played dual roles in muscle cell differentiation (Park et al. 2020) and adipogenesis (Krstic et al. 2018). In addition, MSTRG.262.1 would target MYOZ2 and ANKRD2, while XR_001917125.1 would target MYOZ2. The genes MYOZ2 and ANKRD2 were found to be upregulated in the same Longissimus dorsi muscle tissue of LC goats as those used in the study (Supplementary File 4) and were also reported to positively regulate muscle cell differentiation (Bean et al. 2008;Wang et al. 2018). It was therefore inferred that the three upregulated lncRNAs in LC goats may be responsible for the meat production and intramuscular fat content differences between LC goats and ZB goats.
It was notewor thy that MSTRG.262.1 and XR_001917125.1 also appeared in the interaction network related to meat tenderness. The lncRNAs were co-expressed with myosin light chain 6B (MYL6B) and collagen alpha-1(III) chain (LOC102176755) (Fig. 4C). MYL6B was a crucial muscle structure component, and its nucleotide sequence variations were found to affect beef tenderness (Rosa et al. 2017). LOC102176755 encodes a type of collagen that has a strong negative effect on meat tenderness (Dransfield et al. 2003). The GO analysis result further supported the effect of XR_001917125.1 on meat tenderness as its other two upregulated trans target genes (GREM1 and LOXL2) in Longissimus dorsi muscle of LC goats (Supplementary File 4) were involved in connective tissue development process (P < 0.05) (Supplementary File 6). Connective tissue is mainly composed of collagen, and its content was highly positively correlated with the shear force value of beef (r = 0.95) (Dransfield et al. 2003). In this context, the two lncRNAs may also regulate meat tenderness difference between the two breeds.
Two upregulated lncRNAs (XR_001917386.1, XR_001918614.1) and one downregulated lncRNA The green squares represent differentially expressed miRNAs sponged by these differentially expressed lncRNAs, while the blue circles represent differentially expressed genes targeted by these differentially expressed miRNAs (XR_001297059.2) in LC goats caught our attention as their trans target genes (SMPX and RYR3) were enriched in muscle contraction (P < 0.01) and muscle system process (P < 0.05) (Supplementary File 6). Of the two target genes, upregulated SMPX (Supplementary File 4) would be trans-regulated by upregulated XR_001917386.1 and XR_001918614.1 in LC goats. SMPX was also known as small muscle protein, in that it was a positive regulator of muscle fiber size (Kemp et al. 2001). We therefore speculated that the upregulation of SMPX in skeletal muscle of LC goats may be partly caused by trans-regulation of XR_001917386.1 and XR_001918614.1, eventually resulting in a higher muscle fiber size in LC goats. The lncRNA XR_001297059.2 was downregulated in LC goats with higher intramuscular fat content (Supplementary File 3). This was not surprising as its trans target gene RYR3 played a negative role in adipogenesis (Tsai et al. 2013) and was also downregulated in LC goats (Supplementary File 4).
In interaction network related to muscle development, NR5A2 was common trans target gene of seven differentially expressed lncRNAs (Fig. 4A). NR5A2 has been reported to regulate muscle morphogenesis and glucose metabolism in muscle cells (Bolado-Carrancio et al. 2014;Sheela et al. 2005). Similarly, RORC would be collectively trans-regulated by six differentially expressed lncRNAs in interaction network related to intramuscular fat deposition (Fig. 4B) and meat tenderness (Fig. 4C). RORC was enriched in connective tissue development process related to meat tenderness (P < 0.05) (Supplementary File 6). The variations in RORC were also associated with intramuscular fat and marbling score of cattle (Barendse et al. 2010). These suggest that these multiple lncRNAs may play roles in skeletal muscle development and intramuscular fat deposition in goats by collectively trans-regulating NR5A2 or RORC.
Other differentially expressed lncRNAs involved in the lncRNA-mRNA interaction networks may also regulate the phenotype differences between the two goat breeds in meat yield and meat quality. For example, upregulated MSTRG.12645.1 in Longissimus dorsi muscle of LC goats compared to ZB goats would be co-expressed with KLK7, which was also identified as an upregulated gene in LC goats (Supplementary File 4). KLK7 led to an increase in body fat mass in mice (Kunath et al. 2021). The upregulated XR_001918823.1 in LC goats would target SLC7A8 which was upregulated in Longissimus dorsi muscle of LC goats compared to ZB goats (Supplementary File 4). Increased expression of SLC7A8 in skeletal muscle of yak improved the conversion efficiency of protein intake from total mixed rations, by weaking muscle proteolysis, and enhancing amino acid absorption and self-proteins synthesis in muscle tissue. These proteins are crucial for muscle growth (Liu et al. 2021b). However, 15 cis target genes predicted for differentially expressed lncRNAs obtained in this study were not found in the three lncRNA-mRNA interaction networks. This suggests that these differentially expressed lncRNAs screened in this study play roles in skeletal muscle development and meat quality mainly through a trans-regulatory mechanism.
The lncRNA can also act as a sponge for miRNA to relieve the repression of the target mRNA by miRNA, with an accompanying increase in expression level of mRNAs in mammalian cells (Zhan et al. 2019;Cai et al. 2017). In the study, upregulated XR_001918832.1 in LC goats was predicted to be a sponge of miR-200b and miR-200c to regulate the expression of CRHBP (Fig. 5). The miR-200b was reported to suppress proliferation of C2C12 myoblast (Yao et al. 2013) and differentiation of ovine preadipocytes (Jin et al. 2021), while miR-200c inhibited differentiation of C2C12 myoblast (D'Agostino et al. 2018). It was therefore speculated that upregulated XR_001918832.1 in LC goats may sponge-absorb miR-200b/miR-200c to positively regulate caprine muscle development and intramuscular fat. Similarly, upregulated XR_001917125.1 in LC goats would increase the expression of NCAM2 by relieving the repression effect of miR-34-3p and miR-460-5p (Fig. 5). Given that NCAM2 promoted myogenesis of mice (Ronzoni et al. 2021), upregulation of XR_001917125.1 would result in a higher skeletal muscle mass in LC goats compared to ZB goats. These suggest that ceRNA mechanisms may partly explain the phenotype differences in skeletal muscle mass and intramuscular fat between LC goats and ZB goats. However, the lncRNA-miRNA-mRNA ceRNA networks predicted in the study need to be further corroborated.

Conclusions
In summary, we identified 136 differentially expressed lncRNAs in muscle tissues between LC goats and ZB goats. The target genes of the differentially expressed lncRNAs were related to phenotype differences in meat production, intramuscular fat deposition, and meat tenderness between the two caprine breeds. Several differentially expressed lncRNAs may responsible for the differences via ceRNA mechanism, including XR_001918832.1-miR-200b/miR-200c-CRHBP and XR_001917125.1-miR-34-3p/ miR-460-5p-NCAM2.