Integrative analysis of circRNA, miRNA, and mRNA profiles to reveal ceRNA regulation in chicken muscle development from the embryonic to post-hatching periods

The growth and development of skeletal muscle are regulated by protein-coding genes and non-coding RNA. Circular RNA (circRNA) is a type of non-coding RNA involved in a variety of biological processes, especially in post-transcriptional regulation. To better understand the regulatory mechanism of circRNAs during the development of muscle in chicken, we performed RNA-seq with linear RNA depletion for chicken breast muscle in 12 (E 12) and17 (E 17) day embryos, and 1 (D 1), 14 (D 14), 56 (D 56), and 98 (D 98) days post-hatch. We identified 5755 differentially expressed (DE)-circRNAs during muscle development. We profiled the expression of DE-circRNAs and mRNAs (identified in our previous study) at up to six time points during chicken muscle development and uncovered a significant profile (profile 16) for circRNA upregulation during aging in muscle tissues. To investigate competing endogenous RNA (ceRNA) regulation in muscle and identify muscle-related circRNAs, we constructed a circRNA-miRNA-mRNA regulatory network using the circRNAs and mRNAs from profile 16 and miRNAs identified in our previous study, which included 361 miRNAs, 68 circRNAs, 599 mRNAs, and 31,063 interacting pairs. Functional annotation showed that upregulated circRNAs might contribute to glycolysis/gluconeogenesis, biosynthesis of amino acids, pyruvate metabolism, carbon metabolism, glycogen and sucrose metabolism through the ceRNA network, and thus affected postnatal muscle development by regulating muscle protein deposition. Of them, circRNA225 and circRNA226 from the same host gene might be key circRNAs that could regulate muscle development by interacting with seven common miRNAs and 207 mRNAs. Our experiments also demonstrated that there were interactions among circRNA225, gga-miR-1306-5p, and heat shock protein alpha 8 (HSPA8). Our results suggest that adequate supply of nutrients such as energy and protein after hatching may be a key factor in ensuring chicken yield, and provide several candidate circRNAs for future studies concerning ceRNA regulation during chicken muscle development.


Background
The broiler industry is a key player in contributing to sustainable food sources worldwide. Skeletal muscle is the most important component of the bird and directly correlates with meat production and quality in the broiler industry. Therefore, unveiling the molecular mechanisms underneath skeletal muscle formation and development Open Access *Correspondence: jqsyzslj@163.com 1 Poultry Institute, Shandong Academy of Agricultural Sciences, Ji'nan 250023, China Full list of author information is available at the end of the article is of vital interest. Muscle development is a complex process that can be regulated by several genes, transcription factors, and some non-coding RNAs [1][2][3][4][5].
CircRNAs are a novel type of non-coding RNA that comprise a class of RNA covalently closed-loop structures without 5′-3′ polarities. They are created by RNA splicing events that occur at a characteristic "head to tail" splice junction, where an acceptor splice site at the 5′ end of an exon and a donor site at the 3′ end of a downstream exon are joined [6]. CircRNAs can regulate gene expression at a transcriptional, post-transcriptional, or translational level by the use of microRNA sponges [7]. Increasing evidence has demonstrated that circRNAs modulate gene expression during processes of myogenesis in chicken. CircSVIL promotes both myoblast proliferation and differentiation by sponging miR-203 and increasing the expression of its targets, Jun proto-oncogene (c-JUN) and myocyte enhancer factor 2C (MEF2C) [7]. CircFGFR2 promotes proliferation and differentiation in chicken skeletal muscles by acting as a decoy for miR-133a-5p and miR-29b-1-5p [8]. However, the completed catalog of circRNAs involved in the muscle development of chickens remains mostly unknown. Ouyang et al. (2017) identified 462 differentially expressed (DE) circRNAs during embryonic leg muscle development in chicken [9].
Muscle growth and development go through two important periods: hyperplasia and hypertrophy. Hyperplasia refers to increases in the number of cells or muscle fibers that occur mainly in the embryonic period. Whereas, hypertrophy refers to increases in cell size that occur mainly after birth [10]. Our previous study on a protein-coding gene expression profile in the breast muscle of Shouguang chickens of 12 (E 12) and 17 (E 17) day embryos and at 1 (D 1), 14 (D 14), 56 (D 56), and 98 (D 98) days post-hatch also showed that there were distinct molecular processes that occurred in chicken muscle development between embryonic and post-hatching periods [11]. The systematic identification and characterization of circRNAs from the embryonic to post-hatching periods will contribute to a completed catalog of circR-NAs in chicken muscle development and in uncovering the regulatory mechanisms of muscle development.

Animals and sample collection
One thousand Shouguang chicken eggs were obtained from an experimental farm of the Poultry Institute, Shandong Academy of Agricultural Sciences (PI, SAAS, Ji'nan, China). All eggs were incubated by a normal procedure and chicks were reared in cages under continuous lighting using standard conditions of temperature, humidity, and ventilation on the PI, SAAS, farm. The same diet was fed to all chickens and a three-phase feeding system was used: the starter ration (d 1 to d 28) with 21.0% crude protein and 12.12 MJ/kg, the second phase (d 28 to d 56) with 19.0% crude protein and 12.54 MJ/kg, and the last phase (after d 56) with 16.0% crude protein and 12.96 MJ/kg. Feed and water were provided ad libitum during the experiment. Breast muscles were used at E 12, E 17, D 1, D 14, D 56, and D 98. The embryos, bodies and breast muscles on both sides were weighed. The data were analyzed by one-way ANOVA and multiple comparison using SAS 9.4, and the results were shown as means ± SD. Sub-samples were immediately fixed in 4% paraformaldehyde and held at room temperature and additional samples were snap frozen in liquid nitrogen and held at -80 °C until RNA extraction [11].The sex of the chickens was determined by PCR amplification using sex-specific primers. Chickens with two bands of 450 and 600 bp were females, whereas those with one band of 600 bp were males [12]. Eighteen female chickens from the six developmental stages (E12, E17, D1, D14, D56, and D98) were used for sequencing. Three biological replicates for each stage.

Histology
Fixed tissues were dehydrated through an ascending ethanol series, embedded in paraffin, and sectioned (3-5 μm). After dewaxing in xylene and rehydration using a descending alcohol gradient, mounted muscle sections were stained with hematoxylin and eosin (H&E). Images were captured and processed with Image-Pro Plus 6.0 software [13].

CircRNA library construction and Illumina sequencing
Total RNA was isolated and purified using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's procedure. The amount and purity of the RNA for each sample were quantified using a NanoDrop ND-1000 (NanoDrop, Wilmington, DE, USA). The RNA integrity was assessed by an Agilent 2100 with RIN number > 7.0 [11]. Approximately 5 μg of total RNA was used to deplete ribosomal RNA according to the protocol of a Ribo-Zero rRNA Removal Kit (Illumina, San Diego, USA). Then, RNAs were treated with RNase R (Epicenter Inc, Madison, WI, USA) to remove linear RNAs and to enrich circRNAs. After removing ribosomal and linear RNAs, the enriched circRNAs were fragmented into small pieces using divalent cations under a high temperature. Then the cleaved RNA fragments were reversetranscribed to create the cDNA, which was then used to synthesize U-labeled second-stranded DNAs with E. coli DNA polymerase I, RNase H, and dUTP. An A-base was then added to the blunt ends of each strand to prepare for the ligation to indexed adapters. Each adapter contained a T-base overhang for ligating the adapter to the A-tailed fragmented DNA. Single-or dual-index adapters were ligated to the fragments, and size selection was performed with AMPureXP beads. After the heatlabile UDG enzyme treatment of the U-labeled secondstranded DNAs, the ligated products were amplified by PCR using the following conditions: initial denaturation at 95 °C for 3 min; eight cycles of denaturation at 98 °C for 15 s, annealing at 60 °C for 15 s, and extension at 72 °C for 30 s; and then a final extension at 72 °C for 5 min. The average insert size for the final cDNA library was 300 ± 50 bp. Finally, we performed paired-end sequencing on an Illumina Hiseq 4000 (LC Bio, Hangzhou, China) following the vendor's recommended protocol.

Time series expression profile clustering
Co-expression circRNAs and mRNAs were clustered by Short Time-Series Expression Miner (STEM, version 1.3.11) [21]. Expression profiles of circRNAs and mRNAs were clustered based on their log 2 (normalized data) and their correlation coefficients. The maximum unit change in model profiles between time points was adjusted to two and the maximum number of model profiles to 20. The statistical significance of the number of circRNAs and mRNAs to each profile versus the expected number was computed using the algorithm proposed by Ernst and Bar-Joseph [21].

Construction of the CircRNA-miRNA-mRNA Network
MicroRNA target sites in exons of circRNA loci and target mRNAs of miRNAs were identified using TargetScan and miRanda (http:// www. micro rna. org/ micro rna/ home. do) with a score of 50 or higher and free energy -10 or lower [22]. The circRNA-miRNA-mRNA network was constructed according to the prediction of miRNA binding sites. Cytoscape software was used to construct circRNA-miRNA-mRNA networks [23]. GO annotation and KEGG pathway [24][25][26] analysis were implemented in Omicstudio (https:// www. omics tudio. cn/ index) for mRNAs regulated by circRNAs through ceRNA.

Dual-luciferase reporter assay
For luciferase reporter assay, 293T cells were seeded in 48-well plates and co-transfected with wild type (WT) or mutated (MT) reporter vector and miR-1306-5p mimics or NC duplexes. Firefly and Renilla luciferase activities were measured at 48 h post transfection using a Dual-GLO Luciferase Assay System Kit (Promega, Madison, USA), following the manufacturer's instructions. Luminescence was measured using a Fluorescence/Multi-Detection Microplate Reader (BioTek, Vermont, USA) and firefly luciferase activities were normalized to Renilla luminescence in each well.

Quantitative real-time PCR analysis
The circRNAs were validated with convergent and divergent primers according to a previous study [27]. Details of divergent and convergent primers are in Table S1. PCR products of divergent and convergent primers for cDNA and genomic DNA were analyzed by agarose gel electrophoresis. Back-splicing sites of circRNAs were confirmed by Sanger sequencing at Tsingke Biotech Co. Ltd. (Beijing, China). To check the sensitivity of circRNA to RNaseR, qRT-PCR was also performed using RNA samples with and without RNaseR treatment. The expression of PCR products from divergent primers for each sample was validated using qRT-PCR. The qRT-PCR program was implemented using ABI7500 (Life Technologies, Carlsbad, USA) with SYBR Green (TaKaRa, Dalian, China) in a final volume of 20 μL. Each assay was performed in triplicate using the following cycling conditions: 95 °C/30 s and 40 cycles of 95 °C/5 s, and 60 °C/34 s. The △△ Ct method was used to compare gene expression, with β actin as a reference gene. Three independent replications were used for each assay and data were presented as means ± SD.

Characteristics of the body weight and breast muscle dynamic variation
To investigate the muscle development of chicken, we explored the body weight and breast muscle dynamic variation of Shouguang chickens from the embryonic to post-hatching periods (E12, E17, D1, D14, D56, and D98). The result showed that the body weight (BW) and breast muscle weight (BrW) was increased with growth of chicken, and the BW and BrW were increased markedly after hatching, especially from D14 to D98 (Fig. 1A-B). Moreover, the muscle tissues of chickens were performed using hematoxylin-eosin (H-E) staining to visual the muscle fibers from the embryonic to post-hatching periods (Fig. 1C).

Overview of CircRNA deep sequencing data
To identify the circRNA profile during skeletal muscle development in chickens, breast muscle tissues of three female Shouguang chickens at E 12, E 17, D 1, D 14, D 56, and D 98 were used for RNA sequencing after rRNAdepletion and RNase R treatment. For RNA-seq, a total of 1,310,663,702 reads were generated from 18 samples, and each sample yielded more than 60 million reads. The average GC content was 55.72% (Table 1). Raw data were processed to remove adapter and low quality sequences and then mapped to the chicken reference genome (Gallus-gallus-5.0/galGal5). From the 18 muscle tissues, a total of 30,840 circRNAs were detected. To generate a confident set of circRNAs, we only retained circRNAs expressed in at least three or more samples. These steps resulted in 7,342 circRNAs, and all the circRNAs originated from 2,291 chicken genes (Table S2). These circR-NAs were located on chicken chromosomes 1-28, 30, 32, 33, W, and Z (Fig. 2). We identified 5,755 circRNAs that were DE in pairwise comparisons between the libraries of breast muscle at the six developmental stages (p < 0.05, log 2 (fold change) ≥ 1 or log 2 (fold change) ≤ -1) (Table S3).

STEM analysis of DE-circRNAs and DE-mRNAs
We investigated the DE-mRNAs in breast muscle at E12, E17, D1, D14, D56, and D98 in a previous study [11], and identified 11,975 DE-mRNAs (Table S4). CircRNA can regulate mRNA expression. There may be a regulatory relationship between circRNAs and mRNAs with the same expression pattern. Therefore, we performed trend analysis using STEM to identify DE-circRNAs and DE-mRNAs with the same expression pattern during breast muscle development. We identified three significant profiles (Fig. 3A). Profile 3 with a downregulated pattern contained one circRNA and 5,018 mRNAs ( Fig. 3B and Table S5). Profile 16 with an upregulated pattern contained 69 circRNAs and 1,173 mRNAs ( Fig. 3C and Table  S6). Profile 14, with two circRNAs and 722 mRNAs was the third pattern. There was an increase in regulation from E12 to E17 and a decrease from E17 to D14, and it remained stable from D14 to D56 and finally increased from D56 to D98 (Fig. 3D and Table S7).

Functional analysis of the CircRNAs during muscle development
To investigate the potential functional implication of the circRNAs during breast muscle development, we performed GO and KEGG analysis of the mRNAs in the circRNA-miRNA-mRNA triple network constructed with profile 16 (Table S12). GO functional annotation showed that upregulated circRNAs in profile 16 contributed to the glycolytic process, gluconeogenesis, AMP-activated protein kinase activity, peptidyl-cysteine S-trans-nitrosylation, and oxidoreductase activity (Fig. 4A). The KEGG pathway enrichment displayed that glycolysis/gluconeogenesis, biosynthesis of amino acids, pyruvate metabolism, carbon metabolism, glycogen and sucrose metabolism were the top pathways (Fig. 4B). These results showed that the expression of the circR-NAs that regulated metabolism increased with increased age, and they may play an important role in muscle development.

Validation of CircRNAs
Six circRNAs were randomly selected and the back-site junctions were verified by agarose gel electrophoresis, Sanger sequencing, and quantitative real-time PCR (qRT-PCR) (Fig. 5). We designed divergent and convergent primers for each circRNAs, and used both genomic DNA (gDNA) and complementary DNA (cDNA) as templates for PCR. Divergent primers from each circRNA amplified the expected fragments using cDNA as a template, but they could not amplify PCR products using the gDNA template, suggesting the presence of back-site junctions (Fig. 5A), which were also validated by Sanger sequencing (Fig. 5B). We quantified the six candidate circRNAs with RNaseR treatment to detect the resistance of circRNA to the digestion by RNaseR. The result showed that cir-cRNAs had much more resistant than the linear mRNA (Fig. 5C). The qRT-PCR validation for the five circRNAs was consistent with the trends obtained from circRNA sequencing data (Fig. 5D).

Verification of the interaction among circRNA, miRNA, and mRNA
Our network analysis predicted that there might be ceRNA regulation among circRNA225, gga-miR-1306-5p, and HSPA8 (Fig. 6A). The expression of circRNA225 and HSPA8 were upregulated with age (Table S3, Table S4), while the expression of gga-miR-1306-5p was downregulated with age (Table S8). Therefore, the target relationship of these transcripts was validated using a luciferase reporter gene assay. As demonstrated in Fig. 6B, gga-miR-1306-5p significantly reduced the firefly luciferase activity of the wild type of the circRNA225 reporter compared with negative control, suggesting that circRNA225 directly targets chicken gga-miR-1306-5p. Meanwhile, the target relationship between gga-miR-1306-5p and HSPA8 also showed that gga-miR-1306-5p could combined with the site of wild type reporter, but not the mutant type reporter (Fig. 6C).

Discussion
Studying the mechanism underlying skeletal muscle development would be beneficial to the genetic improvement of the quality and quantity of meat. Therefore, the most critical goal of meat production science is analyzing and understanding the development of muscles. Precisely timed gene expression and alternative splicing patterns play important roles in muscle development. Back-splicing events lead to a large number of specific circRNA species, some of which have important regulatory potential in muscle development [29,30]. In chicks, the circRNA expression profile of embryonic leg muscle (E 11, E 16, and D 1) was uncovered, and 13,377 circRNAs were identified with 462 DE-circRNAs [9]. In the present study, we explored the circRNA expression profile of chicken breast muscle from embryonic (E 12, E 17, and D 1) to post-hatching (D 14, D 56, and D 98) stages and identified 5,755 DE-circRNAs (p-value < 0.05, log 2 (fold change) ≥ 1 or log 2 (fold change) ≤ -1). Our study identified many circRNAs that were highly expressed . Moreover, muscle tissue was also increased rapidly after hatching. Therefore, our results may provide a novel explanation for the phenotype of the differences in muscle development between the embryonic period and the post-hatching period at circRNA expression levels. The main mechanism of circRNA may be in acting as miRNA sponges to modulate post-transcriptional regulation [31]. Our previous studies have explored mRNA and miRNA expression at the same time points as the present research [11,28]. Therefore, we performed circRNA-miRNA-mRNA regulatory network construction using the circRNAs and mRNAs with the same expression pattern and the 587 miRNAs identified in the chicken muscle in our previous studies. As a result, we obtained a circRNA-miRNA-mRNA regulatory network consisting of 68 circRNAs, 599 mRNAs, and 361 miRNAs, and these circRNAs and mRNAs showed a tendency of upregulation with an increase in age. Functional analysis indicated that these up-regulating circRNAs took part in muscle development through regulating genes in muscle metabolism, such as glycolytic processes, gluconeogenesis, AMP-activated protein kinase activity, glycolysis/ gluconeogenesis, biosynthesis of amino acids, and pyruvate metabolism. A previous study also showed that genes associated with the biosynthesis of amino acids, glycolysis/gluconeogenesis, and the tricarboxylic acid (TCA) cycle were upregulated during the goat muscle development [32]. Studies in poultry have shown that muscle mass increased by hypertrophy (increased cellular protein content) after hatching, and this process was controlled by the muscle protein synthesis and degradation [33]. Protein metabolism is often accompanied by energy metabolism, and the glycolysis/gluconeogenesis produces the energy needed for protein turnover during skeletal muscle development [34]. Meanwhile, carbohydrate metabolism can provide a carbon skeleton for nonessential amino acid synthesis. Our ceRNA network analysis also showed that several circRNAs not only play important roles in glycolysis/gluconeogenesis but also contribute to amino acid synthesis through the ceRNA regulation after hatching, which is consistent with the phenotype of rapid increase of muscle tissue and muscle fiber diameter after hatching. Therefore, our study once again demonstrated that the rapid increase in muscle tissue after hatching may be due to the activation of pathways associated with muscle protein deposition, suggesting that ensuring proper energy and protein nutrition after hatching is critical for chicken production (Fig. 7). Interestingly, we observed that two circRNA isoforms (circRNA225 and circRNA226) from the GPD2 gene in this regulatory network. circRNA225 and circRNA226 co-regulated seven miRNAs (gga-miR-106-3p, gga-miR-183, gga-miR-12239-5p, gga-miR-12283-3p, gga-miR-1306-5p, gga-miR-1773-3p, and gga-miR-3594-3p) and 207 mRNAs (e.g., PKM, PGAM1, TPI1, LDHA, UCP3, ATP1A1, HSPA8, and CHAC1) through ceRNA regulation (Fig. S1). The miR-183 family gene cluster plays multiple roles in a wide range of physiological and pathological processes, such as cell proliferation, apoptosis, and metabolism [35]. The overexpression of miR-183-5p inhibits the myogenic differentiation of C2C12 myoblasts [36]. MiR-106b has been found to regulate cellular cholesterol efflux by targeting ABCA1 in macrophages [37], miR-106 plays a role in various dystrophies such as facioscapulohumeral muscular dystrophy, limb-girdle muscular dystrophies, and Miyoshi myopathy [38]. Genes, such as Pyruvate kinase muscle (PKM), Phosphoglycerate mutase 1 (PGAM1), triosephosphate isomerase (TPI1), or lactate dehydrogenase A (LDHA), encoding key enzymes involved in gluconeogenesis/glycolysis and biosynthesis of amino acids KEGG pathway, were regulated by circRNA225 and circRNA226. Pyruvate kinase muscle (PKM) is a key enzyme propelling glycolysis in muscle [39]. Proliferating C2C12 cells exhibit preferential transcription of the PKM2 splice isoform, which is proposed to be essential for the generation of sufficient intermediates in cells for the generation of new macromolecules [40]. Phosphoglycerate mutase 1 (PGAM1) reversibly catalyzes a unique step during glycolysis, controlling the metabolite levels of its substrate in the later stages of the glycolytic pathway [41]. The expression levels of the PGAM1 protein in muscle was positively correlated with chicken aging [42]. Comparison of the longissimus thoracis muscle transcriptome from 15-and 19-month-old Charolais bull calves selected divergently for high or low muscle growth revealed that about two-thirds of the genes (including TPI1 and LDHA) involved in glycolysis were upregulated at 15 and 19 months of age in high muscle growth bull calves. Moreover, other genes that were regulated by circRNA225 and circRNA226 were also involved in muscle metabolism and affected muscle development. Uncoupling proteins (UCPs) can uncouple ATP production from mitochondrial respiration, thereby dissipating Fig. 6 Experimental verification of the circRNA225-gga-miR-1306-5p-HSPA8 regulatory interaction. A Schema of miR-1306-5p binding site in chicken circRNA225 and HSPA8 3'UTR sequence. B Luminescence was measured after co-transfecting the wild-type or mutant sequence of circRNA225 with miR-1306-5p mimics (or mimics negative control (NC)) in 293 T cells. C Luminescence was measured after co-transfecting the wild-type or mutant sequence of HSPA8 3'UTR with miR-1306-5p mimics (or mimics NC) in 293 T cells Lei et al. BMC Genomics (2022) 23:342 Fig. 7 The new circRNA-miRNA-mRNA regulatory network for muscle development of chichken. Blue circles represent circRNAs, red triangles represent miRNAs, and green rectangles represent mRNAs energy as heat and affecting energy metabolism efficiency [43]. UCP3, a fatty acid anion exporter, supports high rates of fatty acid oxidation in muscles [44]. The Na/K-ATPase (NKA) α1 isoform is encoded by the ATPase Na + /K + transporting subunit alpha 1 (ATP1A1) gene. A skeletal muscle-specific ablation of NKA α1 mice had a 35% reduction in skeletal muscle mass and a switch from oxidative to glycolytic fibers [45]. Cation transport regulator-like protein 1 (CHAC1) is a pro-apoptotic protein that has γ-glutamylcyclotransferase activity. The complete knockout of CHAC1 is embryonic lethal. CHAC1heterozygote (het) mice have decreased muscle mass [46]. Heat shock protein alpha 8 (HSPA8) is a molecular chaperone and a member of the heat shock protein family that plays an integral role in amino compound metabolism and lipid homeostasis [47]. The expression levels of HSPA8 in 6 weeks old Pekin Duck were significantly higher than at 2-and 4-weeks-old [48], which was similar with our result. Our experiment results showed that there were target sites of gga-miR-1306-5p in the cir-cRNA225 sequence and HSPA8 mRNA 3´-UTR. Therefore, circRNA225 and circRNA226 may be the potential key factors in the regulation of muscle development by regulation of these miRNAs and mRNAs involved in metabolism.

Conclusion
In summary, our study revealed the expression profiles and potential functions of circRNAs in the breast muscle development of chicken. In total, 5,755 DE-circRNAs were identified during muscle development.
We profiled the expression of DE-circRNAs and DE-mRNAs (identified in our previous study) at up to six time points during chicken muscle development and uncovered a significant profile (profile 16) for circRNA upregulation during aging in muscle tissues. We then constructed a circRNA-miRNA-mRNA regulatory network using the circRNAs and mRNAs in profile 16 and miRNAs identified in our previous study. These circRNAs mainly contribute to metabolism in muscle development by ceRNA regulation. Our study suggested that postnatal nutrient regulation is critical for chicken production and identified many circRNAs that influence muscle development by regulating muscle metabolism.