Transcriptome analysis reveals potential mechanisms and pathways underlying embryonic development with respect to muscle growth and egg production in slow and fast growing chickens

Background: Chicken is one of the important meat sources throughout the globe. Muscle development and egg production are important genetic traits in commercially raising chickens. However, not much information is available in the fast and slow growth of chicken to determine the expression of genes involved in muscle development and egg production in embryo initiation and developmental stages. This study was designed to investigate why improved Aseel (PD4) growing slowly compared with the control broiler (CB), microarray was conducted with the 7 th -day embryo and 18 th -day thigh muscle of improved Aseel (PD4) and control broiler (CL), respectively. Results: In the differential transcripts screening, all the transcripts obtained by microarray of slow and fast growth groups were screened by fold change ≥ 1 and false discovery rate (FDR) <0.05. In total, 19022 transcripts were differentially expressed between the 7 th -day embryo and 18 th -day thigh muscle of improved Aseel compared to the control broiler. Further analysis showed that a high number of transcripts are differentially regulated in the 7 th -day improved Aseel embryo (15382) and fewer transcripts were differentially regulated (3640) in the 18 th -day thigh muscle of improved Aseel compared to control broiler. In the 7 th and 18 th -day improved Aseel embryo, 10127, 2102, 5255, and 1538 transcripts were up and down-regulated, respectively. The commonly up and downregulated transcripts are 545 and 381 between the 7 th and 18 th -day of embryos. In this study, we have selected 18 Gallus gallus candidate reference genes from NCBI and total RNA was isolated from control broiler, improved Aseel embryo tissues, and studied their expression proles by real-time quantitative PCR (qPCR). The best housekeeping gene was identied by using geNorm, NormFinder, BestKeeper, Delta CT, and RefFinder analytical software. The result showed that the TFRC gene is the most stable and further it is used for qPCR data normalization. Further, to validate the differentially expressed genes (DEGs) related to muscle growth, myostatin signaling and development, fatty acid metabolism genes in improved Aseel (PD4) and control broiler embryo tissues by qPCR. Conclusion: Our study identied DEGs that regulate myostatin signaling and differentiation pathway, glycolysis and gluconeogenesis, fatty acid metabolism, Jak-STAT, mTOR, and TGF-β signaling pathways, tryptophan metabolism, PI3K-Akt signaling pathways in improved Aseel. The results revealed that the gene expression architecture is present in the improved Aseel exhibiting embryo growth that will help to improve muscle development, differentiation, egg production, as well as protein synthesis in improved Aseel native chicken. Our ndings may be used as a model for improving the growth in improved Aseel as well as optimizing the growth in the control broiler. study. For fast-growing breast meat yield, functional analysis revealed the favoring metabolic shifts towards alternative catabolic pathways, oxidative stress, inammatory, regeneration, brosis processes, cellular defense, and remodeling [8]. A transcriptome proling analysis was performed in two chicken lines i.e. high (pHu+) and low (pHu-) using an Agilent custom chicken 8×60K microarray. Between two lines, 1436 differentially expressed (DE) genes were found and they are related to biological processes for muscle development and remodeling, carbohydrate, and energy metabolism [9]. A genome-wide association study (GWAS) was conducted for assessing body weight in an F2 chicken population and a microarray expression study was done with the liver of high and low-weight chickens. Also, identied miR-16 as a key regulator, it will suppress the chicken embryo cell proliferation and cellular growth. The mutated miR-16 by inserting 54bp is showed a signicant increase in body weight, bone size, and muscle mass [10]. The comparative transcriptome analysis of grouper sh muscle in hybrid and its parents showed up-regulation of genes related to glycolysis, calcium signaling, and troponin pathways led to enhance muscle growth in the hybrid grouper [11]. Insulin-like growth factor 1 (IGF1) and a cascade of intracellular components (protein kinase B, mTOR, GSK3β, and FoxO) play a major role in the regulation of skeletal muscle growth during development and regeneration [12]. the 7 th and 18 th embryos of the control broiler and improved Aseel were GAPDH (0.371) and TFRC (0.4). BestKeeper calculates the gene expression variation based on Ct values and also calculates the Pearson correlation coecient by pairwise correlation analysis for all reference genes and found the stable genes. According to BestKeeper stability criteria, the most stable genes in the 7 th and 18 th embryos of the control broiler and improved Aseel were HSP70 (1.08) and PPP2CB (1.11). The Delta CT results supported geNorm and NormFinder results and it showed the best stable genes in 7 th and 18 th embryos of control broiler and improved Aseel i.e. TFRC (2.1) and PGK2 (2.13). RefFinder conclusive the calculations using the all above-mentioned algorithms and suggested the stable genes in 7 th and 18 th embryos of control broiler and improved Aseel i.e. TFRC (2.21) and PGK2 (2.85). To overcome different software program limitations, the stability of candidate reference genes was determined based on the consensus ranking for gene expression normalization in 7 th and 18 th -day embryos of control broiler and improved Aseel. Our study identied the most stable genes and indicated that, TFRC and PGK2 for 7 th and 18 th -day embryos of control broiler and improved Aseel. Our observations further strengthen the necessity to analyze the stability of candidate reference genes as suitable references. acid pH line, 3-hydroxymethyl-3-methylglutaryl-CoA lyase (HMGCL) nal leucine ketone acetyl-CoA acetyltransferase-2 (ACAT2) in β-oxidation or degradation of ketogenic amino acids, and nudix hydrolases (NUDT7, NUDT12, NUDT19) hydrolyse a nucleoside di and triphosphates, dinucleoside and diphosphoinositol polyphosphates, nucleotide sugars and RNA caps, were up-regulated. In the present study, 3-hydroxymethyl-3-methylglutaryl-CoA lyase like-1 (HMGCLL1), acetyl-CoA acetyltransferase 2 (ACAT2), nudix type motif 7 (NUDT7), and nudix type motif 21 (NUDT21) were down-regulated and carnitine/palmitoyl-transferase 1 (CPT1) were up-regulated in the 7 th -day embryo of improved Aseel They may regulate the β-oxidation in peroxisomes as well as mitochondria, this is the region may be fats required for initial embryo development and excess fats involved in β-oxidation and nally provide the energy for embryo development. eukaryotic translation initiation factor 4E (EIF4E) key downstream targets for mTORC1 protein synthesis. In low FE birds, o p70S6k higher The muscle tissue of RNAseq transcriptomic higher expression of eukaryotic initiation, elongation, and translation genes in high FE compared to low FE PedM phenotype [142, 156]. In the present study, the late endosomal/lysosomal adaptor MAPK and MTOR activator 3 (LAMTOR3), protein kinase, cAMP-dependent, regulatory, type I, alpha (tissue specic extinguisher 1) (PRKAR1A), solute carrier family 2 (facilitated glucose transporter) member 8 (SLC2A8)/ GLUT8/ GLUTX1, solute carrier family 2 (facilitated glucose transporter), member 3 (SLC2A3), ribosomal protein S6 kinase, 90kDa, polypeptide 1 (RPS6KA1), ribosomal protein S6 kinase, 90kDa, polypeptide 3, ribosomal protein S6 kinase, 52kDa, polypeptide 1, ribosomal protein S6 kinase, 90kDa, polypeptide 6, ribosomal protein S6 kinase-like 1, KIAA1328, KIAA1324, and eukaryotic translation initiation factor 4E The comparative transcriptome study between slow growth improved Aseel and fast growth control broiler revealed the DEGs and their signicantly enriched pathways in slow growth improved Aseel, which inferred that they play an important role in regulating the growth and development of improved Aseel. The transcriptome data provides a theoretical basis for improving the performance of the slow growth improved Aseel as well as how to control the growth performance in fast grown control broiler chicken and provides reference data for revealing the molecular mechanism of slow growth improved Aseel as well as fast growth control broiler chicken. In this study, the mechanistic picture of gene expression data shows the embryo development, muscle development, egg production, plumage development, and energy production in improved Aseel, which would be fostered by a combination of (a) Differential regulation of MSTN, activin-like kinases and up-regulation of SMADs expect SMAD7 in the myostatin signaling pathway combined with down-regulation of caveolin’s (CAV1, CAV2, and CAV3) and differential regulation of insulin-like growth factor binding proteins, (b) Up-regulation of HSP70, NCF1, Map2k2, and down-regulation of MYOD1, and MYOZ2, (c) Up-regulation of fatty acid synthesis and β-oxidation genes (ACACA, ACACB, FASN, and CPT1), (d) Differential MAPK signaling pathway genes (MAP2K2, MKA, NES, SAMSN1, SOS2, and TAB2) (e) Differential regulation of Jak-STAT, mTOR, and TGF-β signaling pathway genes (IGF1, IGF2, IRS1, IRS2, PI3K, Akt1, Akt2, FoxO1, FoxO3, TSC22D1, TSC22D2, RHOA, RHOB, RHOC, RHOF, RHOQ and EIF4EBP1), (f) Differential regulation of mitochondrial genes (ND1, ND2, ND3, ND4, ND4L, ND5, ND6, ATP6, ATP1A2, COX1, COX2, COX3, CYTB

The differentially expressed genes between the 7 th and 18 th -day of the embryo were subjected to hierarchical cluster analysis using the Cluster 3.0 program [30]. We imported the matrix with as many columns as stages and many rows as genes, where each cell contains the log2 transformed fold change value for the gene and individual into the Cluster 3.0 program, normalizing on rows. We applied both rows and columns in hierarchical clustering and exported the resulting dendrogram as an image le.

Functional Characterization
Functional annotation, classi cation, and annotation clustering of selected gene sets were carried out by DAVID Tools 6.7 using biological processes, molecular function gene ontology categories, and KEGG pathways [31,32]. A threshold for signi cance of P<0.05 was considered to choose the signi cant functional categories.

Pathway analyses
Ingenuity Pathway Analysis (IPA; Qiagen, Valencia, CA; http://www.ingenuity.com) software was used for functional annotation, canonical pathway analysis, upstream analysis, and network discovery. The chicken DEGs data sets functionalities are primarily based on mammalian biological mechanisms because IPA is based on bioinformatics in humans. We have attempted to draw possible conclusions based on avian-based literature but biomedical research biases the functional annotations towards human disease Selection of candidate reference genes A total of 24 candidate reference genes were chosen based on their previous use/study in chicken or other avian species and the sequences were downloaded from NCBI (https://www.ncbi.nlm.nih.gov/) and CDS (coding DNA sequence) region was identi ed by using ExPASy translation tool (https://web.expasy.org/translate/) ( Table 1).

Real-Time Quantitative PCR analysis
Microarray expression data was validated using two-step Real-Time Quantitative PCR (qPCR) for speci c con rmation of the differentially expressed genes. First-strand cDNA was synthesized from 2 µg of total RNA using Thermo Scienti c Verso cDNA Synthesis Kit (Thermo Scienti c, USA). Gene specific qPCR primers were designed for 24 housekeeping genes and 83 DEGs using PrimerQuest software (http://eu.idtdna.com; Table S1; Table 2). The qPCR was performed using the BrightGreen 2X qPCR MasterMix-No Dye (Applied Biological Materials Inc., Canada) in Insta Q96™ Real-Time Machine (Himedia Laboratories, India) detection system. The PCR was performed under the following program; 5 min at 95 0 C, followed by 40 cycles of amplification with 15 s of denaturation at 95 0 C, 60 s of annealing/extension at 60 0 C. Three biological replicates were used. To check the specificity of the amplified products melt curve analysis was performed. The 2 -ΔΔCt calculated the relative expression level of each gene, and Transferrin (TFRC; Accession No: X55348.1) from G. gallus was used as a housekeeping gene to normalize the amount of template cDNA added in each reaction.

Statistical Analysis
The stability level of the 18 candidate reference genes from 7 th and 18 th day embryos of control broiler and improved Aseel was determined by using ve statistical algorithms geNorm, NormFinder, BestKeeper, Delta CT and RefFinder [33,34,35,36,37]. Comparison of mean expression values for qRT-PCR between the control broiler and PD4 improved Aseel groups were made using student's t-test and P ≤ 0.05 were considered as statistically signi cant.

Results And Discussion
The understanding of the gene functions, rapid and reliable quanti cation of gene expression, currently quantitative real-time PCR (qPCR) is one of the most sensitive techniques widely used. In quantitative gene expression, normalization is a key factor for accuracy and reliability, for that endogenous reference gene were used. According to the MIQE guidelines, the selection of suitable reference genes may vary for different species, different varieties, experimental conditions, and tissues and have to be validated before gene expression study [38]. Previous and recent studies also described different expression patterns of reference genes and focused on the validation of reference genes applied on different avian tissues [39,40,41]. However, so far validation of genes for their stable expression pattern in different embryo tissues such as 7 th and 18 th -day embryos of control broiler and improved Aseel was not performed.
Source, selection, primer design, veri cation of candidate reference genes In the present study, to identify the suitable reference genes for 7 th and 18 th -day embryos of control broiler and improved Aseel, 24 candidate reference genes with a wide range of biological functions were selected based on previous studies of various avian and non-avian species. They are 18SrRNA, ALB, B2MG, β-ACT, EEF1A1, GAPDH, GUSB, HMBS, HSP10, HSP70, L-LDBC, MRPS27, MRPS30, PGK2, PPP2CB, RPL5, RPL13, RPL14, RPL19, RPL23, SDHA, TBP, TFRC, and DNAJC24. The chicken orthologous genes were obtained from NCBI, CDS region was found and ampli ed with gene-speci c primers (Table S1). For all the primer pairs, melting curve analysis was performed to con rm the speci c ampli cation for each reference gene, which showed presence of single peak indicating speci city of the ampli cation of the genes.
Expression stability and ranking of candidate reference genes Expression levels of all candidate reference genes were measured in the samples collected from 7 th and 18 th -day embryos of the control broiler and improved Aseel. We observed that not all selected reference genes were expressing uniformly across different embryos of control broiler and improved Aseel (Fig. 3).
Differential transcription abundance levels were observed among all 18 genes. The mean Ct values of the selected reference gene ranged from 12.34 to 28.27 irrespective of different embryos. 18S rRNA and DNAJC24 genes showed the most (Ct=12.34) and the least (Ct=33.88) abundant transcripts respectively. Each reference gene had different expression ranges across all sample sets. The results show that the GUSB reference gene had the least variation of expression with mean Ct values ranging from 17.78 to 22.52 whereas the RPL5 gene showed a much higher expression variation with mean Ct values ranging from 13.07 to 32.89 across all sample sets (Fig. 3). It is important to note that there was a wide range of variation among selected reference genes, and it shows that not a single reference gene expresses constantly across the 7 th and 18 th day of embryos of control broiler and improved Aseel in the present study (Fig. 3). Therefore, it is of pivotal importance to choose the most reliable reference gene for expression pro ling of gene/s in different embryos of control broiler and improved Aseel.
For accurate gene expression, to select the best and most reliable reference gene and rank all the candidate reference genes according to their stability value, the most commonly used statistical programs were used i.e. geNorm, NormFinder, BestKeeper, Delta CT, and RefFinder [39,40,41,42,43,44,45,46,47,48]. These algorithms showed some differences in the stability ranking of stable reference genes, which may be due to the differences in each statistical program (Table  1). It is important to rank the stabilities of 18 genes and to con rm the number of reference genes necessary for accurate gene expression pro ling under embryo initiation and developmental stages. The variation among the reference genes determined by geNorm is stability measure (M value) and pairwise comparison expression ratio and provides optima number of genes in a given experiment [49,50]. According to geNorm stability criteria, the most stable genes in the 7 th and 18 th embryos of the control broiler and improved Aseel were PGK2 (0.38) and TFRC (0.38). NormFinder measures the reference gene stability by overall expression variation and across samples variation to reduce sensitivity towards co-regulation [49,51]. According to NormFinder stability criteria, the most stable genes in the 7 th and 18 th embryos of the control broiler and improved Aseel were GAPDH (0.371) and TFRC (0.4). BestKeeper calculates the gene expression variation based on Ct values and also calculates the Pearson correlation coe cient by pairwise correlation analysis for all reference genes and found the stable genes. According to BestKeeper stability criteria, the most stable genes in the 7 th and 18 th embryos of the control broiler and improved Aseel were HSP70 (1.08) and PPP2CB (1.11). The Delta CT results supported geNorm and NormFinder results and it showed the best stable genes in 7 th and 18 th embryos of control broiler and improved Aseel i.e. TFRC (2.1) and PGK2 (2.13). RefFinder conclusive the calculations using the all above-mentioned algorithms and suggested the stable genes in 7 th and 18 th embryos of control broiler and improved Aseel i.e. TFRC (2.21) and PGK2 (2.85).
To overcome different software program limitations, the stability of candidate reference genes was determined based on the consensus ranking for gene expression normalization in 7 th and 18 th -day embryos of control broiler and improved Aseel. Our study identi ed the most stable genes and indicated that, TFRC and PGK2 for 7 th and 18 th -day embryos of control broiler and improved Aseel. Our observations further strengthen the necessity to analyze the stability of candidate reference genes as suitable references.

Differentially expressed transcripts during embryo development stages
To study the effect of muscle development, genome-wide expression analysis was carried out at muscle initiation (7 th -day embryo) and muscle development (18 th -day thigh muscle) stages (Fig. 1a, b). Labeled RNA was hybridized to Affymetrix GeneChip™ Chicken Genome Array. After statistical data analysis, transcripts with FDR adjusted P-value ≤ 0.01 and fold change ≥ 1 were considered as significantly DETs (Fig. 1c). The complete list of the DETs in improved Aseel during embryo development stages as compared to their respective control broiler samples is presented in Supporting Document 1 & 2. In total 19,022 transcripts (16,252 unigenes) which accounted for approximately 58% of the total transcripts present on the GeneChip™ Chicken Genome Array were showed differential expression in improved Aseel at various stages analyzed. The maximum number of transcripts (15,380,46.9% of total DETs) showed differential expression on the 7 th day of the embryo and the least number of transcripts (3642, 11.11% of total DETs) showed differential expression on the 18 th day of the embryo. Commonly up-and down-regulated muscle responsive transcripts were identified among the embryo development stages to find out the degree of overlap (Fig. 1d, e). The maximum number of commonly up-regulated (9582) and down-regulated transcripts (4874) was observed in a 7 th -day embryo. Less number of up-regulated (1557) and down-regulated (1157) transcripts were commonly differentially expressed on the 18 th -day of the thigh muscle. Five commonly differentially expressed (545 up-regulated and 381 down-regulated) transcripts were identified among the embryo development stages, respectively (Fig. 1d, e; Supporting Document1&2).

Cluster analysis of differentially expressed transcripts
To profile the gene expression patterns in response to muscle slow growth during embryo development, the 19022 DETs were classified using hierarchical clustering. The expression patterns were separated into eight major clusters (I-VIII) based on tree branching (Fig. 2). Transcripts present in each stage within each cluster are presented in Supplementary Table S1. Among the eight major clusters, upregulated transcripts were enriched in cluster V, VI and VII and downregulated transcripts were enriched in cluster II, III and IV.

Expression of muscle related genes
The difference in fast and slow growth is generally dependent on the combination of environment and genes. The embryos collected in this experiment are under the same growth environment. In this experiment, the DEGs related to the main cause of growth and development differences mainly include muscle system process, muscle tissue morphogenesis and muscle organ morphogenesis, etc. (Fig. 4). The genes enriched by these entries are mainly muscle-related genes such as TNNC1, TNNT2, MYL3, MYH7, and FBXO32. The contraction of skeletal muscle-related genes is TNNC1, TNNT2, MYL3, and MYH7. In animals, skeletal muscle contraction affects the growth traits like muscle protein, muscle bre diameter, and muscle bre density [52,53]. In the corpus callosum, muscle is the most important component and the sarcomere is composed of thick and thin muscles. Myosin molecules are mainly involved in thick lament formation and the thin lament is a complex containing troponin, actin, and tropomyosin [54]. The skeletal muscle contraction is mainly based on the relative sliding of thick and thin laments [55]. For controlling the poultry muscle, according to previous reports, the troponin family gene is an important candidate gene [52,56,57]. In this experiment, the DEGs obtained related to the troponin family and important components of skeletal muscle are troponin T type 3 (TNNT3), slow muscle troponin T (TNNT1), and troponin I type 1 (TNNI1). The myosin superfamily is a highly conserved family of proteins and is composed of a heavy chain, an alkaline light chain, and a light chain and that is widely present in eukaryotic cells [58]. In striated muscle, smooth muscle, and nonmuscle, the myosin is involved in myo laments as a class II. The MYH gene family is a key subunit in the myosin class II molecule and encoded the myosin heavy chain (MYHC) [59]. The studies revealed that MYH7 gene mutations could cause skeletal muscle disease or skeletal muscle disease with cardiomyopathy [60,61]. In the chicken breast muscle, the MYH7 gene expression is high and it is a hypothesis that it has a regulatory role in muscle tissue growth and development [62]. MYL3 gene is primarily expressed in slow muscle and it is a member of the myosin light chain (MYL) gene family [63]. Inhibition of myosin light chain gene expression in cell lines showed myoblast proliferation [64]. The chicken embryonic leg muscles proteome analysis found that MYL3 protein is closely related to the regulation of muscle contraction and the long non-coding RNAs Inc00037615 and Inc00037619 together regulate the muscle development [65]. In this experiment, the DEGs obtained related to the myosin family are MYH7 and MYL3. In this study, ve genes related to muscle contraction were screened i.e. TNNT3, TNNT1, TNNI1, MYL3, and MYH7. Transcriptome study showed that the expression of these ve genes is high in improved Aseel slow growth compared with control broiler fast growth. Also, these genes simultaneously present in two signi cantly enriched pathways i.e. adrenergic signaling in cardiomyocyte and cardiac muscle contraction. Maybe due to the upregulation of these genes and down-regulation of the Adrenergic signaling in cardiomyocyte and cardiac muscle contraction, improved Aseel will grow slow compare to the control broiler. Previous studies showed that in animals under fasting, FBXO32 gene expression was signi cantly increased and muscles were degraded due to lack of food, and maybe it is associated with muscle atrophy [66,67]. The FBXO32 gene expression was found in chicken's leg muscle, heart, and chest muscle and plays an important role in the 7 th -week growth of chickens [68]. In this experiment, the results of transcriptome data showed that the expression of FBXO32 and FBXO7 genes were signi cantly lower in slow-growth improved Aseel than that in the fast-growth control broiler, which was consistent with previous studies. The DEGs KEGG pathway enrichment analysis revealed three signi cantly enriched adrenergic signaling in cardiomyocytes, Cardiac muscle contraction, and tight junction signaling pathways. In cell junctions, the tight junction is an important component and it acts as a barrier for cells to pass ions and molecules, plasma membrane apical movement regulation, and basal proteins and lipids [69,70]. In all eukaryotic cells, tight junction recognized as the actin cytoskeleton and it is involved in cell division, adhesion, movement, and phagocytosis and also found in chicken leg muscle transcriptome [71,72]. The muscle growth epigenetic transcriptional regulators are differentially regulated in 7 th such as protein arginine N-methyl transferase family (PRMT1, and PRMT3), histone lysine N-methyl transferases (EHMT1, and SETDB1), SWI/SNF chromatin-remodeling enzymes (SmarcB1, and SmarcaA4).

Muscle Development and Myostatin Signaling
The muscle development and differentiation-related genes like MYOD1, MYF6, MYF5, Myoz2, MAP2K6, MAP3K7, CAV1, CAV2, CAV3, HSP70, and NCF2 were differentially regulated in slow-growing improved Aseel then compared to fast-growing control broiler (Fig. 4). In muscle differentiation, MYOD1 promotes the muscle-speci c gene expression and function together with MYF5 and MYOG [73]. MYOD1 combined with transient placeholder protein, prevents the binding of other transcription factors to DNA, and retains the inactive conformation of the DNA [74]. One of the key functions of MYOD is to stop the differentiated myocytes proliferation by enhancing the transcription of p21 and myogenin to remove cells from the cell cycle [75]. Altogether, up-regulation of MYOD1 involved in skeletal muscle phenotype establishment by regulation of precursor cell proliferation and promoting irreversible cell cycle arrest, facilitate differentiation and sarcomere assembly by activation of sarcomeric, and muscle-speci c genes [76]. In this study, transcriptome data shows down-regulation of MYOD1 in the 7 th day improved Aseel embryo, maybe due to this reason the muscle-speci c gene proliferation is slow in improved Aseel. The myozenin is a α-actinin-and γ-lamin-binding protein of Z line skeletal muscle binds to calcineurin, and involved in skeletal muscle myocyte differentiation [77,78]. Therefore, the up-regulation of MYOZ2 and MYOZ3 genes in muscle tissues suggested that they are involved in muscle growth and development and directly in uence meat quality. Besides, MYOZ1 and MYOZ2 genes expressed in mice showed a signi cant reduction of calcineurin gene expression [79,80,81]. Our transcriptome study shows less expression of the MYOZ2 gene and high expression of regulator of calcineurin 1 gene in 7 th day improved Aseel embryo, maybe this is the region the muscle development is slow in improved Aseel. Mitogen-activated protein kinases are major components of pathways controlling embryogenesis, cell differentiation, cell proliferation, cell death, muscle development, and response to environmental stress ( Fig. 8) [82,83,84,85]. The MAP2K6 and MAP3K7 activate MAP kinase and nuclear factor-kappa β (NFκB) and play an important role in its signal transduction pathway, respectively [86,87]. In a proteomic study, predicted the expression of MAP2K1, MAP2K2, and MAP2K4, MAP4K4 gene may inhibit the low feed e ciency phenotype [88]. In our transcriptome data, the mitogen-activated protein kinase family genes are up-regulated in the 7 th -day embryo of improved Aseel. Our results correlating with these results and maybe in this region the improved Aseel has less feed e ciency and slow muscle growth.
In this study, the expression of caveolin family genes like CAV1, CAV2, and CAV3 was downregulated in the 7 th -day embryo of improved Aseel. In cell signaling, the caveolin genes act as sub-cellular structures by assisting the attribute of hormonal signals after binding hormone to the target receptor on the cell surface. The CAV3 is acting as a muscle-speci c isoform for the caveolin protein and mutations or different expressions of CAV3 can result in muscle myopathies [89,90]. In pigs, the CAV3 expression was up-regulated during muscle hyperplasia and it may use as a genetic marker for meat production in swine [91]. In mouse, the high or low expression of CAV3 made muscle cells more susceptible to oxidative stress and reduced survival through PI(3)K/Akt signaling [92]. In the low FE PedM broiler phenotype, the higher expression of CAV3 contributed to higher oxidative stress and enhance muscle development [93]. In high FE breast muscle, the CAV1 protein is involved in insulin signaling [88]. Based on previous reports, down-regulation of caveolin family genes in the 7 th -day embryo of improved Aseel, maybe the improved Aseel muscle development will be reduced. In high FE breast muscle, the up-regulation of HSP70 maintains muscle bre integrity and enhances muscle regeneration and recovery from damage [94]. HSP70 is also responsible for correct folding and assembly of nuclear-encoded proteins an important chaperone for mitochondrial DNA-encoded proteins as components of the mitochondrial electron transport chain targeted for import into the mitochondria [95,96]. In low and high FE phenotype, higher expression of HSP90 and HSPB2 was in response to oxidative stress, respectively [3,88].
In our transcriptome study correlated with previous studies and we observed up-regulation of HSP70 and HSPB1 on the 7 th -day and down-regulation of HSP90 on the 18 th -day of improved Aseel embryo. Maybe the improved Aseel muscle is growing this region slowly compared to the control broiler. In NADPH oxidase 2 (NOX2), the NADPH/NADH oxidase is a critical component and encoded by neutrophil cytosolic factor 2 (NCF2) [97]. In muscle, superoxide was generated by NOX2 in the sarcoplasmic reticulum, and it is a major source of oxidative stress [97,98]. In neutrophils, NADPH generated superoxide during phagocytosis.
The nuclear factor erythroid 2-like 2 (NFE2L2) is a downstream target for NOX2 and activates genes that contain an antioxidant response element in their promoter regions [99,100]. In high FE animals, predicted the NFE2L2 expression should be upregulated [5,88]. In a high FE commercial broiler, the upregulation of NCF2 was associated with muscle remodeling and hypertrophy [5]. In our transcriptome data, the NCF1, NOX1, and NOX3 were up-regulated on the 7 th -day and NCF2 was down-regulated on the 18 th -day of improved Aseel embryo. Maybe due to the differential expression of these genes, the improved Aseel is more resistant to oxidative stress, low FE, and slow growth.
Myostatin is a member of tumor growth factor β (TGF-β) and is known as growth/differentiation factor 8 (GDF-8) [101,102]. In the myostatin (MSTN) signaling pathway, MSTN binds to its receptors ActIIA/ActIIB and activates ALK4 and ALK5 that phosphorylate Smad2/3 leads to its binding with Smad4 and translocation of the complex to the nucleus and where it blocks the transcription of genes responsible for the myogenesis [103,104,105,106,107] (Fig. 8). The myostatin is solely expressed in skeletal muscle during embryogenesis to control the differentiation and proliferation of the myoblasts [101]. However, in the adult stage, it is not only expressed in skeletal muscle but also expressed in other tissues like the heart, adipose tissue, mammary gland [101,108,109,110,111,112]. In turkey satellite cells, MSTN is a strong negative regulator for skeletal muscle growth, differentiation, and proliferation [101,113].
The relation between MSTN and growth performance studies in broilers shows that MSTN is a polymeric gene in which different alleles of the gene can affect performance [114,115,116]. In PedM broiler, the FE differences may be due to different haplotypes of the MSTN gene [107]. Myostatin knockdown by RNAi shows muscle growth enhancement in transgenic sheep and chicken [117,118,20]. In the present study, MSTN has differentially regulated in the 7 th -day improved Aseel embryo. Maybe differential regulation of myostatin is need for myoblasts differentiation and proliferation in initial embryogenesis. Follistatin (FSTN) regulates the MSTN by inhibiting or limit the MSTN activity. Follistatin-like 1 (FSTL1) is a glycoprotein and rich in cysteine (SPARC) family and comprises a secretion signal, a Follistatin-and a Kazal-like domain, two EF-hand domains, and a von Willebrand factor type C domain [119] (http://www.uniprot.org/uniprot/Q12841). In mouse, the FSTL1 were broadly expressed throughout the entire embryo and restricted most of the tissues at the end of gestation but in the adult mouse, it is highly expressed in heart, lung, and subcutaneous white adipose tissue [120,121]. In this present study, FSTL1 was up-regulated and follistatin/kazal down-regulated in the 7 th and 18 th -day of improved Aseel embryo, respectively. Initial up-regulation and later downregulation of FSTL1 may initiate muscle proliferation in the 7 th -day embryo and 18 th -day embryo slowdown the muscle development in improved Aseel. In humans, activin receptor type-1B (ACVR1B) or ALK4 is a protein that acts as a transducer of activin or activin-like ligand signals [122]. ACVR1B forms a complex with ACVR2A/ACVR2B and goes on to recruit the SMAD2/SMAD3 [123]. Besides, ACVR1B transduces nodal, GDF-1, and Vg1 signals combined with other coreceptor molecules like protein cripto [124]. Transforming Growth Factor-β (TGFβ) is a key player in cell proliferation, differentiation, and apoptosis and TGFβ receptors are single-pass serine/threonine kinase receptors and can be eminent by their structural and functional properties [125]. Due to the similar ligand-binding a nities, the transforming growth factor beta receptor I (TGFβR1)/ALK5 and TGFβR2 can be distinguished from each other by peptide mapping only. In mouse, the TGFβ1 mRNA/protein has been present in cartilage, endochondral, membrane bone, and skin and play a role in the growth and differentiation of these tissues [126] (Fig. 8). In the present study, activin A receptor type IB (ALK4) and transforming growth factor beta receptor II (TGFBR2) were up-regulated, and transforming growth factor beta receptor I (TGFBR1/ALK5) was down-regulated in the 7 th -day embryo of improved Aseel. The differential expression of ALK4 and ALK5 may control the myostatin signaling pathway. The SMADs are important for regulating cell development and growth and they have structurally similar proteins and main signal transducers for TGFβ receptors (Fig. 8). The eight SMAD genes are distributed into three sub-types of SMADs, they are R-SMADs, Co-SMADs, and I-SMADs [127,128]. The R-SMADs consist of Smad1, Smad2, Smad3, Smad5, and Smad8/9 and are involved in direct signaling from the TGFβ receptor [129,130]. The Co-SMADs consist of the only SMAD4 and working jointly with R-SMADs to recruit co-regulators to the complex [131]. R/Co-SMADs are primarily located in the cytoplasm, following TGFβ signaling later accumulate in the nucleus, where they can bind to DNA and regulate transcription. I-SMADs consist of SMAD6 and SMAD7 and are predominantly found in the nucleus, and they can act as direct transcriptional regulators. SMAD6 has speci cally associated BMP signaling and SMAD7 is a TGFβ signal inhibitor and suppresses the activity of R-SMADs [132,133,134]. In the present transcriptome study, the SMAD family member 1 (SMAD1), SMAD speci c E3 ubiquitin protein ligase 2, SMAD family member 3 (SMAD3), SMAD family member 5 (SMAD5), TGF-beta signal pathway antagonist Smad7 (SMAD7B) up-regulated in 7 th day of improved Aseel embryo and they may control the myostatin signaling pathway in improved Aseel embryonic stage. Summarizes the initial steps in the MSTN signaling pathway in the present study that would potentially exert a negative effect on muscle differentiation and proliferation in the slow-growing improved Aseel.

Energy Sensing
In humans and animals, the adenosine monophosphate-activated protein kinase (AMPK) gene regulates diverse biological functions [135]. The mammalian 5′ AMPK gene has two α subunits i.e. AMPKα1 and AMPKα2 that are encoded by Prkaa1 and Prkaa2 gene, respectively. The knockout mouse clearly demonstrated that the AMPKα2 controls the homeostasis in skeletal muscle [136,137]. Also observed reduction of bre numbers (~25%) and sizes (~20%) in the soleus muscle of AMPKα1 knockout mice [138]. However, in AMPKα2 knockout mice, both bre size and muscle mass were signi cantly increased, while the muscle bre number has remained similar to WT animals. The muscle mass reduced and increased differentially expressed alternative polyadenylation sites (DE-APSs) were down-regulated in AMPKα1 knockout mice, but up-regulated in AMPKα2 knockout mice, respectively [46]. The ve genes i.e. carbonic anhydrase 3 (Car3), myosin light chain kinase family, member 4 (Mylk4), nebulin (Neb), obscurin (Obscn), and phosphofructokinase, muscle (Pfm) are utilized different APSs showed potential effects on muscle function [46]. The high FE phenotype birds show up-regulation of both AMPKα1 and AMPKα2 [3]. In low energy level conditions the AMPK genes expression increases for ATP production by inhibiting the ATP consuming pathways like fatty acid synthesis, protein synthesis, and gluconeogenesis and stimulates the ATP producing pathways like mitochondrial biogenesis and oxidative phosphorylation, glycolysis, and lipolysis [139,140,141]. In the present study, the AMPK genes like 5'-AMP-activated protein kinase gamma-1 non-catalytic subunit variant 1 (PRKAG1), protein kinase cAMP-dependent regulatory type I alpha (tissue speci c extinguisher 1) (PRKAR1A), protein kinase AMP-activated beta 2 non-catalytic subunit (PRKAB2), protein kinase AMP-activated gamma 2 non-catalytic subunit (PRKAG2), carbonic anhydrase XIII, myosin light chain kinase (MYLK), atrial/embryonic alkali myosin light chain, are down-regulated in the 7 th -day of improved Aseel embryo, maybe due to the down-regulation of energyproducing pathways, the improved Aseel muscle will grow slowly compare to control broiler. Curiously, creatine kinase (muscle isoform, CK-M) and creatine kinase (brain isoform, CK-B) were up-regulated in high and low FE phenotype, respectively [88,142,143]. The reason for this discrepancy may be the high FE phenotype broiler breast muscle has enhanced capabilities for mitochondrial oxidative phosphorylation as well as creatine and phosphorylated creatine shuttle in and out of mitochondria [143]. In the present study, the creatine kinase muscle (CKM), creatine kinase mitochondrial 1A (CKMT1A), creatine kinase brain (CKB), creatine kinase, mitochondrial 2 (sarcomeric) (CKMT2) genes were up-regulated in the 7 th -day embryo of improved Aseel. In skeletal muscle, nitric oxide was synthesized by nitric oxide synthases and it is regulated key homeostatic mechanisms like mitochondrial bioenergetics, network remodeling, mitochondrial unfolded protein response (UPRmt), and autophagy ( Fig. 4; Fig. 5). In mice, nitric oxide synthase de ciency inhibits the Akt-mammalian target of the rapamycin pathway and dysregulated the Akt-FoxO3-mitochondrial E3 ubiquitin-protein ligase 1 (Mul-1) axis [144]. Thus, the mitochondrial biogenesis and body energy balance were controlled by the nitric oxide-cGMP-dependent pathway [145]. In detail, the inhibition of nNOS/NO/cGMP-dependent-protein kinases enhanced the FoxO3 transcriptional activity and increased the Mul-1 expression. The absence of the nitric oxide synthases signi cantly impaired muscle bre growth with muscle force, decreased resistance to fatigue, and degeneration/damage post-exercise. In our study, nitric oxide synthase 2 was up-regulated and, cGMP-dependent protein kinase type I, and FOXO1 was down-regulated in the 7th-day of improved Aseel embryo, maybe this is the region, the improved Aseel muscle strength was high and they are more energetic compared to the control broiler.
The comparative muscle transcriptome analysis between high and less pH chicken showed that most of the glycolysis pathway genes are up-regulated in the less pH chicken [146]. The previous study shows Aseel and broiler chicken's meat do not have any signi cant pH variation but the heavier bird's meat had a signi cantly higher pH [147]. In this study, glycolysis metabolism-related genes are differentially regulated on the 7 th -day of improved Aseel embryo. The glucose-6-phosphate isomerase, fructose bisphosphate aldolase, phosphoglycerate kinase, and enolase were down-regulated and GAPDH, phosphoglycerate mutase, and pyruvate kinase were up-regulated in the 7 th -day embryo of improved Aseel. Fructose bisphosphate aldolase is a key enzyme in glycolysis as well as gluconeogenesis and involved in the stress-response pathway during hypoxia [148]. The high pH chickens have increased oxidative stress, maybe the higher expression of fructose bisphosphate aldolase is linked to its function in the stress-response pathway rather than to its role in ATP biosynthesis [149]. In our study, down-regulation of fructose bisphosphate aldolase enhanced the ATP synthesis, maybe this is the region improved Aseel birds have more energy than control broiler. Noteworthy, the up-regulation of glycolysis pathway genes are increased the pyruvate levels and entering the citric acid cycle, and thus, higher levels of ATP are produced in improved Aseel. The protein phosphatase-1 regulatory subunit 3A (PPP1R3A) binds glycogen with high a nity, and activates glycogen synthase (GYS), and inhibits glycogen phosphorylase kinase (PHK) by dephosphorylation through the protein phosphatase-1 catalytic (PPP1C) subunit. In this study, the glycogen metabolism genes i.e. protein phosphatase-1 regulatory subunit 2 (PPP1R2), glycogenin 1, glycogen phosphorylase, and protein phosphatase-1 catalytic subunit beta (PPP1CB) were down-regulated in the 7 th -day embryo of improved Aseel, maybe this is the region that the improved Aseel muscle have more glycogen. The AMP-activated protein kinase (AMPK) complex is another key regulator of glycogen turnover and it consists of one α catalytic and two non-catalytic subunits, β, and γ. The β subunit binds to glycogen along with α and γ subunits and forms a heterotrimeric AMPK complex. In the muscle cell, the γ subunits of the AMPK complex act as energy sensors and binding to AMP and ATP [150]. In our study, AMP-activated protein kinase beta 2 non-catalytic subunit (PRKAB2), cAMP-dependent protein kinase regulatory type I alpha (tissue-speci c extinguisher 1) (PRKAR1A), AMP-activated protein kinase gamma 2 non-catalytic subunit (PRKAG2) were down-regulated, and 5'-AMP-activated protein kinase gamma-1 noncatalytic subunit variant 1 (PRKAG1), AMP-activated protein kinase gamma 3 non-catalytic subunit (PRKAG3), AMP-activated protein kinase alpha 2 catalytic subunit (PRKAA2) were up-regulated. The down-regulation of β subunits and up-regulation of α and γ subunits may balance the glycogen accumulation and increase the ATP molecules for energy production in improved Aseel muscle, maybe this is the region the improved Aseel has more strong than the control broiler. Apart from these, several other genes indirectly in uence glycogen storage in muscle. The phosphodiesterase 3B (PDE3B) gene activated by insulin and induce antiglycogenolytic effects and the mitochondrial creatine kinase (CKMT2) transfers the high-energy phosphate from mitochondria to creatine. In our study, phosphodiesterase 3A (PDE3A), phosphodiesterase 4D, phosphodiesterase 8A, mitochondrial creatine kinase 2 (CKMT2), mitochondrial creatine kinase 1A (CKMT1A), mitochondrial creatine muscle (CKM) were up-regulated in the 7 th -day of embryo improved Aseel, maybe this is the region, the higher expression of these genes in muscle, improved Aseel birds are more energetic compared to control broiler. To produce energy and compensate for the lack of energy due to carbohydrates and glycolysis, the high pH chickens muscle asks for more intense oxidative pathways, such as lipid β-oxidation and ketogenic amino acid degradation [149]. In the high pH muscle line, the 3-hydroxymethyl-3-methylglutaryl-CoA lyase (HMGCL) catalysis the nal step of leucine metabolism and ketone metabolism, acetyl-CoA acetyltransferase-2 (ACAT2) involved in β-oxidation or degradation of ketogenic amino acids, and nudix hydrolases (NUDT7, NUDT12, NUDT19) hydrolyse a nucleoside di and triphosphates, dinucleoside and diphosphoinositol polyphosphates, nucleotide sugars and RNA caps, were up-regulated. In the present study, 3-hydroxymethyl-3-methylglutaryl-CoA lyase like-1 (HMGCLL1), acetyl-CoA acetyltransferase 2 (ACAT2), nudix type motif 7 (NUDT7), and nudix type motif 21 (NUDT21) were down-regulated and carnitine/palmitoyl-transferase 1 (CPT1) were up-regulated in the 7 th -day embryo of improved Aseel (Fig. 6). They may regulate the β-oxidation in peroxisomes as well as mitochondria, this is the region may be fats required for initial embryo development and excess fats involved in β-oxidation and nally provide the energy for embryo development.

Protein Synthesis
To promote cell growth, the mTORC1 complex increases protein synthesis and lipid metabolism, and autophagy inhibition, and regulates the transcription of several genes [151]. In high FE phenotype, the cDNA microarray data shows higher expression of mTORC1 [152,153]. The mTORC1 complex has two major components i.e. mTOR and RAPTOR [154]. The RAPTOR and mTOR were up and down-regulated in high and low FE birds, respectively, and the up-regulation of RAPTOR could be a positive effect on protein synthesis [154]. In low FE phenotype, the PRKAR1A and GLUT-8 were up-regulated. The p70S6k and eukaryotic translation initiation factor 4E (EIF4E) are the key downstream targets for mTORC1 and involved in enhancing the protein synthesis. In low FE birds, the expression o p70S6k was higher [155]. The muscle tissue of RNAseq transcriptomic data showed higher expression of eukaryotic initiation, elongation, and translation genes in high FE compared to low FE PedM phenotype [142,156]. In the present study, the late endosomal/lysosomal adaptor MAPK and MTOR activator 3 (LAMTOR3), protein kinase, cAMP-dependent, regulatory, type I, alpha (tissue speci c extinguisher 1) (PRKAR1A), solute carrier family 2 (facilitated glucose transporter) member 8 (SLC2A8)/ GLUT8/ GLUTX1, solute carrier family 2 (facilitated glucose transporter), member 3 (SLC2A3), ribosomal protein S6 kinase, 90kDa, polypeptide 1 (RPS6KA1), ribosomal protein S6 kinase, 90kDa, polypeptide 3, ribosomal protein S6 kinase, 52kDa, polypeptide 1, ribosomal protein S6 kinase, 90kDa, polypeptide 6, ribosomal protein S6 kinase-like 1, KIAA1328, KIAA1324, and eukaryotic translation initiation factor 4E binding protein 1 gene were down-regulated in the 7 th -day embryo of improved Aseel compared to control broiler (Fig. 4). May be due to the down-regulation of mTORC1 complex and ribosomal machinery genes, the protein synthesis for muscle growth is less in improved Aseel compared to control broiler.

Insulin Signaling
In chicken, the SHC1 is only activated by nutritional changes and suggesting that insulin signaling in chickens have a tissue-speci c manner [157,158]. When insulin binds to the insulin receptor, both IRS-1 and SHC1 are activated by a phosphoinositide-3 kinase (PI3K) mediated tyrosine phosphorylation activity [3]. The skeletal myoblast is mainly differentiated by two key modulators i.e. insulin-like growth factor 1 (IGF1) and broblast growth factor 2 (FGF2) [159,160]. In L6 and C2C12 myoblasts, a high concentration of IGFs reduces their differentiation, whereas a low concentration enhances their differentiation [161,162]. The IGF binding proteins (IGFBP-1 to IGFBP-6) have highly conserved regions and binds with high a nity to IGF-1 and IGF-2. In the extracellular matrix, IGFBP3 may regulate the interaction of IGFs and it is present in rat soleus muscle (type I muscle bre) [163,164]. In humans, IGFBP3 was also playing a role in myoblasts differentiation [165]. The broblast growth factor 2 (FGF2), transforming growth factor beta (TGFb), and oncogenic Ras also inhibit the skeletal myoblast differentiation [160]. In 23A2 myoblast cell lines, showing the inhibition of 23A2 myoblasts differentiation by IGF1 and FGF2 by stimulating the signaling through mitogen-activated protein kinase (MAPK) kinase (MEK) to MAPK [166]. In our study IGF2, FGF2, IGFBP2, IGF2BP1, IRS-1, IRS-2, INSIG1, PIK3CA, PIK3CD, TGFb, and Ras oncogene family genes are up-regulated in the 7 th -day embryo of improved Aseel (Fig. 7). Maybe this is also one of the regions for improved Aseel muscle is differentiated slowly when compared to the control broiler.

Expression of plumage development genes
In vertebrate coloration, melanin pigmentation is an important component and regulated by strong genetic control [167,168]. In chicken, the plumage coloration development is extremely complex and can be categorised as structural or pigment-based [169]. For animal colouration, melanin is a common component, synthesized in melanocytes, and deposited in various organs as granules [170]. Different pigment patterns are formed based on the presence of melanocytes modulating, arrangement or differentiation, and associated with a series of functional genes [171,172]. The melanogenesis genes such as HOX, CHAC1, GPX3, BMP5, PITX2, RGN, MITF, TYR, KIT, OCA2, ASIP, MCIR, KITLG, IRF4, SLC24A4, SLC45A2, EDN, TYRP1, and TYRP2 are involved in melanin pigmentation [173,174,175,176,177]. The homeobox (HOX) genes are transcription factors and involved in skin appendages development including hair follicles and feathers [178,179,180,181,182,183] (Fig. 8). In black-bone chickens, four HOX genes i.e. HOXB9, HOXC8, HOXA9, and HOXC9 were identi ed for melanin pigmentation [172]. The Wnt signaling is essential for skin organogenesis and its appendages like hairs, feathers, and scales, melanocyte development, and differentiation [184,185,186,187]. The HOXB9 is identi ed as a target gene for Wnt signaling and HOXC8 is expressed in the rst stage of feather morphogenesis like dorsal dermal and epidermal cells [188,189]. In this study , HOXA2, HOXA9, HOXB3, HOXB5, HOXB7, HOXB8, HOXB9, HOXC11,   HOXD1, and HOXD3 are up-regulated in the 7 th -day improved Aseel embryo, maybe this is the region the improved Aseel plumage has multiple colors.
In animals, melanogenesis is regulated by GSH and it is closely associated with melanin deposition in the skin of humans and other mammals [190,191,192,193,194]. The low and high level of GSH indicates eumelanin-type pigmentation and phaeomelanin-producing melanocytes found in skin, respectively [195]. Two feather melanin pigmentation genes were identi ed in black-bone chickens such as ChaC glutathione-speci c gammaglutamylcyclotransferase 1 (CHAC1) and glutathione peroxidase 3 (GPX3) [172]. The CHAC1 cleavage GSH into 5-oxoproline and Cys-Gly dipeptide and GSH over-expressed mammalian cells caused GSH depletion [196,197]. Hence, CHAC1 expression is associated with GSH metabolism and plays an important role in the melanogenesis process. In eumelanin and phaeomelanin synthesis, the hydrogen peroxide reduced by GSH-dependent peroxidase and GPX3 belongs to the GSH peroxidase family and catalyzes the GSH to glutathione disulphide (GSSG) [194,198]. The melanoma cells pigmentation regulated by GSH levels, glutathione peroxidase, and glutathione reductase suggesting that GSH mediated redox process plays an important role in melanogenesis regulation [199].
Hence, the expression of GPX3 plays an active role in chicken feather melanogenesis. In this study, CHAC1, CHAC2, gamma-glutamylcyclotransferase (GGCT), and GPX8 were down-regulated on the 7 th -day of improved Aseel embryo, maybe this is the region the improved Aseel has multiple colors on their feathers. In black-bone chicken, two pathways i.e. TGF-β signaling pathway, and ascorbate and aldarate metabolism were identi ed for plumage melanogenesis [172].
The TGF-β involved in the regulation of chicken retinal pigment epithelial cell proliferation and melanin synthesis [200]. The BMP5 and PITX2 genes are involved in the TGF-β signaling pathway and play a role in chicken melanin synthesis, and the BMP5 and PITX2 were highly expressed in white feather and black feather bulbs, respectively [200]. The BMP3 gene was highly expressed in embryonic and post-embryonic stages of control layer when compare to broiler chicken and BMP4 gene was differentially expressed in juvenile stages of broiler and layer chicken, respectively [201,202]. The regucalcin (RGN) is a calciumbinding protein involved in the ascorbate and aldarate metabolism pathway and plays a crucial role in intracellular calcium homeostasis maintains [203]. In this study, transforming growth factor beta receptor II (TGFBR2), BMP1, BMP1A, BMP4, BMP7, BMPR1A, BMPR2, PITX3 up-regulated and RGN was downregulated in the 7 th -day embryo of improved Aseel, maybe this is the region the improved Aseel plumage has multiple colors.
In melanin synthesis, TYR is a rate-limiting enzyme and involved in different oxidative steps [204,205]. In black vs. white skin chicken, the TYR is highly expressed and it is consistent with sheep coat color studies [173,206,207]. In black-coated vs. white-coated sheep, the TYRP1 gene was highly expressed [207]. KIT is a receptor tyrosine kinase, the mutation in KIT showed piebaldism and auburn hair color in humans and it plays an important role in UVB-induced melanogenesis in the epidermis, and inhibition of KIT expression may result in the lighten human skin color [208,209]. In black skin chicken, KIT is highly expressed and black skin color is due to increased melanin compared to the white skin color [173]. In melanocytes development, microphthalmia-associated transcription factor (MITF) is playing a role and mutations in the MITF gene is responsible for Japanese quail and chicken plumage color and it is supported to alternative splicing of MITF gene in the skin of sheep [210,211]. In ducks, the TYR and MITF expression may involve in the formation of black and white plumage [212]. Melanocortin-1 receptor (MC1R) is binding to melanocyte stimulating hormone (MSH) and initiates the melanogenesis cascade and regulates the mammalian skin pigmentation and hair color [213,214,215,216]. The agouti signaling protein (ASIP) is responsible for the skin color of both white and black-coated sheep and mutations in ASIP cause black and tan pigmentation phenotype in pigs [206,217]. The ASIP binds to MC1R and reduced the MITF and TYR gene expression and nally, the phaeomelanin would be reduced in epidermal tissues. In black-skinned chickens, the expression of ASIP is higher than compared to the white-skinned chicken and it can suppress the MC1R gene expression in black-skinned chickens [173]. The Oculocutaneous albinism type 2 (OCA2) is a common skin pigmentation disorder caused by a mutation in the OCA gene. In black chickens, the OCA2 was up-regulated and it may be related to black skin color [173]. In chickens, the endothelins (EDN1, EDN2, and EDN3) and their receptors (EDNRA, EDNRB, and EDNRB2) are involved in the regulation of pigmentation and plumage [218]. The expression of EDNRB2 was signi cantly different between adult black and non-black chicken [219]. In this study, TYRP1, KITLG, MITF, MC1R, AGRP, YRK, and P56LCK were up-regulated and EDNRA and EDN1 were down-regulated in the 7 th -day of improved Aseel embryo and this is the region the improved Aseel plumage has multiple colors.

Expression of genes related to egg production
In chicken, the reproductive system is regulated by hypothalamic-pituitary-ovarian (HPO) axis hormones, while ovulation, the GnRH-I triggers the pituitary gland to release FSH and LH, and stimulates the secretion of estradiol and progesterone in the ovary [220]. Several reproductive hormone regulation genes were identi ed between high and low egg production chickens, such as HADH, HMGCR, RAB11FIP1, and FAM3D [221]. In the pituitary gland, lipid metabolic process, prolactin and, MAPK signaling pathway genes i.e. HMGCR, HMGCS1, NFKB1, VAV3, SOS1, IL1R1, MEF2C, STK3 were highly expressed. In chicken, the anterior pituitary gland synthesized and released prolactin and involved in reproduction, laying eggs, metabolism, development, and hypothalamic-pituitary-gonadal axis regulation [222,223]. In chicken, the HMGCR gene variants (G-789-A, C-937-G, and A-2316-C) and high and low concentrations of VLDL showed higher and lower egg production, respectively [224]. In laying chickens, the APOB is a primary organizing protein for chylomicrons and VLDL and responsible for the transport of lipoprotein and circulating in the plasma, and stored into the oocytes to form an egg yolk [225,226]. In our study, GNRHR, HADHB, HMGCS1, HMGCS2, RAB11FIP2, RAB11FIP3, RAB11FIP4, NFKB2, VAV2, SOS2, MEF2D, STK3, PRL, PRLR, PRLH, and PRLHR2 genes were up-regulated, and FSHR, VAV3, IL1RL1, and IL1RAPL2 genes are down-regulated, and family with sequence similarity genes and apolipoprotein B were differentially regulated in the 7 th -day embryo of improved Aseel. Maybe this is the region the improved Aseel has less egg production than the commercial chicken. In avian species, the genes SPP1, BPIFB3, and EDIL3 are mainly involved in egg and oviposition, development of reproduction system, and vesicle-mediated eggshell calcication, respectively [227,228,229,230]. In this study, secreted phosphoprotein 1 (SPP1) and secreted phosphoprotein 2 (SPP2) genes are down-regulated, and EDIL3 is up-regulated in the 7 th -day of improved Aseel embryo, due to this, egg and oviposition are less and eggshell calcication is more in improved Aseel. In nandanyao chicken, FN1, FGF7, SOX2, ALDOB, HSPA2 genes are expressed in the ovary, and UQCRH, COX5A, FN1, TGFB, ACTN1 genes are expressed in the uterus and involved in egg production [231]. In this study, FN1, FNDC3A, FGFR1, FGFR3, FRS3, FRS2, FGFR2, FGFRL1, FGF8, FGF18, FGF3, FGF12, SOX2, SOX3, SOX4, SOX5, SOX7, SOX8, SOX9, SOX11, SOX17, HSPA2, HSPA4, COX1, COX2, COX3, TGFB4, and ACTN1 were up-regulated and ALDOB, HSPA5, HSPA8, HSPA9, HSP12A, UQCRFS1, UQCRB, were down-regulated in the 7 th -day improved Aseel embryo. The differential expression of the ovary and uterus related genes are differentially expressed in the 7 th -day of improved Aseel embryo, due to this region the egg production is less in improved Aseel. The DEGs related to the pituitary gland between high and low egg production chickens are mainly involved in mTOR and Jak-STAT signaling pathways [221]. In mouse, the mTOR signaling pathway will regulate the granulosa cell proliferation and differentiation [232]. In this study, the mTOR and Jak-STAT signaling pathways were upregulated in the 7 th -day improved Aseel embryo.
In stressful conditions, peripheral and brain tryptophan levels can alter by stimulating the immune system and activating the hypothalamic-pituitary-adrenal axis [233,234]. In this study, tryptophan metabolism was down-regulated in the 7 th -day of the embryo and up-regulated in the 18 th -day thigh muscle, this is the region the improved Aseel has less egg production. In high egg productions, the hypothalamus genes are highly expressed such as EXFABP, SNRNP25, FAM114A1, and SIX1 [221]. In the hypothalamus, nerve growth factor response, lipid metabolism, canonical Wnt signaling pathway genes were highly expressed i.e. SIX1, RPS15, IGFBP7 thus play a role in chicken egg production. In laying hens, the dietary corticosterone treatment shows low levels of extracellular fatty acid-binding protein (EXFABP) and suggested that the egg white proteins synthesis and secretion may effect by environmental stress [235].
Many studies reported that ovarian follicular development is stimulated by IGFBPs and plays a role in the ovary to FSH action [236,237]. In chicken adipose tissue, the lipid metabolism gene like insulin-like growth factor binding protein 7 (IGFBP7) was highly expressed and it is correlated with egg production [238,239]. In this study, FAM114A, FAM116A, FAM116B, FAM117A, FAM117B, SIX1, SIX2, RPS13, RPS24, IGFBP2, IGFBP3, and IGFBP5 genes are up-regulated and FABP1, FABP2, FABP3, FABP5, SNRPB, SNRPB2, IGFBP1, and IGFBP7 genes are down-regulated in the 7 th -day embryo of improved Aseel, these differentially expressed genes may cause less egg production in improved Aseel (Fig. 7). The cuticle or organic matrix of the eggshell-related genes i.e. MEPE, BPIFB3, RARRES1, and WAP are highly expressed in oviposition [240,241,242]. In this study, RARA, RARB, POSTN, CDH4, CDH13, CDH8, CDH11, and CDH20 were upregulated, and RARRES1, and CDH1, were down-regulated in the 7 th -day improved Aseel embryo. Due to the differential expression of these genes may be the improved Aseel eggshell thickness is more than commercial laying eggs.

Conclusion
The comparative transcriptome study between slow growth improved Aseel and fast growth control broiler revealed the DEGs and their signi cantly enriched pathways in slow growth improved Aseel, which inferred that they play an important role in regulating the growth and development of improved Aseel. The transcriptome data provides a theoretical basis for improving the performance of the slow growth improved Aseel as well as how to control the growth performance in fast grown control broiler chicken and provides reference data for revealing the molecular mechanism of slow growth improved Aseel as well as fast growth control broiler chicken. In this study, the mechanistic picture of gene expression data shows the embryo development, muscle development, egg production, plumage development, and energy production in improved Aseel, which would be fostered by a combination of (a) Differential regulation of MSTN, activin-like kinases and up-regulation of SMADs expect SMAD7 in the myostatin signaling pathway combined with down-regulation of caveolin's (CAV1, CAV2, and CAV3) and differential regulation of insulin-like growth factor binding proteins, (b) Up-regulation of HSP70, NCF1, Map2k2, and down-regulation of MYOD1, and MYOZ2, (c) Up-regulation of fatty acid synthesis and β-oxidation genes (ACACA, ACACB, FASN, and CPT1), (d) Differential MAPK signaling pathway genes (MAP2K2, MKA, NES, SAMSN1, SOS2, and TAB2) (e) Differential regulation of Jak-STAT, mTOR, and TGF-β signaling pathway genes (IGF1,   IGF2, IRS1, IRS2, PI3K     Hierarchical clustering of differentially expressed genes in 7th and 18th day of embryo (Aseel vs. control broiler). Hierarchical clustering of differentially expressed genes led to the formation of eight distinct clusters: I, II, III, IV, V, VI, VII and VIII, which include genes up and down regulated in muscle imitation and elongation stages of Aseel compared with control broiler, de nes the speci c molecular regulation of muscle growth. Each row represents the expression pattern of a single gene, and each column corresponds to a single sample: column 1, 7th day of embryo and column 2, 18th day of embryo. The expression levels are represented by color chat, with red representing the up regulation, green representing the down regulation and black representing the missing values or no change.

Figure 3
Identi cation of stable reference genes and variability of Ct values in the 7th-day embryo and 18th-day thigh muscle of improved Aseel and control broiler. Total RNA was isolated from embryo and thigh muscle tissues and converted to cDNA, and gene-speci c primers were used in qPCR to determine Ct values.
Mean±SD of Ct values are shown.

Figure 4
Relative expression of DEGs that are involved in muscle development, myostatin signaling, muscle metabolism (energy sensing and storage), and protein synthesis. Y-axis represents relative mRNA expression level and X-axis represents tissue samples used for qPCR study ( Control broiler, PD4)  Diagrammatic representation of DEGs analysis in the 7th-day embryo and 18th-day thigh muscle of aseel as compared to respective controls of control broiler. The red and green arrows represents the up and down-regulated genes (P<0.05), respectively. → the activation of the process, ┤inhibition of the process. 1.
Canonical pathway of Smad activation. Myostatin binds to ActRIIB and induces its assembly with activin type I receptor, subsequent phosphorylation of Smad2/3 leads to its binding with Smad4 and translocation of the complex to the nucleus where it blocks the transcription of genes responsible for the myogenesis. Smad6 and Smad7 compete for the binding with activin type I receptor and Smad7 can also prevent the formation of the