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 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 RNA-seq (mRNA and lncRNA sequencing) and small RNA-seq (miRNA sequencing).

For long RNA-seq, 9 μg of total RNA in each sample was used to remove ribosomal RNA (rRNA) using the Ribo-Zero™ rRNA Removal Kit (Epicentre, Madison, WI, United States). The total RNA of nine samples (TC_E10-1, TC_E10-2, TC_E10-3, TC_E14-1, TC_E14-2, TC_E14-3, TC_E18-1, TC_E18-2, and TC_E18-3) from which rRNA was removed was used as the input material. According to the NEBNext^® Ultra™ Directional RNA Library Prep Kit for Illumina (NEB, United States), nine cDNA libraries were constructed from RNAs obtained in the previous step, and three biological replicates per developmental stage were prepared. Finally, the cDNA libraries were sequenced on Illumina HiSeq 2500 (Majorbio Bio-pharm Technology Co., Ltd., Shanghai, China), and 125-bp paired-end reads were generated. For small RNA-seq, approximately 3 μg of total RNA per sample (nine samples) was used for preparing the libraries following the manufacturer’s recommendations in the SMARTer smRNA-Seq Kit for Illumina (Clontech, United States). Subsequently, ligation of the RNA 3′- and 5′-adapters, reverse transcription of the synthetic first strand, PCR amplification, and PAGE gel purification were performed. The cDNA libraries for small RNA-seq were constructed. Finally, the libraries were also sequenced on Illumina HiSeq 2500, and 50-bp single-end reads were generated.

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.

For long RNA-seq and small RNA-seq, quality-qualified clean data (reads) were mapped to the chicken reference genome (the resource of the reference genome: http://asia.ensembl.org/Gallus_gallus/Info/Index, the version of the reference genome: GRCg6a) using HISAT2 (http://ccb.jhu.edu/software/hisat2/index.shtml) (Kim et al., 2019) and Bowtie2 (v2.2.9) (Langmead and Salzberg, 2012), respectively. Furthermore, the possibility of clean reads that span exons was considered, and many clean reads that failed to align were segmented for re-mapping using TopHat2 (http://ccb.jhu.edu/software/tophat/index.shtml). The mapped reads were assembled using Cufflinks (http://cole-trapnell-lab.github.io/cufflinks/) and StringTie (http://ccb.jhu.edu/software/stringtie/) software applications (Pertea et al., 2015).

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.

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.

Differentially expressed mRNAs, lncRNAs, and miRNAs were analyzed using DESeq2, DEGseq, and edgeR software programs, for which the filter criteria are P _adj [corrected p-value by BH (fdr correction with Benjamini/Hochberg)] < 0.05, P _adj < 0.001, and P _adj < 0.05, respectively, and |log₂(fold change)| ≥ 1 (Love et al., 2014).

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 trans-acting 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.

CeRNAs are important and novel transcriptional regulators. The ceRNA network is one of the most important methods for post-transcriptional 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).

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.

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.

Three comparisons (E10 vs. E14, E10 vs. E18, and E14 vs. E18) between E10, E14, and E18 were performed during the skeletal muscle growth and development of Tibetan chickens. In E10 vs. E14, 1,600 DE mRNAs (Figure 2A), 210 DE lncRNAs (Figure 2B), and 52 DE miRNAs (Figure 2C) were found, among which 931 DE mRNAs, 148 DE lncRNAs, and 21 DE miRNAs were upregulated and 669 DE mRNAs, 62 DE lncRNAs, and 31 DE miRNAs were downregulated (Supplementary Data S5). The top 20 upregulated and downregulated mRNAs, lncRNAs, and miRNAs based on log₂(fold change) in E10 vs. E14 are summarized in Table 1 and Supplementary Data S6.

GO and KEGG functional enrichment analyses were carried out with DE mRNAs and target mRNAs of DE lncRNAs and DE miRNAs between E10 and E14. Significantly enriched GO terms and KEGG pathways (P _adj ≤ 0.05) of the top 20 are displayed in Figures 2A–C. In total, 1,600 DE mRNAs are mainly enriched in GO terms such as sarcomere organization, muscle structure development, embryonic limb morphogenesis, muscle contraction, and muscle fiber development. The KEGG enrichment analysis revealed the involvement of DE mRNAs in extracellular matrix (ECM)–receptor interaction, PI3K-Akt signaling pathway, and other pathways (Figure 2A).

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).

A total of 4,610 DE mRNAs (Figure 3A), 573 DE lncRNAs (Figure 3B), and 137 DE miRNAs (Figure 3C) were identified in the comparison group E10 vs. E18. Of these, 2,390 DE mRNAs, 393 DE lncRNAs, and 69 DE miRNAs were upregulated, and 2,220 DE mRNAs, 180 DE lncRNAs, and 68 DE miRNAs were downregulated (Supplementary Data S7). The top 20 upregulated and downregulated mRNAs, lncRNAs, and miRNAs based on log₂(fold change) in E10 vs. E18 are summarized in Table 2 and Supplementary Data S8.

For 4,610 DE mRNAs, significantly enriched GO terms (P _adj ≤ 0.05) in the top 20 mainly contained striated muscle contraction, embryonic forelimb morphogenesis, intermediate filament organization, sarcomere organization, and ATP biosynthetic process. The significantly enriched KEGG pathways were the citrate cycle (TCA cycle), DNA replication, valine, leucine, and isoleucine degradation, glycolysis/gluconeogenesis, and some pathways associated with diseases (Figure 3A).

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 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).

In comparison between E14 and E18, we screened 2,166 DE mRNAs (1,255 upregulation and 911 downregulation), 234 DE lncRNAs (126 upregulation and 108 downregulation), and 33 DE miRNAs (22 upregulation and 11 downregulation) (Figures 4A–C; Supplementary Data S9). The top 20 upregulated and downregulated mRNAs, lncRNAs, and miRNAs based on log₂(fold change) in E14 vs. E18 are summarized in 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 collagen-containing 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).

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 comparisons (E10 vs. E14, E10 vs. E18, and E14 vs. E18) were screened by Venn analysis. Therefore, 531 common DE mRNAs, 16 common DE lncRNAs, and eight common DE miRNAs were obtained (Figures 5A–C; Supplementary Data S11). Subsequently, we performed GO and KEGG enrichment analyses for these common DE mRNAs and target DE mRNAs during leg muscle growth and development in Tibetan chickens. For common DE mRNAs, target mRNAs of common DE lncRNAs, and common DE miRNAs, a total of 129, 118, and 3 GO terms and 8, 9, and 16 KEGG pathways were significantly enriched, respectively (Supplementary Data S12).

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 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).

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.

The ceRNA regulatory network related to muscle growth and development was constructed from 50 common DE mRNAs, 13 DE miRNAs (gga-miR-1458, gga-let-7b, gga-miR-1769-3p, gga-miR-145-5p, gga-let-7j-5p, gga-miR-18b-3p, gga-let-7a-5p, gga-miR-29b-3p, gga-miR-22-3p, gga-miR-302d, gga-miR-2184a-5p, gga-miR-449d-5p, and 1_4798), and six common DE lncRNAs (ENSGALT00000097778, NONGGAT000245.2, MSTRG.5174.3, ENSGALT00000096019, MSTRG.1173.4, and MSTRG.7875.25) that yielded 67 pairs of candidate ceRNAs (lncRNA–miRNA–mRNA) (Figure 6; Supplementary Data S16).

The relative expression levels of six DE mRNAs (MYL3, CSRP3, LDB3, MAT1A, KLHL31, and MYH1B) (Figure 7A), six DE lncRNAs (ENSGALT00000102361, ENSGALT00000104754, NONGGAT000245.2, ENSGALT00000107014, MSTRG.22805.1, and MSTRG.6624.1) (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.