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 Genomics volume 18, Article number: 82 (2017)
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.
Angiotensin I converting enzyme
Angiotensin II type 1 receptor
Adenosine monophosphate-activated protein kinase
Ankyrin repeat domain 1
Androgen response elements
Aryl hydrocarbon receptor nuclear translocator-like
ATPase Na+/K+ transporting beta 1 polypeptide
ATPase Ca++ transporting cardiac muscle slow twitch 2
B-cell CLL/lymphoma 6
Barred Plymouth Rock
BTG family member 2
Carbonic anhydrase 3
Chemokine C-C motif ligand 5
Cysteine dioxygenase type 1
Cellular FBJ murine osteosarcoma viral oncogene homolog
UDP glycosyltransferase 8
Cellular Jun protooncogene
Creatine kinase mitochondrial 1A
Collagen type I alpha 2
Cysteine and glycine-rich protein 3
Catenin alpha 1
Dual specificity phosphatase 1
Elongation factor 2 RNA polymerase II
The European Molecular Biology Laboratory
Fatty acid binding protein 4
Dickkopf WNT signaling pathway inhibitor 2
False discovery rate
Fibrinogen beta chain
Four and a half LIM domains 2
Follistatin like 1
Follistatin like 3
Gamma-aminobutyric acid (GABA) A receptor delta
Glyceraldehyde 3-phosphate dehydrogenase
Glycinamide ribonucleotide transformylase; Phosphoribosylglycinamide formyltransferase
Gene Expression Omnibus
Glutathione S-transferase alpha 3
- HO-1 or HMOX1:
Heme oxygenase 1
Institutional Animal Care and Use Committee
Joining chain of multimeric IgA and IgM
Insulin like growth factor 2
Interleukin 23 receptor
Ingenuity Pathway Analysis
Inter-alpha-trypsin inhibitor heavy chain family member 5
Lysine-specific demethylase 2B
Kruppel-like factor 2
Kruppel-like factor 9
Kelch-like family member 38
LIM and cysteine-rich domains 1
V-maf avian musculoaponeurotic fibrosarcoma oncogene homolog
Mitogen activated protein kinase kinase 4
Mitogen-activated protein kinase kinase kinase 1
MAPK/ERK kinase kinase
Mesenchyme homeobox 2
Microsomal glutathione S-transferase 1
Mitogen activated protein kinase 4
Musculoskeletal embryonic nuclear protein 1
Myosin binding protein C 1 slow type
Myosin heavy chain 15
Nuclear factor, erythroid 2-like 2
Nicotinamide riboside kinase 2
NADPH oxidase organizer 1
Nuclear related factor 2
Nuclear receptor interacting protein 1
Proprotein convertase subtilisin/kexin type 1
Platelet derived growth factor
Platelet derived growth factor receptor
Phosphorylase kinase gamma 1
Protein kinase C
Phosphotyrosine interaction domain containing 1
Phospholipase C gamma
- PPARGC1A or PGC1α:
Peroxisome proliferator-activated receptor gamma coactivator 1 alpha
Peroxisome proliferator-activated receptor gamma, coactivator 1 beta
Pancreatic progenitor cell differentiation and proliferation factor
Protein kinase AMP-activated gamma 3 non-catalytic subunit
Protein interacting with cyclin A1
Quantitative reverse transcription polymerase chain reaction
Receptor (G protein-coupled) activity modifying protein 1
Residual feed intake
Regulator of cell cycle
Ras homolog family member F
Reactive oxygen species
Reads per kilobase per million
Reads per million
Ras-related associated with diabetes
Sodium channel voltage gated type V alpha subunit
Serpin peptidase inhibitor clade H heat shock protein 47 member 1
Solute carrier family 38 member 4
Solute carrier family 7 member 2
Single nucleotide polymorphisms
Superoxide reductase 3
ST8 alpha-N-acetyl-neuraminide alpha-2,8-sialyltransferase 4
Signal transducer and activator of transcription 1
Signal transducer and activator of transcription 3
Succinate-CoA ligase ADP-forming beta subunit
Transforming growth factor beta
Troponin I 1
Troponin T type 2 cardiac
Tubulin polymerization-promoting protein family member 3
- TRXR1 or TXNRD1:
Thioredoxin reductase 1
UDP-glucose pyrophosphorylase 2
Vascular cell adhesion molecule 1
VCP-interacting membrane protein
Tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein gamma
Niemann H, Kuhla B, Flachowsky G. Perspectives for feed-efficient animal production. J Anim Sci. 2011;89(12):4344–63.
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.
Bottje W, Kong BW. Cell Biology Symposium: feed efficiency: mitochondrial function to global gene expression. J Anim Sci. 2013;91(4):1582–93.
Zhou N, Lee WR, Abasht B. Messenger RNA sequencing and pathway analysis provide novel insights into the biological basis of chickens' feed efficiency. BMC Genomics. 2015;16:195-015-1364-0.
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.
Dorshorst BJ, Ashwell CM. Genetic mapping of the sex-linked barring gene in the chicken. Poult Sci. 2009;88(9):1811–7.
Lopez G, de Lange K, Leeson S. Partitioning of retained energy in broilers and birds with intermediate growth rate. Poult Sci. 2007;86(10):2162–71.
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.
Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008;5(7):621–8.
Hudson NJ, Dalrymple BP, Reverter A. Beyond differential expression: the quest for causal mutations and effector molecules. BMC Genomics. 2012;13:356.
Hudson NJ, Reverter A, Wang Y, Greenwood PL, Dalrymple BP. Inferring the transcriptional landscape of bovine skeletal muscle by integrating co-expression networks. PLoS One. 2009;4(10):e7249.
Hudson NJ. Symmorphosis and livestock bioenergetics: production animal muscle has low mitochondrial volume fractions. J Anim Physiol Anim Nutr. 2009;93(1):1–6.
Mishra P, Varuzhanyan G, Pham AH, Chan DC. Mitochondrial Dynamics Is a Distinguishing Feature of Skeletal Muscle Fiber Types and Regulates Organellar Compartmentalization. Cell Metab. 2015;22(6):1033–44.
Kiessling KH. Muscle structure and function in the goose, quail, pheasant, guinea hen, and chicken. Comp Biochem Physiol B. 1977;57(4):287–92.
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.
Sihvo HK, Immonen K, Puolanne E. Myodegeneration with fibrosis and regeneration in the pectoralis major muscle of broilers. Vet Pathol. 2014;51(3):619–23.
Seidah NG, Khatib AM, Prat A. The proprotein convertases and their implication in sterol and/or lipid metabolism. Biol Chem. 2006;387(7):871–7.
Seidah NG, Prat A. The biology and therapeutic targeting of the proprotein convertases. Nat Rev Drug Discov. 2012;11(5):367–83.
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.
Cabello-Verrugio C, Cordova G, Salas JD. Angiotensin II: role in skeletal muscle atrophy. Curr Protein Pept Sci. 2012;13(6):560–9.
Cook DR, Doumit ME, Merkel RA. Transforming growth factor-beta, basic fibroblast growth factor, and platelet-derived growth factor-BB interact to affect proliferation of clonally derived porcine satellite cells. J Cell Physiol. 1993;157(2):307–12.
Doumit ME, Cook DR, Merkel RA. Fibroblast growth factor, epidermal growth factor, insulin-like growth factors, and platelet-derived growth factor-BB stimulate proliferation of clonally derived porcine myogenic satellite cells. J Cell Physiol. 1993;157(2):326–32.
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.
Yu Y, Sweeney M, Zhang S, Platoshyn O, Landsberg J, Rothman A, Yuan JX. PDGF stimulates pulmonary vascular smooth muscle cell proliferation by upregulating TRPC6 expression. Am J Physiol Cell Physiol. 2003;284(2):C316–30.
Hayes JD, McMahon M. The double-edged sword of Nrf2: subversion of redox homeostasis during the evolution of cancer. Mol Cell. 2006;21(6):732–4.
McPherron AC, Lawler AM, Lee SJ. Regulation of skeletal muscle mass in mice by a new TGF-beta superfamily member. Nature. 1997;387(6628):83–90.
Gilson H, Schakman O, Combaret L, Lause P, Grobet L, Attaix D, Ketelslegers JM, Thissen JP. Myostatin gene deletion prevents glucocorticoid-induced muscle atrophy. Endocrinology. 2007;148(1):452–60.
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.
Li ZB, Kollias HD, Wagner KR. Myostatin directly regulates skeletal muscle fibrosis. J Biol Chem. 2008;283(28):19371–8.
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.
Lee SJ. Quadrupling muscle mass in mice by targeting TGF-beta signaling pathways. PLoS One. 2007;2(8):e789.
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.
Sasajima H, Nakagawa K, Yokosawa H. Antiproliferative proteins of the BTG/Tob family are degraded by the ubiquitin-proteasome system. Eur J Biochem. 2002;269(14):3596–604.
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.
Rath N, Wang Z, Lu MM, Morrisey EE. LMCD1/Dyxin is a novel transcriptional cofactor that restricts GATA6 function by inhibiting DNA binding. Mol Cell Biol. 2005;25(20):8864–73.
Liu L, Askari A. Beta-subunit of cardiac Na + −K + −ATPase dictates the concentration of the functional enzyme in caveolae. Am J Physiol Cell Physiol. 2006;291(4):C569–78.
Nicholson B, Manner CK, Kleeman J, MacLeod CL. Sustained nitric oxide production in macrophages requires the arginine transporter CAT2. J Biol Chem. 2001;276(19):15881–5.
Hannafon BN, Sebastiani P, de las Morenas A, Lu J, Rosenberg CL. Expression of microRNA and their gene targets are dysregulated in preinvasive breast cancer. Breast Cancer Res. 2011;13(2):R24.
Wu X, Johansen JV, Helin K. Fbxl10/Kdm2b recruits polycomb repressive complex 1 to CpG islands and regulates H2A ubiquitylation. Mol Cell. 2013;49(6):1134–46.
Nilsson EC, Long YC, Martinsson S, Glund S, Garcia-Roves P, Svensson LT, Andersson L, Zierath JR, Mahlapuu M. Opposite transcriptional regulation in skeletal muscle of AMP-activated protein kinase gamma3 R225Q transgenic versus knock-out mice. J Biol Chem. 2006;281(11):7244–52.
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.
Donsmark M, Langfort J, Holm C, Ploug T, Galbo H. Hormone-sensitive lipase as mediator of lipolysis in contracting skeletal muscle. Exerc Sport Sci Rev. 2005;33(3):127–33.
Anthonsen MW, Degerman E, Holm C. Partial purification and identification of hormone-sensitive lipase from chicken adipose tissue. Biochem Biophys Res Commun. 1997;236(1):94–9.
Claessens F, Denayer S, Van Tilborgh N, Kerkhofs S, Helsen C, Haelens A. Diverse roles of androgen receptor (AR) domains in AR-mediated signaling. Nucl Recept Signal. 2008;6:e008.
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.
Callewaert F, Sinnesael M, Gielen E, Boonen S, Vanderschueren D. Skeletal sexual dimorphism: relative contribution of sex steroids, GH-IGF1, and mechanical loading. J Endocrinol. 2010;207(2):127–34.
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.
Brasaemle DL, Dolios G, Shapiro L, Wang R. Proteomic analysis of proteins associated with lipid droplets of basal and lipolytically stimulated 3 T3-L1 adipocytes. J Biol Chem. 2004;279(45):46835–42.
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.
Conley CA. Leiomodin and tropomodulin in smooth muscle. Am J Physiol Cell Physiol. 2001;280(6):C1645–56.
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.
Sly WS, Hu PY. Human carbonic anhydrases and carbonic anhydrase deficiencies. Annu Rev Biochem. 1995;64:375–401.
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.
Liu C, Gersch RP, Hawke TJ, Hadjiargyrou M. Silencing of Mustn1 inhibits myogenic fusion and differentiation. Am J Physiol Cell Physiol. 2010;298(5):C1100–8.
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.
Franck E, Madsen O, van Rheede T, Ricard G, Huynen MA, de Jong WW. Evolutionary diversity of vertebrate small heat shock proteins. J Mol Evol. 2004;59(6):792–805.
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.
Banfi B, Clark RA, Steger K, Krause KH. Two novel proteins activate superoxide generation by the NADPH oxidase NOX1. J Biol Chem. 2003;278(6):3510–3.
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.
Kranias EG, Hajjar RJ. Modulation of cardiac contractility by the phospholamban/SERCA2a regulatome. Circ Res. 2012;110(12):1646–60.
Winegrad S. Cardiac myosin binding protein C. Circ Res. 1999;84(10):1117–26.
Chen Z, Zhao TJ, Li J, Gao YS, Meng FG, Yan YB, Zhou HM. Slow skeletal muscle myosin-binding protein-C (MyBPC1) mediates recruitment of muscle-type creatine kinase (CK) to myosin. Biochem J. 2011;436(2):437–45.
O'Sullivan LA, Liongue C, Lewis RS, Stephenson SE, Ward AC. Cytokine receptor signaling through the Jak-Stat-Socs pathway in disease. Mol Immunol. 2007;44(10):2497–506.
Lagenaur C, Kunemund V, Fischer G, Fushiki S, Schachner M. Monoclonal M6 antibody interferes with neurite extension of cultured neurons. J Neurobiol. 1992;23(1):71–88.
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.
Razzaque MS, Lanske B. The emerging role of the fibroblast growth factor-23-klotho axis in renal regulation of phosphate homeostasis. J Endocrinol. 2007;194(1):1–10.
Olauson H, Lindberg K, Amin R, Jia T, Wernerson A, Andersson G, Larsson TE. Targeted deletion of Klotho in kidney distal tubule disrupts mineral metabolism. J Am Soc Nephrol. 2012;23(10):1641–51.
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.
Mao B, Wu W, Li Y, Hoppe D, Stannek P, Glinka A, Niehrs C. LDL-receptor-related protein 6 is a receptor for Dickkopf proteins. Nature. 2001;411(6835):321–5.
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.
Mellon SH, Griffin LD. Neurosteroids: biochemistry and clinical significance. Trends Endocrinol Metab. 2002;13(1):35–43.
de Moerloose P, Casini A, Neerman-Arbez M. Congenital fibrinogen disorders: an update. Semin Thromb Hemost. 2013;39(6):585–95.
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
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.
Primers used for qRT-PCR. The first column indicates primer names and the second column shows sequences for the forward and reverse primers. (XLSX 9 kb)
List of entire 2464 DE genes in the comparison of PeM and BPR. (XLSX 141 kb)
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 DE genes in TCA cycles and oxidative phosphorylation. Gene symbols, full gene name and fold change values were displayed. (XLSX 11 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)
About this article
Cite this article
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 Genomics 18, 82 (2017). https://doi.org/10.1186/s12864-016-3471-y
- Androgen Receptor
- Differentially Express
- Ingenuity Pathway Analysis
- Differentially Express Gene
- Feed Efficiency