In this study, Chengkou mountain chicken embryos were used as experimental animals and were purchased from Chongqing Xuanpeng Agricultural Development Co., Ltd. A total of 200 eggs were incubated following the conventional incubation procedure (37.8°C, 55% humidity) and embryos were harvested on Day 12, 16, 19, and 21 of incubation. In order to maintain uniformity across samples from the four periods, we opted for leg muscles as the primary sample due to their ease of collection during the embryonic period. We ensured comprehensive collection of the leg muscles, obtaining three biological replicates per period and stored them at −80°C. Twelve embryonic skeletal muscle samples were fixed with 4% paraformaldehyde and stored at 4°C for histological observation.

After paraffin embedding, the tissue samples were stained with hematoxylin and eosin, and the muscle fibers were imaged using an OLYMPUS microscope imaging system (Olympus, Tokyo, Japan). For each picture, 30 muscle fibers were randomly selected (19).

Total RNA was extracted from the muscles of the four periods according to the instruction method of Trizol reagent (TaKaRa, Dalian, China), and the concentration and purity of RNA were determined using a NanoDrop microspectrophotometer (Thermo Fisher Scientific, MA, United States). After removal of ribosomal RNA (rRNA) (Epicentre, United States), linear RNA was digested with RNase R enzyme (QIAGEN, Germany) to obtain circular RNA. The fragmented RNA was then treated with a fragmentation buffer to generate short fragments. Random hexamers were used to synthesize the first strand of cDNA, and a buffer, dNTPs, RNase H, and DNA polymerase I (QIAGEN, Germany) were added to synthesize the second strand of circular RNA using the linear RNA as a template. Purification was performed using a QiaQuick PCR kit (QIAGEN, Germany), followed by elution with EB buffer, end repair, sequential addition of base A and sequencing adapter, and recovery of the target fragment by agarose gel electrophoresis. The target fragments were then amplified by PCR to complete the library preparation. The samples were subsequently sent to Gene Denovo Biotechnology Co., Ltd. (Guangzhou, China) for sequencing using Illumina HiSeqTM 2,500 (Illumina, CA, United States).

To ensure data quality, raw reads containing adapter sequences with more than 10% N were removed. The remaining reads were further filtered to remove low-quality reads, where the number of bases with a quality value of Q ≤ 10 accounted for more than 50% of the entire read. The resulting high-quality (HQ) clean reads were then matched to the ribosome database using Bowtie2 in the Hisat2 tool (20, 21). The ribosome-mapped reads were removed, leaving only the unmapped reads for subsequent analysis. These unmapped reads were then aligned to the reference genome. The Find_circ software was used to identify circRNAs, and the resulting circRNA identification results were filtered to obtain highly credible circRNAs (22). A difference in circRNA with a | log2 (FC) | > 1 and FDR (false discovery rate) < 0.05 is defined as significant.

Using DAVID25 software, we mapped all of the source genes to the Gene Ontology database.¹ We then calculated the number of genes associated with each GO entry through hypergeometric inspection, defining the background of source genes and determining significant enrichment of GO terms relative to the genome. KEGG pathway enrichment of source genes was performed using KOBAS v2.0.

Utilizing the Short Time-Series Expression Miner (STEM) software,² time-series clustering analysis of differentially expressed circRNAs was performed to elucidate various gene expression patterns during the embryonic stage. The polygenic screening required a minimum of 2 variants and allowed for a maximum of 20 trends. Data normalization was carried out using log2 (RPM), with a p < 0.05 considered indicative of a reliable trend.

The WGCNA R package was used to cluster genes with similar expression patterns, investigate the association between modules and specific phenotypes or traits, and examine the expression patterns of multiple genes (23). The soft threshold was set to 2, with all other settings remaining at their default values. A p-value of <0.05 was considered indicative of module association.

Utilizing the Cytoscape software,³ we drew ceRNA networks and predicted potential relationships among circRNAs, mRNAs, and miRNAs using the miRDB⁴ online database. The construction of the ceRNA network followed the Cytoscape network construction template, incorporating mRNA and miRNA data from our team’s previous research results (19, 24–26). The screening criteria in miRDB included only functional miRNAs, while excluding targets with a prediction score of less than 80, and miRNAs with more than 2,000 predicted targets in the genome.

Here, we selected six circRNAs to verify the sequencing results via RT-qPCR. Primers were designed using Primer Premier, as shown in Supplementary Table S3. The specific PCR method was based on previous experiments conducted by our team (19).

The data were processed using SPSS 20.0 (SPSS Inc., United States) software and presented as mean ± standard deviation (mean ± SD). A comparison between the two groups was conducted using the unpaired t-test and Duncan multiple range test. Graphs were generated using GraphPad Prism 9 (San Diego, CA, United States). p < 0.05 was considered statistically significant.

We stained sections of Chengkou mountain chicken embryonic muscles at four stages and compared the muscle fibers to determine muscle development at the respective stages (Figure 1). Hematoxylin-Eosin sections showed that on Day 12 (E12), the muscle fibers were compounded. By Day 16 (E16), the muscle fibers’ outline had emerged and their diameter had notably expanded. At Day 19 (E19), the differentiation of myofibers became more pronounced, and their diameter continued to increase. The basic structure of the muscle fibers was visible, becoming clearest by Day 21 (E21) when the fiber structure had fully matured and could be clearly identified. The analysis revealed a gradual increase in muscle fiber diameter over the course of time (Supplementary Figure S1). The HE results showed the basic rules of muscle growth in Chengkou mountain chickens, providing a preliminary basis for subsequent analyses.

The 12 sequenced 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) generated 1,334,509,224 high-quality clean reads after filtering out reads with adapter sequences, >10% N ratio, and low-quality bases (Supplementary Table S1). After filtering, the base quality of all samples (Q20 values) was greater than 98% and the GC content was approximately 45%, with a balanced base composition clean reads (Supplementary Figure S2). Approximately 93% of the reads were mapped to the chicken reference genome (Table 1). The suspected circRNA reads (anchor reads) were subsequently analyzed by comparing them with the reference genome using bowtie to identify circRNAs (Table 2).

A total of 9,904 novel circRNAs were identified from 3,273 genes, with 1,945 source genes producing only one circRNA and the most, 35, coming from the same gene (Figure 2A). CircRNAs were predominantly derived from exons, accounting for 78.27%, followed by introns (13.05%), intergenic regions (4.48%), and antisense circRNAs (4.20%) (Figure 2B). The identified circRNAs had a range of lengths, from 73 bp to 97,968 bp (Figure 2C). Most of these circRNAs were derived from chicken chromosome one (NC_006088.5) (Figure 2D).

A principal component analysis (PCA) showed that the 12 samples clustered according to the biological repeats of different time nodes (Figure 3A), indicating the biological repeatability and significance of between-group differences. There were 444 circRNAs expressed only at E12, 1,195 only at E16, 719 only at E19, and 980 only at E21 (Figure 3B). A total of 1,983 common circRNA-derived genes were identified across the four periods (Figure 3C), which were mainly associated with protein binding, intracellular parts, and cellular protein modification processes (Figure 3D). Pathway analysis revealed that these genes were mainly implicated in focal adhesion, the GnRH signaling pathway, and the MAPK signaling pathway (Figure 3E). A total of 917 novel circRNAs exhibiting differential expression (|log2(FC)| > 1, FDR < 0.05) were identified through pairwise comparison of the four periods (Figure 4A), including 240 in E12 vs. E16, 330 in E12 vs. E19, 450 in E12 vs. E21, 174 in E16 vs. E19, 320 in E16 vs. E21 and 118 in E19 vs. E21 (Supplementary Table S1). In comparison to E12, E16, and E19, E21 had 121, 150, and 196 up-regulated circRNAs, respectively, and 119, 180, and 254 down-regulated circRNAs, respectively, (Figure 4B). The differentially expressed circRNAs were analyzed for expression trend (the expression of circRNA was expressed in Reads of exon model per million mapped reads(RPM), and the data were normalized by the method of log2) to reveal the potential effect mode of circRNA on embryonic muscle development. The differential circRNAs were enriched in 20 trends, of which six were significant (p < 0.05) (Figure 4C). The circRNA expression trend timeline is shown in Figure 4D. The significant trends included two upward trends and three downward trends. In total, 225 circRNAs were enriched in the upward trend and 196 circRNAs were enriched in the downward trend (Supplementary Data Sheet 2).

The circRNA expression network was constructed to identify potential factors related to muscle development. Therefore, two soft thresholds were chosen to ensure that the weighted gene co-expression network analysis (WGCNA) module conformed to a scale-free distribution (Figure 5A). Thus, we marked the 17 WGCNA modules (ME) identified in this study with different colors (Figure 5B). There were correlations between genes within modules, as well as connections between different modules (Figure 5C). The modules corresponding to different embryonic muscle development stages were significantly different, indicating that circRNAs play different functions at different stages. ME Salmon had the most significant correlations with muscle development (r = 0.86, p < 0.05) in E12, ME Red in E16 (r = 0.60, p < 0.05), ME Grey60 in E19 (r = 0.73, p < 0.05), and ME MidnightBlue in E21 (r = 0.85, p < 0.05) (Figure 5D). Among the circRNAs of the most significantly associated WGCNA modules, the ten most enriched items in E12 ME Salmon were the cell leading edge, the protein serine/threonine phosphatase complex, the phosphatase complex, the lamellipodium, the Rad17 RFC-like complex, the extrinsic component of the neuronal dense core vesicle membrane, the extrinsic components of the dense core granule membrane, cytosol, intracellular, and intracellular parts (Figure 5E). A KEGG pathway analysis also revealed that pathways related to cell physiology and organism development were significantly enriched, including the MAPK signaling pathway, autophagy, the VEGF signaling pathway, the hedgehog signaling pathway, and the FoxO signaling pathway (Figure 5F).

Using Miranda and miRDB, we predicted the associations between circRNAs and microRNAs (miRNAs), as well as between miRNAs and mRNAs. After filtering out irrelevant circRNAs, a competitive endogenous RNA (ceRNA) network was constructed using 754 circRNAs, 68 miRNAs, 2,180 mRNAs, and 3,002 interaction nodes (Supplementary Figure S3). To further explore the effects of circRNAs on skeletal muscle development, this study established the interaction between circRNA and four muscle-related genes, namely MYH1D, MYH1B, CAPZA1, and PERM1 (Figure 6). Novel_circ_002810, novel_circ_005018, novel_circ_005654, novel_circ_005015, and novel_circ_003935 regulate pathways of nine miRNAs, including gga-miR-23b-3p, gga-miR-9-5p, gga-miR-1706, gga-miR-15c-5p, gga-miR-1644, gga-miR-218-5p, gga-miR-1808, gga-miR-1626-5p, and gga-miR-1773-3p.

Six randomly selected circRNAs were used for RT-qPCR, including novel_circ_002016, novel_circ_004526, novel_circ_005116, novel_circ_007014, novel_circ_001058, and novel_circ_005196. The consistency in expression patterns of circRNAs between circRNA-seq and RT-qPCR results confirms the accuracy of the sequencing (p < 0.05) (Figure 7).