Abstract Background: Muscle development, egg production, and plumage colors are different between native and broiler chickens. The study was designed to investigate why improved Aseel (PD4) is colorful, stronger, and grew slowly compared with the control broiler (CB). Methods: A microarray was conducted using the 7th-day embryo (7EB) and 18th-day thigh muscle (18TM) of improved Aseel and broiler, respectively. Also, we have selected 24 Gallus gallus candidate reference genes from NCBI, and total RNA was isolated from the broiler, improved Aseel embryo tissues, and their expression profiles were studied by real-time quantitative PCR (qPCR). Furthermore, microarray data were validated with qPCR using improved Aseel and broiler embryo tissues. 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, 8,069 transcripts were differentially expressed between the 7EB and 18TM of PD4 compared to the CB. A further analysis showed that a high number of transcripts are differentially regulated in the 7EB of PD4 (6,896) and fewer transcripts are differentially regulated (1,173) in the 18TM of PD4 compared to the CB. On the 7th^- and 18th-day PD4 embryos, 3,890, 3,006, 745, and 428 transcripts were up- and downregulated, respectively. The commonly up- and downregulated transcripts are 91 and 44 between the 7th- and 18th-day of embryos. In addition, the best housekeeping gene was identified. Furthermore, we validated the differentially expressed genes (DEGs) related to muscle growth, myostatin signaling and development, and fatty acid metabolism genes in PD4 and CB embryo tissues by qPCR, and the results correlated with microarray expression data. Conclusion: Our study identified DEGs that regulate the myostatin signaling and differentiation pathway; glycolysis and gluconeogenesis; fatty acid metabolism; Jak-STAT, mTOR, and TGF-β signaling pathways; tryptophan metabolism; and PI3K-Akt signaling pathways in PD4. The results revealed that the gene expression architecture is present in the improved Aseel exhibiting embryo growth that will help improve muscle development, differentiation, egg production, protein synthesis, and plumage formation in PD4 native chickens. Our findings may be used as a model for improving the growth in Aseel as well as optimizing the growth in the broiler. Keywords: fast and slow growth chicken, 7th- and 18th-day embryo tissues, microarray, reference gene, quantitative real-time PCR Introduction Animal agriculture production is essential for supplying protein nutrition to the increasing global human population. The broiler chickens are genetically selected with highly improved production efficiency through rapid growth and high feed efficiency compared to improved Aseel native chicken birds. Therefore, understanding mechanisms regulating rapid muscle growth and high feed efficiency between control broiler and improved Aseel may improve the quality of improved Aseel animal production systems ([32]Niemann et al., 2011). In high and low production efficiency breast muscle phenotypes, male pedigree broiler breeder chickens were used for a global gene expression cDNA microarray study ([33]Kong et al., 2011; [34]Bottje et al., 2012; [35]Bottje and Kong, 2013). Also, RNAseq global gene expression studies have been performed with breast muscle and duodenal tissue in commercial broilers and low and high residual feed intake broilers, respectively ([36]Lee et al., 2015; [37]Zhou et al., 2015). Global gene expression studies mostly showed that production ability could also be related to different cellular mechanisms such as mitochondrial oxidative stress, inflammatory response, protein degradation, stress responses, growth hormone signaling, cell cycle, apoptosis, and fatty acid transportation. A recent transcriptome study reported that differentially expressed genes are enriched in myogenic growth and differentiate on the 6th and 21st day of breast muscle in modern pedigree broiler chickens compared with legacy chicken lines ([38]Davis et al., 2015). A transcriptome analysis was performed with the pectoralis major muscles of slow- and fast-growing chickens (n = 8) to understand the myopathies related to structural changes and molecular pathways using an 8 × 60 K Agilent chicken microarray histological study. For fast-growing breast meat yield, a functional analysis revealed the favoring of metabolic shifts toward alternative catabolic pathways, oxidative stress, inflammation, regeneration, fibrosis processes, cellular defense, and remodeling ([39]Pampouille et al., 2019). A transcriptome profiling analysis was performed in two chicken lines, that is, high (pHu+) and low (pHu-), using an Agilent custom chicken 8 × 60 K microarray. Between these two lines, 1,436 differentially expressed (DE) genes were found, and they were related to biological processes for muscle development and remodeling and carbohydrate and energy metabolism ([40]Beauclercq et al., 2017). A genome-wide association study (GWAS) was conducted to assess body weight in an F2 chicken population, and a microarray expression study was conducted with the liver of high and low-weight chickens. Also, we identified miR-16 as a critical regulator that will suppress chicken embryo cell proliferation and cellular growth. The mutated miR-16 by inserting 54bp showed a significant increase in body weight, bone size, and muscle mass ([41]Jia et al., 2016). The comparative transcriptome analysis of grouper fish muscle in the hybrid and its parents showed up-regulation of genes related to glycolysis, calcium signaling, and troponin pathways that enhanced muscle growth in the hybrid grouper ([42]Sun et al., 2016). In addition, insulin-like growth factor 1 (IGF1) and a cascade of intracellular components (protein kinase B, mTOR, GSK3β, and FoxO) play a significant role in the regulation of skeletal muscle growth during development and regeneration ([43]Schiaffino and Mammucari, 2011). Muscle growth contains complex network combinations, that is, cell proliferation, differentiation, and metabolism ([44]Fuentes et al., 2013). In mammals, during the embryonic period, the total skeletal muscle fiber number is initiated, and after birth, only muscle hypertrophy occurs ([45]Rowe and Goldspink, 1969; [46]Timson and Dudenhoeffer, 1990; [47]Rehfeldt et al., 2000). In teleosts, hypertrophic and hyperplastic muscle growth can happen simultaneously during the entire life ([48]Stickland, 1983; [49]Weatherley et al., 1988). However, chickens’ total skeletal muscle fiber number is initiated and fixed during the embryonic period ([50]Bhattacharya et al., 2015; [51]Bhattacharya et al., 2019). As a consequence, chicken muscle mass accounts for a better proportion of body weight, and it is an excellent experimental model for studying fundamental growth regulatory mechanisms in vertebrates ([52]Weatherley and Gill, 1985). Therefore, for controlling muscle growth, understanding the mechanisms in chicken is necessary to optimize poultry. In America, during the mid-19th century, dual-purpose chicken, that is, Barred Plymouth Rock (BPR), a foundational or heritage breed of the modern commercial broilers, was developed by crosses with Black Java, Black Cochin, and Dominique breed with alternating white and black bars of feather pigmentation ([53]Lopez et al., 2007; [54]Dorshorst and Ashwell, 2009). For egg production, multiple gene interactions in various organs regulate energy metabolism, protein synthesis, and storage ([55]Silversides and Villeneuve, 1999). Previous genomic and transcriptomic reports identified genes associated with reproduction traits ([56]Ciacciariello and Gous, 2005). In total, 26 differentially expressed genes (DEGs) were identified in ovaries between pre-laying and egg-laying periods ([57]Kang et al., 2009). The 12 genes identified were related to reproduction regulation pathways such as GnRH, G protein-coupled receptor, calcium-signaling pathways, biosynthesis of steroid hormones, oocyte meiosis, and progesterone-mediated oocyte maturation ([58]Luan et al., 2014). In chickens, nine transcripts related to high egg production were identified in the hypothalamus and the pituitary gland ([59]Shiue et al., 2006). Recently, a comparative transcriptome study was conducted between the hypothalamus and the pituitary gland in Chinese dagu chickens ([60]Wang and Ma, 2019). However, no studies reported how to regulate the genes in embryos for oogenesis and egg development in chickens. The genetic and developmental foundation of morphological complexity is one of the most significant questions in evolutionary biology, and because avian feathers come in a variety of shapes, they make a great model system for research on the evolution and development of unique morphological features ([61]Losos et al., 2013; [62]Chen B. et al., 2015). Feathers are an excellent model to study the molecular basis of phenotypic variation of an important structure in a single species because they have evolved to have different forms in color, morphology, and mechanical properties not only among different bird species but also in different body regions of an individual bird. The structure and shape of a body feather vary dynamically from the distal end to the proximal end, with the distal end forming before the proximal end. A body feather’s barbs change from being mostly pennaceous at the proximal end to plumulaceous at the distal end ([63]Ng et al., 2015). A great example of exaptation is the feather, which may have originally been developed to regulate body temperature but was later appropriated for show and flight. These and other evolutionary innovations most likely resulted from altered feather-related gene expression patterns. The morphological novelties of feathers are the result of the origin and evolution of plesiomorphic molecular signaling modules ([64]Prum, 2005). A novel technological platform made possible by systems biology research can disclose the molecular expression profiles connected to various morphological processes. The identification of genes linked to changes in feather and scale will be aided by transcriptome investigations and bioinformatic analysis ([65]Ozsolak and Milos, 2011; [66]Chang et al., 2015). A transcriptomic study was conducted on two feather types at different times during their regeneration after plucking and then compared the gene expression patterns in different types of feathers and different portions of a feather and identified morphotype-specific gene expression patterns ([67]Ng et al., 2015). Recently, transcriptome data from yellow and white feather follicles from 7- to 11-week-old F3 chickens were generated to screen for genes involved in the production of pheomelanin particles ([68]Zheng et al., 2020). However, it is still completely unknown what causes feather variance genetically, particularly during embryo development. Understanding the molecular dynamics of embryonic development concerning the process of feather growth can help us better comprehend how different feather shapes have evolved through time. Thus, the objectives of the present study are to explore mechanisms and pathways underlying embryonic development with respect to muscle growth, egg production, and plumage formation in slow-growing native and fast-growing broiler chickens. Materials and methods Animals The study was conducted on the fast-growing broiler pure line (developed from the synthetic population) and improved Aseel (PD4) chicken lines maintained at the institute farm, ICAR-Directorate of Poultry Research, Hyderabad, India ([69]Figure 1A). The improved Aseel (PD-4) has been developed from the Indian native Aseel breed of chicken by imposing selection for body weights at 8 weeks of age for the last 10 generations. The body weight of these birds at 8 weeks during the S-10 generation was 551.0 ± 3.60 g. The control broiler birds are random-bred broilers, and there is no selection imposed on this population. The body weight of the control broiler line at 6 weeks was 951.0 ± 1.20 g. The birds of both populations were maintained under an intensive management system. A total of 60 fertile eggs were kept for hatching (30 for each group) in the incubator (Global Incubators, Hyderabad, India) at 100.3°F temperature and 79.2°F humidity. After the 7 and 18 days of incubation, eggs were harvested (15 for each group), and embryos were collected and stored at −80°C up to total RNA isolation. FIGURE 1. FIGURE 1 [70]Open in a new tab Differentially expressed transcripts/genes (DEGs) in improved Aseel (PD4) during embryo development stages (7th- and 18th-day embryos). (A) Fast- (control broiler) and slow-growth (improved Aseel, PD4) birds are used in this study. (B) Embryo tissues (complete embryos from the 7th day and thigh muscles from the 18th day) are used for the transcriptome analysis. (C) Number of up- and downregulated transcripts (DETs) (p value ≤0.05 and fold change ≥1) in Aseel as compared to control broiler samples. (D) and (E) Number of common and differentially expressed transcripts among the initial (7EB) and developed (18TM) embryo stages. RNA extraction and evaluation For RNA isolation, the complete embryo from the 7th-day and thigh muscle from the 18th-day embryo were used from the control broiler and PD4 lines. The tissue samples were collected from three independent embryos during each time point to isolate RNA and consider each replicate as one biological replicate during each period. Total RNA was isolated using Trizin RNA extraction reagent (GCC Biotech (India) Pvt. Ltd.), according to the manufacturer’s protocol, and RNA was puried by DNase treatment (DNase I solution, HiMedia, India) to remove a trace amount of DNA. The purity and quantity were monitored on 1.2% denatured agarose gels and the NanoDrop 1,000 Spectrophotometer (Thermo Scientic, United States). The quality of total RNA was assessed by checking 200–300 ng of total RNA on an RNA nano chip using an Agilent Bioanalyzer 2100 (Agilent Technologies, United States), according to the manufacturer’s instructions. RNA labeling, amplification, and hybridization The Agilent Quick Amp Kit (Part number: 5190–0442) was used for sample labeling. In addition, 500 ng of total RNA was reverse transcribed using an oligo dT primer tagged to a T7 promoter sequence, and in the same reaction, the cDNA thus obtained was converted to double-stranded cDNA. Labeled cRNA preparation and hybridization on GeneChip and scanning were done following Affymetrix protocols ([71]http://www.affymetrix.com). In the in vitro transcription step, the cDNA was converted to cRNA using the T7 RNA polymerase enzyme, Cy3 dye was added to the reaction mix and incorporated into the newly synthesized strands, and the obtained cRNA was cleaned up using Qiagen RNeasy columns (Qiagen, Cat No: 74106). The concentration and amount of dye incorporated were determined using a NanoDrop 1,000 Spectrophotometer (Thermo Fisher Scientic, United States). The QC-passed samples for specific activities were taken for hybridization. Then 600 ng of labeled cRNA was hybridized on the array using the Gene Expression Hybridization Kit (Part Number 5190–0404; Agilent Technologies, United States) in Sure Hybridization Chambers (Agilent technologies, United States) at 65°C for 16 h. Agilent Gene Expression Wash Buffers (Part No: 5188–5327) were used for washing the hybridized slides, and then the slides were scanned on a G2505C scanner (Agilent Technologies, United States). Microarray data analysis After scanning, DAT, CEL, CHP, XML, and JPEG image were generated for each array with Feature Extraction Software (Version-10.7, Agilent Technologies, United States). The CELles containing estimated probe intensity values were further analyzed with GeneSpring GX-11.0 software (Agilent Technologies, United States). Normalization of the data was performed in GeneSpring GX using the 75th percentile shift, and this normalization takes each column in an experiment independently and computes the n ^th percentile of the expression values for this array across all spots; fold change was calculated concerning specific control samples. Genes were up- and down-regulated showing one-fold and above within the samples concerning the control sample were identified, and for the replicates, a Student’s t-test p-value was calculated. The expression data obtained have been submitted to the Gene Expression Omnibus database (GEO, [72]http://www.ncbi.nlm.nih.gov/geo) at the National Center for Biotechnology Information with the accession numbers [73]GSE62443-[74]GSE62445. Hierarchical clustering analysis The differentially expressed genes between the 7th^- and 18th-day of the embryo were subjected to a hierarchical cluster analysis using the Cluster 3.0 program ([75]Eisen et al., 1998). We imported the matrix with as many columns as stages and 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. To demonstrate the expression pattern and tree diagram of DETs, we applied rows and columns to the cluster 3.0 software and carried out hierarchical clustering using the complete linkage approach with the Euclidean distance. The resulting dendrogram was then exported as an image file. Functional characterization Biological data and analysis tools are combined in the Database for Annotation Visualization and Integrated Discovery (DAVID; [76]https://david.ncifcrf.gov/), which provides systematic functional annotation for a large number of genes and proteins. GO annotation linked with biological processes, molecular function, and KEGG pathway enrichment analyses were carried out using the online DAVID tool version 6.7 (11,12) to examine the possible functions of discovered DEGs ([77]Huang et al., 2009a,[78]b). The FDR p-value of < 0.05 and fold change of > 1 were considered to be significant. Pathway analyses Ingenuity Pathway Analysis (IPA; Qiagen, Valencia, CA; [79]http://www.ingenuity.com) software was used for functional annotation, canonical pathway analysis, upstream analysis, and network discovery. The chicken DEGs data set functionalities are primarily based on mammalian biological mechanisms because IPA is based on human bioinformatics. We have attempted to draw possible conclusions based on avian-based literature, but biomedical research biases the functional annotations toward 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; the sequences were downloaded from NCBI ([80]https://www.ncbi.nlm.nih.gov/), and the CDS (coding DNA sequence) region was identified by using the ExPASy translation tool ([81]https://web.expasy.org/translate/) ([82]Supplementary Table S1). Real-time quantitative PCR analysis Microarray expression data were validated using two-step real-time quantitative PCR (qPCR) for specific confirmation of the differentially expressed genes. First-strand cDNA was synthesized from 2 µg of total RNA using the Thermo Fisher Scientific Verso cDNA Synthesis Kit (Thermo Scientific, United States). Gene speciac qPCR primers were designed for 24 housekeeping genes and 83 DEGs using PrimerQuest software ([83]http://eu.idtdna.com; [84]Supplementary Table S1; [85]Table 1). The qPCR was performed using the BrightGreen 2X qPCR MasterMix-No Dye (Applied Biological Materials Inc. Canada) in the Insta Q96™ Real-Time Machine (HiMedia Laboratories, India) detection system. The PCR was performed under the following program: 5 min at 95°C followed by 40 cycles of amplication with 15 s of denaturation at 95°C and 60 s of annealing/extension at 60°C. A total of three biological replicates were used. A melt curve analysis was performed to check the special city of the amplied products. The 2^−ΔΔCt calculated the relative expression level of each gene, and Transferrin (TFRC; Accession No: [86]X55348.1) from G. gallus was used as a housekeeping gene to normalize the amount of template cDNA added in each reaction. TABLE 1. List of primers used for qPCR to validate microarray data S. No Gene Name Accession No. Primer Sequence Tm GC% Amplicon Size (bp) 1 Destrobrevin alpha (D alpha) [87]CR733292.1 5′-CAA​CCC​TTT​GTG​GAG​GAA​AGA-3′ 62 47.6 114 5′-GAA​CCT​CCC​GCA​GAA​ACA​A-3′ 62 52.6 2 Uncharacterized protein (UP5) [88]ES605836.1 5′-GAA​CCA​AAT​GCT​GGC​AGA​AG-3′ 62 50 112 5′-AAA​TAC​TCT​CTG​GGT​GAA​CAG​G-3′ 62 45.5 3 Toll-interacting protein (TOLLIP) [89]NM_001006471.1 5′-GTG​TAA​CGA​AGA​GGA​CCT​GAA​A-3′ 62 45.5 95 5′-TGT​TCC​CTC​TCT​GAG​CTT​CTA-3′ 62 47.6 4 Asw [90]CN225783.1 5′-GGC​AAC​ACG​TGA​AAT​CCA​TTC-3′ 62 47.6 119 5′-GCG​CAC​GTC​TCT​GTA​TTC​T-3′ 62 52.6 5 Chain A, fibrinogen alpha subunit (Chain A FAS) [91]BX935039.1 5′-TGA​CGA​CAC​AGA​CCA​GAA​TTA​C-3′ 62 45.5 106 5′-GGT​TTC​CAC​AAT​TAC​CCG​ATT​G-3′ 62 45.5 6 Hypothetical protein (HP29) [92]AW198329.1 5′-CCC​AGA​TGA​CAG​AAG​AAC​AAA​TAA​AG-3′ 62 38.5 106 5′-CCC​TCT​TCT​CCA​AAG​CAT​GTA​T-3′ 62 45.5 7 Fibrinogen gamma chain precursor (FGCP) [93]BG642009.1 5′-CTG​GTC​ACC​TCA​ATG​GAC​AAT​A-3′ 62 45.5 106 5′-CAT​CGG​TCA​CGC​CAT​GTT-3′ 62 55.6 8 Apolipoprotein B precursor (ALPBP) [94]NM_001044633.1 5′-CTT​GAG​GCC​AAC​TCC​AAA​GTA-3′ 62 47.6 102 5′-GTG​CTC​CCA​GAC​TGC​ATA​AA-3′ 62 50 9 Maestro heat-like repeat–containing protein family member 2B (MHCRCP2B) [95]CR406681.1 5′-CTG​GAA​CAC​ACC​ACA​GAC​TT-3′ 62 50 130 5′-CCC​GAT​AGA​TGT​CCT​TTC​CAT​AC-3′ 62 47.8 10 Activin A receptor, type IB (AARIB) [96]XM_001231300 5′-GCA​CGG​ATC​TCT​CTT​TGA​CTA​C-3′ 62 50 120 5′-TGA​GTA​CCC​ACG​ATC​TCC​AT-3′ 62 50 11 cAMP responsive element modulator (CREM) [97]NM_204387 5′-CAA​GAG​AGA​GCT​GCG​ACT​TAT​G-3′ 62 50 102 5′-AGC​ACA​GCC​ACA​CGA​TTT-3′ 62 50 12 Caveolin 1, caveolae protein, 22kDa (CAV1) [98]NM_001105664 5′-CAT​TCC​CAT​GGC​ACT​CAT​CT-3′ 62 50 106 5′-GCA​CTG​GAT​C′CAA​TCA​GGT​AG-3′ 62 50 13 Caveolin 2 (CAV2) [99]NM_001007086 5′-TGC​TGT​ACA​AGC​TGC​TGA​G-3′ 62 52.6 140 5′-CAC​TGA​AGG​CAA​GAC​CAT​GA-3′ 62 50 14 Follistatin-like 1 (FSTL1) [100]NM_204638 5′-CGA​TGA​CAT​GTG​AAG​GGA​AGA-3′ 62 47.6 105 5′-TCT​GCA​GCT​CCT​GAA​CAT​ATC-3′ 62 47.6 15 WAP, follistatin/kazal, immunoglobulin, kunitz, and netrin domain containing 1 (WAPFK) NP9672441 5′-GAG​GGC​AAC​AAC​AAC​AAC​TTC 62 47.6 109 5′-TCA​GCA​CCA​TCT​TGC​TCT​TC-3′ 62 50 16 Glucose-6-phosphate isomerase (GPI) [101]NM_001006128 5′-CAC​TTC​TGC​CCT​ATG​ACC​AAT​A-3′ 62 45.5 110 5′-GTA​GTC​CAC​ACG​AGA​TCC​TTT​C-3′ 62 50 17 Solute carrier family 2 (facilitated glucose transporter), member 3 (SLC2A3) [102]NM_205511 5′-GTG​GTA​CAC​AGG​ATG​TAT​CTC​AAG-3′ 62 45.8 112 5′-CGA​TAG​TTT​GGA​GAG​CGG​AAT​AG-3′ 62 47.8 18 Heat shock 60kDa protein 1 (chaperonin) (HSPD1), nuclear gene encoding mitochondrial protein [103]NM_001012916 5′-GGT​GAG​AAG​GCT​CAG​ATT​GAA-3′ 62 47.6 122 5′-GCT​ACT​CCG​TCA​GAT​AGT​TTG​G-3′ 62 50 19 Heat shock 70kDa protein 5 (glucose-regulated protein, 78kDa) (HSPA5) [104]NM_205491 5′-TGA​GAC​AGT​TGG​AGG​TGT​AAT​G-3′ 62 45.5 103 5′-AGT​GGG​CTG​ATT​GTC​AGA​AG-3′ 62 50 20 Heat shock 70kDa protein 8 (HSPA8) [105]NM_205003 5′-AGT​TTG​AGC​TGA​CCG​GTA​TTC-3′ 62 47.6 122 5′-CTC​CTT​GCC​AGT​GCT​CTT​ATC-3′ 62 52.4 21 Heat shock factor binding protein 1 (HSBP1) [106]NM_001112809 5′-ATG​CAG​GAC​AAA​TTT​CAA​ACC​A-3′ 62 36.4 118 5′-CTA​CTC​CCG​CTT​GTG​TCA​TC-3′ 62 55 22 Heat shock transcription factor 1 (HSTF 1) [107]BM440477 5′-GCA​GCA​GAA​GGT​GGT​CAA​TA-3′ 62 50 146 5′-AGT​ACT​GGC​GGC​TGT​ATT​TC-3′ 62 50 23 Partial mRNA for heat shock protein 70 (hsp70 gene) [108]AJ301880 5′-CCC​AGT​AAG​TGC​GGG​TCA​TAA-3′ 64 52.4 85 5′-CGC​TCC​GCC​AGT​CAC​TT-3′ 64 64.7 24 Homeobox C9 (HBC9) [109]BX950823 5′-AGA​TGT​CCG​TAC​ACA​AAG​TAT​CA-3′ 62 39.1 105 5′-GTT​TAG​GAC​TCG​GGC​TAC​TTC-3′ 62 52.4 25 Insulin-like growth factor 1 receptor (IGF1R) [110]NM_205032 5′-TGT​GTA​CGT​TCC​AGA​CGA​ATG-3′ 62 47.6 104 5′-CCT​TGG​CTA​TTC​CCT​CAT​ACA​C-3′ 62 50 26 Insulin-like growth factor binding protein 1 (IGFBP1) [111]NM_001001294 5′-CAG​GAC​CAG​ATG​CTG​AAC​TAT​C-3′ 62 50 134 5′-CCC​TGT​TCT​TTC​CAT​TTC​TTG​TG-3′ 62 43.5 27 Mitogen-activated protein kinase 8 interacting protein 3 (MAPK8IP3) [112]XM_424591 5′-GTG​ATG​ACA​ACA​GCG​ACA​AAT​C-3′ 62 45.5 118 5′-CCA​GGC​ACA​GAG​ACA​AAG​AA-3′ 62 50 28 Mitogen-activated protein kinase kinase kinase 4 (MAPKKK4) [113]CR523470 5′-AGT​GGA​TGA​ACT​ACG​TGC​TAA​C-3′ 62 45.5 120 5′-CCG​GGA​GAG​CCG​AAA​TAA​AT-3′ 62 50 29 Mitogen-activated protein kinase-activated protein kinase 3 (MAPKAPK3) [114]XM_414262 5′-CTG​AAG​ACT​GAC​CCA​ACA​GAG-3′ 62 52.4 139 5′-CTT​CAT​CCC​AGT​GGT​CCT​TAT​C-3′ 62 50 30 Myozenin 2 (MZ2) [115]BX930590 5′-GAA​ACA​ACA​AGC​ATC​AGC​CAT​TA-3′ 62 39.1 121 5′-GCT​GAG​TGT​TGA​TAG​TTC​CTC​TAC-3′ 62 45.8 31 Angiotensin II receptor, type 1 (AGTR1) [116]NM_205157 5′-TTC​CTG​GAT​TCC​TCA​TCA​AGT​G-3′ 62 45.5 103 5′-GGG​CAT​AGC​TGT​ATC​CAC​AAT​A-3′ 62 45.5 32 Angiotensin II receptor–associated protein (AGTRAP) [117]BX930324 5′-CTT​CAA​CAT​AGG​TCT​CAA​CCG​T-3′ 62 45.5 106 5′-CTG​AGC​TGC​CTT​GCT​TGA-3′ 62 55.6 33 CD9 molecule (CD9) [118]NM_204762 5′-TAC​TAC​AAT​GCC​ATG​CCC​TAA​A-3′ 62 40.9 134 5′-TAG​CAC​AGC​AAA​GAA​CCA​TAC​T-3′ 62 40.9 34 Dickkopf homolog 2 (DKK2) [119]XM_420494 5′-CGC​AAC​AAG​AAG​AAC​AGT​CAT​TAT-3′ 62 37.5 105 5′-GGG​ATC​ACC​TTC​ATG​TCC​TTT​A-3′ 62 45.5 35 Glycoprotein M6A (GPM6A) [120]NM_001012579 5′-GGA​TCT​TCG​CCA​GTA​TGG​TAT​T-3′ 62 45.5 97 5′-TAG​CTC​ATT​CGA​GTC​ACA​CAT​C-3′ 62 45.5 36 Glycoprotein M6B (GPM6B) [121]NM_001012545 5′-GAA​CAT​CTG​CAA​CAC​GAA​TGA​G-3′ 62 45.5 124 5′-GGC​CCA​GTT​AGA​AGA​CAG​TAT​C-3′ 62 50 37 Janus kinase 1 (JAK1) [122]NM_204870 5′-CAA​GGA​ACT​AGC​TGA​CCT​GAT​G-3′ 62 50 98 5′-CCT​CCA​GTT​TGT​TGA​TGT​CTC​T-3′ 62 45.5 38 Janus kinase 2 (JAK2) [123]NM_001030538 5′-GAT​GGA​TGC​CCT​GAT​GAG​ATT-3′ 62 47.6 92 5′-CGC​TGA​GCA​AGA​TCC​CTA​AA-3′ 62 50 39 Janus kinase and microtubule interacting protein 2 (JAKMIP2) [124]CR390426 5′-GAC​TGC​ATC​AGT​TCA​TCA​TTT​CTC-3’ 62 41.7 130 5’-ACA​GGA​ACA​CAT​TGC​TGG​T-3′ 62 47.4 40 Janus kinase and microtubule interacting protein 3 (JAKMIP3) [125]XM_426548 5′-TAT​CAA​CTT​CCA​CCA​CGT​TCC-3′ 62 47.6 100 5′-CAT​CAG​CTC​TGC​CAC​TAC​TAT​G-3′ 62 50 41 Leiomodin 3 (fetal) (LMOD3) [126]BX935813 5′-GAG​AAT​GAC​TGC​AGA​GGA​GAT​G-3′ 62 50 97 5′-TTT​GTA​GTG​CCG​CTC​CTT​C-3′ 62 52.6 42 Musculoskeletal, embryonic nuclear protein 1 (MUSTN1) [127]NM_213580 5′-CCA​AGT​CAT​GAA​GCA​GTG​TGA-3′ 63 47.6 94 5′-TGA​CTT​CTC​AAA​GAC​CGT​TTC​G-3′ 63 45.5 43 Myosin binding protein C, fast type (MYBPC2) [128]NM_001044659 5′-CTG​ATG​GAG​CGC​AAG​AAG​AA-3′ 62 50 105 5′-GAA​GAC​GCC​CTC​GAT​CAT​TT-3′ 62 50 44 Myosin binding protein C, slow type (MYBPC1) [129]BX935207 5′-CCT​GAA​ACG​TAG​GGA​GGT​TAA​A-3′ 62 45.5 131 5′-TGC​CTC​TCA​GGT​CAG​TGA​TA-3′ 62 50 45 Perilipin 1 (PLIN1) [130]NM_001127439 5′-CCA​GAA​GAG​GAG​GAG​GAA​GAT-3′ 62 52.4 100 5′-TAG​CAC​TGT​GAG​CCC​TGT​A-3′ 62 52.6 46 Phospholamban (PLN) [131]NM_205410 5′-CGA​TAG​CAG​GGT​TTC​CAT​ACT​T-3′ 62 45.5 117 5′-TGT​CAG​CTC​TCT​CCA​GTA​GAA-3′ 62 47.6 47 RCD1 required for cell differentiation1 homolog (S. pombe) (RS. pombe 001006521 5′-TGA​TTG​GAG​CCT​TGG​TGA​AA-3′ 62 45 105 5′-GTT​CAC​TGC​CAG​ACT​CCA​TAA​T-3′ 62 45.5 48 Slow muscle troponin T (TNNT1) [132]NM_205114 5′-CCC​TCC​ACA​TTG​AGC​ACA​T-3′ 62 52.6 104 5′-CTC​CAT​CAG​GTC​GAA​CTT​CTC-3′ 62 52.4 49 Troponin T type 3 (skeletal, fast) (TNNT3) [133]NM_204922 5′-GAA​GCA​AAC​AGC​TAG​AGA​GAC​A-3′ 62 45.5 125 5′-GGT​ATA​ACC​AGT​CCC​ACA​GTT​C-3′ 62 50 50 Troponin I type 1 (skeletal, slow) (TNNI1) [134]BX931462 5′-TCT​CTT​CGT​CCA​CAA​TCT​CAA​C-3′ 62 45.5 128 5′-ACA​GTC​GGA​GAA​GGA​GAG​ATA​C-3′ 62 50 51 Myostatin (MSTN) [135]NM_001001461.1 5′-GGA​TGG​GAC​TGG​ATT​ATA​GCA​C-3′ 62 50 97 5′-GGT​GAG​TGT​GCG​GGT​ATT​T-3′ 62 52.6 52 Follistatin (FST) [136]NM_205200.1 5′-ACA​ACT​TAT​CCG​AGC​GAG​TG-3′ 62 50 112 5′-CTT​CCT​CTG​GGT​CTT​CGT​TAA​T-3′ 62 45.5 53 Activin A receptor type 2A (ACVR2A) [137]NM_205367.1 5′-GCA​AGA​ATG​TGC​TGC​TGA​AA-3′ 62 45 109 5′-CCA​ACC​TGT​CCA​TGT​GTA​TCT-3′ 62 47.6 54 Activin A receptor type 2B (ACVR2B) [138]NM_204317.1 5′-GAA​GTG​TTA​GAG​GGA​GCA​ATC​A-3′ 62 47.5 118 5′-CTG​GAC​CAT​CAA​CTG​CTC​TAC-3′ 62 52.4 55 SMAD family member 2-Z (SMAD2Z) [139]NM_204561.1 5′-GGG​AGT​GCG​TCT​CTA​TTA​CAT​C-3′ 62 50 110 5′-CAG​GAT​GCC​AGC​CAT​ATC​TT-3′ 62 50 56 Activin A receptor type 1B (ACVR1B) [140]XM_015300267.2 5′-GCA​CGG​ATC​TCT​CTT​TGA​CTA​C-3′ 62 50 120 5′-TGA​GTA​CCC​ACG​ATC​TCC​AT-3′ 62 50 57 Transforming growth factor beta receptor 1 (TGFBR1) [141]NM_204246.1 5′-TCG​TGT​GCC​AAG​TGA​AGA​AG-3′ 62 50 102 5′-CCA​GAG​CCT​GAA​GTT​GTC​ATA​TC-3′ 63 47.8 58 Myogenin (MYOG) [142]NM_204184.1 5’-GGC​TGA​AGA​AGG​TGA​ACG​AA-3’ 62 50 116 5’-GCG​CTC​GAT​GTA​CTG​GAT​G-3’ 62 57.9 59 Mitogen-activated protein kinase kinase 6 (MAP2K6) [143]XM_003642348.2 5′-CTC​AGC​AGA​G′TCG​TCG​ATT​T-3′ 62 47.6 101 5′-GCA​GGG​TGA​AGA​AAG​GAT​GT-3′ 62 50 60 Mitogen-activated protein kinase kinase kinase 7 (MAP3K7) [144]XM_015284683.2 5′-CAG​CCC​TTG​TTT​CAG​GAG​AAG-3′ 63 52.4 101 5′-GCC​TCG​TTT​AGG​CTT​GGA​ATA​G-3′ 63 50 61 Caveolin 3 (CAV3) [145]NM_204370.2 5′-GCT​TTG​ATG​GTG​TCT​GGA​AAG-3′ 61 47.6 142 5′-ATG​TGG​CAG​AAG​GAG​ATG​AG-3′ 61 50 62 Protein kinase AMP-activated catalytic subunit alpha 1 (PRKAA1) [146]NM_001039603.1 5′-CTT​GAC​GAT​CAC​CAT​CTG​TCT​C-3′ 62 50 140 5′-TGC​CAC​TTC​GCT​CTT​CTT​AC-3′ 62 50 63 Protein kinase AMP-activated catalytic subunit alpha 2 (PRKAA2) [147]NM_001039605.1 5′-GGA​GGT​CTG​TGA​GAA​GTT​TGA​G-3′ 62 50 124 5′-GTT​CAT​GAT​CCT​CCG​GTT​GT-3′ 62 50 64 Creatine kinase, M-type (CKM) [148]NM_205507.1 5′-CGA​CCA​CTT​CCT​GTT​CGA​TAA-3′ 62 47.6 109 5′-GAA​CGT​CTT​GTT​GTC​GTT​GTG-3′ 62 47.6 65 Mechanistic target of rapamycin (serine/threonine kinase) (MTOR) [149]XM_417614.5 5′-AAG​GTT​TCT​TCC​GGT​CCA​TAT​C-3′ 62 45.5 98 5′-ATC​AGG​CCA​GTG​ACC​ATA​ATC-3′ 62 47.6 66 Ribosomal protein S6 kinase A1 (RPS6KA1) [150]NM_001109771.2 5′-GGA​ACC​CAG​CCA​ACA​GAT​TA-3′ 62 50 104 5′-TTC​CCT​TCG​GTA​CAG​CTT​ATT​C-3′ 62 45.5 67 Carnitine palmitoyltransferase I (CPT1) [151]DQ314726.1 5′-GCC​TTC​GTG​CGC​AGT​AT-3′ 62 58.8 146 5′-ACG​TAG​AGG​CAG​AAG​AGG​T-3′ 62 52.6 68 Acyl-CoA synthetase long-chain family member 1 (ACSL1) [152]NM_001012578.1 5′-GCC​AGT​ACG​TAG​GCA​TCT​TT-3′ 62 50 116 5′-TGC​TTC​AGT​TCC​CAG​TGT​ATC-3′ 62 47.6 69 Enoyl-CoA hydratase, short chain 1 (ECHS1) [153]NM_001277395.1 5′-CAG​GTG​GGA​GCT​ATT​GTC​ATC-3′ 62 52.4 97 5′-CAT​AGC​ACT​CCT​GGA​AGG​TTT-3′ 62 47.6 70 Hydroxyacyl-CoA dehydrogenase (HADH) [154]NM_001277897.1 5′-GCT​ATC​CCA​TGG​GTC​CAT​TT-3′ 62 50 100 5′-AGA​GGA​TTG​TTG​GGC​TCT​ATT​G-3′ 62 45.5 71 Acyl-CoA oxidase 2 (ACOX2) [155]XM_015293306.2 5′-TGC​CAC​CAT​CTG​TCA​CTT​ATC-3′ 62 47.6 141 5′-TAG​CTG​CTG​TGC​TGC​TTA​TC-3′ 62 50 72 Sterol regulatory element binding protein 1 (SREBP1) [156]AJ310768.1 5′-CAT​GGA​GGT​GGC​GAA​GG-3′ 62 64.7 134 5′-TGT​CAG​GCT​CGG​AGT​CA-3′ 62 58.8 73 Fibroblast growth factor 2 (FGF2) [157]NM_205433.1 5′-TTC​GAG​CGC​TTG​GAA​TCT​AAT​A-3′ 62 40.9 94 5′-GCT​TGT​ACT​GTC​CAG​TCC​TTT-3′ 62 47.6 74 Fibroblast growth factor receptor 1 (FGFR1) [158]NM_205510.1 5′-CGT​CAC​CAA​AGT​GGC​TGT​A-3′ 62 52.6 98 5′-TGC​CGA​TCA​TCT​TCA​TCA​TCT​C-3′ 62 45.5 75 DNA methyltransferase 3 alpha (DNMT3A) [159]NM_001024832.1 5′-CCT​TCT​TCT​GGC​TCT​TTG​AGA​A-3′ 62 45.5 111 5′-CAG​ACA​CCT​CTT​TGG​CAT​CA-3′ 62 50 76 Forkhead box O3 (FOXO3) [160]MK861853.1 5′-CTC​TCA​GGC​TCC​TCT​TTG​TAT​TC-3′ 62 47.8 109 5′-CAC​ACT​CCA​AGC​TCC​CAT​T-3′ 62 52.6 77 Peroxisome proliferator–activated receptor gamma (PPARγ) [161]AF163811.1 5′-CCC​AAG​TTT​GAG​TTT​GCT​GTG-3′ 62 47.6 99 5′-TGG​GCG​ATC​TCC​ACT​TAG​TA-3′ 62 50 78 Myogenic factor 6 (MYF6) [162]FJ882409.1 5′-GCT​GGA​TCA​GCA​GGA​CAA​A-3′ 62 52.6 100 5′-GCA​GGT​GCT​CAG​GAA​GTC-3′ 62 61.1 79 Acetyl-CoA carboxylase α(ACACA) [163]NM_205505.1 5′-CAG​ATT​TGT​TGT​CAT​GGT​GAC-3′ 60 42.9 162 5′-ACA​GCC​TGC​ACT​GGA​ATG​C-3′ 60 57.9 80 Acetyl-CoA carboxylase β(ACACB) [164]XM_025155692.1 5′-GCT​CCT​GCT​GCC​CAT​ATA​TTA-3′ 60 47.6 94 5′-GTC​CGT​GAT​GAC​ACC​TTT​CT-3′ 60 50 81 Fatty acid synthase (FASN) [165]NM_205155.3 5′-GTT​CTC​TGT​ACA​GAG​AAT​GTG-3′ 60 42.9 168 5′-CCA​TGT​TTG​ACT​TGG​TTG​ATC-3′ 60 42.9 [166]Open in a new tab The statistical analysis To assess the expression variation of the candidate reference genes, all the samples were divided into three broad categories: the combination of 7EB and 18TM samples of control broiler and improved Aseel, 7EB and 18TM samples of control broiler alone, and 7EB and 18TM samples of improved Aseel alone. The qRT-PCR machine–generated Ct values for each of the cDNA samples were then used to determine the degree of data variability between the samples. The stability level of the 24 candidate reference genes from the 7EB and 18TM of control broiler and improved Aseel was determined using five statistical algorithms: geNorm, NormFinder, BestKeeper, Delta CT, and RefFinder ([167]Vandesompele et al., 2002; [168]Andersen et al., 2004; [169]Pfaffl et al., 2004; [170]Silver et al., 2006; [171]Xie et al., 2012). GeNorm and NormFinder, which use relative expression values as input data and convert Ct values to linear scale expression quantities using the 2^−∆Ct method, as well as BestKeeper, which uses the raw Ct value directly and the comparative Ct method, were used to determine the expression stability level. RefFinder, a program that combines the most important computational tools currently available (geNorm, Normfinder, BestKeeper, and the comparative 2^−∆Ct method) to compare and rank the stability of candidate reference genes with the geometric mean of individual gene appropriate weight, was used to determine the overall final ranking of reference genes across all samples. Comparison of mean expression values for qRT-PCR between the control broiler and PD4 improved Aseel groups were made using the Student’s t-test and p ≤ 0.05 was considered statistically significant. Results Source, selection, primer design, and verification of candidate reference genes In the present study, to identify the suitable reference genes for the 7th- and 18th-day embryos of the 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. These are 18S rRNA, 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, and the CDS region was found and amplified with gene-specific primers ([172]Supplementary Table S1). For all the primer pairs, the melting curve analysis was performed to confirm the specific amplification for each reference gene, a single peak with no visible primer-dimer formation and genomic DNA contamination was observed, and no signals were detected in the non-template controls. Expression stability and ranking of candidate reference genes Expression levels of all candidate reference genes were measured in the samples collected from the 7EB and 18TM of the control broiler and improved Aseel. Each reference gene had different expression ranges across all sample sets, and the 18S rRNA and DNAJC24 genes showed the most (Ct = 12.34) and the least (Ct = 33.88) abundant transcripts, respectively. In the combined analysis, we observed that not all selected reference genes were expressed uniformly across 7EB and 18TM of the control broiler and improved Aseel. The genes with the lowest global variability were GUSB, PP2CB, and HSP70 ([173]Figure 2A). The results show that the GUSB reference gene had the least variation in 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 ([174]Figure 2A). In control broiler samples, PP2CB, ALB, and GUSB were the top three genes with the lowest variation ([175]Figure 3B); whereas HSP70, GUSB, and β-2MG showed little variation in improved Aseel ([176]Figure 3C). 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 was expressed constantly across the 7EB and 18TM of control broiler and improved Aseel in the present study. Therefore, it is essential to choose the most reliable reference gene for expression profiling gene/s in different embryos of the control broiler and improved Aseel. The most popular statistical tools geNorm, NormFinder, BestKeeper, Delta CT, and RefFinder were used for the analysis to choose the best and most trustworthy reference gene and rank all the potential reference genes according to their stability values for accurate gene expression ([177]Table 2). The variation among the reference genes determined by geNorm is stability measure (M value) and pairwise comparison expression ratio and provides an optimal number of genes in a given experiment. NormFinder measures the reference gene stability by overall expression variation and across samples variation to reduce sensitivity toward co-regulation. BestKeeper calculates the gene expression variation based on Ct values, calculates the Pearson correlation coefficient by pairwise correlation analysis for all reference genes, and finds the stable genes. The Delta CT method directly used the raw Ct values and found the best stable genes. RefFinder is conclusive of the calculations using the aforementioned algorithms and suggested stable genes. FIGURE 2. FIGURE 2 [178]Open in a new tab Relative level of mRNA expression of the 24 candidate reference genes in the 7EB and 18TM of both control broiler and improved Aseel. The data are presented as mean cycle threshold (Ct) values and shown as box and whisker plots. The boxes represent the interquartile range of the mean Ct values, whereas the middle, up, and lower bars represent the mid-hinge, maximum, and minimum Ct values, respectively. The X-axis represents the gene names, and Y-axis represents the Ct values of all the tissues. (A) 7EB and 18TM samples of control broiler and improved Aseel; (B) 7EB and 18TM samples of control broiler; (C) 7EB and 18TM samples of PD4. FIGURE 3. FIGURE 3 [179]Open in a new tab Hierarchical clustering of DEGs on the 7EB and 18TM of the 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 downregulated, defining the specific molecular regulation of Aseel growth. Each row represents the expression pattern of a single gene, and each column corresponds to a single sample: column 1, 7EB; column 2, 18TM. The expression levels are represented by a color chat (p value ≤0.05 and fold change ≥1), with red representing upregulation, green representing downregulation, and black representing the missing values or no change. TABLE 2. Ranking of the candidate reference genes according to their stability value per indicated software Gene Name geNorm NormFinder BestKeeper ΔCT Comprctensive 7EB & 18TM of CB & PD4 combined analysis 7EB & 18TM of CB analysis 7EB & 18TM of PD4 analysis 7EB & 18TM of CB & PD4 combined analysis 7EB & 18TM of CB analysis 7EB & 18TM of PD4 analysis 7EB & 18TM of CB & PD4 combined analysis 7EB & 18TM of CB analysis 7EB & 18TM of PD4 analysis 7EB & 18TM of CB & PD4 combined analysis 7EB & 18TM of CB analysis 7EB & 18TM of PD4 analysis 7EB & 18TM of CB & PD4 combined analysis 7EB & 18TM of CB analysis 7EB & 18TM of PD4 analysis M R M R M R SV R SV R SV R SD R SD R SD R SD R SD R SD R GM R GM R GM R 18S rRNA 2.05 18 2.27 18 1.70 18 3.02 19 3.25 18 2.26 19 1.83 6 2.13 7 1.35 5 3.75 20 4.09 19 2.96 21 14.43 18 14.6 18 13.95 18 ALB 3.00 23 2.48 19 2.60 23 5.99 24 3.22 17 7.95 24 2.54 10 1.13 2 4.56 24 6.32 24 3.86 18 8.08 24 19.28 23 10.52 12 24 24 B2MG 2.70 22 3.27 23 1.82 19 4.02 23 5.48 24 2.27 20 1.66 4 2.25 8 0.99 3 4.4 23 5.63 24 2.93 20 14.85 19 18.24 22 12.45 16 βActin 1.48 13 1.64 14 1.04 11 2.38 16 2.74 16 2.01 16 4.12 23 4.52 22 3.72 22 2.88 16 3.18 16 2.44 14 16.94 21 17.05 21 15.59 20 DNAJC24 1.60 15 1.72 15 1.16 13 2.71 17 3.25 19 2.28 21 4.37 24 4.98 23 3.77 23 3.13 17 3.59 17 2.63 17 18.25 22 18.57 23 18.41 22 EEF1A1 1.54 14 1.83 16 0.48 2 0.87 4 1.25 6 0.33 2 1.73 5 1.39 4 2.08 8 2.43 9 2.85 13 1.86 3 7.21 6 8.53 7 3.46 3 GAPDH 1.18 8 1.23 8 1.10 12 0.37 1 0.06 1 0.57 4 1.91 7 1.97 6 1.85 7 2.28 5 2.52 4 2.01 7 4.21 3 3.83 3 7.1 6 GUSB 2.38 20 2.86 21 1.99 21 3.07 20 3.83 21 2.40 22 1.14 3 1.29 3 0.89 2 3.68 19 4.28 21 3.02 22 12.44 15 13.06 16 12.08 15 HMBS 0.67 2 0.54 2 0.91 8 1.02 5 0.87 4 1.20 9 3.01 14 3.21 14 2.82 16 2.24 3 2.37 3 2.08 10 5.01 4 4.74 4 10.67 12 HSP10 1.01 6 0.97 5 0.99 10 1.86 14 2.0 14 1.81 14 3.79 21 4 20 3.58 21 2.54 12 2.73 10 2.29 13 12.54 16 11.38 13 14.32 19 HSP70 2.54 21 3.06 22 1.91 20 3.55 21 4.76 23 2.23 18 1.08 1 1.62 5 0.42 1 4.04 22 5.01 23 2.92 19 10.04 10 15.71 19 9.21 8 L-LDBC 1.32 10 1.42 10 1.28 15 1.09 6 1.31 7 0.95 7 2.45 9 2.74 9 2.16 9 2.45 10 2.72 8 2.19 12 8.78 8 8.63 8 10.49 11 MRPS27 1.10 7 1.32 9 0.95 9 1.36 10 1.21 5 1.57 13 3.08 16 2.94 13 3.22 19 2.43 8 2.62 6 2.18 11 10.06 11 7.9 6 12.84 17 MRPS30 0.92 5 1.03 6 0.77 5 1.48 12 1.80 13 1.20 8 3.26 18 3.43 17 3.10 18 2.42 7 2.74 11 1.96 6 9.76 9 11.42 14 8.49 7 PGK2 0.38 1 0.35 1 0.41 1 0.52 3 0.51 3 0.48 3 2.59 11 2.8 10 2.38 11 2.13 2 2.35 2 1.82 2 2.85 2 2.78 2 2.85 2 PPP2CB 2.23 19 2.69 20 1.56 17 2.91 18 3.80 20 1.86 15 1.11 2 1.12 1 1.1 4 3.55 18 4.25 20 2.67 18 10.67 13 9.57 11 11.81 14 RPL13 0.83 4 0.75 3 0.82 6 1.15 7 1.47 9 0.82 5 3.06 15 3.39 16 2.72 14 2.26 4 2.54 5 1.89 4 6.77 5 7.33 5 6.65 5 RPL14 1.66 16 1.13 7 1.41 16 1.45 11 1.46 8 1.54 12 2.15 8 2.86 11 1.42 6 2.66 13 2.72 9 2.55 15 11.81 14 8.92 9 11.64 13 RPL19 1.25 9 1.53 12 0.62 3 1.33 8 1.76 12 0.87 6 2.92 13 3.24 15 2.60 13 2.46 11 2.9 14 1.91 5 10.34 12 13.45 17 6.28 4 RPL23 1.43 12 1.59 13 0.70 4 2.09 15 2.60 15 1.22 11 3.64 20 4.46 21 2.82 15 2.77 15 3.09 15 2.04 9 15.55 20 16.04 20 9.28 9 RPL5 1.83 17 2.03 17 2.10 22 3.59 22 4.52 22 2.83 23 3.93 22 5.32 24 2.54 12 4.04 21 4.75 22 3.43 23 20.68 24 21.38 24 19.55 23 SDHA 1.37 11 1.48 11 1.22 14 1.84 13 1.68 11 2.11 17 3.45 19 3.43 18 3.47 20 2.7 14 2.82 12 2.60 16 14.27 17 12.99 15 16.90 21 TBP 0.78 3 0.88 4 0.87 7 1.35 9 1.61 10 1.21 10 3.23 17 3.57 19 2.89 17 2.35 6 2.66 7 2.03 8 7.78 7 9.03 10 10.21 10 TFRC 0.38 1 0.35 1 0.41 1 0.4 2 0.49 2 0.27 1 2.59 12 2.88 12 2.31 10 2.1 1 2.33 1 1.79 1 2.21 1 2.21 1 1.78 1 [180]Open in a new tab M, gene expression stability measure; SD, standard deviation value; SV, stability value; GM, geomean value; and R, ranking Differentially expressed transcripts during embryo development stages To study the effect of muscle development, genome-wide expression analysis was carried out at muscle initiation (7th-day embryo) and muscle development (18th-day thigh muscle) stages ([181]Figures 1A,B). Labeled RNA was hybridized to the Affymetrix GeneChip™ Chicken Genome Array. After statistical data analysis, transcripts with an FDR-adjusted p-value ≤ 0.05 and a fold change ≥ 1 were considered as differentially expressed transcripts (DETs) ([182]Figure 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 [183]Supplementary Material S1. In total, 8,069 transcripts, which accounted for approximately 24% of the total transcripts present on the GeneChip™ Chicken Genome Array, showed differential expression in improved Aseel at various stages analyzed. The maximum number of transcripts (6,896, 21% of total DETs) showed differential expression on the 7th day of the embryo, and the least number of transcripts (1,173, 3.5% of total DETs) showed differential expression on the 18th day of the embryo. Commonly up- and downregulated muscle-responsive transcripts were identied among the embryo development stages to nd out the degree of overlap ([184]Figures 1D,E). The maximum unique number of upregulated (3,799) and downregulated transcripts (2962) was observed in a 7th-day embryo. A small number of upregulated (654) and downregulated (384) transcripts were uniquely differentially expressed on the 18th-day of the thigh muscle. The commonly differentially expressed (91 upregulated and 44 downregulated) transcripts were identified among the embryo development stages, respectively ([185]Figures 1D,E; [186]Supplementary Material S1). Cluster analysis of differentially expressed transcripts To profile the gene expression patterns in response to muscle slow growth and egg production during embryo development, the 8,069 DETs were classied using hierarchical clustering. The expression patterns were separated into eight major clusters (I–VIII) based on tree branching ([187]Figure 3). Transcripts and involved pathways present in each stage within each cluster are presented in [188]Supplementary Material S2. Among the eight major clusters, upregulated transcripts were enriched in clusters I, VI, and VII, and downregulated transcripts were enriched in clusters II, III, and IV. GO annotation and pathway enrichment analysis of differentially expressed genes DAVID 6.7 was used to annotate and enrich the DEGs related to GO annotations and pathways between the 7th-day up and the 18th-day down and vice versa. The results showed that coenzyme metabolism, cell division and chromosome partitioning, outer membrane, and transcription/cell division and chromosome partitioning functions were upregulated in the 7th-day embryo and downregulated in the 18th-day thigh muscle, whereas chaperones, lipid metabolism, outer membrane/carbohydrate transport and metabolism, protein turnover, and posttranslational modification functions were downregulated in the 7th-day embryo and upregulated in the 18th-day thigh muscle. Cell envelope biogenesis, cytoskeleton, and inorganic ion transport and metabolism functions were differentially regulated between the 7th-day embryo and 18th-day thigh muscle compared to the respective controls ([189]Supplementary Material S3). The KEGG pathway enrichment analysis of DEGs was performed using the IPA tool. The results showed that differential regulation of pathways between 7th-day embryo and 18th-day thigh muscle of PD4 compared to their respective controls, that is, Cell cycle, Cell adhesion molecules (CAMs), SNARE interactions in vesicular transport, Oocyte meiosis, Endocytosis, Apoptosis, ABC transporters, Calcium signaling pathway, MAPK signaling pathway, Wnt signaling pathway, Jak-STAT signaling pathway, Toll-like receptor signaling pathway, TGF-βsignaling pathway, cytokine–cytokine receptor interaction, Basal transcription factors, Focal adhesion, Tight junction, Regulation of actin cytoskeleton, Cardiac muscle contraction, Vascular smooth muscle contraction, Insulin signaling pathway, Oxidative phosphorylation, Glutathione metabolism, Glycolysis/Gluconeogenesis, Citrate cycle (TCA cycle), Pentose phosphate pathway, Pyruvate metabolism, Fatty acid biosynthesis, Fatty acid metabolism, Glycerophospholipid metabolism, Heparan sulfate biosynthesis, N-Glycan biosynthesis, Purine metabolism, Pyrimidine metabolism, Tryptophan metabolism, Serine and threonine metabolism, and Valine, leucine and isoleucine degradation ([190]Supplementary Materials S2, S3). Validation of DEGs by qPCR In this study, the expression of DEGs between the fast (CB) and slow growth (PD4) chickens of the 7th-day embryo and 18th-day thigh muscle was verified by qPCR ([191]Supplementary Table S2). The verified transcripts were divided into four groups: i. muscle development, myostatin signaling, muscle metabolism, and protein synthesis ([192]Figure 4), ii. Embryo development ([193]Figure 5), iii. Fatty acid metabolism ([194]Figure 6), and iv. Cell signaling and egg production ([195]Figure 7). The results showed that the expression trend of the DEGs between the fast and slow-growing chickens is consistent in qPCR results, and this attests to the reliability of the microarray data. FIGURE 4. [196]FIGURE 4 [197]Open in a new tab Relative expression of DEGs that are involved in muscle development, myostatin signaling, muscle metabolism (energy sensing and storage), and protein synthesis. The Y-axis represents relative mRNA expression level, and the X-axis represents tissue samples used for the qPCR study ( Control broiler, PD4). The p values have been stated on the comparative bars; * indicates the p ≤ 0.0 5, and NS indicates the non-significant. Standard error was used for error bars. FIGURE 5. [198]FIGURE 5 [199]Open in a new tab Relative expression of DEGs that are involved in embryo development. The Y-axis represents relative mRNA expression level, and the X-axis represents tissue samples used for the qPCR study ( Control broiler, PD4). The p values have been stated on the comparative bars; * indicates the p ≤ 0.05, and NS indicates the non-significant. Standard error was used for error bars. FIGURE 6. [200]FIGURE 6 [201]Open in a new tab Relative expression of DEGs that are involved in fatty acid metabolism. The Y-axis represents relative mRNA expression level, and the X-axis represents tissue samples used for the qPCR study ( Control broiler, PD4). The p values have been stated on the comparative bars; * indicates the p ≤ 0.05, and NS indicates the non-significant. Standard error was used for error bars. FIGURE 7. [202]FIGURE 7 [203]Open in a new tab Relative expression of DEGs that are involved in cell signaling and egg production. The Y-axis represents relative mRNA expression level, and the X-axis represents tissue samples used for the qPCR study ( Control broiler, PD4). The p values have been stated on the comparative bars; * indicates the p ≤ 0.05, and NS indicates the non-significant. Standard error was used for error bars. Discussion Globally, chicken is one of the most protein-rich meat sources. Muscle development and egg production are essential genetic traits in commercially grown chickens. However, not much information is available on genes involved in muscle development and egg production in slow and fast-growing chickens. In this study, we selected fast (CB) and slow-growth (PD4) chickens to determine the expression of genes related to embryo initiation and developmental stages. Microarray was conducted with the 7th-day embryo (7EB) and 18th-day thigh muscle (18TM) of PD4 and CB, respectively. According to the MIQE guidelines, selecting suitable reference genes may vary for different species, varieties, experimental conditions, and tissues and has to be validated before gene expression study ([204]Bustin et al., 2009). Previous and recent studies also described different expression patterns of reference genes and focused on validating reference genes applied to different avian tissues ([205]Olias et al., 2014; [206]Bages et al., 2015; [207]Nascimento et al., 2015). However, so far, validation of genes for their stable expression patterns in different embryo tissues, such as 7th and 18th-day embryos of control broiler and improved Aseel has not been performed. Candidate reference genes validation 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, that is, geNorm, NormFinder, BestKeeper, Delta CT, and RefFinder ([208]De Boever et al., 2008; [209]Yue et al., 2010; [210]Yang F. et al., 2013; [211]Olias et al., 2014; [212]Bages et al., 2015; [213]Nascimento et al., 2015; [214]Borowska et al., 2016; [215]Mitra et al., 2016; [216]Khan et al., 2017; [217]Zhang et al., 2018; [218]Mogilicherla et al., 2022). These algorithms showed some differences in the stability ranking of stable reference genes, which may be due to the differences in each statistical program ([219]Table 2). Ranking the stability of 24 genes is crucial, as is confirming the number of reference genes required for precise gene expression profiling in various embryonic tissues. For each gene, geNorm generates a stability measure (M value), allowing for ranking based on expression stability (with the lower value indicating increased gene stability across samples). To assess the value of including more references, it