Skip to main content
Download PDF

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

BMC Genomics volume 23, Article number: 342 (2022) Cite this article

Abstract

Background

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.

Results

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

Conclusions

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.

Peer Review reports

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 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 circRNAs in chicken muscle development and in uncovering the regulatory mechanisms of muscle development.

Methods

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 reverse-transcribed 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 heat-labile UDG enzyme treatment of the U-labeled second-stranded 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.

RNA-seq data analysis

Cutadapt [14] was used to remove the reads that contained adaptor contamination, low quality, and undetermined bases. Then, sequence quality was verified using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). We used Bowtie2 [15] and Tophat2 [16] to map reads to the chicken reference genome (Gallus-gallus-5.0/galGal5). The remaining reads (unmapped reads) were mapped to the genome using TopHat-Fusion [17]. CIRCExplorer [18, 19] was used to de novo assemble the mapped reads to circRNAs. Then, back-splicing reads were identified in unmapped reads by TopHat-Fusion and CIRCExplorer. All samples generated unique circRNAs. The DE circRNAs were selected with log2(fold change) > 1 or log2 (fold change) < -1 and with statistical significance (p-value < 0.05) by R package-edge R [20].

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 log2(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.microrna.org/microrna/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.omicstudio.cn/index) for mRNAs regulated by circRNAs through ceRNA.

Cell culture

293T (human embryonic kidney) cells obtained from ATCC (Cell Systems & cGMP Biorepository, Gaithersburg, MD, USA) and were cultured in DMEM (Gibco, Gaithersburg, MD, USA) supplemented with 10% fetal bovine serum (Hyclone, Logan, UT, USA), 1% penicillin/streptomycin (Invitrogen, Carlsbad, CA, USA). All cells were cultured at 37 °C in a 5% CO2 humidified atmosphere.

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.

Results

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

Fig. 1
figure 1

Changes in body weights, muscle weights and fibers at different ages of Shouguang chickens. A Changes in body weight with age. B Changes in breast muscle with age. C Muscle fibers of the breast muscle stained with hematoxylin and eosin (H&E) (bar, 100 μm). Different uppercase superscripts in (A) and (B) indicate significant differences in body weight and breast muscle weight at different ages (P < 0.01)

Full size image

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 rRNA-depletion 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 circRNAs 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, log2(fold change) ≥ 1 or log2(fold change) ≤ -1) (Table S3).

Table 1 Summary of sequencing results in chicken breast muscle at six stages
Full size table
Fig. 2
figure 2

Distribution of circRNAs on chicken chromosomes

Full size image

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

Fig. 3
figure 3

STEM analysis of transcript profiles. A Each box corresponds to a type expression profile and only colored profiles are significantly different. B Profile 3. C Profile 16. D Profile 14

Full size image

Exploration of the CircRNA-mediated CeRNA regulatory network

Increasing evidence indicates that circRNAs can sequester relevant miRNAs through MREs to post-transcriptionally regulate gene expression. We identified 587 miRNAs in breast muscle at E12, E17, D1, D14, D56, and D98 in a previous study [28] (Table S8). To investigate ceRNA regulation in muscle and identify muscle-related circRNAs, we identified putative miRNA-circRNA and miRNA-mRNA interactions using the 587 miRNAs and transcripts identified in profiles 3, 14, and 16. Due to the small number of circRNAs in profiles 3 (one circRNA) and 14 (two circRNAs), a ceRNA network was not constructed with these profiles. In profile16, we obtained 823 miRNA-circRNA and 15,474 miRNA-mRNA interaction pairs (Table S9) and constructed a circRNA-miRNA-mRNA triple network, which included 361 miRNAs, 68 circRNAs, 599 mRNAs, and 31,063 interaction pairs (Table S10). Interestingly, in profile 16 we observed that two circRNA isoforms (circRNA225, circRNA226) had the same host gene, glycerol-3-phosphate dehydrogenase 2 (GPD2). The regulatory network showed that circRNA225 and circRNA226 co-regulated seven miRNAs (gga-miR-106-3p, gga-miR-12239-5p, gga-miR-12283-3p, gga-miR-1306-5p, gga-miR-1773-3p, gga-miR-183, and gga-miR-3594-3p) and 207 mRNAs through ceRNA regulation (Table S11).

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 circRNAs that regulated metabolism increased with increased age, and they may play an important role in muscle development.

Fig. 4
figure 4

Function annotation of the circRNAs in profile 16. A Top 10 GO terms of the putative target mRNAs in the ceRNA network using transcripts in profile 16. B Top 15 pathways of the putative target mRNAs in the ceRNA network using transcripts in profile 16

Full size image

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 circRNAs 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).

Fig. 5
figure 5

Experimental validation of circular RNAs. A Divergent primers amplify circRNAs in cDNA, but not genomic DNA (gDNA). Black triangles represent convergent primers and white triangles represent divergent primers. The gels were cropped to improve the clarity and conciseness of the presentation (B) Sanger sequencing confirmed the back-splicing junction sequence of circRNAs. C qRT-PCR showing resistance of circRNAs to RNaseR digestion. (D) The qRT-PCR validation of six differentially expressed circRNAs in six development time points

Full size image

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

Fig. 6
figure 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

Full size image

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, log2(fold change) ≥ 1 or log2(fold change) ≤ -1). Our study identified many circRNAs that were highly expressed after birth and found many differences in circRNA expression at different stages before and after birth. For example, there are fewer DE-circRNAs within hatching periods (e.g., E 17 versus E 12 or D 1 versus E 17) or post-hatching periods (e.g., D 98 versus D 56 or D 56 versus D 14), but abundant DE-circRNAs between the hatching period and the post-hatching period (e.g., D 98 versus E 12 or D 14 versus E 17). 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).

Fig. 7
figure 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

Full size image

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 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. CHAC1-heterozygote (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 circRNA225 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.

Availability of data and materials

The raw sequence data reported in this paper have been deposited in the Genome Sequence Archive in BIG Data Center, Beijing Institute of Genomics (BIG), Chinese Academy of Sciences and is publicly accessible at http://bigd.big.ac.cn/gsa (accession no CRA002573 (circRNA), CRA002587 (miRNA), CRA001773 (mRNA)).

Abbreviations

ATP1A1:

ATPase Na + /K + transporting subunit alpha 1

CHAC1:

Cation transport regulator-like protein 1

GPD2:

Glycerol-3-phosphate dehydrogenase 2

HSPA8:

Heat shock protein alpha 8

LDHA:

Lactate dehydrogenase A

PGAM1:

Phosphoglycerate mutase 1

PKM:

Pyruvate kinase muscle

TPI1:

Triosephosphate isomerase

UCP3:

Uncoupling protein 3

References

  1. Luo W, Nie Q, Zhang X. MicroRNAs involved in skeletal muscle differentiation. J Genet Genomics. 2013;40(3):107–16.

    Article  CAS  PubMed  Google Scholar 

  2. Berkes CA, Tapscott SJ. MyoD and the transcriptional control of myogenesis. Semin Cell Dev Biol: 2005: Elsevier. 2005;16:585–95.

    CAS  Google Scholar 

  3. Nie M, Deng Z-L, Liu J, Wang D-Z. Noncoding RNAs, emerging regulators of skeletal muscle development and diseases. Biomed Res Int. 2015;2015:676575.

    PubMed  PubMed Central  Google Scholar 

  4. McDaneld TG, Smith TP, Doumit ME, Miles JR, Coutinho LL, Sonstegard TS, Matukumalli LK, Nonneman DJ, Wiedmann RT. MicroRNA transcriptome profiles during swine skeletal muscle development. BMC Genomics. 2009;10(1):77.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Cassar-Malek I, Passelaigue F, Bernard C, Léger J, Hocquette J-F. Target genes of myostatin loss-of-function in muscles of late bovine fetuses. BMC Genomics. 2007;8(1):63.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Qu S, Yang X, Li X, Wang J, Gao Y, Shang R, Sun W, Dou K, Li H. Circular RNA: a new star of noncoding RNAs. Cancer Lett. 2015;365(2):141–8.

    Article  CAS  PubMed  Google Scholar 

  7. Ouyang H, Chen X, Li W, Li Z, Nie Q, Zhang X. Circular RNA circSVIL promotes myoblast proliferation and differentiation by sponging miR-203 in chicken. Front Genet. 2018;9:172.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Chen X, Ouyang H, Wang Z, Chen B, Nie Q. A novel circular RNA generated by FGFR2 gene promotes myoblast proliferation and differentiation by sponging miR-133a-5p and miR-29b-1-5p. Cells. 2018;7(11):199.

    Article  CAS  PubMed Central  Google Scholar 

  9. Ouyang H, Chen X, Wang Z, Yu J, Jia X, Li Z, Luo W, Abdalla BA, Jebessa E, Nie Q. Circular RNAs are abundant and dynamically expressed during embryonic muscle development in chickens. DNA Res. 2017;25(1):71–86.

    Article  CAS  PubMed Central  Google Scholar 

  10. Ylihärsilä H, Kajantie E, Osmond C, Forsen T, Barker DJ, Eriksson JG. Birth size, adult body composition and muscle strength in later life. Int J Obes. 2007;31(9):1392–9.

    Article  Google Scholar 

  11. Liu J, Lei Q, Li F, Zhou Y, Gao J, Liu W, Han H, Cao D. Dynamic transcriptomic analysis of breast muscle development from the embryonic to post-hatching periods in chickens. Front Genet. 2020;10(1308):1–12.

  12. Fridolfsson, A. K., and Ellegren, H. A simple and universal method for molecular sexing of non-ratite birds. J Avian Biol. 1999;30(1):116–21.

  13. Liu R, Wang H, Liu J, Wang J, Zheng M, Tan X, Xing S, Cui H, Li Q, Zhao G. Uncovering the embryonic development-related proteome and metabolome signatures in breast muscle and intramuscular fat of fast-and slow-growing chickens. BMC Genomics. 2017;18(1):816–816.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 2011;17(1):10–2.

    Article  Google Scholar 

  15. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9(4):357.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):R36.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Kim D, Salzberg SL. TopHat-Fusion: an algorithm for discovery of novel fusion transcripts. Genome Biol. 2011;12(8):R72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Zhang X-O, Dong R, Zhang Y, Zhang J-L, Luo Z, Zhang J, Chen L-L, Yang L. Diverse alternative back-splicing and alternative splicing landscape of circular RNAs. Genome Res. 2016;26(9):1277–87.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Zhang X-O, Wang H-B, Zhang Y, Lu X, Chen L-L, Yang L. Complementary sequence-mediated exon circularization. Cell. 2014;159(1):134–47.

    Article  CAS  PubMed  Google Scholar 

  20. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.

    Article  CAS  PubMed  Google Scholar 

  21. Ernst J, Bar-Joseph Z. STEM: a tool for the analysis of short time series gene expression data. BMC Bioinformatics. 2006;7(1):191.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Agarwal V, Bell GW, Nam J-W, Bartel DP. Predicting effective microRNA target sites in mammalian mRNAs. Elife. 2015;4:e05005.

    Article  PubMed Central  Google Scholar 

  23. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28(11):1947–51.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Kanehisa M, Furumichi M, Sato Y, Ishiguro-Watanabe M, Tanabe M. KEGG: integrating viruses and cellular organisms. Nucleic Acids Res. 2021;49(D1):D545–51.

    Article  CAS  PubMed  Google Scholar 

  27. Jeck WR, Sorrentino JA, Wang K, Slevin MK, Burd CE, Liu J, Marzluff WF, Sharpless NE. Circular RNAs are abundant, conserved, and associated with ALU repeats. RNA. 2013;19(2):141–57.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Liu J, Li F, Hu X, Cao D, Liu W, Han H, Zhou Y, Lei Q. Deciphering the miRNA transcriptome of breast muscle from the embryonic to post-hatching periods in chickens. BMC Genomics. 2021;22:64.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Hong L, Gu T, He Y, Zhou C, Hu Q, Xingwang W, Zheng E, Huang S, Xu Z, Yang J. Genome-wide analysis of circular RNAs mediated ceRNA regulation in porcine embryonic muscle development. Front Cell Dev Biol. 2019;7:289.

    Article  PubMed  PubMed Central  Google Scholar 

  30. Abdelmohsen K, Panda AC, De S, Grammatikakis I, Kim J, Ding J, Noh JH, Kim KM, Mattison JA, de Cabo R. Circular RNAs in monkey muscle: age-dependent changes. Aging (Albany NY). 2015;7(11):903.

    Article  CAS  Google Scholar 

  31. Liu M, Wang Q, Shen J, Yang BB, Ding X. Circbank: a comprehensive database for circRNA with standard nomenclature. RNA Biol. 2019;16(7):899–905.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Zhan S, Zhao W, Song T, Dong Y, Guo J, Cao J, Zhong T, Wang L, Li L, Zhang H. Dynamic transcriptomic analysis in hircine longissimus dorsi muscle from fetal to neonatal development stages. Funct Integr Genomics. 2018;18(1):43–54.

    Article  CAS  PubMed  Google Scholar 

  33. Braun T, Gautel M. Transcriptional mechanisms regulating skeletal muscle differentiation, growth and homeostasis. Nat Rev Mol Cell Biol. 2011;12(6):349.

    Article  CAS  PubMed  Google Scholar 

  34. Duan Y, Li F, Li Y, Tang Y, Kong X, Feng Z, Anthony TG, Watford M, Hou Y, Wu G. The role of leucine and its metabolites in protein and energy metabolism. Amino Acids. 2016;48(1):41–51.

    Article  CAS  PubMed  Google Scholar 

  35. Dambal S, Shah M, Mihelich B, Nonn L. The microRNA-183 cluster: the family that plays together stays together. Nucleic Acids Res. 2015;43(15):7173–88.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Nguyen MT, Min K-H, Lee W. MiR-183-5p induced by saturated fatty acids regulates the myogenic differentiation by directly targeting FHL1 in C2C12 myoblasts. BMB Rep. 2020;53(11):605.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Kim J, Yoon H, Ramírez CM, Lee S-M, Hoe H-S, Fernández-Hernando C, Kim J. MiR-106b impairs cholesterol efflux and increases Aβ levels by repressing ABCA1 expression. Exp Neurol. 2012;235(2):476–83.

    Article  CAS  PubMed  Google Scholar 

  38. Eisenberg I, Eran A, Nishino I, Moggio M, Lamperti C, Amato AA, Lidov HG, Kang PB, North KN, Mitrani-Rosenbaum S. Distinctive patterns of microRNA expression in primary muscular disorders. Proc Natl Acad Sci. 2007;104(43):17016–21.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Liu M, Wang Y, Ruan Y, Bai C, Qiu L, Cui Y, Ying G, Li B. PKM2 promotes reductive glutamine metabolism. Cancer Biol Med. 2018;15(4):389.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Koopman R, Ly CH, Ryall JG. A metabolic link to skeletal muscle wasting and regeneration. Front Physiol. 2014;5:32.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Feng Y, Zhang X, Zhang S, Xu S, Chen X, Zhou C, Wang X, Xie X, Lu W. Paclitaxel resistance in ovarian cancers relies on a PGAM1 mediated glycolytic metabolism. 2020.

  42. Teltathum T, Mekchay S. Proteome changes in Thai indigenous chicken muscle during growth period. Int J Biol Sci. 2009;5(7):679.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Schrauwen P, Hesselink M. UCP2 and UCP3 in muscle controlling body metabolism. J Exp Biol. 2002;205(15):2275–85.

    Article  CAS  PubMed  Google Scholar 

  44. Bezaire V, Spriet L, Campbell S, Sabet N, Gerrits M, Bonen A, Harper ME. Constitutive UCP3 overexpression at physiological levels increases mouse skeletal muscle capacity for fatty acid transport and oxidation. FASEB J. 2005;19(8):977–9.

    Article  CAS  PubMed  Google Scholar 

  45. Kutz LC: Role of ATP1A1 in Skeletal Muscle Growth and Metabolism. 2020.

  46. Crawford RR: Regulation and function of the unfolded protein response gene CHAC1. Louisiana State University Health Sciences Center; 2016.

  47. Scimone C, Alibrandi S, Donato L, Esposito T, Sidoti A, D’Angelo R. VARIANTS OF THE MOLECULAR CHAPERONE HSPA8 AND HSPA1A GENES IN TRIMETHYLAMINURIA: A PILOT STUDY. 2020.

  48. Xu T, Huang W, Zhang X, Ye B, Zhou H, Hou S. Identification and characterization of genes related to the development of breast muscles in Pekin duck. Mol Biol Rep. 2012;39(7):7647–55.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We thank International Science Editing (http://www.internationalscienceediting.com) for editing this manuscript.

Funding

This research was funded by the Natural Science Foundation of Shandong province (ZR2019BC077), the Agricultural Scientific and Technological Innovation Project of Shandong Academy of Agricultural Sciences (CXGC2021A13). China Agriculture Research System of MOF and MARA (CARS-41). Agricultural Breed Project of Shandong Province(2020LZGC013). Shandong Provincial Natural Science Foundation (ZR2020MC169); Collection, Protection and Accurate Identification of Livestock Germplasm Resources (2019LZGC019); China Agriculture Research System of MOF and MARA(CARS-40); Major Scientific and Technological Innovation Project (MSTIP): the Research and Demonstration on Key Technologies of Precision Breeding and Management of Laying Hens(2019JZZY020611).

Author information

Authors and Affiliations

  1. Poultry Institute, Shandong Academy of Agricultural Sciences, Ji’nan, 250023, China

    Qiuxia Lei, Xin Hu, Haixia Han, Jie Wang, Wei Liu, Yan Zhou, Dingguo Cao, Fuwei Li & Jie Liu

  2. Poultry Breeding Engineering Technology Center of Shandong Province, Ji’nan, 250023, China

    Qiuxia Lei & Jie Liu

  3. Institute of Animal Sciences, Chinese Academy of Agricultural Sciences, Beijing, 100193, China

    Xin Hu

  4. Molecular and Cellular Biology, Gembloux Agro-Bio Tech, University of Liège, 5030, Gembloux, Belgium

    Xin Hu

Authors
  1. Qiuxia Lei
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Xin Hu
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Haixia Han
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Jie Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Wei Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Yan Zhou
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Dingguo Cao
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Fuwei Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Jie Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

QL performed experiments and data analysis and draft writing. XH, HH, WL and YZ contributed to animal experiments, JW, DC and FL contributed to data analysis. JL designed experiments and supervised and coordinated the study. All authors reviewed the manuscript. The author(s) read and approved the final manuscript.

Corresponding author

Correspondence to Jie Liu.

Ethics declarations

Ethics approval and consent to participate

All of the animal experiments were conducted in accordance with the Guidelines for Experimental Animals, established by the Ministry of Science and Technology (Beijing, China) and the study complies with the ARRIVE guidelines (https://arriveguidelines.org). Animal experiments were approved by the Science Research Department of the Shandong Academy of Agricultural Sciences (SAAS) (Ji’nan, China). Ethical approval for animal survival was given by the animal ethics committee of SAAS (No. SAAS-2019–029).

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1:

Table S1. Detailed information of primers used in this study.

Additional file 2:

Table S2. The expression of all identified circRNAs. 

Additional file 3:

Table S3. The differentially expressed circRNAs.

Additional file 4:

Table S4. The DE-mRNAs identified in our previous study.

Additional file 5:

Table S5. The circRNAs and mRNAs identified in profile 3.

Additional file 6:

Table S6. The circRNAs and mRNAs identified in profile 16.

Additional file 7:

Table S7. The circRNAs and mRNAs identified in profile 14.

Additional file 8:

Table S8. The miRNAs identified in our previous study.

Additional file 9:

Table S9. The miRNA-circRNA interaction pairs and miRNA-mRNA interaction pairs in profile 16.

Additional file 10:

Table S10. The circRNA-miRNA-mRNA network of profile 16.

Additional file 11:

Table S11. The circRNA-miRNA-mRNA network for circRNA225 and circRNA226.

Additional file 12:

Table S12. Function annotation of the circRNAs in profile 16.

Additional file 13:

Figure S1. Predicted biomathematical circRNA-miRNA-mRNA network for circRNA225 and circRNA226. Only the top 50 mRNAs are shown in the network. Yellow ellipses represent circRNAs, green diamonds represent miRNAs, and red triangles represent mRNAs.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Lei, Q., Hu, X., Han, H. et al. Integrative analysis of circRNA, miRNA, and mRNA profiles to reveal ceRNA regulation in chicken muscle development from the embryonic to post-hatching periods. BMC Genomics 23, 342 (2022). https://doi.org/10.1186/s12864-022-08525-5

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-022-08525-5

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us