Transcriptome analysis of mRNAs, lncRNAs, and miRNAs in the skeletal muscle of Tibetan chickens at different developmental stages

Introduction: As a valuable genetic resource, native birds can contribute to the sustainable development of animal production. Tibetan chickens, known for their special flavor, are one of the important local poultry breeds in the Qinghai–Tibet Plateau. However, Tibetan chickens have a slow growth rate and poor carcass traits compared with broilers. Although most of the research on Tibetan chickens focused on their hypoxic adaptation, there were fewer studies related to skeletal muscle development. Methods: Here, we performed the transcriptional sequencing of leg muscles from Tibetan chicken embryos at E (embryonic)10, E14, and E18. Results: In total, 1,600, 4,610, and 2,166 DE (differentially expressed) mRNAs, 210, 573, and 234 DE lncRNAs (long non-coding RNAs), and 52, 137, and 33 DE miRNAs (microRNAs) were detected between E10 and E14, E10 and E18, and E14 and E18, respectively. Functional prediction showed several DE mRNAs and the target mRNAs of DE lncRNAs and DE miRNAs were significantly enriched in sarcomere organization, actin cytoskeleton organization, myofibril, muscle fiber development, and other terms and pathways related to muscle growth and development. Finally, a lncRNA–miRNA–mRNA ceRNA (competing endogenous RNA) network associated with muscle growth and development, which contained 6 DE lncRNAs, 13 DE miRNAs, and 50 DE mRNAs, was constructed based on the screened DE RNAs by Gene Ontology (GO) enrichment. These DE RNAs may play a critical regulatory role in the skeletal muscle development of chickens. Discussion: The results provide a genomic resource for mRNAs, lncRNAs, and miRNAs potentially involved in the skeletal muscle development of chickens, which lay the foundation for further studies of the molecular mechanisms underlying skeletal muscle growth and development in Tibetan chickens.

Introduction: As a valuable genetic resource, native birds can contribute to the sustainable development of animal production. Tibetan chickens, known for their special flavor, are one of the important local poultry breeds in the Qinghai-Tibet Plateau. However, Tibetan chickens have a slow growth rate and poor carcass traits compared with broilers. Although most of the research on Tibetan chickens focused on their hypoxic adaptation, there were fewer studies related to skeletal muscle development.
Methods: Here, we performed the transcriptional sequencing of leg muscles from Tibetan chicken embryos at E (embryonic)10, E14, and E18.
Results: In total, 1, 600,4,610,and 2,166 DE (differentially expressed) mRNAs,210,573,, and 52, 137,and 33 DE miRNAs (microRNAs) were detected between E10 and E14, E10 and E18, and E14 and E18, respectively. Functional prediction showed several DE mRNAs and the target mRNAs of DE lncRNAs and DE miRNAs were significantly enriched in sarcomere organization, actin cytoskeleton organization, myofibril, muscle fiber development, and other terms and pathways related to muscle growth and development.
Finally, a lncRNA-miRNA-mRNA ceRNA (competing endogenous RNA) network associated with muscle growth and development, which contained 6 DE lncRNAs, 13 DE miRNAs, and 50 DE mRNAs, was constructed based on the screened DE RNAs by Gene Ontology (GO) enrichment. These DE RNAs may play a critical regulatory role in the skeletal muscle development of chickens.

Introduction
The Tibetan chicken is a unique Chinese native breed that is mainly distributed in the Qinghai-Tibet Plateau. It has excellent meat quality and flavor due to the rich content of essential amino acids and flavor amino acids. Tibetan chickens exhibit characteristics such as high adaptability, resistance to extensive management, and strong disease resistance, but with the challenges of slow growth rate and low meat yield (Tang et al., 2015;Chi et al., 2020). The yield and quality of poultry mainly depend on the development of skeletal muscle (Dransfield and Sosnicki, 1999). The number of myofibers in chickens is basically fixed at the embryonic stages, and the differentiation of primary and secondary fibers is completed by 3/4 of the incubation period (Chen et al., 2013). Therefore, it is indispensable to understand global gene expression during embryonic skeletal muscle growth and development for improving the growth rate and muscle mass of chickens.
Skeletal muscle growth and development are extremely complex processes. In addition to mRNAs, plenty of ncRNAs (non-coding RNAs) were shown to be involved in skeletal muscle growth and development in multiple species, including pigs, sheep, chickens, and cattle (Cai et al., 2017;Wang et al., 2022a;Zhan et al., 2022;Zhang et al., 2022). miRNAs can inhibit post-transcriptional gene expression by binding to specific mRNAs. For example, miR-2954 prevented myoblast differentiation into multinucleated myotubes via suppressing the expression of YY1 (YY1 transcription factor) gene in the process of chicken skeletal muscle development at the embryonic stage . lncRNA is highly abundant ncRNA (Guttman and Rinn, 2012) that has complex biological functions such as interfering with gene expression, interacting with proteins, acting as ceRNAs (competing endogenous RNAs), and encoding peptides (Matsumoto et al., 2017;Li et al., 2020;Wang et al., 2020;Zhang et al., 2020). In recent years, several studies have confirmed that lncRNAs are widely involved in many stages of muscle growth and development in livestock and poultry . lncRNA-Six1 affected muscle growth and development via encoding a micropeptide (Cai et al., 2017). lncRNAs work as ceRNAs, which is the most common regulatory function. For example, lncRNA MEG3, as the "sponge" of miRNA-135, promoted bovine skeletal muscle differentiation via interacting with miRNA-135 and MEF2C (myocyte enhancer factor 2C) (Liu et al., 2019), and this lncRNA was also found in pig skeletal muscle (Yu et al., 2018). lncRNA-125b promoted skeletal muscle satellite cell differentiation by functioning as a ceRNA for miR-125b to positively regulate IGF2 (insulin-like growth factor 2) expression (Zhan et al., 2019). To explore the role of coding RNAs and ncRNAs in the development of chicken skeletal muscle, we performed a systematic transcriptome analysis of mRNAs, lncRNAs, and miRNAs in the skeletal muscle of chicken embryos at various developmental stages.
For poultry, skeletal muscle development during the embryonic period plays a fateful role in the potential of post-incubation muscle development. The primary and secondary muscle fibers were formed in approximately 7-12 days, independent myotube formation and differentiation were initiated mainly from 12 to 18 days, and muscle fiber hypertrophy was observed during the post 1/4 of embryo incubation (Stockdale and Miller, 1987;Stickland et al., 2004;Picard et al., 2010). Therefore, in this study, the leg muscles from Tibetan chickens at E10-, E14-, and E18-day were selected for sequencing analysis. DE mRNAs, DE lncRNAs, and DE miRNAs were screened from three comparisons (E10 vs. E14, E10 vs. E18, and E14 vs. E18). Common DE RNAs, which were present in all three comparisons, were obtained through Venn analysis. Subsequently, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses were performed on common DE mRNAs, target genes of common DE lncRNAs, and common DE miRNAs. The ceRNA network contributes to exploring the functions and regulatory mechanisms of genes and facilitates a deeper understanding of many biological phenomena. So, we constructed the ceRNA regulatory networks associated with skeletal muscle growth and development. Some key factors were expected to be revealed for chicken muscle development. These findings will help us understand the regulation of coding and ncRNAs in chicken muscle growth and development.

Animal preparation and sample collection
In this experiment, we collected eggs from Tibetan chickens (TCs) obtained from Mao Xian Jiuding Original Ecological Livestock and Poultry Breeding Co., Ltd. (Aba Tibetan and Qiang Autonomous, Sichuan Province, China). The eggs were incubated at a humidity of 55% and a temperature of 37.8°C. At the target time, we dissected the chicken embryos and determined gender by observing the gonadal features of the fetus. All leg muscle samples from nine Tibetan chickens were collected from healthy male embryos at three developmental stages (E10, E14, and E18), with three samples in each stage. All collected tissues were immediately frozen in liquid nitrogen and subsequently stored at −80°C for RNA extraction and sequencing. All treatment and sample collection procedures were approved by the IACUC (Institutional Animal Care and Use Committee) of Southwest Minzu University (Sichuan, China) with approval number 2020MDLS44.

Total RNA isolation, library construction, and sequencing
Total RNA was isolated from nine leg muscle samples at three differential development stages using the TRIzol reagent (Vazyme, Nanjing, China). The quality of the RNA samples was detected by 1% agarose gel electrophoresis. The purity and concentration of the RNA samples were assessed using the NanoPhotometer ® N60 (Implen, Munich, Germany), and the integrity and RNA integrity number (RIN) values of the isolated RNAs were evaluated using the Agilent Bioanalyzer 2100 system (Agilent Technologies, Massy, France). The RIN values of all nine samples were greater than 8, and all could be used for subsequent analysis. The isolated RNAs that passed these tests were divided into three parts. One part was used for quantitative real-time PCR (RT-qPCR) by reverse transcription into cDNA strands using the PrimeScriptTM RT Reagent Kit with gDNA Eraser (Perfect Real Time) (Takara, Dalian, China), and the other two parts were used for long RNAseq (mRNA and lncRNA sequencing) and small RNA-seq (miRNA sequencing).

Data quality control and read mapping
Large volumes of raw data stored in the FASTQ format were obtained after sequencing. In order to ensure the reliability of the subsequent bioinformatics analysis, the raw data were first quality-controlled. Clean data were obtained from raw data using SeqPrep (https://github.com/jstjohn/SeqPrep) and Sickle (https://github.com/najoshi/sickle) software tools. Specific operations include removing the reads containing adapter, removing the reads containing poly-N (N rate >10%), and removing the low-quality reads (cutting the bases with quality (Q) < 20 at the end of the sequence (3′-end); if there are still bases with Q < 10 in the remaining sequence, the whole sequence will be eliminated; otherwise, it will be retained). For small RNA-seq, in addition to the aforementioned quality control operations, reads shorter than 18 nt and longer than 32 nt were filtered out. After the data quality control, the quality assessment was carried out, and the Q20, Q30, and GC content of clean data were calculated. All further analyses were carried out based on the high-quality, clean data.

Identification of lncRNAs
lncRNA sequences that have been previously reported were downloaded from lncRNA databases, mainly including NONCODE (http://www.noncode.org/index.php), Ensembl, and NCBI databases. Based on these databases, known lncRNAs were identified. For the identification of novel lncRNAs, the transcripts with class code "x (exonic overlap with the reference on the opposite strand)," "i (transfrag falling entirely within a reference intron)," "j (potentially novel transcript)," "u (unknown, intergenic transcript)," and "o (generic exonic overlap with a reference transcript)" were obtained using gffcompare software. On this basis, the transcripts whose length ≥200 bp, exon number ≥2, and ORF length ≤300 bp were screened as candidate lncRNAs. The encoding potential of candidate lncRNAs was screened using the analysis methods of coding potential calculator (CPC), coding-non-coding index (CNCI), coding potential assessment tool (CPAT), and Pfam protein domain analysis tools (Kong et al., 2007;Sun et al., 2013;Wang et al., 2013;Robert et al., 2014), and the screening conditions of these four methods are CPC score <0.5, CNCI score <0, CPAT score <0.5, and Pfam database without comments. The intersection of the results obtained from these four analysis methods represents the set of novel lncRNAs.
Quantitative analysis of lncRNA expression levels was performed using RSEM software (http://deweylab.github.io/ RSEM/). The expression of lncRNAs was measured by TPM (transcripts per million reads). Based on the expression details, a correlation analysis between biological replicate samples for long RNA-seq was performed using the Pearson correlation algorithm.

Identification of miRNAs
Known miRNAs were identified by aligning clean reads to pre-miRNAs (miRNA precursors) and mature miRNAs in the miRbase database (http://www.mirbase.org/). Simultaneously, ribosomal RNAs (rRNAs), transfer RNAs (tRNAs), small nuclear RNAs (snRNAs), small nucleolar RNAs (snoRNAs), and other ncRNAs and repeats were removed, and unannotated reads were obtained. According to the particular hairpin structure of pre-miRNAs, novel miRNAs were predicted using miRDeep2 software (https://www.mdc-berlin.de/content/ mirdeep2-documentation) (Friedländer et al., 2012). First, the unannotated reads were aligned with the reference genome, and the surrounding sequences were intercepted for secondary structure prediction. Next, according to the prediction results, the unannotated reads were filtered using information of the Dicer cutting site, energy values, and other features to identify novel miRNAs. The expression levels of miRNAs were also measured via TPM. Based on the expression details, a correlation analysis between biological replicate samples for small RNA-seq was performed using the Pearson correlation algorithm.
The lncRNA target genes of cis-acting were predicted by satisfying the requirement that the lncRNA affects the expression of surrounding genes (within 100 kb). When the expression of lncRNAs is significantly positively or negatively correlated with the expression of some genes, the lncRNA target genes of transacting were predicted by the correlation analysis of expression between lncRNAs and protein-encoding genes. The prediction of miRNA target genes was carried out using miRanda, TargetScan, and RNAhybrid software applications (Lewis et al., 2003;Betel et al., 2008).
GO (http://www.geneontology.org/) analyses were carried out for DE mRNAs, the target genes of DE lncRNAs, and DE miRNAs using GOATOOLS software. Furthermore, KEGG (http://www. genome.jp/kegg/) pathway enrichment analyses were performed for DE mRNAs, the target genes of DE lncRNAs, and DE miRNAs by Majorbio Bio-pharm Technology Co., Ltd. The terms and pathways with P adj < 0.05 were regarded as significant enrichments.

Construction of the ceRNA (lncRNA-miRNA-mRNA) network
CeRNAs are important and novel transcriptional regulators. The ceRNA network is one of the most important methods for posttranscriptional regulation of non-coding RNAs (Braga et al., 2020). The interactions of DE miRNA-DE mRNA and DE miRNA-DE lncRNA related to skeletal muscle growth and development were predicted using miRanda, TargetScan, and RNAhybrid software applications. Based on the ceRNA network theory, Cytoscape software (v3.6.0) was used to construct and visualize the ceRNA (lncRNA-miRNA-mRNA) network (Poliseno et al., 2010;Salmena et al., 2011).

Validation of DE RNAs by RT-qPCR
To verify the accuracy and reliability of differentially expressed RNAs (DE RNAs) screened by RNA-seq, six DE mRNAs, six DE lncRNAs, and six DE miRNAs were randomly selected for RT-qPCR. We used the PrimeScript ™ RT Reagent Kit with gDNA Eraser (Perfect Real Time) (Takara, Dalian, China) to convert RNAs to cDNAs. The random hexamers were used for DE mRNAs and DE lncRNAs, and stem-loop RT primers were used for DE miRNAs.
RT-qPCR was performed using the TB Green ® Premix Ex Taq ™ Ⅱ Kit (Takara, Dalian, China) according to the manufacturer's instructions. GAPDH for mRNAs and lncRNAs and U6 for miRNAs were used as endogenous controls. All primers used in RT-qPCR are designed using Primer Premier 5 software and exhibited in Supplementary Data S1. RT-qPCR was carried out in the LightCycler480 system (Roche, Swiss Confederation) with three replicates for each sample. The data were expressed as values relative to the E10 group. The relative expression levels of DE mRNAs, DE lncRNAs, and DE miRNAs were calculated using the 2 −ΔΔCT method, and one-way ANOVA in SPSS 20.0 was used to analyze the data. Differences were considered significant at p < 0.05.
All tools, databases, and related information used in this study are shown in Supplementary Data S2.

Overview of RNA-seq data of Tibetan chicken leg muscle at different embryonic stages
After the quality control for raw reads, a total of 1,051,297,502 clean reads for long RNA-seq and 97,422,413 clean reads for small RNA-seq were obtained from all samples. The percentages of Q20 and Q30 were above 97.8% and 93.91%, respectively. The average GC content of all samples was approximately 45.71%, and the average base error rate was below 0.0227% (Supplementary Data S3). The clean reads were mapped to the chicken reference genome with a range of 91.24%-94.68% for long RNA-seq and 96.32%-99.09% for small RNA-seq (Supplementary Data S4).
At the three development stages, there are 29,615 mRNAs, 12,399 lncRNAs (containing 11,014 known and 1,385 novel lncRNAs), and 1,216 miRNAs (containing 944 known and 272 novel miRNAs). The length of lncRNAs was mainly longer than 3,000 bp ( Figure 1A), and the majority of the lncRNAs were intergenic-lncRNAs (76.0%), followed by antisense-lncRNAs (12.3%), sense-exon-overlap lncRNAs (6.4%), bidirectional-lncRNAs (4.5%), and sense-intron-overlap lncRNAs (0.6%) ( Figure 1B). According to the characteristics of miRNAs, the length of miRNAs ranged from 17 to 28 bp, and most of the miRNAs (known and novel miRNAs) were between 20 and 24 bp long ( Figure 1C). The distribution of TPM values for mRNAs, lncRNAs, and miRNAs was found to be similar across the three different developmental stages in Tibetan chicken leg muscles ( Figures 1D-F). The expression levels of mRNAs, lncRNAs, and miRNAs were concentrated, and the results of inter-sample correlation analysis for long RNA-seq ( Figure 1G) and small RNA-seq ( Figure 1H) were valid (the correlation coefficients between the biological replicates were above 0.98), which suggested that these RNA-seq results are reliable for further analysis.
A total of 1,173 cis-and trans-target mRNAs were predicted for 210 DE lncRNAs between E10 and E14. The target mRNAs are mainly enriched in GO terms such as myofibril assembly, muscle fiber development, M band, structural constituent of muscle, sarcomere organization, and muscle cell development. Enriched KEGG pathways also included ECM-receptor interaction and PI3K-Akt signaling pathway; in addition, there are many target mRNAs enriched in the cGMP-PKG signaling pathway ( Figure 2B). For 52 DE miRNAs, we predicted 997 target mRNAs that were enriched in the top 20 GO terms including collagen-containing extracellular matrix, collagen trimer, and regulation of embryonic development. A total of 13 significantly enriched KEGG pathways were obtained from the target mRNAs of DE miRNAs, which mainly included ECM-receptor interaction, PI3K-Akt signaling pathway, and other pathways related to the regulation of diseases ( Figure 2C). Frontiers in Physiology frontiersin.org 05 Li et al. 10.3389/fphys.2023.1225349
A total of 3,474 cis-and trans-target mRNAs of 573 DE lncRNAs were mainly enriched in the top 20 GO terms, including DNA replication initiation, cellular respiration, forelimb morphogenesis, cytochrome complex, DNA origin binding, NADH dehydrogenase activity, and ATP biosynthetic process. In addition, the enriched KEGG pathways were mainly TCA cycle, cell cycle, ECM-receptor interaction, glycolysis/gluconeogenesis, steroid biosynthesis, and p53 signaling pathway ( Figure 3B). For 137 DE miRNAs, we predicted 4,678 target mRNAs that were mainly involved in GO terms such as myofibril assembly, regulation of vasculature development, regulation of embryonic development, and intermediate filament organization. These target mRNAs were Frontiers in Physiology frontiersin.org significantly enriched in the KEGG pathways of ECM-receptor interaction, TCA cycle, arginine biosynthesis, valine, leucine, and isoleucine degradation, steroid biosynthesis, and glycolysis/ gluconeogenesis ( Figure 3C).

Functional prediction of differentially expressed RNAs between E14 and E18
In comparison between E14 and E18, we screened 2,166 DE mRNAs ( Table 3; Supplementary Data S10.
Similarly, GO and KEGG functional analyses were performed. The top 20 significantly enriched GO terms and KEGG pathways of DE mRNAs are shown in Figure 4A. The GO terms mainly included ATP biosynthetic process, striated muscle cell development, muscle cell development, sarcomere organization, muscle contraction, and other terms related to cell respiration. The KEGG pathways of the regulation of diseases, TCA cycle, regulation of cardiac muscle, glycolysis/gluconeogenesis, and calcium signaling pathway were significantly enriched.
There were 1,561 cis-and trans-target mRNAs identified for 234 DE lncRNAs. Significantly enriched top 20 GO terms for the target mRNAs of DE lncRNAs primarily included sarcomere organization, ATP biosynthetic process, reactive oxygen species metabolic process, and other terms related to cell respiration. Significantly enriched top 20 KEGG pathways were primarily associated with the regulation of diseases, amino acid metabolism, TCA cycle, glycolysis/gluconeogenesis, and steroid biosynthesis ( Figure 4B). We predicted 625 target mRNAs for 33 DE miRNAs. The significantly enriched GO terms were collagencontaining extracellular matrix, sarcolemma, Z disc, actin binding, and vascular endothelial growth factor receptor Ⅰ binding. Furthermore, the KEGG pathways were ECM-receptor interaction, focal adhesion, protein digestion and absorption, oxidative phosphorylation, and PI3K-Akt signaling pathway ( Figure 4C).

Selection and functional prediction of common differentially expressed RNAs
To further explore the regulation mechanisms of coding and non-coding RNAs in chicken leg muscle growth and development, the common DE RNAs that were differentially expressed in all three The GO enrichment analysis for common DE mRNAs suggested that significantly enriched top 20 GO terms mainly included sarcomere organization, actomyosin structure organization, striated muscle cell development, actin cytoskeleton organization, muscle system process, muscle fiber development, Z disc, myofibril, and muscle contraction. The KEGG enrichment analysis revealed only eight significant KEGG pathways, which mainly included focal adhesion, ECM-receptor interaction, PI3K-Akt signaling pathway, and other pathways related to cardiac muscle ( Figure 5A).
The target mRNAs of common DE lncRNAs and DE miRNAs were predicted from the three comparisons. Significantly enriched top 20 GO terms for target mRNAs of common DE lncRNAs were primarily involved in sarcomere organization, muscle cell development, actomyosin structure organization, muscle fiber development, muscle contraction, regulation of skeletal muscle contraction, actin filament-based process, and sarcolemma. Furthermore, the significantly enriched KEGG pathways were similar to the pathways of common DE mRNAs ( Figure 5B). The target mRNAs of common DE miRNAs significantly enriched only three GO terms: sarcolemma, postsynaptic membrane, and cell junction. Significantly enriched KEGG pathways included Frontiers in Physiology frontiersin.org regulation of cardiac muscle, regulation of some diseases and cancer, ECM-receptor interaction, dorso-ventral axis formation, PI3K-Akt signaling pathway, progesterone-mediated oocyte maturation, taurine and hypotaurine metabolism, and oocyte meiosis ( Figure 5C).
3.6 Construction of the ceRNA regulatory network related to muscle growth and development DE mRNAs were selected from the significantly enriched GO terms associated with muscle growth and development. DE lncRNAs and DE miRNAs related to muscle growth and development were obtained through targeted relationships. Based on the three comparisons (E10 vs. E14, E10 vs. E18, and E14 vs. E18), we constructed three ceRNA networks related to muscle growth and development. In the comparison group E10 vs. E14, 205 DE mRNAs, 41 DE miRNAs, and 100 DE lncRNAs together constitute 370 lncRNA-miRNA-mRNA interaction pairs (Supplementary Figure S1; Supplementary Data S13). A total of 2,594 lncRNA-miRNA-mRNA interaction pairs were obtained in the ceRNA network of the comparison group E10 vs. E18 constructed with 936 DE mRNAs, 131 DE miRNAs, and 375 DE lncRNAs (Supplementary Figure S2; Supplementary Data S14). In the comparison group E14 vs. E18, the ceRNA network was constructed from 49 DE mRNAs, 19 DE miRNAs, and 53 DE mRNAs that yielded 78 lncRNA-miRNA-mRNA interaction pairs (Supplementary Figure S3; Supplementary Data S15). Considering that these ceRNA networks contain a large amount of information and each relationship cannot be shown in these figures, we constructed a mini-ceRNA network of common DE RNAs.

Validation of RNA-seq results by RT-qPCR
The relative expression levels of six DE mRNAs (MYL3, CSRP3, LDB3, MAT1A, KLHL31, and MYH1B) ( Figure 7A) (Figure 7B), and six DE miRNAs (miR-120a-5p, miR-1a-3p, miR-184-3p, miR-133a-3p, miR-1662, and miR-133a-5p) ( Figure 7C) were measured using RT-qPCR to validate the accuracy and reliability of the RNA-seq results. The results suggested that the expression trend of RT-qPCR for 18 DE RNAs at three developmental stages was generally consistent with the RNA-seq results. The aforementioned analysis indicated that the results of RNA-seq were reliable.

Discussion
Over the past 50 years, chickens have become the most consumed meat in the world (Petracci and Cavani, 2012).
Skeletal muscle is one of the most important components of the biological body, accounting for approximately 40% of the body weight (Güller and Russell, 2010). The growth and development of skeletal muscle have a significant impact on the yield of livestock and poultry meat products (Berri et al., 2007). Therefore, it is of great significance to conduct research on skeletal muscle growth and development in livestock and poultry. In this study, we selected three stages (E10, E14, and E18), which were important for chicken skeletal muscle, to explore the molecular mechanisms of muscle growth and development in chickens. The expression profiles of mRNAs, lncRNAs, and miRNAs at the three stages were detected by RNA-seq and bioinformatics analysis. The potential functions of DE mRNAs and target mRNAs of DE lncRNAs and DE miRNAs in the three comparisons (E10 vs. E14, E10 vs. E18, and E14 vs. E18) were revealed. Then, we constructed the ceRNA regulatory network related to chicken skeletal muscle growth and development using the common DE mRNAs, DE lncRNAs, and DE miRNAs. In summary, our data provided a comprehensive understanding of the molecular regulatory mechanisms underlying skeletal muscle growth and development in chickens. For the GO enrichment analysis of three comparisons, there were many GO terms related to skeletal muscle growth and development and limb formation in the comparison group E10 vs. E14, such as forelimb morphogenesis, sarcomere organization, muscle structure development, striated muscle cell development, and embryonic limb morphogenesis. In addition, many target mRNAs of DE lncRNAs were enriched in GO terms such as myofibril assembly, muscle fiber development, striated muscle contraction, M band, and muscle cell development. Many target mRNAs of DE miRNAs were enriched in the GO terms related to morphogenesis; the reason could be that during embryonic muscle development, the myogenic progenitor cells derived from the myotome continue to migrate and proliferate to the limbs, thereby promoting limb morphogenesis (Shi et al., 2022). These results suggested that the muscle had changed dramatically at E14 compared to E10 after continuous growth and development. Meanwhile, lncRNAs and miRNAs were extensively involved in muscle growth and development. Compared with the comparison group E10 vs. E14, the DE mRNAs and the target mRNAs of DE lncRNAs in the comparison groups E10 vs. E18 and E14 vs. E18 were mainly enriched in the GO terms associated with energy metabolism and mitochondrial electron transport chain besides the GO terms related to muscle growth and development. At the late stage of incubation, chicken embryos had frequent physiological activities and huge energy consumption. The mitochondrial electron transport chain is a metabolic pathway for organic molecules in organisms that helps them obtain energy (Kogot-Levin and Saada, 2014). Approximately 95% of the energy required for cellular activity comes from mitochondria (Chance et al., 1979). The growth and development of skeletal muscle were closely related to energy metabolism (Krischek et al., 2016). Some studies reported that insufficient energy could lead to weight loss in chickens and affect growth performance and muscle production in broilers (Nyo and Uni, 2010;Lamot et al., 2014).
For the KEGG enrichment analysis of three comparisons, ECM-receptor interaction, PI3K-Akt signaling pathway, TCA cycle, and glycolysis/gluconeogenesis were the most frequent pathways among all significantly enriched KEGG analyses. The main components of the ECM mainly included collagen, proteoglycan, and elastin (Yue, 2014). Apart from maintaining tissue structure, laminin in ECM was essential for early embryonic development and organogenesis (Rasmussen and Karsdal, 2019;Di Caprio and Bellas, 2020). The ECM had a beneficial effect on muscle cell differentiation, and ECM proteins Frontiers in Physiology frontiersin.org (specifically, laminin and entactin) could enhance myotube formation (Grefte et al., 2016). It was shown that in the absence of an ECM-cellular receptor interaction signal, the expression of myogenin could not drive skeletal muscle differentiation (Osses and Brandan, 2002). In addition, the ECM also promoted myogenic differentiation and maturation (Yang et al., 2016). The PI3K-Akt pathway was involved in various life activities of the organism, such as cell cycle and cell apoptosis, lipid metabolism, glucose homeostasis, and protein synthesis (Wang et al., 2016;Zhu et al., 2021). Several studies found that the PI3K-Akt pathway played a crucial role in myoblast proliferation and differentiation by mediating transduction of many growth factors (Ling et al., 2021;Lyu et al., 2021;Oh et al., 2021). For example, miR-9-5p inhibited the proliferation and differentiation of chicken skeletal muscle satellite cells by targeting IGF2BP3 (insulin-like growth factor-2 mRNA-binding protein 3) via the IGF2-PI3K/Akt signaling pathway (Yin et al., 2020). miR-485-5p had also been indicated to promote sheep myoblast proliferation by targeting PPP1R13B (protein phosphatase 1 regulatory subunit 13B) via the PI3K-Akt signaling pathway . Except for the mitochondrial electron transport chain, the pathways of TCA cycle and glycolysis/ gluconeogenesis are also the energy sources for skeletal muscle growth, development, and activity. They were present in the E10 vs. E18 and E14 vs. E18 comparison groups, which suggested the energy requirements increased as the skeletal muscle continued to develop and mature. Interestingly, many pathways related to disease occurrence were enriched, which may help us better understand the connections between coding and noncoding RNAs and these diseases. GO and KEGG enrichment analyses demonstrated that the skeletal muscle of Tibetan chickens showed significantly different degrees of development at three different stages. These results showed that there were many Frontiers in Physiology frontiersin.org coding and non-coding RNAs involved in skeletal muscle growth and development and provided a reference for understanding the overall view of the skeletal muscle growth and development process in Tibetan chickens.
To further explore the RNAs associated with skeletal muscle growth and development, we screened 531 common DE mRNAs, 16 common DE lncRNAs, and eight common DE miRNAs that were differentially expressed in all three comparisons. Subsequently, GO and KEGG enrichment analyses were performed for common DE mRNAs, target mRNAs of common DE lncRNAs, and DE miRNAs. As expected, plenty of GO terms related to muscle growth and development were significantly enriched. Similar to the aforementioned KEGG enrichment results, ECM-receptor interaction, PI3K-Akt signaling pathway, and several pathways linked with disease occurrence were mainly enriched. Therefore, numerous common DE mRNAs, common DE lncRNAs, and DE miRNAs regulated the growth and development of skeletal muscle. The ceRNA regulatory network is centered on miRNAs, linking mRNAs with lncRNAs/circRNAs to form a network of mutual influence and regulation (Poliseno et al., 2010;Salmena et al., 2011). Recently, numerous studies have shown that lncRNA-miRNA-mRNA is involved in the regulation of skeletal muscle growth and development (Cui et al., 2022;Shen et al., 2023). Therefore, we constructed a ceRNA regulatory network associated with skeletal muscle growth and development to better clarify the molecular mechanisms underlying skeletal muscle growth and development in chickens.

FIGURE 6
Analysis of the ceRNA (lncRNA-miRNA-mRNA) regulatory network. The blue circle, red V, and green rectangle notes represent the common DE mRNAs, DE miRNAs, and common DE lncRNAs, respectively.
Frontiers in Physiology frontiersin.org poultry (Liu et al., 2018;Ren et al., 2021). Similarly, CACNA1S played a role in the regulation of skeletal muscle growth and development, such as muscle contraction, muscle system process, T-tubule, and I band, in this study. So, in the ceRNA network, the lnc-45.2-gga-let-7-TPM2 and lnc-45.2-gga-let-7-CACNA1S interaction pairs may be involved in the formation of myofibril and sarcomere and regulate muscle contraction. lnc-778 showed a trend of upregulation and then downregulation in the three developmental stages of chickens, which indicated that skeletal muscle growth and development were extremely complex Validation of the RNA-seq data using RT-qPCR. The expression levels of six DE mRNAs (A), six DE lncRNAs (B), and six DE miRNAs (C) were validated with RT-qPCR in Tibetan chicken embryo leg muscles at three developmental stages.
Frontiers in Physiology frontiersin.org 14 processes. It could bind to multiple miRNAs (gga-miR-1458, gga-miR-18b-3p, gga-miR-29b-3p, and gga-miR-22-3p). It had been reported that gga-miR-22-3p inhibited proliferation and promoted differentiation of skeletal muscle cells in porcine (Dang et al., 2020), sheep (Wang et al., 2022b), and C2C12 cells . The ceRNA network analysis showed that lnc-778 acted as a sponge for gga-miR-22-3p to increase the expression of CASQ (calsequestrin 1) and AR (androgen receptor) and suppress the expression of RHOBTB1 (Rho-related BTB domain containing 1). Some studies suggested that CASQ was involved in muscle growth and development (Hou et al., 2021;Yang et al., 2021). AR affected the growth and development of sarcomere myofibrillar organization (Ghaibour et al., 2023). RHOBTB1 has a function in actin filament system organization and cell proliferation regulation (Xiao et al., 2017;Silva et al., 2020). Therefore, the regulation of skeletal muscle cell proliferation and differentiation may be the result of the co-regulation of three interaction pairs composed of lnc-778, gga-miR-22-3p, CASQ, AR, and RHOBTB1 in chickens. In our study, these miRNAs targeted many genes that were involved in the GO terms associated with sarcomere organization, muscle contraction, myofibril assembly, and regulation of the MAPK cascade.
For upregulated lncRNAs, lnc74.3 acted as a sponge for gga-miR-1769-3p to increase the expression of ENSGALG00000042388 and DUSP10 (dual specificity phosphatase 10), among which the DUSP10 gene was a member of the dual-specificity MAPK phosphatase (DUSP) enzyme family. Various DUSP enzymes were vital in the regulation of mitogen-activated protein kinase (MAPK) signaling pathways; the activation states of MAPKs were primarily regulated by a family of DUSPs (Gém et al., 2021). Many studies found that MAPK signaling pathways contributed to muscle growth and development (Catterall, 2000;Dulhunty, 2006). Therefore, DUSP10 may be involved in the regulation of skeletal muscle growth and development via the MAPK signaling pathway. lnc74.3 was also linked with NEXN (nexilin F-actin binding protein), MYOCD (myocardin), CLIP4 (CAP-Gly domain containing linker protein family member 4), ASPN (asporin), and CACNA1S by targeting gga-miR-145-5p. NEXN was mainly related to actin binding, MYOCD was involved in the regulation of muscle contraction and muscle system processes, CLIP4 was related to cytoskeletal development, and ASPN was involved in calcium ion binding and tissue development. lnc-019 could bind to gga-miR-302d, upregulating CAV2 (caveolin 2), MYBPC1 (myosin-binding protein C1), and ASPN and downregulating HSPH1 (heat shock protein family H (Hsp110) member 1), MYH15 (myosin heavy chain 15), and LAMA1 (laminin subunit alpha 1). These mRNAs were all associated with muscle growth and development. For example, the myosin-binding protein C translated by MYBPC1 was present in the A band of the sarcomere, regulating the cycling of actomyosin cross-bridges during muscle fiber contraction by interacting with the thick and thin filaments of the sarcomere (Geist and Kontrogianni, 2016). It is suggested that MYBPC1 affects skeletal muscle satellite cell proliferation and differentiation (Velleman et al., 2021). Therefore, lnc-019 may regulate skeletal muscle growth and development through lnc-019-gga-miR-302d-MYBPC1. It is worth mentioning that other DE lncRNAs, including lnc-73.4, lnc-778, and lnc-75.25, which were involved in the ceRNA regulatory network, also played important roles in the regulation of skeletal muscle growth and development, except for the previously mentioned ceRNA regulation modalities. As a ceRNA for gga-miR-2184a-5p and gga-miR-449d-5p, lnc-73.4 regulated the expression of MATN2 (matrilin 2) and TNNI2 (troponin I2, fast skeletal type). MATN2 played an important role in the formation of primary and secondary myofibers . Nihashi et al. (2019) selected TNNI2 as a candidate gene regulating myoblast proliferation and differentiation in chickens. Later, Yoshimoto et al. (2020) found that the expression of TNNI2 gradually increased with muscle regeneration. Finally, lnc-75.25 acted as a "sponge" for the novel miRNA 1_4798, which could target 14 genes. According to the GO enrichment analysis, these genes were mainly related to the regulation of actomyosin, striated muscle, myofibril development, and regulation of calcium ion transport and binding, and most genes were upregulated. We speculated that this novel miRNA may be involved in skeletal muscle growth and development by regulating gene expression.

Conclusion
In this study, we used RNA-seq to systematically analyze the expression of mRNAs, lncRNAs, and miRNAs in the embryonic leg muscles of Tibetan chickens at three developmental stages. The functions of DE mRNAs, DE lncRNAs, and DE miRNAs in three comparisons were annotated by GO and KEGG enrichment analyses. The expression levels of various RNAs during skeletal muscle growth and development in chicken embryos were comprehensively revealed, which provided a genomic resource for a better understanding of the molecular mechanisms underlying skeletal muscle growth and development in the embryonic stage of chickens. In addition, we constructed a ceRNA network related to skeletal muscle growth and development and revealed several candidate lncRNA-miRNA-mRNA interaction pairs that may regulate skeletal muscle growth in chickens. However, the specific role of these candidate RNAs needs to be further validated. These results provide a theoretical basis for further studies on the mechanisms underlying skeletal muscle growth and development in chickens.

Data availability statement
The datasets presented in this study can be found in online repositories. The name of the repository/repositories and accession Frontiers in Physiology frontiersin.org