Modern broiler chickens exhibit very rapid growth and high feed efficiency compared to unselected chicken breeds. The improved production efficiency in modern broiler chickens was achieved by the intensive genetic selection for meat production. This study was designed to investigate the genetic alterations accumulated in modern broiler breeder lines during selective breeding conducted over several decades.
To identify genes important in determining muscle growth and feed efficiency in broilers, RNA sequencing (RNAseq) was conducted with breast muscle in modern pedigree male (PeM) broilers (n = 6 per group), and with an unselected foundation broiler line (Barred Plymouth Rock; BPR). The RNAseq analysis was carried out using Ilumina Hiseq (2 x 100 bp paired end read) and raw reads were assembled with the galgal4 reference chicken genome. With normalized RPM values, genes showing >10 average read counts were chosen and genes showing <0.05 p-value and >1.3 fold change were considered as differentially expressed (DE) between PeM and BPR. DE genes were subjected to Ingenuity Pathway Analysis (IPA) for bioinformatic functional interpretation.
The results indicate that 2,464 DE genes were identified in the comparison between PeM and BPR. Interestingly, the expression of genes encoding mitochondrial proteins in chicken are significantly biased towards the BPR group, suggesting a lowered mitochondrial content in PeM chicken muscles compared to BPR chicken. This result is inconsistent with more slow muscle fibers bearing a lower mitochondrial content in the PeM. The molecular, cellular and physiological functions of DE genes in the comparison between PeM and BPR include organismal injury, carbohydrate metabolism, cell growth/proliferation, and skeletal muscle system development, indicating that cellular mechanisms in modern broiler lines are tightly associated with rapid growth and differential muscle fiber contents compared to the unselected BPR line. Particularly, PDGF (platelet derived growth factor) signaling and NFE2L2 (nuclear factor, erythroid 2-like 2; also known as NRF2) mediated oxidative stress response pathways appear to be activated in modern broiler compared to the foundational BPR line. Upstream and network analyses revealed that the MSTN (myostatin) –FST (follistatin) interactions and inhibition of AR (androgen receptor) were predicted to be effective regulatory factors for DE genes in modern broiler line. PRKAG3 (protein kinase, AMP-activated, gamma 3 non-catalytic subunit) and LIPE (lipase E) are predicted as core regulatory factors for myogenic development, nutrient and lipid metabolism.
The highly upregulated genes in PeM may represent phenotypes of subclinical myopathy commonly observed in the commercial broiler breast tissue, that can lead to muscle hardening, named as woody breast. By investigating global gene expression in a highly selected pedigree broiler line and a foundational breed (Barred Plymouth Rock), the results provide insight into cellular mechanisms that regulate muscle growth, fiber composition and feed efficiency.
Production efficiency in animal agriculture continues to be of vital importance to meet the protein nutrition needs of an ever increasing global human population. Genetically selected modern broiler chickens acquired highly improved production efficiency through rapid growth and high feed efficiency compared to unselected chicken breeds. Thus, understanding of mechanisms regulating rapid muscle growth and high feed efficiency may play an increasingly important role in developing more efficient, and therefore sustainable, animal production systems .
A global gene expression study using cDNA microarray assay was conducted on breast muscle phenotyped by high and low production efficiency in male pedigree broiler breeder chickens [2–4]. Additional global gene expression studies have been performed with breast muscle tissue using RNAseq in commercial broilers  and duodenal tissue using RNAseq in broilers selected for low and high residual feed intake (RFI) . Several global expression studies showed that production efficiency may be associated with various cellular mechanisms including mitochondrial oxidative stress, inflammatory response, protein degradation, stress responses, growth hormone signaling, cell cycle, apoptosis and fatty acid transportation. A recent study on the transcriptome analyses of breast muscles in modern pedigree broiler chicken and legacy chicken lines reported that the transcriptional profile of differentially enriched genes on day 6 and day 21 in modern broiler muscles includes genes for enhanced myogenic growth and differentiation .
The Barred Plymouth Rock (BPR) chicken, with its characteristic pattern of alternating white and black bars of feather pigmentation were developed by crosses involving the Black Java, Black Cochin, and Dominique breeds in America during the mid-19th century  and is a foundational or heritage breed of the modern commercial broilers. The BPR breed, originally developed for the dual purposes of both meat and egg production, has much slower growth rate (589 g empty body weight at 42 days) compared to commercial broilers (2276 g at 42 days) .
With distinct differences in growth characteristics, the objective of the present study was to conduct global gene expression analyses on breast muscle tissue obtained from modern pedigree male (PeM) broilers (rapid growth and muscle development) compared with the foundational BPR chickens (slow growth and lower efficiency) . Recent popular RNA-sequencing method using massively parallel nucleotide sequencing techniques were used for the comparison of global gene expression in this study. From this comparison, DE genes showing greater differences between the two groups may provide more definitive insights into cellular regulatory mechanisms in muscle growth and feed efficiency.
The present study was conducted in accordance with the recommendations in the guide for the care and use of laboratory animals of the National Institutes of Health. All procedures for animal care complied with the University of Arkansas Institutional Animal Care and Use Committee (IACUC): Protocol #14012.
Animals, breast muscle tissue and RNA extraction
Breast muscle tissue analyzed in this study were obtained from immature (≤8 weeks old) pedigree broiler males (PeM) and BPR. The PeM samples had been obtained previously  and from an inbred male line highly selected for growth and feed efficiency . Briefly, immature PeM and BPR chickens (<8 weeks old, n = 6 per breed) were killed by an overdose of sodium pentobarbital (i.v.) and breast muscle tissue obtained and flash frozen in liquid nitrogen. RNA was extracted from breast muscle using TRIzol reagent (Invitrogen Life Technologies, Thermo-Fisher Scientific, Carlsbad, CA) following the manufacturer’s instructions. After the first RNA extraction, the RNA samples were treated with DNase I (New England Biolabs Inc., Ipswich, MA) and extracted again with TRIzol reagent. The RNA quality was then assessed by Agilent 2200 TapeStation instrument (Santa Clara, CA). The RNA samples were used for RNA sequencing (RNAseq) analysis.
RNAseq and data analysis
RNAseq library preparation for individual samples with barcoding and Illumina sequencing were performed by the Research Technology Support Facility at Michigan State University (East Lansing, MI). Illumina HiSeq system 2x100 bp paired end read technology was used for RNAseq. The RNA sequence FASTA files were mapped to the chicken genome (galgal4) retrieved from Ensembl using CLC Genomics Workbench software (Qiagen, Valencia, CA; licensed to CSIRO) following the software manual (http://resources.qiagenbioinformatics.com/manuals/clcgenomicsworkbench/752/User_Manual.pdf; http://www.clcbio.com/wp-content/uploads/2012/08/RNA-Seq_analysis_part_I.pdf) and the RNAseq analytical pipeline recommended in . Briefly, raw reads were processed by reads’ quality control and adapter trimming with default setting. Filtered reads were subjected to reference based-RNAseq pipeline of CLC software with default settings. The EMBL annotation files of reference chicken genome were downloaded from Ensembl (ftp://ftp.ensembl.org/pub/release-87/embl/gallus_gallus/) and tracks for sequence and genes were generated with track conversion. The parameters of initial RNAseq mapping are as followings: Use reference genome annotated with genes and transcripts; Eukaryote origin; Map to gene region only; mapping options with 2 mismatch costs, 3 insertion costs, 3 deletion costs, 0.7 length fraction, 0.8 similarity fraction, both strand specific, 10 maximum number of hits for a read; total counts for expression value. Total mapped counts were transformed to log2 values of the number of reads per million (RPM) to stabilize the variance and then a further quantile normalization was performed . Normalized values were subjected to further statistical analyses performed by JMP Genomics 6 (SAS Institute Inc., Cary, NC). Genes with less than 10 average raw read counts in both comparison groups were removed and not considered for further analyses. The t-statistics was used to compare between PeM and BPR. Genes showing >1.3 fold differences and <0.05 p-value in the comparison between PeM and BPR were considered differentially expressed (DE) genes. The p-value correction (FDR calculation) by multiple tests was not applied in this study since we used a less stringent approach on a gene by gene basis, and this allowed us to import much more informative lists and assess the genome-wide data using various alternative approaches, such as the hyper-geometric stats used by functional enrichments for pathway analyses, the co-expression patterns in the hierarchical heatmaps, and the binomial distributions in the mitoproteome plot [12, 13].
The complete human mitoproteome was downloaded from http://mitominer.mrc-mbu.cam.ac.uk/release-3.1/begin.do generating a list of 1,046 nuclear and mitochondrial genes encoding mitochondrial proteins. Protein names were converted to gene names and 699 genes matched the chicken RNAseq data in this study. The expression of mitoproteome genes was then plotted in MA format, and all those with higher values in the PeM group were color coded red and those with lower values in the PeM color coded blue. The skew in distribution away from the null expectation of 50:50 was quantified by binomial statistics. This approach formalizes the extent to which the mitoproteome data is biased to one or the other of the two groups.
Hierarchical clustering analysis
The genes with greatest DE (20 up- and 20 down-regulated) between PeM and BPR were subjected to hierarchical cluster analysis using the JMP Genomics program. We imported the matrix with as many columns as birds (12) and as many rows as genes (40), where each cell contains the log2 transformed fold change value for that gene and individual into JMP Genomics program, normalizing on rows. We applied hierarchical clustering on both rows and columns and exported the resulting dendrogram as an image file.
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. Because IPA is based on the bioinformatics in humans, functionalities for differentially expressed genes in our chicken datasets are based primarily on mammalian biological mechanisms. As the investment in biomedical research biases the functional annotations towards human disease, we have attempted to draw plausible conclusions based on avian based literature.
One microgram of total RNA was obtained from 6 muscle samples each for PeM and BPR for general validation of the RNAseq results and for specific confirmation of the top 10 genes, respectively. RNA was converted into cDNA with qScript™ cDNA SuperMix (Quanta Biosciences, Gaithersburg, MD) following the manufacturer’s instructions. The cDNA samples were diluted in a 1:10 ratio and a portion (2 μL) of the cDNA was subjected to a real-time quantitative PCR (qPCR) reaction using an ABI prism 7500HT system (ThermoFisher Scientific, Waltham, MA) with PowerUp™ SYBR® Green Master Mix (ThermoFisher Scientific, Waltham, MA). The specific oligonucleotide primers were designed by the PRIMER3 program (http://frodo.wi.mit.edu) and the list of primers were listed in Additional file 1. The conditions of real-time qPCR amplification were 1 cycle at 95 °C for 2 min, 40 cycles at 95 °C for 30 s, 65 °C for 30 s. The chicken glyceraldehyde 3-phosphate dehydrogenase (GAPDH) gene was used as the internal control. Dissociation curves were performed at the end of amplification for validating data quality. All qPCR reactions were performed 3 times and the values of average cycle threshold (Ct) were determined for each sample, and 2−ΔΔCt values for the comparison of PeM and BPR were used for relative quantification by fold-change and statistical significance.
Results and discussion
Twelve RNAseq libraries were constructed using RNA samples of breast muscles from 6 chickens each for PeM and BPR and sequenced with 2x 100 bp paired end reads on 4 lanes. Totally ~800 million 100 bp sequence reads were obtained with an average of ~67 million reads per sample and ~80% of sequence reads were mapped to the chicken reference genome (galgal4). After quantile normalization with log2 transformation of the number of reads per million (RPM) and filtering out very low expression genes (raw read counts <10), 9,263 genes remained and the MA plot showed that up- and down-regulated genes were evenly distributed in the comparison between PeM and BPR (Fig. 1a). Finally, 2,464 DE genes showing <0.05 p-value and >1.3 fold change were identified in the comparison between PeM and BPR (Additional File 2). To validate RNAseq results, a subset of 20 genes selected from the DE list were subjected to qPCR. These were genes associated with muscle structures [slow twitch subunits (MYH15, TPM3, MYOZ2, TNNI1, MYBPC1), slow muscle associated proteins (MB, CA3), muscle fat content (PLN, FABP4), muscle growth/differentiation (MSTN)], mitochondrial biogenesis and proteins (PPARGC1A, PPARGC1B, CKMT1A, NOXO1), and five out of the top 10 downregulated genes (IL23R, GPM6A, PPDPF, ACE, CD9). Results indicated that expression levels for 16 genes out of 20 tested were well-matched between RNAseq and qPCR analyses in terms of the direction and magnitude of change (Table 1). Four genes were not matched between RNAseq and qPCR for reasons that are not clear to us, but one possibility relates to the fact that the two methods use quite different normalization approaches. Those transcripts where the two methods make similar DE predictions can be considered the most robust of our reported data. Additionally, hierarchical clustering showed clear discrimination the 12 birds into the correct group of origin (Fig. 2).
Of the 699 genes that were detected that encode proteins localized to the mitochondrion, 194 and 447 were numerically higher and lower, respectively, in PeM compared to BPR, while 58 genes were not changed (Fig. 1b). This is a very significant skew (binomial P value < 0.000001) implying the PeM birds have a lowered mitochondrial content. Recently, a breast muscle transcriptomic analyses using RNAseq method was reported on the comparison between the broiler pedigree Ross and a legacy Illinois chicken line . When raw sequencing data, which is available on Gene Expression Omnibus (GEO), were reanalyzed for gene expression, the 608 genes encoding proteins localized to the mitochondrion were differentially expressed as 251 and 357 were up- and downregulated, respectively, in the 21 day old Ross chickens compared to Illinois chickens (data not shown; Seo et al., manuscript in preparation). Similar findings of these two independent studies showing a lower number of upregulated genes encoding mitochondrial proteins in modern broiler compared to unselected chicken breeds may strongly support reduced numbers of mitochondrial content in modern breeds selected for growth and feed efficiency.
This mitochondrial content finding is in line with the economic physiological design argument proposed by Hudson  in the context of production animal feed efficiency. This argument proposes that an economically designed feed converter should have enough, but not too much, fuel burning capacity. This is because spare capacity bears costs of construction (mitochondrial biogenesis), maintenance (trans-membrane proton gradient) and load (space that could be used for other organelles) that should be avoided if they are unnecessary.
Hierarchical clustering on extreme DE
The top 20 up- and 20 downregulated genes were subjected to hierarchical clustering analysis (Fig. 2). The results clearly discriminate the 12 birds into the correct group of origin. Functional groups contained in the 20 upregulated genes in PeM include cytoskeleton and muscle contraction [LMOD2 (leiomodin 2, cardiac), CA3 (carbonic anhydrase III), MUSTN1 (musculoskeletal, embryonic nuclear protein 1), TNNT2 (troponin T type 2, cardiac), MYBPC1 (myosin binding protein C, slow type), MB (myoglobin), MYH15 (myosin heavy chain 15)], extracellular matrix [or secreting proteins; TNC (tenascin C), ANGPTL5 (angiopoietin-like 5), NOV (nephroblastoma overexpressed), IGJ (joining chain of multimeric IgA and IgM)], transporter [PLN (phospholamban), MB], protease [PCSK1 (proprotein convertase subtilisin/kexin type 1)], and ion channel [SCN5A (sodium channel, voltage gated, type V alpha subunit)]. Functional groups contained in the 20 downregulated genes include transcription factors [FHL2 (four and a half LIM domains 2), KLF9 (Kruppel-like factor 9)], calcium signaling [PROCA1 (protein interacting with cyclin A1), GPM6A (glycoprotein M6A)], synthetic enzymes [GART (phosphoribosylglycinamide formyltransferase), CGT (UDP glycosyltransferase 8), RHOF (ras homolog family member F), GSTA3 (glutathione S-transferase alpha 3)], protein ubiquitination [KLHL38 (kelch-like family member 38)], membrane bound [SLC38A4 (solute carrier family 38, member 4), CD9 (CD9 molecule), ACE (angiotensin I converting enzyme), GPM6A, IL23R (interleukin 23 receptor)], and secreting protein [DKK2 (dickkopf WNT signaling pathway inhibitor 2), KL (klotho)]. More details on the function of a narrower list for the 10 most DE genes are described below.
The 10 most DE genes
The full name of the 10 most DE genes between PeM and BPR muscle and their functional characteristics are shown in Tables 2 and 3. Interestingly, highly upregulated genes in PeM muscle are closely associated with expression of gene and phenotypes of subclinical myopathy commonly observed in the commercial broiler breast tissue, as discussed below. Five of the 10 most upregulated genes (LMOD2, CA3, MUSTN1, TNNT2, and MYBPC1) in PeM are mainly expressed in skeletal or cardiac muscles (Tables 2 and 3), particularly in slow twitch fibers (type 1), except for MUSTN1 which is not expressed in slow twitch fibers. This observation indicates that breast muscle in PeM chickens could have higher slow fiber composition compared to BPR muscle. When considering the mitoproteome and mitochondrial contents discussed earlier, the possible higher slow fiber composition in PeM is counter to the lowered gene expression for mitoproteome, which are typically associated with fast twitch muscle fibers. Generally, slow-twitch or oxidative fibers have high mitochondrial content, slow contraction rates, increased reliance on oxidative phosphorylation (OXPHOS), high resistance to fatigue, and high representation in postural muscles. In contrast, fast-twitch or glycolytic fibers have lower mitochondrial content, rapid contractions, decreased reliance on OXPHOS, low resistance to fatigue, and high representation in muscle groups used for directional movement . In general, we know that chickens have a white breast muscle with a very high preponderance of fast twitch glycolytic fibers bearing a low mitochondrial content (e.g. ). Based on the DE genes encoding mitochondrial proteins (Fig. 1b), modern PeM breast muscles may contain an even higher content of fast twitch fibers compared to BPR muscle. However, the most upregulated genes showed the possibility of muscle composition, that contain higher gene expression for the slow muscle fiber composition in modern PeM compared to BPR. This indicates that chicken breast muscle mostly consists of fast fibers, more rapidly growing modern PeM broiler muscle contains increased slow fiber compositions, and the increased mass of slow fibers may be regulated independently of mitochondrial content and function. It also suggests that the mitochondria in avian muscle may function to regulate muscle mass growth by modulating oxidative stress, but may not regulate fiber composition. These results may represent a fundamentally new finding about muscle physiology in birds. In an earlier report by Mutryn et al. , wooden breast myopathy muscle (which in commercial broilers showed a dramatic upregulation of TNNI1, MB, MYBPC1, and MYOZ2 genes compared to unaffected muscle), suggested fiber type switching from fast- to slow type fibers. In addition to the 10 most upregulated genes, expression of MB (10.85 fold) and MYOZ2 (5.90 fold) were also highly upregulated in PeM compared to BPR (Additional file 2). Taken together, PeM muscles may contain increased slow fiber compositions compared to BPR.
Regarding upregulation of PLIN1 (perilipin 1) in PeM, mammalian PLIN1 is known to be expressed in adipose tissue, but not in skeletal muscle . In this study, PLIN1 from each bird showed a very low expression level (range of 10 – 66 raw read counts) in BPR breast muscle while relatively higher expression (range of 260 – 1,700 raw read counts) was observed in PeM breast muscles. Possibly, these data indicate that a potentially new pathway related to hormone sensitive lipase is activated in the large muscle mass of modern broilers compared to the leaner mass in the unselected foundation BPR line. Additionally, increased expression of adipose genes showed that PeM muscles may retain increased mass of adipose tissues, where the wooden breast myopathy showed increased occurrence of irregular adipose tissue throughout the muscle . Higher expression of NOXO1 (NADPH oxidase organizer 1) indicates that oxidative stress and scavenging antioxidants are closely associated with feed efficiency, which has been reported previously [2, 3, 5]. PCSK1 (proprotein convertase subtilisin/kexin type 1) is not restricted to the expression in muscle tissues, and its functions in glucose and lipid metabolism [20, 21] implies that PeM muscle has more active cellular pathways for energy metabolism. Active energy metabolism has been known to be closely connected to oxidative stress in muscle .
In contrast to the 10 most up-regulated genes, expression of the 10 most down-regulated genes in PeM muscles are not restricted to muscle, but are expressed in most tissues and organs (Tables 2 and 3). Of the 10 most down-regulated genes, six genes are classified to membrane bound channels and secreting factors [IL23R, KL, GPM6A, GABRD (gamma-aminobutyric acid (GABA) A receptor, delta), FGB (fibrinogen beta chain), PPDPF (pancreatic progenitor cell differentiation and proliferation factor) and CD9], indicating that PeM muscles may retain characteristics of lower amounts of secreting factors compared to BPR muscles. The expression of ACE, in combination with angiotensin II and AGTR1 (angiotensin II type 1 receptor), can produce muscle atrophy , indicating that down-regulation of ACE in PeM muscle may contribute to muscle hypertrophy compared to BPR chickens. Davis et al.  also observed elevated ACE in a heritage (Illinois) broiler compared to modern broiler (Ross) and hypothesized that this would play a role in lowering protein synthesis and muscle development in the slower growing Illinois broiler.
Overview of IPA analysis
For more stringent pathway analyses, DE genes in muscle from PeM and BPR were subjected to biological functional analyses using Ingenuity Pathway Analysis (IPA, Qiagen, Valencia, CA) and the analytic settings were limited to publications on muscle tissue only. The summary of IPA analysis in diseases, molecular, cellular and physiological functions for DE genes are presented in Table 4. The most relevant biological functions in skeletal muscles comparison between PeM and BPR included organismal injury/abnormalities, carbohydrate metabolism, cellular growth/proliferation, small molecule biochemistry, abnormal morphology of muscle, angiogenesis of skeletal muscle, and myogenesis of skeletal muscles. Further functional analyses including major functions of DE genes, canonical pathways, biological functions, upstream analyses, and molecular networks that can influence muscle growth and feed efficiency are discussed below.
Downregulation of mitochondrial function in PeM muscle
A total of 168 differentially expressed mitoproteome genes were independently analyzed by IPA to characterize functional alterations of mitochondrial energy expenditure and cellular respiration. As expected, mitochondrial functions including oxidative phosphorylation, and TCA cycle activity were downregulated in PeM muscles compared to BPR muscles (Additional file 3: Figure S1 and Additional file 4). Mitochondrial complex 2, 3, 4, and 5 were downregulated in PeM (with no change in complex 1) and the genes responsible for oxidative phosphorylation were decreased in PeM compared to BPR. Likewise, the responsible genes for TCA cycles were decreased in PeM compared to BPR (Additional file 4).
IPA analysis generated 380 canonical pathways with DE genes and most of them were associated with general signal transduction pathways which are not limited to muscle growth/mass or feed efficiency. Thus, the two most relevant pathways for rapid growth and greater mass of muscle (up-regulated PDGF signaling pathway; [24, 25]) and feed efficiency (Up-regulation of NFE2L2 mediated oxidative stress pathway; [5, 26]) are discussed below. Other potentially important canonical pathways will be discussed in future study.
Up-regulated PDGF signaling pathway
The PDGF signaling pathway appeared to be activated in PeM compared to BPR chickens. PDGFR α and β subunits, and their downstream signal transduction molecules, including JAK, SRC, STAT1 (signal transducer and activator of transcription 1), STAT3, RAS (ras oncogene), PLCγ (phospholipase C, gamma), PKC (protein kinase C), cJUN (cellular Jun protooncogene), cFOS (cellular FBJ murine osteosarcoma viral oncogene homolog), MEKK1 (MAPK/ERK kinase kinase, a.k.a. MAP3K1, mitogen-activated protein kinase kinase kinase 1), MKK4 (mitogen activated protein kinase 4, a.k.a. MAP2K4, mitogen activated protein kinase kinase 4), and JNK1 were upregulated or predicted to be activated (Fig. 3). Since skeletal muscle cell itself does not proliferate and PDGF is known to stimulate the proliferation of muscle satellite cells [24, 25], increased activity of PDGF signaling pathways to STAT1/STAT3 complex  and to cJUN/cFOS complex  may influence satellite cell formation and its proliferation, resulting in greater muscle cell growth in PeM chickens compared with BPR. Additionally, upregulation of PDGF signaling pathways in PeM muscle may represent a possible subclinical condition of monocyte infiltration of the regenerating muscle tissues, which has been observed as part of wooden breast myopathy in modern commercial broilers [17, 19].
Up-regulation of NFE2L2 mediated oxidative stress pathway
Recently, Zhou et al. , reported that with commercial broiler chickens, NFE2L2 (nuclear factor, erythroid 2-like 2; also known as NRF2; Fig. 4) mediated oxidative stress pathway was predicted to be activated in breast muscles in PeM compared to low feed efficiency muscle in commercial broilers. In addition, our recent proteomics results showed that NFE2L2 was also predicted to be activated in the high feed efficiency pedigree male . NFE2L2 as a key transcription factor induces a number of antioxidant genes in response to oxidative stress as well as tissue injury and inflammation . In the present study, NFE2L2 was predicted to be activated in PeM compared to BPR based on the increased expression of downstream targets including TRXR1 (thioredoxin reductase 1 or TXNRD1 in Additional file 2), SOD (superoxide reductase or SOD3 in Additional file 2), CAT (catalase), HO-1 (heme oxygenase 1 or HMOX1 in Additional file 2) and MAF (v-maf avian musculoaponeurotic fibrosarcoma oncogene homolog or MAFA/MAFF in Additional file 2) (Fig. 4). However, the expression of other molecules that respond to NFE2L2 were variable (Table 5). In addition, many other genes [CAT, CCL5 (chemokine C-C motif ligand 5), DUSP1 (dual specificity phosphatase 1), FOS, JUN, MGST1 (microsomal glutathione S-transferase 1), PRDX4 (peroxiredoxin 4), SOD3, STAT3, VCAM1 (vascular cell adhesion molecule 1), and XDH (xanthine dehydrogenase)], along with several antioxidant genes associated with oxidative stress, were upregulated in PeM muscle compared to BPR (Table 6). Earlier biochemical studies for mitochondrial functions and oxidative respirations with breast muscle tissues of high- and low feed efficiency chickens reported lowered ROS levels in high feed efficiency chickens [10, 22], indicating high feed efficiency muscles (or rapidly growing muscles) retain more efficient scavenging mechanisms to reduce ROS levels compared to low feed efficiency muscle. Thus, rather than the possibilities that higher levels of oxidative stress may play a role for PeM, the fact that maintaining low ROS levels caused by higher antioxidant functions is necessary for rapid growth would be a better explanation for the biological functions of the DE genes. Taken together, NFE2L2 mediated antioxidant/detoxifying pathway appears to be more active in muscle of modern PeM broilers compared to BPR chickens.
Upstream regulators and network analyses
Upstream regulators and gene network analyses, which represent the intermolecular connections among interacting genes based on functional knowledge inputs, were performed on the DE genes using the IPA program. IPA analyses selected 11 potentially important upstream regulators (p-value <0.05) and 9 meaningful molecular networks (data not shown). Interestingly most networks were centered with the upstream regulators. The upstream regulators centered in the molecular networks, that are possibly important for muscle growth and feed efficiency are presented below and gene information including full gene name and fold change values for focus molecules in each network are listed in Additional file 5.
The first molecule of interest was MSTN (myostatin), a member of TGF-β family, which is an important negative regulator of skeletal muscle growth . MSTN expression was higher in PeM muscle compared to BPR muscle (Additional file 2). Although the DE of a few downstream target molecules [e.g. DES (desmin) and IGF2 (insulin like growth factor 2)] was indicated as being inconsistent with previous reports [31, 32], the downstream molecules of MSTN including increased expression of FN1 (fibronectin 1), PID1 (phosphotyrosine interaction domain containing 1), SERPINH1 [serpin peptidase inhibitor, clade H (heat shock protein 47), member 1, (collagen binding protein 1)], CTSV (cathepsin V) and decreased expression of PPARGC1A [peroxisome proliferator-activated receptor gamma, coactivator 1 alpha; PGC1α] [31–34] indicated that the functional MSTN activity or expression was enhanced in PeM compared to BPR (Fig. 5). Since PeM chickens grow much faster than BPR, MSTN functions in PeM were expected to be less active than BPR, but the result was opposite. The qPCR of MSTN was consistent with RNAseq result (Table 1) and no single nucleotide polymorphisms (SNP) altering an amino acid residue were found in MSTN cDNA sequences in both PeM and BPR chickens (data not shown). According to a recent RNAseq report comparing a modern pedigree broiler (ROSS) muscle and to that in heritage chicken (Illinois) muscles, MSTN expression was very similar in both chicken lines and showed relatively high expression levels (~30 RPKMs) in both lines, though slightly higher expression was observed in pedigree broiler muscles than heritage chicken muscles . Similarly, the MSTN expression in this study showed relatively abundant expression levels in both BPR (range of 380 – 1000 raw read counts) and in PeM (range of 1000 – 2100 raw read counts) breast muscles (data not shown). Thus, assuming no mutations were involved in the MSTN gene in both lines, high expression levels of MSTN could be hypothesized to prevent extreme muscle growth in both PeM and BPR. The regulatory factors that modulate the function of TGF-β family including MSTN, secreted protein follistatin (FST), together with its associated proteins (follistatin like 1, FSTL1 and follistatin like 3, FSTL3) are known as potent antagonists of myostatin that take advantage of its ability to hinder access to signaling receptors on skeletal muscle [35, 36]. The DE genes in this study showed that FST and FSTL1 expression were upregulated by 3.86 and 1.73 fold changes, respectively, in PeM compared with BPR muscles (Additional file 2), while FSTL3 expression was not detected. Function of abundantly expressed MSTN appears to be attenuated by the autocrine or paracrine effects of upregulated FST and its associated FSTL1 in PeM muscles compared to BPR muscles, resulting in the greater muscle growth in PeM chickens. Therefore, the MSTN functions on the muscle growth may be mediated by the interactions between MSTN and its regulatory factors (e.g. FST and FSTLs), not by the differential expression of MSTN alone. The functional roles of MSTN, its downstream effectors, and FST regulation in PeM compared with BPR need to be investigated more in future studies.
The regulations centered with PRKAG3 (protein kinase, AMP-activated, gamma 3 non-catalytic subunit) were closely associated with functions in organismal injury/abnormalities, and cellular assembly/organization (Fig. 6). Although there was no expression difference observed between PeM and BPR, mRNA sequence analysis of PRKAG3 in BPR muscle revealed an allele specific expression at nucleotide position (c.943G > A) that induced a nonsynonymous amino acid change (p. V315M) in 5 out of 6 BPR chickens; however, this was not observed in PeM sequence results (data not shown). Three and two BPR samples showed 50% (heterozygous GA genotype) and 100% mutations (homozygous AA genotype), respectively, for PRKAG3. Therefore, DE genes associated with PRKAG3 appeared to be caused by the differential protein/kinase function. Expression for TPPP3 (tubulin polymerization-promoting protein family member 3), ATP1B1 (ATPase, Na+/K+ transporting, beta 1 polypeptide), BTG2 (BTG family, member 2), KDM2B (lysine-specific demethylase 2B), RGCC (regulator of cell cycle), KLF2 (Kruppel-like factor 2), VIMP (VCP-interacting membrane protein), CD93 (CD93 molecule), SORD (sorbitol dehydrogenase), LMCD1 (LIM and cysteine-rich domains 1), ITIH5 (inter-alpha-trypsin inhibitor heavy chain family, member 5), and CTNNA1 [catenin (cadherin-associated protein), alpha 1, 102 kDa] were increased in PeM compared to BPR, while expression of SLC7A2 [solute carrier family 7 (cationic amino acid transporter, y + system), member 2], ELL2 (elongation factor, RNA polymerase II, 2), NRIP1 (nuclear receptor interacting protein 1) and UGP2 (UDP-glucose pyrophosphorylase 2) were decreased by PRKAG3 in PeM. All of these downstream molecules of PRKAG3 are involved in either a) glucose metabolism (SORD , UGP2 , VIMP ), b) cytoskeletal formation (TPPP3 , CD93 , CTNNA1 ), c) cell cycle regulation (BTG2 , CD93 , RGCC ), d) transcription (ELL2 , KLF2 , LMCD1 ), e) transporters in cell/mitochondrial membranes (ATP1B1 , SLC7A2 ), f) hormone dependent nuclear receptor (NRIP1 ), g) ubiquitination (KDM2B , BTG2 ), or h) inflammatory response (VIMP ). The PRKAG3 associated interactions were derived from a list of DE genes identified by microarray assay that was performed with skeletal muscle tissues of PRKAG3 transgenic mice including missense point mutation and knock-out genotypes . PRKAG3 is the muscle-specific predominant γ isoform of adenosine monophosphate (AMP)-activated protein kinase (AMPK) and has been shown to play an important role in glucose uptake, glycogen synthesis, and fat oxidation in white skeletal muscle (fast twitch fiber) found in pig and mouse models carrying a mutation of the gene [53–56]. Differential expression of various downstream effector molecules of PRKAG3 suggests that muscle specific AMPK pathway may play important roles in metabolism of glucose and lipids with additional cellular functions in the muscle growth, regeneration of the degenerated muscle tissues, differential muscle fiber composition, and feed efficiency phenotypic expression in chicken. This concept is supported by the fact that the most relevant molecular and cellular function identified by IPA analysis was carbohydrate metabolism that involved with 16 molecules (Table 4). AMPK has been known to play important roles to regulate feed efficiency related to mitochondrial functions in chicken muscle [2–4]. AMPK functions in monitoring cell energy status and adjusts ATP production to meet cellular needs. The AMPK is stimulated by increased AMP to ATP ratio and enhances energy production by stimulating mitochondrial biogenesis.
The regulation of LIPE (lipase E, hormone sensitive), which is predicted to be activated in PeM compared to BPR, is involved in cell morphology, connective tissue development/function, and tissue morphology (Fig. 7). LIPE associated regulation was derived from an earlier report which identified DE genes and proteins in soleus muscle tissues of LIPE null mice compared with wild-type littermates using microarray and two dimensional gel electrophoresis based mass-spectrometry proteomics assays . Expressions for RRAD (Ras-related associated with diabetes), ST8SIA4 (ST8 alpha-N-acetyl-neuraminide alpha-2,8-sialyltransferase 4), CDO1 (cysteine dioxygenase type 1), COL1A2 (collagen, type I, alpha 2), YWHAG (tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein, gamma), CSRP3 [cysteine and glycine-rich protein 3 (cardiac LIM protein)], MEOX2 (mesenchyme homeobox 2), BCL6 (B-cell CLL/lymphoma 6), and ANKRD1 [ankyrin repeat domain 1 (cardiac muscle)] were up-regulated in PeM compared to BPR, while expression of RAMP1 [receptor (G protein-coupled) activity modifying protein 1], AQP1 (aquaporin 4), PHKG1 [phosphorylase kinase, gamma 1 (muscle)], ARNTL (aryl hydrocarbon receptor nuclear translocator-like) and SUCLA2 (succinate-CoA ligase, ADP-forming, beta subunit) were downregulated. Hormone-sensitive lipase E is a key enzyme in fatty acid mobilization in adipocytes and muscle LIPE has been shown to be activated by adrenaline-mediated protein kinase A phosphorylation in mammals  as well as by a contraction-induced mechanism, which is independent of protein kinase A . Although LIPE may have critical functions in mammalian muscle, to date, the avian LIPE gene has not been identified in the genome, transcriptome or proteome. Previously, only a partial 86 kDa protein for the LIPE enzyme was purified from chicken adipose tissue and lipase activity was confirmed . Since mammalian LIPE genes showed high sequence similarity, we tried to identify muscle LIPE gene expression in chicken based on mammalian LIPE gene information but this was unfortunately not successful (data not shown). Therefore, detailed functional discussion on avian muscle LIPE will have to wait until the avian LIPE gene is identified.
Androgen receptor (AR)
A steroid hormone receptor, AR (androgen receptor) was predicted to be inhibited in PeM muscle (Fig. 8). AR expression observed from the re-analysis of muscle transcriptome on Ross and Illinois chickens showed a significant downregulation (−2.0 fold, p-value = 0.0003) in Ross chicken muscle at 21 day post hatch compared to Illinois chicken . These results support the inhibition of AR functions in modern broiler muscle tissues compared to unselected chicken breeds. AR is a ligand (androgen)-inducible transcription factor that binds to specific DNA sequences called androgen response elements (AREs) via its DNA-binding domain (DBD) and recruits coactivators that facilitate the transcription of target genes . AR and androgen activate muscle-specific genes encoding structural proteins and myogenic factors . Androgens and AR may activate the estrogen receptor after aromatization  and play a role in anabolic action in skeletal muscle through crosstalk with other signaling molecules . Although MSTN and NMRK2 (nicotinamide riboside kinase 2), which are known to be transcriptionally activated by AR, were not consistent with earlier findings [65, 66], expression of IGF2, CSRP3 [cysteine and glycine-rich protein 3 (cardiac LIM protein)], ATP2A2 (ATPase, Ca++ transporting, cardiac muscle, slow twitch 2), and TNNT2 [troponin T type 2 (cardiac)] were upregulated in AR receptor knock-out mice [65, 66], suggesting that increased expression of these genes in PeM muscle may be caused by the inhibition of AR functions compared to BPR. Predicted inhibition of AR indicated that steroid receptor- mediated transcriptions were down-regulated in modern broiler PeM compared to progenitor BPR muscles.
Summary and conclusion
In this study, we have demonstrated DE genes in breast muscle of modern PeM (rapid growth, high feed efficiency) broiler breeder compared to BPR unselected foundational breed (slow growth, low feed efficacy). Thousands of DE genes involved a variety of cellular and physiological functions in muscle growth and energy metabolism. In addition, many differences reflecting differential muscle fiber composition were characterized in the two chicken breeds. The highly upregulated genes in PeM may represent phenotypes of subclinical myopathy commonly observed in the commercial broiler breast tissue, that can lead to muscle hardening, named as woody breast. Importantly, this study suggests that the fiber composition of avian breast muscle is not consistent with mammalian white muscle in regard to differential mitochondrial content and function. The results of this study provide insights into the genetic alterations accumulated in modern broiler breeder line during selective breeding conducted over several decades.
Kong BW, Song JJ, Lee JY, Hargis BM, Wing T, Lassiter K, Bottje W. Gene expression in breast muscle associated with feed efficiency in a single male broiler line using a chicken 44 K oligo microarray. I. Top differentially expressed genes. Poult Sci. 2011;90(11):2535–47.
Bottje WG, Kong BW, Song JJ, Lee JY, Hargis BM, Lassiter K, Wing T, Hardiman J. Gene expression in breast muscle associated with feed efficiency in a single male broiler line using a chicken 44 K microarray. II. Differentially expressed focus genes. Poult Sci. 2012;91(10):2576–87.
Lee J, Karnuah AB, Rekaya R, Anthony NB, Aggrey SE. Transcriptomic analysis to elucidate the molecular mechanisms that underlie feed efficiency in meat-type chickens. Mol Genet Genomics. 2015;290(5):1673–82.
Davis RV, Lamont SJ, Rothschild MF, Persia ME, Ashwell CM, Schmidt CJ. Transcriptome analysis of post-hatch breast muscle in legacy and modern broiler chickens reveals enrichment of several regulators of myogenic growth. PLoS One. 2015;10(3):e0122525.
Bottje W, Iqbal M, Tang ZX, Cawthon D, Okimoto R, Wing T, Cooper M. Association of mitochondrial function with feed efficiency within a single genetic line of male broilers. Poult Sci. 2002;81(4):546–55.
Mutryn MF, Brannick EM, Fu W, Lee WR, Abasht B. Characterization of a novel chicken muscle disorder through differential gene expression and pathway analysis using RNA-sequencing. BMC Genomics. 2015;16:399-015-1623-0.
Macpherson RE, Vandenboom R, Roy BD, Peters SJ. Skeletal muscle PLIN3 and PLIN5 are serine phosphorylated at rest and following lipolysis during adrenergic or contractile stimulation. Physiol Rep. 2013;1(4):e00084.
Bottje W, Brand MD, Ojano-Dirain C, Lassiter K, Toyomizu M, Wing T. Mitochondrial proton leak kinetics and relationship with feed efficiency within a single genetic line of male broilers. Poult Sci. 2009;88(8):1683–93.
Kong BW, Lassiter K, Piekarski-Welsher A, Dridi S, Reverter-Gomez A, Hudson NJ, Bottje WG. Proteomics of Breast Muscle Tissue Associated with the Phenotypic Expression of Feed Efficiency within a Pedigree Male Broiler Line: I, Highlight on Mitochondria. PLoS One. 2016;11(5):e0155679.
Catania A, Iavarone C, Carlomagno SM, Chiariello M. Selective transcription and cellular proliferation induced by PDGF require histone deacetylase activity. Biochem Biophys Res Commun. 2006;343(2):544–54.
Durieux AC, Amirouche A, Banzet S, Koulmann N, Bonnefoy R, Pasdeloup M, Mouret C, Bigard X, Peinnequin A, Freyssenet D. Ectopic expression of myostatin induces atrophy of adult skeletal muscle by decreasing muscle gene expression. Endocrinology. 2007;148(7):3140–7.
Bonala S, McFarlane C, Ang J, Lim R, Lee M, Chua H, Lokireddy S, Sreekanth P, Leow MK, Meng KC, Shyong TE, Lee YS, Gluckman PD, Sharma M, Kambadur R. Pid1 induces insulin resistance in both human and mouse skeletal muscle during obesity. Mol Endocrinol. 2013;27(9):1518–35.
Hill JJ, Davies MV, Pearson AA, Wang JH, Hewick RM, Wolfman NM, Qiu Y. The myostatin propeptide and the follistatin-related gene are inhibitory binding proteins of myostatin in normal serum. J Biol Chem. 2002;277(43):40735–41.
Lee AY, Chung SK, Chung SS. Demonstration that polyol accumulation is responsible for diabetic cataract by the use of transgenic mice expressing the aldose reductase gene in the lens. Proc Natl Acad Sci U S A. 1995;92(7):2780–4.
Reynolds 4th TH, Pak Y, Harris TE, Manchester J, Barrett EJ, Lawrence Jr JC. Effects of insulin and transgenic overexpression of UDP-glucose pyrophosphorylase on UDP-glucose and glycogen accumulation in skeletal muscle fibers. J Biol Chem. 2005;280(7):5510–5.
Gao Y, Walder K, Sunderland T, Kantham L, Feng HC, Quick M, Bishara N, de Silva A, Augert G, Tenne-Brown J, Collier GR. Elevation in Tanis expression alters glucose metabolism and insulin sensitivity in H4IIE cells. Diabetes. 2003;52(4):929–34.
Vincze O, Tokesi N, Olah J, Hlavanda E, Zotter A, Horvath I, Lehotzky A, Tirian L, Medzihradszky KF, Kovacs J, Orosz F, Ovadi J. Tubulin polymerization promoting proteins (TPPPs): members of a new family with distinct structures and functions. Biochemistry. 2006;45(46):13818–26.
Dean YD, McGreal EP, Akatsu H, Gasque P. Molecular and cellular properties of the rat AA4 antigen, a C-type lectin-like receptor with structural homology to thrombomodulin. J Biol Chem. 2000;275(44):34382–92.
Janssens B, Goossens S, Staes K, Gilbert B, van Hengel J, Colpaert C, Bruyneel E, Mareel M, van Roy F. AlphaT-catenin: a novel tissue-specific beta-catenin-binding protein mediating strong cell-cell adhesion. J Cell Sci. 2001;114(Pt 17):3177–88.
Badea T, Niculescu F, Soane L, Fosbrink M, Sorana H, Rus V, Shin ML, Rus H. RGC-32 increases p34CDC2 kinase activity and entry of aortic smooth muscle cells into S-phase. J Biol Chem. 2002;277(1):502–8.
Bartholomeeusen K, Xiang Y, Fujinaga K, Peterlin BM. Bromodomain and extra-terminal (BET) bromodomain inhibition activate transcription via transient release of positive transcription elongation factor b (P-TEFb) from 7SK small nuclear ribonucleoprotein. J Biol Chem. 2012;287(43):36609–16.
Pei L, Leblanc M, Barish G, Atkins A, Nofsinger R, Whyte J, Gold D, He M, Kawamura K, Li HR, Downes M, Yu RT, Powell HC, Lingrel JB, Evans RM. Thyroid hormone receptor repression is linked to type I pneumocyte-associated respiratory distress syndrome. Nat Med. 2011;17(11):1466–72.
Milan D, Jeon JT, Looft C, Amarger V, Robic A, Thelander M, Rogel-Gaillard C, Paul S, Iannuccelli N, Rask L, Ronne H, Lundstrom K, Reinsch N, Gellin J, Kalm E, Roy PL, Chardon P, Andersson L. A mutation in PRKAG3 associated with excess glycogen content in pig skeletal muscle. Science. 2000;288(5469):1248–51.
Lebret B, Le Roy P, Monin G, Lefaucheur L, Caritez JC, Talmant A, Elsen JM, Sellier P. Influence of the three RN genotypes on chemical composition, enzyme activities, and myofiber characteristics of porcine skeletal muscle. J Anim Sci. 1999;77(6):1482–9.
Ciobanu D, Bastiaansen J, Malek M, Helm J, Woollard J, Plastow G, Rothschild M. Evidence for new alleles in the protein kinase adenosine monophosphate-activated gamma(3)-subunit gene associated with low glycogen content in pig skeletal muscle and improved meat quality. Genetics. 2001;159(3):1151–62.
Barnes BR, Marklund S, Steiler TL, Walter M, Hjalm G, Amarger V, Mahlapuu M, Leng Y, Johansson C, Galuska D, Lindgren K, Abrink M, Stapleton D, Zierath JR, Andersson L. The 5'-AMP-activated protein kinase gamma3 isoform has a key role in carbohydrate and lipid metabolism in glycolytic skeletal muscle. J Biol Chem. 2004;279(37):38441–7.
Hansson O, Donsmark M, Ling C, Nevsten P, Danfelter M, Andersen JL, Galbo H, Holm C. Transcriptome and proteome analysis of soleus muscle of hormone-sensitive lipase-null mice. J Lipid Res. 2005;46(12):2614–23.
Langfort J, Ploug T, Ihlemann J, Baranczuk E, Donsmark M, Gorski J, Galbo H. Additivity of adrenaline and contractions on hormone-sensitive lipase, but not on glycogen phosphorylase, in rat muscle. Acta Physiol Scand. 2003;178(1):51–60.
Estrada M, Liberona JL, Miranda M, Jaimovich E. Aldosterone- and testosterone-mediated intracellular calcium response in skeletal muscle cell cultures. Am J Physiol Endocrinol Metab. 2000;279(1):E132–9.
Dubois V, Laurent M, Boonen S, Vanderschueren D, Claessens F. Androgens and skeletal muscle: cellular and molecular action mechanisms underlying the anabolic actions. Cell Mol Life Sci. 2012;69(10):1651–67.
Dubois V, Laurent MR, Sinnesael M, Cielen N, Helsen C, Clinckemalie L, Spans L, Gayan-Ramirez G, Deldicque L, Hespel P, Carmeliet G, Vanderschueren D, Claessens F. A satellite cell-specific knockout of the androgen receptor reveals myostatin as a direct androgen target in skeletal muscle. FASEB J. 2014;28(7):2979–94.
MacLean HE, Chiu WS, Notini AJ, Axell AM, Davey RA, McManus JF, Ma C, Plant DR, Lynch GS, Zajac JD. Impaired skeletal muscle development and function in male, but not female, genomic androgen receptor knockout mice. FASEB J. 2008;22(8):2676–89.
Wang H, Hu L, Dalen K, Dorward H, Marcinkiewicz A, Russell D, Gong D, Londos C, Yamaguchi T, Holm C, Rizzo MA, Brasaemle D, Sztalryd C. Activation of hormone-sensitive lipase requires two steps, protein phosphorylation and binding to the PAT-1 domain of lipid droplet coat proteins. J Biol Chem. 2009;284(46):32116–25.
Tsukada T, Pappas CT, Moroz N, Antin PB, Kostyukova AS, Gregorio CC. Leiomodin-2 is an antagonist of tropomodulin-1 at the pointed end of the thin filaments in cardiac muscle. J Cell Sci. 2010;123(Pt 18):3136–45.
Li J, Chen Y, Wang YG, Zhao XL, Gilbert ER, Liu YP, Wang Y, Hu YD, Zhu Q. MUSTN1 mRNA Abundance and Protein Localization is Greatest in Muscle Tissues of Chinese Meat-Quality Chickens. Int J Mol Sci. 2013;14(3):5545–59.
Tobacman LS, Nihli M, Butters C, Heller M, Hatch V, Craig R, Lehman W, Homsher E. The troponin tail domain promotes a conformational state of the thin filament that suppresses myosin activity. J Biol Chem. 2002;277(31):27636–42.
Vos MJ, Zijlstra MP, Kanon B, van Waarde-Verhagen MA, Brunt ER, Oosterveld-Hut HM, Carra S, Sibon OC, Kampinga HH. HSPB7 is the most potent polyQ aggregation suppressor within the HSPB family of molecular chaperones. Hum Mol Genet. 2010;19(23):4677–93.
Matsuno K, Yamada H, Iwata K, Jin D, Katsuyama M, Matsuki M, Takai S, Yamanishi K, Miyazaki M, Matsubara H, Yabe-Nishimura C. Nox1 is involved in angiotensin II-mediated hypertension: a study in Nox1-deficient mice. Circulation. 2005;112(17):2677–85.
Alfonso J, Fernandez ME, Cooper B, Flugge G, Frasch AC. The stress-regulated protein M6a is a key modulator for neurite outgrowth and filopodium/spine formation. Proc Natl Acad Sci U S A. 2005;102(47):17196–201.
Wang X, Ying W, Dunlap KA, Lin G, Satterfield MC, Burghardt RC, Wu G, Bazer FW. Arginine decarboxylase and agmatinase: an alternative pathway for de novo biosynthesis of polyamines for development of mammalian conceptuses. Biol Reprod. 2014;90(4):84.
Li X, Liu P, Liu W, Maye P, Zhang J, Zhang Y, Hurley M, Guo C, Boskey A, Sun L, Harris SE, Rowe DW, Ke HZ, Wu D. Dkk2 has a role in terminal osteoblast differentiation and mineralized matrix formation. Nat Genet. 2005;37(9):945–52.
Jiang Z, Song J, Qi F, Xiao A, An X, Liu NA, Zhu Z, Zhang B, Lin S. Exdpf is a key regulator of exocrine pancreas development controlled by retinoic acid and ptf1a in zebrafish. PLoS Biol. 2008;6(11):e293.
Nakamura Y, Handa K, Iwamoto R, Tsukamoto T, Takahasi M, Mekada E. Immunohistochemical distribution of CD9, heparin binding epidermal growth factor-like growth factor, and integrin alpha3beta1 in normal human tissues. J Histochem Cytochem. 2001;49(4):439–44.
We thank CSIRO and Cell and Molecular Biology Graduate Program of University of Arkansas to support the licenses for CLC Genomics Workbench and JMP Genomics 6, respectively. We also thank Toni Reverter at CSIRO for statistical support.
USDA-NIFA (#2013–01953), Arkansas Biosciences Institute (Little Rock, AR) and the Agricultural Experiment Station (University of Arkansas, Fayetteville).
Availability of data and materials
Datasets supporting the results of this article are also included in the Additional files 1, 2, 4, and 5, Additional file 3: Figure S1.
BWK and WB designed the experiments, performed the experiments, analyzed the data, bioinformatic analyses and wrote the manuscript. NH analyzed RNAseq data and contribute to edit the manuscript preparation. SL analyzed RNAseq data. DS and BK analyzed qPCR assay. KL, DC, and AP contribute to analyze sample collection, perform experiments, qPCR and bioinformatic analyses. SD and NA care the animal handling, sample collection, and data analyses. All authors read and approved the final manuscript.
BWK, SL, DS, BK, KL, DC, AP, SD, NA and WB: Department of Poultry Science, Center of Excellence for Poultry Science, University of Arkansas, Fayetteville, Arkansas, NH: School of Agriculture and Food Science, University of Queensland St. Lucia, AU.
The authors declare that they have no competing interests.
Consent for publication
Ethics approval and consent to participate
All procedures for animal care complied with the University of Arkansas Institutional Animal Care and Use Committee (IACUC): Protocol #14012.
Authors and Affiliations
Department of Poultry Science, Center of Excellence for Poultry Science, University of Arkansas, Fayetteville, Arkansas, USA
Byung-Whi Kong, Dongwon Seo, Seok Lee, Bhuwan Khatri, Kentu Lassiter, Devin Cook, Alissa Piekarski, Sami Dridi, Nicholas Anthony & Walter Bottje
School of Agriculture and Food Science, University of Queensland, Gatton, Australia
The canonical pathway of oxidative phosphorylation. The differentially expressed genes in breast muscle associated with the electron transport chain on (Complex I, II, III, IV, and V) that were downregulated (outlined in green) in the PeM. Pumping of hydrogen ions (H+) creates a proton motive force between the inner (IMM) and outer (OMM) mitochondrial membranes that is used to drive ATP synthesis. (PPTX 254 kb)
List of focus molecules in gene networks and upstream analysis. Gene symbols, full gene name and fold change values were displayed for the illustrations of network and upstream analyses. (XLSX 14 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Kong, BW., Hudson, N., Seo, D. et al. RNA sequencing for global gene expression associated with muscle growth in a single male modern broiler line compared to a foundational Barred Plymouth Rock chicken line.
BMC Genomics18, 82 (2017). https://doi.org/10.1186/s12864-016-3471-y