Dynamic Transcriptomic Analysis of Breast Muscle Development From the Embryonic to Post-hatching Periods in Chickens

Skeletal muscle development and growth are closely associated with efficiency of poultry meat production and its quality. We performed whole transcriptome profiling based on RNA sequencing of breast muscle tissue obtained from Shouguang chickens at embryonic days (E) 12 and 17 to post-hatching days (D) 1, 14, 56, and 98. A total of 9,447 differentially expressed genes (DEGs) were filtered (Q < 0.01, fold change > 2). Time series expression profile clustering analysis identified five significantly different expression profiles that were divided into three clusters. DEGs from cluster I with downregulated pattern were significantly enriched in cell proliferation processes such as cell cycle, mitotic nuclear division, and DNA replication. DEGs from cluster II with upregulated pattern were significantly enriched in metabolic processes such as glycolysis/gluconeogenesis, insulin signaling pathway, calcium signaling pathway, and biosynthesis of amino acids. DEGs from cluster III, with a pattern that increased from E17 to D1 and then decreased from D1 to D14, mainly contributed to lipid metabolism. Therefore, this study may help us explain the mechanisms underlying the phenotype that myofiber hyperplasia occurs predominantly during embryogenesis and hypertrophy occurs mainly after birth at the transcriptional level. Moreover, lipid metabolism may contribute to the early muscle development and growth. These findings add to our knowledge of muscle development in chickens.


INTRODUCTION
In chicken production, skeletal muscle development is closely associated with the amount of meat production and its quality, ultimately affecting the economic benefits. Therefore, elucidating the molecular mechanisms underlying chicken skeletal muscle development is of vital interest. The muscle mass is determined by cell numbers and unit cell size. Hyperplasia refers to the increases in cell number or muscle fiber number that occur mainly in the embryonic period, as the number of muscle fibers is fixed by the day of hatching. However, hypertrophy refers to the increase in cell size that occurs mainly after birth (Ylihärsilä et al., 2007;Liu et al., 2017b;Ouyang et al., 2017). Therefore, there may be distinct molecular processes that occur in chicken muscle development between the embryonic and posthatching periods.
Over the past few years, there has been much progress in exploring the molecular mechanisms underlying muscle growth and development in chickens, but most of the studies focused on the embryonic or post-hatching period. Davis et al. (2015) characterized the transcriptome of Ross 708 chicken breast muscle at specified time points from 6 to 21 days after hatching. Li et al. (2019) explored the messenger RNA (mRNA) and microRNA (miRNA) profiles of Gushi chicken muscle tissues in the late postnatal stage (6, 14, 22, and 30 weeks). Our previous study examined the protein expression profiles in the breast muscle of Beijing-You chickens at ages 1, 56, 98, and 140 days using isobaric tags for relative and absolute quantification (Liu et al., 2016). Li et al. (2017) and Ouyang et al. (2017) explored the transcriptome and protein expression profiles in leg muscle tissues of Xinghua chicken at embryonic days (E) 11 and E16 and post-hatching day (D) 1, respectively. However, few studies have paid attention to the whole muscle development from embryonic to post-hatching periods in chickens. Only Liu et al. (2017b) investigated the proteomes of breast muscle in Cobb and Beijing-You chickens at E12, E17, D1, and D14.
Shouguang chickens that have been breed in China for 2,000 years are a dual-purpose breed with large bodies (Gao et al., 2008), which may be excellent material for studying muscle development. Therefore, we chose the critical breast muscle developmental stages in the embryonic to post-hatching periods (E12 and E17 and D1, D14, D56, and D98) of Shouguang chickens for quantitative analysis of the gene expression profile of breast muscle, which may help us explore the development-related genes expression signatures in breast muscle and its distinction between embryonic and posthatching periods.

Animals
Shouguang chicken eggs were obtained from the experimental farm of the Poultry Institute (PS), Shandong Academy of Agricultural Sciences (SAAS, Jinan, China). All eggs were incubated with the normal procedure and chicks were reared in cages using standard conditions of temperature, humidity, and ventilation at the farm of the PS, SAAS. The same diet was fed to all chickens and a three-phase feeding system was used: starter ration (days 1-28) with 21.0% crude protein and 12.12 MJ/kg; second phase (days 28-56) with 19.0% crude protein and 12.54 MJ/kg; and final phase (after day 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 E12, E17, D1, D14, D56, and D98. All fresh breast muscle tissue samples were collected, frozen in liquid nitrogen, and stored at −80°C until RNA extraction. The sex of chicken embryos was identified by polymerase chain reaction (PCR) of the CHD1 gene (Fridolfsson and Ellegren, 1999). Chickens with two bands of 600 and 450 bp were born as female, and those with one band of 600 bp were born as male.

RNA Extraction, cRNA Library Construction, and Sequencing
Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). Three female chickens at each stage (except E17) were used for further experiments. Total RNA quantity and purity were analyzed using a Bioanalyer 2100 (Agilent, Santa Clara, CA, USA) with RNA integrity number >7.0. Approximately 10 mg total RNA was used to deplete rRNA using the Epicentre Ribo-Zero Gold Kit (Illumina, San Diego, CA, USA). Following purification, the poly(A)− or poly(A)+ RNA fraction was fragmented into small pieces using divalent cations under elevated temperature. The cleaved RNA fragments were reverse-transcribed to create the final complementary DNA (cDNA) library in accordance with the protocol for the RNA sequencing (RNA-Seq) sample preparation kit (Illumina). The average insect size for the paired-end libraries was 300 ± 50 bp. We performed paired-end sequencing on an Illumina Hiseq 4000 at LC-Bio, China.

RNA-Seq Reads Mapping and DEG Analysis
We aligned reads to the genome of Gallus_gallus 5.0 (GCA_000002315.3) using HISAT package (Kim et al., 2015), which initially removed reads based on quality information accompanying each read and then mapped the reads to the reference genome. The mapped reads of each sample were assembled using StringTie (Pertea et al., 2015). All transcriptomes from samples were merged to reconstruct a comprehensive transcriptome using perl scripts. After the final transcriptome was generated, StringTie and edgeR (Robinson et al., 2010) were used to estimate the expression levels of all transcripts. StringTie was used to perform expression level for mRNAs by calculating fragments per kilobase of transcript per million fragments mapped (FPKM). Differentially expressed genes (DEGs) were selected with log2 (fold change) > 1 or log2 (fold change) less than −1 with statistical significance (Q < 0.01) by R package. 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. CRA001773).

Time Series Expression Profile Clustering
The non-parametric clustering algorithm of STEM (Short Time-Series Expression Miner, version 1.3.11) (Ernst and Bar-Joseph, 2006) was used to cluster and visualize the expression patterns of DEGs. Expression profiles of DEGs were clustered based on their log2 (FPKM values) and their correlation coefficients. The maximum unit change in model profiles between time points was adjusted to 2 and the maximum number of model profiles to 50. The statistical significance of the number of DEGs to each profile versus the expected number was computed using the algorithm proposed by Ernst and Bar-Joseph (2006).

Functional Annotation
Functional analysis of DEGs was performed using the DAVID (Database for Annotation, Visualization and Integrated Discovery) tool (http://david.abcc.ncifcrf.gov/) (Dennis et al., 2003). The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was performed using KOBAS version 3.0 (Xie et al., 2011). Gene Ontology (GO) terms and KEGG pathways with P < 0.05 were considered significantly enriched groups of genes possibly contributing to muscle development.

qRT-PCR Confirmation
To confirm our differential expression results, we conducted quantitative reverse transcription PCR (qRT-PCR) for six selected genes (MYOG, MYH11, TNNI2, TNNT3, TNNC2, and TPM2). The total RNA was used for first-strand cDNA synthesis using a commercial kit (TaKaRa, Dalian, China). cDNA was subsequently used for qRT-PCR analyses with an ABI 7500 Detection System (Applied Biosystems, Foster City, CA, USA) and primers designed using Primer Premier version 5.0 (PREMIER Biosoft, Palo Alto, CA, USA), as listed in Table S1. mRNA abundance of candidate genes was determined using the KAPA SYBR ® FAST qPCR Master Mix (2×) Universal Cocktail (KAPA Biosystems, Boston, MA, USA).
qRT-PCR was performed following the instructions of ABI 7500 with default parameters. The 2 −DDCt method (Livak and Schmittgen, 2001) was used to calculate the relative mRNA abundance. The beta actin gene (ACTB) was used as the housekeeping gene. Three independent replications were used for each assay and data were presented as means ± SD.

Overall Assessment for Sequencing Data Mapping Statistics
To identify mRNA expressed in breast muscle tissue development of chickens, we constructed 17 cDNA libraries (E12_1, E12_2, E12_3, E17_1, E17_2, D1_1, D1_2, D1_3, D14_1, D14_2, D14_3, D56_1, D56_2, D56_3, D98_1, D98_2, and D98_3) from breast muscle samples at six developmental stages. As shown in Table 1, 68,431,306-103,358,850 raw reads were generated in the 17 libraries, and 65,384,366-100,882,250 clean reads were obtained after discarding adaptor sequences and low-quality reads. We mapped clean reads to chicken reference genome Gallus_gallus 5.0 and found that 84.41-90.10% of the clean reads in the libraries were mapped to the chicken reference genome ( Table 1).

Differential Expression Analysis of Genes
In pairwise comparisons between the libraries of breast muscle at the six developmental stages, a total of 9,447 genes were differentially expressed (Q < 0.01, fold change > 2) ( Figure 1 and Table S2). There were 2,502, 4,582, 4,394, 3,689, and 4,607 DEGs in E17, D1, D14, D56, and D98 compared to E12. Comparing successive ages within each region, 2,502, 2,429, 1,839, 262, and 144 DEGs were found in E17 versus E12, D1 versus E17, D14 versus D1, D56 versus D14, and D98 versus D56, respectively. The numbers of DEGs were greatest in E17 versus E12 and lowest in D98 versus D56, which indicated that regional differences in gene expression were greatest during the earlier stages of embryo development.

STEM Analysis of DEG Profiles
As our data were collected at different time points, STEM was used to cluster and visualize possible changes in the profiles of 9,447 DEGs at six time points of breast muscle development. Within the 50 model profiles, five expression profiles containing 5,269 genes were statistically significant (P < 0.05; Figure 2A and Table S3). Of these, profiles 8 and 12 with downregulated patterns contained 3,233 and 693 DEGs, respectively ( Figure  2B and Table S3), while profiles 39 and 49 with upregulated patterns contained 380 and 156 DEGs, respectively ( Figure 2C and Table S3). Profile 25 with 717 genes as the third pattern showed an increase from E17 to D1 and reached a peak at D1, then decreased from D1 to D14 and remained stable from D14 to D98 ( Figure 2D and Table S3). Thus, the expression pattern of DEGs can be divided into three clusters: cluster I (profiles 8 and 12, total of 3,926 DEGs) with downregulated pattern; cluster II (profiles 39 and 49, total of 536 DEGs) with upregulated pattern; and cluster III (profile 25, total of 717 DEGs). The results provide new information related to further characterization of novel molecules associated with skeletal muscle development in chickens.

GO Enrichment Analysis
To explore the biological function of DEGs, GO enrichment analysis was performed based on cluster analysis. The genes in cluster I (profiles 8 and 12) were significantly enriched in 138 GO terms (68 under biological process, 38 under cellular component, and 32 under molecular function) ( Table S4). Within the biological process category, the most abundant GO terms consisted of DNA replication, cell division, ATP-dependent chromatin remodeling, mitotic nuclear division, and DNA repair ( Figure 3A). The genes in cluster II (profiles 39 and 49) were significantly enriched in 34 GO terms (16 under biological process, six under cellular component, and 12 under molecular function) ( Table S5). Within the biological process category, the most abundant GO terms consisted of carbohydrate metabolic process, xanthine catabolic process, gluconeogenesis, and positive regulation of interferon-g production ( Figure 3B). Skeletal muscle contraction was also included in this category ( Figure 3B). We also analyzed the biological function of genes in cluster III (profile 25). Eighteen GO terms (10 under biological process, three under cellular component, and five under molecular function) were significantly enriched (Table S6). Within the biological process category, the most abundant GO terms consisted of fatty acid beta-oxidation using acyl-CoA dehydrogenase, fatty acid beta-oxidation, positive regulation of focal adhesion assembly, lipid homeostasis, and regulation of stress fiber assembly ( Figure 3C).

KEGG Enrichment Analysis
We used KEGG pathway analysis to explore the signaling pathways of DEGs involved in cluster I (profiles 8 and 12), cluster II (profiles 39 and 49), and cluster III (profile 25). For cluster I, the DEGs were significantly enriched in 13 pathways ( Figure 4A and Table S7), and half of these pathways are involved in cell division, such as cell cycle, DNA replication, nucleotide excision repair, mismatch repair, oocyte meiosis, and spliceosome. As for cluster II, the DEGs were significantly enriched in 16 pathways ( Figure 4B and Table S8), and it was interesting that all of these pathways were directly or indirectly involved in the main metabolic processes of the organism, such   Figure 4C and Table S9), and most of these pathways were associated with metabolism. Among these metabolic pathways, the most enriched pathways were those related to lipid metabolism, such Frontiers in Genetics | www.frontiersin.org January 2020 | Volume 10 | Article 1308 as fatty acid degradation, propanoate metabolism, butanoate metabolism, fatty acid elongation, synthesis and degradation of ketone bodies, and fatty acid metabolism.

Validation of DEGs by qRT-PCR
The qRT-PCR assays were conducted to validate six selected DEGs from RNA-Seq: MYOG, MYH11, TNNI2, TNNT3, TNNC2, and TPM2. Relative expression changes of qRT-PCR data were highly (r = 0.83-0.99) correlated with RNA-Seq data ( Figure 5), suggesting the reliability of the RNA-Seq approach.

DISCUSSION
Skeletal muscle growth and development includes a series of closely regulated changes in gene expression level from embryo to adult, and uncovering the gene expression patterns underneath chicken skeletal muscle development contributes to meat production. Previous transcriptome analysis of chicken muscle only concentrated on the embryonic period or adult stage, and few studies have systematically examined the transcriptome of chicken skeletal muscle development from the embryonic period to the growing period. To investigate the mechanisms of skeletal muscle development systematically, we used RNA-Seq to generate extensive cDNA libraries for six developmental stages of chickens from E12 to D98. As shown in Figure 1 and Table S2, a total of 9,447 DEGs were identified in pairwise comparisons between the libraries of breast muscle at the six developmental stages, and the regional differences in gene expression were greatest during the earlier stages of embryo development (E17 versus E12) than the late postnatal stage (D98 versus D56). Previous studies have demonstrated the increase in cell number or muscle fiber number which occurs mainly in embryonic periods as the numbers of muscle fibers were fixed by the day of hatching and then impacted on the postnatal accretion of muscle mass (Smith, 1963;Sporer et al., 2011). Thus, embryonic periods are the critical periods for muscle development and more genes are active in these periods. Six selected DEGs involved in muscle development were validated using qRT-PCR, and the results were consistent with those from RNA-Seq, suggesting reliability of the identified DEGs through RNA-Seq ( Figure 5).
Since muscle development was accompanied by the differential expression of related genes in different growth periods and our data were also collected at different time points, we used STEM software, which is widely used to study dynamic biological processes (Guo et al., 2016;Ma et al., 2018;Zhan et al., 2018), to investigate the dynamic genetic changes during breast muscle development. We experimented with various numbers of profiles and found that five profiles (combining with three clusters) best captured the expression patterns of DEGs (Figure 2 and Table S3). The genes in a cluster have similar temporal expression patterns and may be involved in the same biological process. Therefore, we performed GO and KEGG analyses to explore the function of the DEGs with similar temporal expression patterns.
For cluster I with downregulated pattern ( Figure 2B and Table S3), 3,926 genes were significantly clustered, which was more than 40% of the total DEGs. These genes were more highly expressed in the early periods of muscle development than the late stages of growth, which further confirmed that the early development might play key roles in muscle growth. GO functional annotation and KEGG analysis both showed that the downregulated genes were significantly enriched in cell proliferation, including DNA replication, cell cycle, cell division, mitotic nuclear division, and DNA repair ( Figures 3A  and 4A), which was similar to previous research on goat muscle development from gestation to birth which showed that genes with downregulated patterns were also involved in cell proliferation processes (Zhan et al., 2018). Therefore, these results further support the hypothesis that the total number of skeletal myofibers is defined by hyperplasia during embryogenesis. Among these downregulated genes, CCNA2, CCNB2, and CDK1, which encode the cyclins and their cognate cyclin-dependent protein kinases, were not only significantly enriched in the biological process of cell division but also significantly enriched in cell cycle pathway (Tables S4 and S7). Cyclin A2 possesses a unique role in its two-point control of the cell cycle, first by interacting with CDK2 in controlling the G1/S transition into DNA synthesis and then by interacting with CDKs 1 and 2 to control the G2/M entry into mitosis (Li et al., 1998). Previous study has been demonstrated that constitutive expression of cyclin A2 in a transgenic mouse yields robust postnatal cardiomyocyte mitosis and hyperplasia (Chaudhry et al., 2004). Cyclin B2 was also demonstrated to have a regulatory role in chicken breast muscle development (Li et al., 2019). Moreover, CDK1 and CDK2 play integral roles in reducing MyoD activity during myoblast proliferation by phosphorylating MyoD (Kitzmann et al., 1999). These results suggest that the genes with a downregulated pattern of expression play regulatory roles in chicken breast muscle development through the processes involved in the early stages of cell proliferation, and genes related to cyclins and their cognate cyclin-dependent protein kinases may be critical factors in regulating cell proliferation.
For cluster II with an upregulated pattern of expression ( Figure 2C and Table S3), the functional annotation and pathway analysis both showed that these genes were significantly enriched in metabolism such as carbohydrate metabolism, glycolysis/gluconeogenesis, calcium signaling pathway, insulin signaling pathway, and biosynthesis of amino acids ( Figures 3B and 4B). A previous study of goat muscle development also found that genes with upregulated patterns of expression were related to metabolic pathways, such as biosynthesis of amino acids, glycolysis/gluconeogenesis, and the TCA cycle (Zhan et al., 2018). Among these genes, PHKG1, PPP1R3C, and FBP2 were significantly enriched not only in glycogen biosynthetic process and gluconeogenesis but also in insulin signaling pathway and calcium signaling pathway (Tables S5 and S8). PHKG1, as a key factor in insulin signaling and calcium signaling pathways, encodes the catalytic subunit of phosphorylase kinase, which functions in the cascade activation of glycogen breakdown in muscle tissue (Ma et al., 2014). Fructose-1,6-bisphosphatase encoded by FBP2 catalyzes the hydrolysis of fructose-1,6-bisphosphate to fructose-6-phosphate and inorganic phosphate, which plays a regulatory role in the synthesis of glycogen/glucose. Previous findings point to FBP2 as an important link between calcium-induced muscle contractive and metabolic (glycolytic) activity, mitochondrial function, and cell survival (Pirog et al., 2014). PPP1R3C was also called protein targeting to glycogen (PTG) and regulated glycogen metabolism (Ji et al., 2019). Moreover, glycolytic process and biosynthesis of amino acids pathway were both significantly enriched in PKLR, PGK2, and TPI1 (Tables S5 and S8); thus, these genes may be important regulatory switches for protein and energy conversion and ultimately influence muscle development. A number of studies have demonstrated that muscle mass increased by hypertrophy (increased cellular protein content) after hatching and was controlled by synthesis of muscle proteins or their degradation (Braun and Gautel, 2011). Protein and energy metabolism are tightly coupled, and the energy from glycolysis/ gluconeogenesis is needed for protein turnover during skeletal muscle development (Duan et al., 2016;Liu et al., 2016). Therefore, these results show that genes involved in metabolism may be critical for postnatal myofiber growth, muscle hypertrophy, and muscle regeneration, and the protein synthesis and energy metabolism of skeletal muscle regulated by insulin signaling pathway and calcium signaling pathway may be important for coordinating muscle development.
Cluster III (profile 25) showed an increase from E17 to D1 and reached a peak at D1, then decreased from D1 to D14 and remained stable from D14 to D98 ( Figure 2D and Table S3). The functional annotation showed that processes and pathways involved in lipid metabolism were significantly enriched, such as fatty acid b-oxidation, lipid homeostasis, fatty acid degradation, and propanoate metabolism ( Figures 3C and 4C). Previous studies have demonstrated that lipids stored in the adipocytes during embryonic life are transferred to the muscle fibers and used for growth and energy requirements at the early stage, while muscle again stores lipids in later life (Chartrin et al., 2007;Liu et al., 2016). Interestingly, ACADL, ACAT1, HADHA, ACADS, ECHS1, and AUH, which are significantly enriched in fatty acid beta-oxidation, were active during embryonic life in the present study (Table S6), which further demonstrated that the lipids were the important energy source for muscle development and growth at the early stage. Moreover, some of these genes were also the key regulatory molecules for intramuscular fat (IMF) deposition. For example, ACADL and HADHA have been identified as candidate biomarkers for IMF deposition in Cobb and Beijing-You chickens (Liu et al., 2017a), and it was interesting that their expression patterns in Shouguang chickens were similar to Cobb and Beijing-You chickens. Moreover, ACAT1 expression was significantly lower in muscle of AA chickens with low IMF content than in Beijing-You chickens with abundant IMF (Liu et al., 2018), suggesting that ACAT1 may contribute to IMF deposition. These results suggest that genes involved in lipid metabolism, and especially those related to fatty acid beta-oxidation, play important roles in early muscle development and deposition of IMF.

CONCLUSION
In the present study, we systematically identified DEGs and investigated their temporal expression profiles during chicken breast muscle development from E12 to D98. A total of 9,447 DEGs were identified in chicken breast muscle and showed three significantly different expression patterns. Functional enrichment analysis suggests that genes with downregulated patterns contribute to cell proliferation processes, while genes with upregulated patterns are mainly involved in metabolism. Genes related to lipid metabolism change dramatically around the time of birth, which may play important roles in early muscle development and deposition of IMF. In summary, our study will facilitate understanding of the mechanisms underlying the phenotype that myofiber hyperplasia occurs predominantly during embryogenesis and hypertrophy occurs mainly after birth at the transcriptional level. These findings elucidate the regulatory mechanisms involved in chicken breast muscle development.

DATA AVAILABILITY STATEMENT
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 CRA001773).

ETHICS STATEMENT
All experiments were approved by the Animal Care Committee of the Shandong Academy of Agricultural Sciences (Ji'nan, China). The experimental procedures with chickens were performed according to the Guidelines for Experimental Animals established by the Ministry of Science and Technology (Beijing, China).

AUTHOR CONTRIBUTIONS
JL and QL performed experiments and data analysis and draft writing. FL, YZ, JG, and WL contributed to animal experiments and data analysis. HH and DC designed experiments and supervised and coordinated the study. All authors reviewed the manuscript.