Genome-wide characteristics and potential functions of circular RNAs from the embryo muscle development in Chengkou mountain chicken

The Chengkou mountain chicken, a native Chinese poultry breed, holds significant importance in the country’s poultry sector due to its delectable meat and robust stress tolerance. Muscle growth and development are pivotal characteristics in poultry breeding, with muscle fiber development during the embryonic period crucial for determining inherent muscle growth potential. Extensive evidence indicates that non-coding RNAs (ncRNAs) play a regulatory role in muscle growth and development. Among ncRNAs, circular RNAs (circRNAs), characterized by a closed-loop structure, have been shown to modulate biological processes through the regulation of microRNAs (miRNAs). This study seeks to identify and characterize the spatiotemporal-specific expression of circRNAs during embryonic muscle development in Chengkou mountain chicken, and to construct the potential regulatory network of circRNAs-miRNA-mRNAs. The muscle fibers of HE-stained sections became more distinct, and their boundaries were more defined over time. Subsequent RNA sequencing of 12 samples from four periods generated 9,904 novel circRNAs, including 917 differentially expressed circRNAs. The weighted gene co-expression network analysis (WGCNA)-identified circRNA source genes significantly enriched pathways related to cell fraction, cell growth, and muscle fiber growth regulation. Furthermore, a competitive endogenous RNA (ceRNA) network constructed using combined data of present and previous differentially expressed circRNAs, miRNA, and mRNA revealed that several circRNA transcripts regulate MYH1D, MYH1B, CAPZA1, and PERM1 proteins. These findings provide insight into the potential pathways and mechanisms through which circRNAs regulate embryonic muscle development in poultry, a theoretical support for trait improvement in domestic chickens.


Introduction
Chicken meat is a prominent consumer livestock product globally, with its consumption steadily rising (1).In China, an increasing number of consumers favor local meat breeds.However, local breeds such as the Chengkou mountain chicken are encountering challenges in meeting the demand for high-quality poultry due to their slow growth and development.Consequently, the focus of China's poultry industry has always been on improving chicken production and ensuring high-quality meat (2).
Skeletal muscle, as the predominant constituent of animal meat products, significantly influences poultry meat yield by virtue of its growth and development (3).The activation of skeletal myogenesis is governed by the myokines Myf5, which are present in cells located in the dorsomedial portion of the somites (4).Upon interaction with neural crest cells carrying Wnt1, Myf5 triggers expression in the dermis, leading to the downregulation of Pax3 expression and the subsequent formation of the primary myotome (5).The process of muscle fiber growth and development is intricate and can be categorized into three stages: early embryonic, late embryonic, and postnatal (6).During the early embryonic stage, monocytic muscle cells, originating from Pax3 + dermal cell progenitors, align along the entire cranio-caudal length of the somite to form the initial sarcomere (7).Subsequently, in the late embryonic stage, embryonic myoblasts fuse into myotubes, and myogenic progenitors differentiate into multinucleated myofibers under the regulation of MRF4 (8).As the developmental process unfolds, muscle cells produce primary and secondary muscle fibers, with satellite cells gradually forming to supply nuclei for the growing muscle fibers (9).After birth, the muscle fibers mature, and satellite cells enter a quiescent state, reactivating only to repair damaged muscle fibers (10).Enhanced insight into the intricate developmental regulatory network of embryonic skeletal muscle can offer more precise guidance in poultry breeding, thereby improving the economic benefits of breeding.
Circular RNAs (circRNAs) discovered in 1976 are single-stranded circular transcripts that are important in muscle development (11,12).Trans-splicing without a 5′ cap and 3′poly (A) tail produces circRNAs, which have a longer half-life than linear mRNA (13).Ouyang et al. (14) showed that circRNAs are abundant and dynamically expressed during chick embryonic muscle development.CircRNAs affect cell proliferation in skeletal muscle development (15).An in vivo experiment showed that circZfp609, a gene marker of ZNF609, inhibits myogenic differentiation by sponging miR-194-5p (16).The isolated circZfp609 inhibits BCLAF1 and affects the expression of monoclonal antibody to Myf5 and MyoG.Chen et al. (17) showed that overexpressing circCLTH promotes muscle cell differentiation and fusion in buffalo.Furthermore, Liu et al. (18) showed that circARID1A regulates skeletal muscle regeneration in mice by acting as a sponge for miR-6368.Undoubtedly, circRNAs play a pivotal role in muscle development, and investigating their regulatory mechanism is imperative for enhancing poultry meat production performance.
The identification and regulatory functions of circRNAs in Chengkou mountain chicken remain poorly understood.In a previous study, we examined the mechanism of embryonic muscle development in Chengkou mountain chicken using transcriptomes (19).In this study, we aim to elucidate the specific functions of circRNAs in muscle development and explore potential regulatory pathways by analyzing differentially expressed circRNAs at various stages of embryonic muscle development and predicting associated regulatory pathways.

Experimental animals and materials
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.

Phenotypic identification of different stages of muscle development
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).

cDNA library construction and circRNA sequencing
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).

Identification and statistics of circRNAs
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.

Function enrichment analysis of differentially expressed circRNAs
Using DAVID25 software, we mapped all of the source genes to the Gene Ontology database. 1We 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.

Time series analysis
Utilizing the Short Time-Series Expression Miner (STEM) software,2 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.

WGCNA analysis
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.

The ceRNA network construction
Utilizing the Cytoscape software, 3 we drew ceRNA networks and predicted potential relationships among circRNAs, mRNAs, and miRNAs using the miRDB4 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][25][26].The screening criteria in miRDB included only functional miRNAs, while excluding targets with a

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

Statistical analysis
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.

Phenotypes at different stages of muscle development
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.

Differential expression analysis of circRNAs
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 downregulated 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

Weighted gene co-expression network analysis
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 Frontiers in Veterinary Science 07 frontiersin.orgdistribution (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).
and RT-qPCR results confirms the accuracy of the sequencing (p < 0.05) (Figure 7).

Discussion
Limited research exists on the regulatory mechanisms of growth and development in the local Chengkou mountain chicken breed in China.While our previous study analyzed the role of mRNA in muscle development using transcriptomics, the investigation of circRNAs remains insufficient (19).This study aims to address this gap by focusing on the role and regulation of circRNAs in embryonic muscle development of Chengkou mountain chickens.
The stained muscle sections indicated that the embryonic muscle development pattern of the Chengkou mountain chicken closely resembled that of the Tibetan chicken and other local breeds in China (27,28).During E12, the muscle fibers were intermixed, and their distinct boundaries were indiscernible.However, by E19, the muscle fibers displayed clear boundaries and were easily identifiable.This insight holds relevance for examining of muscle development in poultry breeds in other regions.Similar to other animal studies, this study identified over 9,000 novel circRNAs, including 917 differentially expressed circRNAs at different stages (29).Our findings are similar to those of Yuan et al. (30), who suggest that chromosome 1 is crucial in skeletal muscle development, as most new circRNAs are derived from it.The comparison of differentially expressed circRNAs between different periods also indicated that different circRNAs regulate myogenesis in different stages.Notably, some of the circRNAs were derived from genes that are highly associated with muscle development, such as MEF2C.The MEF2C gene affects skeletal muscle growth and development by regulating calcium-mediated carbonic anhydrase III (CAIII) expression (31).It is evident that circRNAs are involved in the process of muscle development.
The ME salmon module, the principal WGCNA module linked to muscle development, exhibited the strongest correlation with this trait.Hence, we chose it for analysis of muscle-related pathways and gene enrichment.CircRNAs were found to be enriched in GO and KEGG pathways associated with muscle development, such as significant enrichment in MAPK and Wnt signaling pathways.The MAPK signaling pathway is linked to cell  (36)(37)(38).In this study, we hypothesized that circRNA is important for controlling the growth and development of embryonic muscles.Indeed, the 754 identified circRNAs target 68 known miRNAs (24).These 68 miRNAs and the 2,180 differentially express mRNAs formed the ceRNA network.All transcripts in the ceRNA regulatory network were differentially expressed in the circRNA-miRNA-mRNA pathway, affirming the reliability of the results.Thus, a more comprehensive understanding of the entire regulatory pathway of muscle development can be attained through multi-omics analysis (39).Furthermore, the miRNAs in the network contain regulatory factors related to muscle development.For instance, miR-9-5p stimulates myogenic differentiation through Dlx3/Myf5, while miR-15c-5p can regulate muscle generation by activating the IGF1-PI3K/AKT signaling pathway (40,41).The results show that the detected circRNAs may influence muscle growth and development during the embryonic stage by controlling these miRNAs.
The miRDB prediction criteria (scores >80 and miRNA length <2,000) for targeted gene prediction revealed four candidate circRNAregulated genes related to muscle development, including MYH1B, MYH1D, CAPZA1, and PERM1 (42).MYH1B and MYH1D are members of the Myosin Heavy Chain (MyHC) gene family, which plays an important role in skeletal muscle growth and development (43,44).Yu et al. (45) demonstrated the regulatory influence of LncRNA-FKBP1C, a long non-coding RNA, on MYH1B, thereby facilitating myoblast differentiation in poultry growth and development.Additionally, during embryonic development, MYH1D is significantly enriched in the broiler genome and MYH1F is essential for leg muscle growth in Chengkou mountain chicken (19,46).In this study, the MyHC family members were up-regulated over time, indicating that this family is important for muscle fiber growth.
Huang et al. (47) discovered that CAPZA1 impedes hepatocellular carcinoma (HCC) cell metastasis by managing actin cytoskeleton remodeling in HCC cells, mainly through epithelial-mesenchymal transition.CAPZA1 encodes the α1 subunit of the actin-conjugated CapZ, and its function is related to the assembly of actin filaments (48).In this experiment, the ceRNA network revealed that miR-9-5p suppresses 104 circRNAs, up-regulating CAPZA1.This pattern could remodel the actin skeleton, enabling the normal growth and development of muscle (49).The ceRNA analysis demonstrated that decreasing circRNAs up-regulate miR-338-3p and down-regulate PERM1.PERM1 is a muscle-specific regulator that regulates mitochondrial biogenesis and oxidation and enhances the spare respiratory capacity in muscles by enhancing mitochondrial function and vascular formation in skeletal muscles (50).However, whether CAPZA1 influences skeletal muscle growth through EMT remains unknown.In addition, PERM1 does not increase the total fiber mass, and the underlying reasons for its regulatory behavior remain unknown.Therefore, the analysis of the ceRNA network clarified that circRNAs are able to influence the activity of genes important for skeletal muscle development by regulating miRNAs, which may have an impact on normal skeletal muscle development.

Conclusion
This study generated 9,904 circRNAs, including 7,102 "exon" type circRNAs, from 12 cDNA libraries.CircRNA expression is time-specific, and key circRNA-derived genes significantly enriched GO terms and KEGG pathways for cell composition, cell development regulation, and muscle system regulation.This study used sequencing, screening, and circRNA identification in the skeletal muscles of developing Chengkou mountain chicken embryos to identify the macro effects and action pathways of circRNAs on muscle development.Therefore, the results help explore the specific biological functions and mechanisms of key circRNAs in subsequent studies of chicken embryonic muscle development.The numerous differentially expressed circRNAs provide data for subsequent improvement of Chengkou mountain chicken breeding.The discovery of signaling pathways and genes also provides theoretical support for trait improvement of domestic chickens.

FIGURE 2
FIGURE 2 Identification of circRNAs.(A) Source genes of the identified circRNAs.(B) Type identification of all novel circRNAs.(C) Length distribution of all novel circRNAs.(D) The distribution of identified circRNAs in different chromosomes.

FIGURE 3
FIGURE 3 Differential expression analysis of circRNAs.(A) Principal component analysis of samples.(B) Distribution of circRNA number in different periods.(C) Upset plot of circRNA-derived genes in different periods.(D) GO analysis of circRNA-derived genes at the same intersection in four periods.(E) KEGG analysis of circRNA-derived genes at the same intersection in four periods.

FIGURE 4
FIGURE 4 Analysis of differential circRNAs in four periods.(A) The number of differentially expressed circRNAs was compared between the four periods.(B) Upand down-regulated circRNA analysis.(C) Distribution trend of differential circRNA expression, color means significant difference (p < 0.05), gray means not significant (p > 0.05).(D) Time line of differential circRNAs.

FIGURE 5
FIGURE 5 Weighted gene co-expression network analysis of circRNAs.(A) The power value curve.(B) Module eigenvalue clustering.(C) Module gene correlation analysis.(D) Correlation analysis of traits.(E) Top 20 enriched entries in GO pathway analysis.(F) Top 20 enriched entries in KEGG pathway analysis.

FIGURE 6
FIGURE 6 Predicted interactions between circRNA, miRNA and mRNA.(A) Networks associated with MYH1B and MYH1D.(B) Networks associated with PERM1.(C) Networks associated with CAPZA1.Orange triangle represents circRNA and blue circle represents gene.The green diamond shape represents the miRNA.

TABLE 2 Table
of comparison between anchor reads and the chicken reference genome.

TABLE 1
Statistical table of the alignment of total reads to the reference genome.