Transcriptome RNA Sequencing Reveals That Circular RNAs Are Abundantly Expressed in Embryonic Breast Muscle of Duck

Simple Summary Research into the genetic mechanisms of skeletal muscle development in ducks is of great significance to the breeding and industrialization of indigenous ducks in China. This study examined the crucial candidates during two developmental stages of duck embryonic breast muscles by circular RNA sequencing. There were 16,622 circular RNAs detected between embryonic day 13 (E13) and embryonic day 19 (E19), of which 260 showed differential expression. Additionally, we examined the protein–protein interaction (PPI) with the parental genes of the differentially expressed circular RNAs, as well as the enriched GO terms and KEGG pathways. In addition, we selected one circular transcript of the growth arrest-specific gene 2 (GAS2), or circGAS2-2, to perform functional validations in duck embryonic myoblasts. The results showed that circGAS2-2 accelerated cell cycling and promoted the proliferation of myoblasts. Our results supplement the duck circular RNA database and demonstrate that circular RNAs are abundantly and differentially expressed during duck skeletal muscle development. Abstract Circular RNAs are widespread in various species and have important roles in myogenesis. However, the circular RNAs involved in breast muscle development in ducks have not yet been studied. Here, to identify circular RNAs during duck skeletal muscle development, three pectorales from Shan Ma ducks at E13 and E19, which represent undifferentiated and differentiated myoblasts, respectively, were collected and subjected to RNA sequencing. A total of 16,622 circular RNAs were identified, of which approximately 80% were exonic circular RNAs and 260 were markedly differentially expressed between E19 and E13. The parental genes of the differentially expressed circular RNAs were significantly enriched in muscle-related biological processes. Moreover, we found that the overexpression of circGAS2-2 promoted cell cycle progression and increased the proliferation viability of duck primary myoblasts; conversely, knockdown of circGAS2-2 retarded the cell cycle and reduced the proliferation viability of myoblasts. Taken together, our results demonstrate that circular RNAs are widespread and variously expressed during the development of duck skeletal muscle and that circGAS2-2 is involved in the regulation of myogenesis.


Introduction
China is the biggest duck meat producer in the world. The consumption of poultry meat in China has significantly increased as a result of African swine fever (http://data.stats.gov.cn/, accessed on 1 September 2022). The development of poultry skeletal muscle is influenced by many factors, including genetics, diet, and the environment [1,2]. Myogenesis can be divided into the following several stages: somite differentiation into myoblasts; myoblasts' migration and expression of particular myogenic transcription factors (MTFs); myoblasts' differentiation into myotubes under the influence enrich the circRNAs. As a result of fragmentation and reverse transcription of the enriched circRNAs, cDNA was synthesized, and U-labeled second-stranded DNAs were synthesized with DNA polymerase I, RNase H, and dUTP. To prepare the blunt ends of each strand for ligation to the indexed adapters, an A-base was added. The fragments were ligated to single-or dual-index adapters, and AMPureXP beads were used for the size selection. A PCR amplification of the UDG-treated second-stranded DNAs was carried out after the products had been heat-labilely treated with the enzyme. In the final cDNA library, the average insert size was 300 ± 50 bp. In the final step, paired-end sequencing was performed on an Illumina HiSeq 4000 platform (LC Bio, Hangzhou, Zhejiang, China) using the vendor's recommended protocols.

Annotation of Duck circRNAs and Differential Expression Analysis between E13 and E19
The reads with adaptors, low-quality bases, or undetermined bases were removed with Cutadapt [33]. The quality of the sequence was verified using FastQC (http://www. bioinformatics.babraham.ac.uk/projects/fastqc/, accessed on 1 February 2020). The sequencing reads were mapped to the genome of Anas platyrhynchos (http://asia.ensembl. org/Anas_platyrhynchos/Info/Index, accessed on 3 February 2020) using Bowtie 2 and TopHat 2 [34,35]. The mapping of the remaining reads was performed using TopHat-Fusion [36]. The mapped reads were assembled into circular RNAs using CIRCExplorer2 and CIRI. TopHat-Fusion was then used to identify the back-splicing reads among the unmapped reads. All of the samples produced distinct circular RNAs [37][38][39].
The relative expression levels of the circRNAs were represented by spliced reads per billion mapping (SRPBM). The circular RNA expression levels from the different groups (E13 and E19) were calculated using in-house scripts. The R package was used for the comparisons, with a |log 2 fold change| of ≥ 1 and a p-value of ≤ 0.05 regarded as indicating differential expression [40].

Enrichment Analysis of the Differentially Expressed circRNAs
Using the DAVID functional annotation tool [41], we performed a pathway enrichment for the parental genes of the differentially expressed circRNAs using Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG). The parental genes were also subjected to gene network analyses using STRING v11.0 and Cytoscape v3.7.2 [42].

cDNA Synthesis and qRT-PCR
The Monad MonScript TM All-in-One Kit (Biopro, Shanghai, China) with DNase was used to reverse-transcribe the total RNA obtained from the animal tissues or cell lines. The synthesized cDNA samples were then diluted in nuclease-free water at a ratio of 1:4 for qPCR. The relative mRNA or circular RNA expression were determined by qRT-PCR in a final volume of 20 µL with 2 × T5 Fast qPCR Mix from TsingKe Bio (Beijing, China). The internal control was the duck GAPDH gene. In order to identify the circRNAs, qRT-PCR was performed on the RNA samples with and without digestion by RNase R. The qRT-PCR program was run by an ABI Q5 system (Thermo Fisher, Waltham, MA, USA) using the following procedure: 95 • C for 3 min; 40 cycles of 95 • C for 10 s; annealing temperature for 1 min; collection of fluorescence from 65 to 95 • C. Each sample was examined with three technical repeats, and the expression levels were calculated using the comparative 2 −∆∆Ct method and presented as means ± S.E.M. A Student's t-test and an one-way ANOVA test were used to compare the expression levels among the different groups. Table S1 lists all of the primers utilized in the current study.

Validation of circRNAs
The circRNAs were validated by RT-PCR with divergent and convergent primers using cDNA or gDNA as the templates. The products were gel purified and cloned with a pClone007 vector kit (TsingKe). The DNASTAR (SeqMan) software (http://www.dnastar.com/, accessed Vet. Sci. 2023, 10, 75 4 of 15 on 3 January 2021) was used to align and analyze the results of positive clones sent to TsingKe for Sanger sequencing.

Overexpression Plasmid Construction and siRNA Synthesis
For the construction of the circGAS2-2 overexpression vector, exons 5, 6, and 7 of the duck GAS2 gene were amplified using cDNA from the pectorales and inserted into a pK25ssAAV-ciR vector (Geneseed, Guangzhou, China) between the EcoRI and BamHI restriction sites. Geneseed was used to design and synthesize the circGAS2-2 siRNA based on the sequence listed in Table S1.

Cell Isolation, Culture, and Electroporation
Duck primary myoblasts were extracted from the pectorales of E13 Shan Ma duck embryos (Tianyun). The breasts of E13 ducks were collected and washed in pre-cooled phosphate-buffered saline with 0.5% penicillin/streptomycin (Invitrogen). The muscle samples were cut into pieces using scissors and then trypsinized (Gibco, Grand Island, NY, USA) at 37 • C for 40 min after the skin and bones had been removed. The digestion was terminated with fetal bovine serum (FBS, Gibco, Grand Island, NY, USA). The mixture was filtered using a 40-µm cell strainer (Sangon, Shanghai, China) and centrifuged at 1100 rpm for 5 min. After serial plating, the cells were grown in complete Dulbecco's modified Eagle medium (Gibco, Grand Island, NY, USA) with 15% FBS and 0.2% penicillin/streptomycin and incubated at 37 • C in a 5% CO 2 humidified atmosphere. The myoblasts were induced to differentiate using 0.5% FBS. A total of 5 µg of DNA plasmid and 5 nmol of siRNA or siRNA NC were mixed with 5 × 10 6 myoblasts in 200 µL of Entranster-E electroporation solution (Engreen, Beijing, China), and cell electroporation was performed on a BTX ECM2001 (BTX, San Diego, CA, USA).

Flow Cytometry of Cell Cycle Analysis
A total of 5 × 10 6 myoblasts were collected and electroporated with an overexpression plasmid or vector control and siRNA or siRNA NC. Then, the cells were seeded in a six-well plate. After 48 h, the myoblasts were collected in 70% ethanol and stored overnight at −20 • C, followed by incubation with 50 mg/mL of propidium iodide (Sigma, St Louis, MO, USA) containing 10 mg/mL of RNase A (Takara, Otsu, Shiga, Japan) and 0.2% (v/v) of Triton X-100 (Sigma) at 4 • C for 30 min. Based on the DNA content in the different phases (G1, S, and G2), the results of the BD AccuriC6 flow cytometer (BD, San Jose, CA, USA) were analyzed using Flow-Jo7.6 (https://www.flowjo.com/, accessed on 10 July 2020). The results were presented as means ± S.E.M. A Student's t-test was used to compare the cell numbers among the different phases.

Cell-Counting Assay
After the cell electroporation, the myoblasts were seeded in a 96-well plate. Cell proliferation was detected at 12 h, 24 h, 36 h, 48 h, and 60 h using a TransDetect CCK kit (TransGen, Beijing, China), following the manufacturer's instructions. Briefly, 10 µL of CCK solution was added to each well and incubated for 1.5 h at 37 • C with 5% CO 2 . Then, an Infinite 200 PRO plate reader (Tecan, Männedorf, Switzerland) was used to measure the absorbance at a wavelength of 450 nm.

Overview of RNA Sequencing Data and Identification of Duck Circular RNA
To identify the circular RNAs in the duck muscle samples, we performed genome-wide circular RNA sequencing with RNase R digestion, followed by functional analysis in duck primary myoblasts. HPF/HPR primers were used to amplify the CHD1 and distinguish the sex of the embryos ( Figure S1). Breast muscle tissues from three female Shan Ma duck embryos at E13 and E19 were subjected to circular RNA sequencing after ribosomal RNA (rRNA) depletion and RNase R digestion.  (Table S2). After eliminating the low-quality reads and adaptors, 159,352,398 clean reads were mapped to the duck genome (Table S2). The unique mapped reads, multi-mapped reads, pair-end mapped reads, reads mapped to sense strands, and reads mapped to antisense strands are shown in Table S2. Then, the circular RNAs were identified using the following criteria: a mismatch of ≤ 2, back-spliced junction reads of ≥ 1, and a distance between two splice sites of < 100 kb. Using these criteria, a total of 19,886 circular RNAs, generated by 10,882 genes, were found among 4,284,753 candidate back-spliced junction reads (Table S3). Of these circular RNAs, 13,034 were circRNAs, 1583 were ciRNAs, and 2005 originated from intergenic sequences ( Figure 1B). They were mostly found on chromosomes 1, 2, 3, 4, 5, and Z ( Figure 1A). The predominant type of circular RNA in the duck skeletal muscle was circRNA (approximately 80%), while the ciRNAs and intergenic circular RNAs were the minor types ( Figure 1C). A total of 1175 circular RNAs were detected in all of the six samples; 1663 circular RNAs were expressed in three of the E13 samples, and 1890 circular RNAs were found in all of the E19 samples ( Figures 1D and S2). Overall, the RNA sequencing data showed that circular RNAs are copiously expressed in duck embryonic skeletal muscle.
To identify the circular RNAs in the duck muscle samples, we performed genomewide circular RNA sequencing with RNase R digestion, followed by functional analysis in duck primary myoblasts. HPF/HPR primers were used to amplify the CHD1 and distinguish the sex of the embryos ( Figure S1). Breast muscle tissues from three female Shan Ma duck embryos at E13 and E19 were subjected to circular RNA sequencing after ribosomal RNA (rRNA) depletion and RNase R digestion.
A total of 290,131,458 reads were sequenced from all of the six samples (Table S2). After eliminating the low-quality reads and adaptors, 159,352,398 clean reads were mapped to the duck genome (Table S2). The unique mapped reads, multi-mapped reads, pair-end mapped reads, reads mapped to sense strands, and reads mapped to antisense strands are shown in Table S2. Then, the circular RNAs were identified using the following criteria: a mismatch of ≤ 2, back-spliced junction reads of ≥ 1, and a distance between two splice sites of < 100 kb. Using these criteria, a total of 19,886 circular RNAs, generated by 10,882 genes, were found among 4,284,753 candidate back-spliced junction reads (Table S3). Of these circular RNAs, 13,034 were circRNAs, 1583 were ciRNAs, and 2005 originated from intergenic sequences ( Figure 1B). They were mostly found on chromosomes 1, 2, 3, 4, 5, and Z ( Figure 1A). The predominant type of circular RNA in the duck skeletal muscle was circRNA (approximately 80%), while the ciRNAs and intergenic circular RNAs were the minor types ( Figure 1C). A total of 1175 circular RNAs were detected in all of the six samples; 1663 circular RNAs were expressed in three of the E13 samples, and 1890 circular RNAs were found in all of the E19 samples ( Figures 1D and S2). Overall, the RNA sequencing data showed that circular RNAs are copiously expressed in duck embryonic skeletal muscle.

Identification of Differentially Expressed Circular RNAs
To explore the potential functions of the circular RNAs expressed in the duck skeletal muscle, we performed a differential expression analysis using the following criteria: a |log 2 fold change| of ≥ 1 and a p-value of ≤ 0.05 (Figure 2A). A total of 260 circular RNAs, of which 155 were upregulated and 105 were downregulated, were differentially expressed between the E19 group and the E13 group ( Figure 2B). The top 100 differentially expressed circular RNAs were then clustered. The data were normalized by SRPBM. The results showed that the upregulated and downregulated circular RNAs were quite distinct between the E13 and E19 groups ( Figure 2C). A total of 170 circular RNAs were expressed in both of the groups, of which 57 were specifically expressed in the E19 group Vet. Sci. 2023, 10, 75 6 of 15 and 33 were uniquely expressed in the E13 group ( Figure 2D). The results demonstrate that duck embryonic breast muscle development is characterized by the differential expression of many circular RNAs.
pressed between the E19 group and the E13 group ( Figure 2B). The top 100 differentially expressed circular RNAs were then clustered. The data were normalized by SRPBM. The results showed that the upregulated and downregulated circular RNAs were quite distinct between the E13 and E19 groups ( Figure 2C). A total of 170 circular RNAs were expressed in both of the groups, of which 57 were specifically expressed in the E19 group and 33 were uniquely expressed in the E13 group ( Figure 2D). The results demonstrate that duck embryonic breast muscle development is characterized by the differential expression of many circular RNAs.

Enrichment Analysis of the Differentially Expressed Circular RNAs
Both circular RNAs and mRNAs originate from pre-mRNA splicing, and circular RNAs can regulate the expression of other transcripts of their parental genes. To predict the function of circular RNAs in embryonic skeletal muscle development, a Gene Ontology (GO) enrichment was performed for the parental genes of the differentially expressed circular RNAs in terms of cellular components, biological processes, and molecular functions. The parental genes were associated with 182 GO terms. Protein binding, ATP binding, protein kinase activity, nucleus, membrane, integral component of membrane, protein phosphorylation, intracellular signal transduction, and G protein-coupled receptor signaling were the most highly enriched GO terms ( Figure 3A). The parental genes were also enriched in muscle-related GO terms, including regulation of skeletal muscle cell differentiation (p < 0.01) and cardiac muscle tissue regeneration (p < 0.01) (Table S4). In addition, a Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis showed that the parental genes of the differentially expressed circular RNAs were significantly enriched in 16 signaling pathways, including gap junctions, the calcium signaling pathway, purine metabolism, progesterone-mediated oocyte maturation, adrenergic signaling in cardiomyocytes, vascular smooth muscle contraction, the GnRH signaling pathway, and cardiac muscle contractions ( Figure 3B, Table S5). Furthermore, we performed a protein-protein interaction (PPI) analysis of the parental genes using STRING. A total of 36 proteins encoded by the parental genes of the differentially expressed circular RNAs interacted with other parental gene products, including some muscle-related proteins, such as the fibroblast growth factor receptor 2 (FGFR2), members of the adenylate cyclase and EPH families, and erythropoietin-producing hepatoma (RBFOX1) ( Figure 3C). These results indicate that differentially expressed circular RNAs are associated with muscle development.

Verification of the Differential Expression Circular RNAs
To validate the differential expression of the circular RNAs identified in this study, seven circRNAs were randomly selected for qRT-PCR detection. To increase the accuracy of the qRT-PCR, primers were designed with the forward primer crossing the back-splicing junction site. The expression patterns of the selected circRNAs according to the qPCR verification were consistent with the RNA sequencing results, except in the case of circSETD3-1 ( Figure 4). CircGAS2-2, circFGFR2-2, circPEPD-1, and circGLI3-1 were significantly downregulated, whereas circSTK39 and circMAPKBP1-1 were significantly upregulated. These results indicate that the sequencing analysis was reliable.

Experimental Validation and Spatiotemporal Expression of circGAS2-2
Owing to the abundant expression of downregulated circular RNAs, circGAS2-2 was selected as the candidate circRNA for the subsequent functional validation experiments in duck primary breast myoblasts. The genomic structure of the duck GAS2 gene is shown in Figure 5A. The GAS2 gene is located on duck chromosome 5 and contains eight exons. Interestingly, the seven circular RNAs derived from the GAS2 gene based on our sequencing data were all exonic circular RNAs. Of these circular RNAs, circGAS2-2 had the highest expression level. A divergent and convergent primer set based on the duck GAS2 genomic sequence was designed to confirm the sequence and junction site of circGAS2-2. The PCR procedure with the divergent and convergent primers was run using the cDNA and genomic DNA as templates. The electrophoresis results revealed that the convergent primer PCR products were amplified by both templates as expected, whereas no divergent primer PCR product was found with the genomic DNA template ( Figure 5B). The RNase R treatment also confirmed the existence of circGAS2-2. The qRT-PCR data results demonstrated that RNase R had little effect on circGAS2-2; however, the RNA levels of GAPDH, which was used as a linear RNA control, were significantly decreased ( Figure 5C). Furthermore, an analysis of the PCR products with the divergent primers by Sanger sequencing showed that the 5 site of exon five and the 3 site of exon seven were linked, indicating that circGAS2-2 is generated from the fifth, sixth, and seventh exons of GAS2 ( Figure 5D). The expression levels of circGAS2-2 on E10, E13, E16, E19, E22, and 1 day post-hatching (P1) were detected by qRT-PCR. The expression trend showed a gradual decline from E13 to P1 ( Figure 5E,F). These findings demonstrate that circGAS2-2 exists and is expressed differently throughout embryonic skeletal muscle development, suggesting that circGAS2-2 may have key roles in duck muscle development.

Experimental Validation and Spatiotemporal Expression of circGAS2-2
Owing to the abundant expression of downregulated circular RNAs, circGAS2-2 was selected as the candidate circRNA for the subsequent functional validation experiments in duck primary breast myoblasts. The genomic structure of the duck GAS2 gene is shown in Figure 5A. The GAS2 gene is located on duck chromosome 5 and contains eight exons. Interestingly, the seven circular RNAs derived from the GAS2 gene based on our sequencing data were all exonic circular RNAs. Of these circular RNAs, circGAS2-2 had the highest expression level. A divergent and convergent primer set based on the duck GAS2 genomic sequence was designed to confirm the sequence and junction site of circGAS2-2. The PCR procedure with the divergent and convergent primers was run using the cDNA and genomic DNA as templates. The electrophoresis results revealed that the convergent primer PCR products were amplified by both templates as expected, whereas no divergent primer PCR product was found with the genomic DNA template ( Figure 5B). The RNase R treatment also confirmed the existence of circGAS2-2. The qRT-PCR data results demonstrated that RNase R had little effect on circGAS2-2; however, the RNA levels of GAPDH, which was used as a linear RNA control, were significantly decreased ( Figure  5C). Furthermore, an analysis of the PCR products with the divergent primers by Sanger sequencing showed that the 5′ site of exon five and the 3′ site of exon seven were linked, indicating that circGAS2-2 is generated from the fifth, sixth, and seventh exons of GAS2 ( Figure 5D). The expression levels of circGAS2-2 on E10, E13, E16, E19, E22, and 1 day

circGAS2-2 Promotes the Proliferation of Duck Breast Primary Myoblasts
Given that circGAS2-2 was differentially expressed at several phases of the duck breast muscle development, we hypothesized that it might regulate the proliferation and differentiation of breast primary myoblasts. To confirm our hypothesis, duck primary myoblasts were isolated, cultured, and induced to differentiate. To investigate the function of circGAS2-2, an overexpression vector and short interfering RNAs (siRNAs) were constructed and electroporated into myoblasts. The group with the overexpression vector expressed circGAS2-2 with high efficiency compared with the group electroporated with pK25ssAAV-ciR, and the siRNA significantly knocked down the level of circGAS2-2 compared with the siRNA negative control (NC) group ( Figure 6A). The cell cycle in the primary myoblasts from the ducks was analyzed using flow cytometry after the electroporation. The results revealed that the overexpression of circGAS2-2 significantly enhanced cell cycling by decreasing the number of cells in the G1 phase and increasing the number in the G2 phase ( Figure 6B). By contrast, knockdown of circGAS2-2 markedly reduced the cell population in the G2 stage and retarded the cell cycle of the primary myoblasts ( Figure 6C). The duck primary myoblasts were tested for proliferation vitality using the cell-counting assay. The groups electroporated with the overexpression plasmid had markedly higher proliferation vitality than the pK25ssAAV-ciR group at 36 h, 48 h, and 60 h after the electroporation ( Figure 6D). Conversely, the electroporation with si-circGAS2-2 significantly suppressed the cell growth rate compared with the siRNA NC ( Figure 6E). Taken together, these results suggest that circGAS2-2 promotes the proliferation of duck breast primary myoblasts. We also investigated the function of circGAS2-2 in myoblast differentiation. After electroporation with pK25ssAAV-ciR-circGAS2-2/pK25ssAAV-ciR and si-circGAS2-2/siRNA NC, the expression levels of the myoblast differentiation marker genes, myogenic factor 5 (MYF5), MRF4, myogenin (MYOG), and myogenic differentiation 1 (MYOD1), were determined. The qRT-PCR results showed that neither ectopic expression nor knockdown of circGAS2-2 significantly changed the expression levels of MYOD1, MRF4, MYF5, or MYOG ( Figure S3), indicating that circGAS2-2 had no regulatory effect on myoblast differentiation.  2/pK25ssAAV-ciR and si-circGAS2-2/siRNA NC, the expression levels of the myoblast differentiation marker genes, myogenic factor 5 (MYF5), MRF4, myogenin (MYOG), and myogenic differentiation 1 (MYOD1), were determined. The qRT-PCR results showed that neither ectopic expression nor knockdown of circGAS2-2 significantly changed the expression levels of MYOD1, MRF4, MYF5, or MYOG ( Figure S3), indicating that circGAS2-2 had no regulatory effect on myoblast differentiation.

Discussion
Circular RNAs have been found to be abundant and dynamically expressed in different tissues of various species. Non-canonical RNA splicing was detected in large transcriptome sequencing data, indicating that circular RNAs are universal in the transcriptome [43]. When circular RNAs in mice were identified, their biological function became a focus of speculation [44]. Circular RNA sequencing databases for various domestic animals, including pigs, cattle, and sheep, have now been constructed [27,45,46]. However, few studies have considered circular RNAs in poultry [28,47,48]. In this study, we investigated the expression of circular RNAs in two stages of duck breast muscle development. A total of 16,622 circular RNAs were identified in the E13 and E19 groups, which was greater than the number detected in duck follicles (4204) and chicken embryonic leg muscles (13,377) [28,48]. These results demonstrate that circular RNAs are universally abundant and dynamically expressed in different tissues. The circular RNAs identified in this study were mainly derived from exons, consistent with studies in chicken and pig muscles but not with those in duck follicles, indicating that back-splicing of pre-mRNA is tissue-specific [27,28,48].
We identified 260 circular RNAs that were differentially expressed between the E13 and E19 groups. These differentially expressed circular RNAs were significantly enriched in GO terms associated with skeletal muscle cell differentiation and muscle tissue regeneration. A KEGG analysis showed that the circular RNAs were significantly enriched in the following: gap junctions, the calcium signaling pathway, purine metabolism, vascular smooth muscle contraction, the GnRH signaling pathway, the MAPK signaling pathway, and cardiac muscle contraction. Gap junction proteins are associated with cell communication in cardiac muscle cells [49]. Calcium homeostasis has important roles in controlling the fate of muscle tissue through multiple pathways, and calcium signaling has a vital regulatory function in duck breast muscle development [50,51]. Moreover, the MAPK and GnRH signaling pathways are associated with animal body weight and contribute to muscle growth [52,53]. According to the PPI analysis, many proteins encoded by the parental genes of the differentially expressed circular RNAs had interactions with other parental proteins. FGFR2 is an important activator of satellite cells in skeletal muscle growth, and circFGFR2 promotes the proliferation and differentiation of primary myoblasts [2,54]. Taken together, these results indicate that differentially expressed circular RNAs may participate in skeletal muscle growth and development via regulation of their parental genes. Circular RNAs, which play a major role in myogenesis by acting as sponges for miR-NAs, have not yet been fully understood [29,55]. As a scaffold for miR-378a-3p, circLMO7 increases histone deacetylase 4 (HDAC4) expression, decreases myocyte enhancer factor 2A (MEF2A) expression, and promotes the differentiation of myoblasts [56]. Circular RNAs also regulate myogenesis by producing functional peptides [14]. Genes containing multiple exons usually generate different circular transcripts by alternative back-splicing [37]. In this study, many duck genes were found to transcribe several circular RNAs, most of which were circRNAs. A total of seven circRNAs were produced by the GAS2 gene. We selected a circRNA derived from the duck GAS2 gene for validation experiments in myoblasts, owing to its high expression levels among the downregulated circular RNAs. Before the myoblast experiments, we used divergent-convergent primer PCR, Sanger sequencing, and RNase R treatments to validate the presence of circGAS2-2. The results demonstrated that circGAS2-2 is abundantly expressed in duck breast muscle.
The growth arrest-specific (GAS) gene family has been linked to reversible growth arrest. The actin cytoskeleton and cell shape are swiftly reset by GAS2-regulated microfilament changes in response to growth arrest triggered by environmental stimuli, such as apoptosis, distinct proliferative stimulation, and varied growth factors [57]. The GAS2 protein is an apoptotic regulator that interacts with caspase-3 and caspase-7 and promotes the proliferation of T-cell acute lymphoblastic leukemia by activating the Wnt/β-catenin pathway [58][59][60]. The function of GAS2 in muscle development has not yet been reported. In this study, we found that circGAS2-2 could promote the proliferation of duck breast primary myoblasts by facilitating the progression of the cell cycle. However, neither overexpression nor knockdown of circGAS2-2 had any significant effect on the expression levels of the myoblast differentiation-determining factors (MYOD1, MYOG, MYF5, and MRF4), indicating that circGAS2-2 does not regulate duck myoblast differentiation. Thus, further investigation is needed to determine the underlying mechanisms by which circGAS2-2 promotes myoblast proliferation.

Conclusions
In conclusion, our results showed that circular RNAs were abundantly and differentially expressed during the process of duck skeletal muscle development. The circular RNA circGAS2-2, which is derived from GAS2, could promote myoblast proliferation in ducks.