Histological characteristics of muscle
To assess the muscle development regulation in the embryo of the Chengkou Mountain Chicken, we obtained embryonal muscle data at multiple time points. Muscle fibers have stage characteristics in the embryonic period, and myofiber morphology is a significant difference at different stages of embryonic development. On the 12th day of the embryo development, muscle fiber's complete shape had not been formed, and the outline of muscle fiber was not clear (Fig. 1A). With time, the cross-section of myofiber revealed a complete structure. The intervals between muscle bundles were gradually clear and distinct. At the same time, the structure of myofiber tended to mature (Fig. 1B- D). The muscle fiber surface of the embryo at E19 (8.01 ± 0.59) was significantly larger (P < 0.01) than the embryo at E16 (6.27 ± 0.50), and E21 (11.17 ± 0.87) was significantly larger (P < 0.01) than E19 (Fig. 1E). The cross-section area of muscle fibers presented a trend similar to diameter (Fig. 1F).
Overview Of Rna-sequencing
To obtain complete and accurate mRNA transcripts of the chicken embryo, we constructed 12 cDNA libraries (E12-1, E12-2, E12-3, E16-1, E16-2, E16-3, E19-1, E19-2, E19-3, E21-1, E21-2, and E21-3) from embryo leg muscle. A total of 1,337,535,812 raw reads were generated from 12 cDNA libraries. Clean reads totaling 1,334,509,224 were obtained after filtering out adaptor, N ratio greater than 10% reads, base reads, and low-quality reads. The percentage of clean reads for each duplicate was greater than 99% (Supplementary Table 1). With an error rate of 0.1%, more than 94% of bases were accurately identified (Supplementary Table 2). A comparison of reference area statistics showed that approximately 50%-60% of reads match the exon region (Supplementary Table 3). Mapping the sequence of the chicken's reference genome, which was about 5% did not match the genome sequence (Supplementary Table 4). About 80% of transcripts had high gene coverage (Supplementary Figure S1). All samples were distributed randomly and uniformly, and the number of genes showed trends to saturation (Supplementary Figure S2).
Using RNA-seq, 16,330 mRNAs transcripts were detected, including 109 novel mRNAs transcripts. Transcript expression was presented by FPKM (fragments per kilo-base of exon per million fragments mapped) value. The FPKM distribution of mRNAs is shown in Fig. 2A, whereas the expression of different samples is shown as a violin chart (Fig. 2B). To effectively find the most "main" element and structure in the data, the complex sample composition relationship was reflected on the two characteristic values of the horizontal and vertical coordinates. This aided in exploring the distance relationship between samples. The 12 samples were divided into four parts, which showed satisfactory repeatability (Fig. 2C). Then, we established a relationship cluster diagram to reflect the relationship between samples (Fig. 2D) intuitively. Sequences showed a reliable clustering effect, which ensured the veracity of the subsequent analysis.
Analysis Of Differentially Expressed Genes (degs)
We used p-value < 0.05 and fold change > 2 as the criteria to screen for differential genes by comparing pairwise differences at four time points during embryo muscle development. A total of 6,726 differentially expressed mRNAs were identified, including 2,251 in E12vsE16, 4,324 in E12vsE19, 5224 in E12vsE21, 1,274 in E16vsE19, 2,735 in E16vsE21, and 857 in E19vsE21 (Figure S3). And the number of DEGs at different time points was summarized in Fig. 3A. Through cluster analysis, we further revealed the differential expression of genes in different periods (Fig. 3B). To identify genes that play a key role in muscle development throughout the embryonic period, we performed Venn on genes at different stages. A total of 50 key genes were found in the intersection, generated from the Venn of DEGs (Fig. 3C).
Sample Time Series Of Degs
To comprehensively reveal muscle development status at different time points, we analyzed the expression trend of differential genes and selected more biologically meaningful target genes with P < 0.05 as the screening condition. Differential genes were enriched into 20 trends, of which six trends appeared significant (P < 0.05); notably, most of the genes were significantly enriched (Fig. 4A). The time-series line of differential gene expression is shown in Fig. 4B. The significant enrichment trend contains most of the differentially expressed genes. The overall gene expression trend was classified as either rising or falling. The significantly enriched trends included 4 down-regulation and 2 up-regulation trends. A total of 3,963 DEGs were significantly enriched in down-regulation trends (Profile 0, profile 2, profile 7, and profile 1), whereas 1,353 DEGs were enriched in up-regulation trends (profile 19 and 17). These findings effectively revealed the gene expression status of muscle development in the middle and late embryo stages.
Function Enrichment Analysis Of The Degs In Descending Trend
Using the GO enrichment analysis, we explored the function of the target genes. The top 20 enrichment terms in the three sections (Cellular Component, Molecular Function, Biological Process) are displayed in Fig. 5A-C. Cellular components contained 91 significance terms (P < 0.05), such as non-membrane-bounded organelle, intracellular non-membrane-bounded organelle, microtubule cytoskeleton, cell, and cell part. In total, 53 significant terms were enriched, for example, c DNA-dependent ATPase activity, motor activity, cytoskeletal protein binding, and helicase activity. And biological processes involved 248 significance terms; the top 5 terms included cellular component organization or biogenesis, organelle organization, cellular component organization, DNA metabolic process, and cell cycle.
The KEGG enrichment analysis of the DEGs is shown in Fig. 5D, including the top 20 KEGG pathways. A total of 183 KEGG pathway terms were enriched, including 17 significant terms, for instance, Cell cycle, ECM-receptor interaction, Necroptosis. Numerous muscle development pathways were reported, including the Wnt signaling pathway (P < 0.05), MAPK signaling pathway, TGF-beta signaling pathway, and mTOR signaling pathway.
We comprehensively analyzed 50 genes selected for development at different times to identify the key mechanism by which they contribute to development. Multiple genes were enriched in muscle development-related pathways, including regulation of muscle system process, regulation of muscle contraction, and muscle system process. Furthermore, KEGG pathway analysis of 50 genes revealed that eight genes were enriched in 10 pathways, among which H2A was found to be associated with three biological pathways (Table 1).
Table 1
50 Key genes KEGG pathways
Pathway | Pathway ID | K_ id | Differentially expressed genes |
Primary immunodeficiency | ko05340 | K03648 | UNG, UDG; uracil-DNA glycosylase |
MicroRNAs in cancer | ko05206 | K04381 | STMN1; stathmin |
Homologous recombination | ko03440 | K10877 | RAD54B; DNA repair and recombination protein RAD54B |
Estrogen signaling pathway | ko04915 | K09571 | FKBP4_5; FK506-binding protein 4/5 |
Base excision repair | ko03410 | K03648 | UNG, UDG; uracil-DNA glycosylase |
p53 signaling pathway | ko04115 | K10129 | GTSE1, B99; G-2 and S-phase expressed protein 1 |
Alcoholism | ko05034 | K11251 | H2A; histone H2A |
Necroptosis | ko04217 | K11251 | H2A; histone H2A |
Systemic lupus erythematosus | ko05322 | K11251 | H2A; histone H2A |
MAPK signaling pathway | ko04010 | K04381 | STMN1; stathmin |
Validation Of Candidate Genes
To reveal the key genes associated with embryonic muscle development, we screened several genes with higher expression levels among the 50 key differentially expressed genes, including MYH1F, SLC25A12 (up-regulation), STMN1, VASH2, and TUBAL3 (down-regulation). Similarly, HADHB was picked out from the 109 novel genes. Upon conducting q-PCR verification on the selected differential genes, we generated consistent findings as with RNA-sEq. Then, the candidate genes were verified via qRT-PCR (Fig. 6). Notably, similar results were reported as those obtained through sequencing, which confirmed the reliability of the sequencing data.