BACKGROUND
in a resource-constrained world supporting a rapidly growing human population there is great interest in enhancing the production efficiency of our major animal and plant food industries.
feed is the single largest commercial cost in animal production  <cit> .
any inefficiency not only affects the bottom line but also negatively impacts resource usage  and waste production .
consequently, there is considerable economic and environmental interest in increasing animal feed conversion efficiency.
small increases in feed efficiency across numerous animals can have a large industrial and environmental impact.
broiler chickens are the most feed efficient of all the vertebrates, with some companies boasting a conversion of just ~ <dig>  units dry weight intake per unit wet weight gain  <cit> .
this makes broilers particularly valuable as biological models for understanding the mechanistic basis of feed efficiency.
along with other avian meat producing species, cf.
turkeys, pheasants, partridge, grouse and quails, chickens are members of the phasianidae taxonomic clade, the largest branch of the galliformes.
the members of this group tend to be sedentary, resident ground-dwelling birds that use short, burst flights to escape predators  <cit> .
this behaviour explains the functionally unusual breast muscle metabolism of the ancestral birds such as the red junglefowl  progenitor of modern domestic chickens, which is dominated by explosive, fast twitch contractile isoforms and a relatively low oxidative capacity metabolism.
the birds we chose for this analysis are pedigree broiler males  from a population of highly inbred, commercial animals subject to very strong historical selection for high muscle mass and elevated feed efficiency   <cit> .
the significant difference in fe  between the two treatment groups is in the order of  <dig>  fold in absolute terms.
some of the variation in the fe phenotype presumably derives from biological modifications in the non-muscle tissues.
this logic, coupled with the prior application of microarray technology on the same individuals  <cit> , gave us the expectation that any molecular differences we may detect in the breast muscle between the two groups would be fairly subtle.
nevertheless, it has been established for some time that commercial broilers possess a muscle structure that is, histologically speaking, composed almost exclusively of low mitochondrial content type iib glycolytic fibres  <cit> .
this fibre homogeneity is advantageous as it provides a very stable transcriptional background against which to perform the various molecular analyses.
here, we have used illumina rna sequencing technology to screen the breast muscle transcriptomes of the high fe  and low fe  groups.
we aimed to both test existing hypotheses as well as generate new ones.
we have previously argued for a role of both isolated mitochondrial physiology  and overall tissue mitochondrial content  <cit>  in driving variation in production animal feed efficiency.
therefore, as part of our analyses we have specifically harvested the data for genes encoding mitochondrial proteins, finding expression was skewed towards the hfe birds.
furthermore, we have also used a combination of basic differential expression  and a recently developed differential connectivity  algorithm called regulatory impact factor  analysis to screen the data in an unbiased fashion .
the rif algorithm detects regulatory molecules that are transcriptionally rewired between the two states, irrespective of whether that molecule is de itself.
this is a potent approach for the identification of transcription factors that tend to be lowly and stably expressed at the mrna level  <cit> , but whose functional activity is controlled at the protein level through changes in cellular localisation, ligand binding or post-translational modifications.
this analysis points to progesterone signalling through the progesterone receptor as a key molecular driver.
methods
animal resources and phenotypes
we used breast skeletal muscle samples from  <dig> male broilers divergent in feed efficiency as previously described  <cit> .
male broilers were raised under standard industry practices and euthanized using carbon dioxide asphyxiation.
cobb-vantress inc. adheres to internal industrial standard guidelines for animal welfare and has standard operating procedures for euthanasia, dissection and other experimental techniques.
in brief, the  <dig> broilers were derived from a single genetic line of a commercial cobb-vantress inc. population.
both the hfe  and lfe  birds were sub-sampled from the extremes of an original group of  <dig> birds.
in turn, those  <dig> birds were selected as the most efficient from a larger population of  <dig> birds.
the high fe group had greater body weight gain from  <dig> to  <dig> weeks but consumed the same amount of feed as the lfe group.
for the purposes of this manuscript we have defined fe as the gain in bird mass relative to feed intake, such that a higher number equates to a more efficient bird.
in our case, the hfe birds yielded efficiencies of  <dig>  ±  <dig>  compared to  <dig>  ±  <dig>   for the lfe birds.
this animal resource is noteworthy in the sense that the two groups are not divergently selected.
this is an important distinction compared to other research populations in this area.
for example, the observed de in divergent selection experiments could be a consequence of functionally irrelevant founder effects that have accumulated over subsequent matings.
in our case the birds are from the same line and generation, but separated on the basis of 1 week of fcr testing in a controlled experimental environment.
mapping, counting and normalising the mrna reads
after phenol chloroform rna extraction we submitted our samples to the research support facility at michigan state university  for illumina hiseq  <dig> base pair paired end read sequencing.
in brief, we used clc genomics workbench  <dig> to map the reads to gallus gallus genome assembly version  <dig>  the clc software follows the analytical pipeline recommended in  <cit> .
we log <dig> transformed the rpm data to stabilise the variance and then performed a further quantile normalisation.
the impact of the two levels of normalisation on whether the  <dig> individual samples could discriminate into the treatment groups of origin was confirmed by principal component analysis .
differential expression
to compare the  <dig> treatment groups we plotted the ma  plot comparing hfe birds to lfe birds .
the x axis, a, represents the average abundance in the  <dig> groups, and the y axis, m, is hfe average transcript abundance minus lfe average transcript abundance.
a subset of functionally relevant outlier de genes was chosen for annotation on the plot.fig.
<dig>
a the slow contractile subunits highlighted are myh <dig>  tpm <dig>  myoz <dig>  tnni <dig>  myl <dig>  tnnt <dig> and csrp <dig> consistent with a shift in fibre composition towards the fast glycolytic and away from the slow aerobic fibres in hfe birds.
b the mrna encoding the mitoproteome is significantly skewed towards the hfe consistent with a slight increase in mitochondrial content.
c the extreme 5% pif used as input for the rif analysis have been highlighted.
these de genes are distributed evenly across the ma plot, closely tracking the overall shape of the plot
mrna encoding the mitoproteome
we downloaded the complete human mitoproteome from http://mitominer.mrc-mbu.cam.ac.uk/release- <dig> /begin.do generating a list of  <dig> nuclear and mitochondrial genes encoding mitochondrial proteins.
we converted the protein names to gene names, and found  <dig> matches in our chicken rnaseq data.
the mitoproteome was then plotted in ma format, and all those with higher values in the hfe group were colour coded blue and those with lower values in the hfe colour coded red .
the skew in distribution away from the null expectation of 50: <dig> was quantified by binomial statistics.
this approach formalises the extent to which the mitoproteome data is biased to one or the other of the two groups.
phenotypic impact factor
we next computed a modified de called pif  which we have found to have a number of appealing characteristics, both numerical and biological.
from a numerical perspective it can be used to establish extreme de genes in a way that accounts for the structure of the distribution of the data .
it de-emphasises lowly abundant genes that are inherently noisy as they approach the detection limit of the technology.
it also draws together the  <dig> major sources of variation in gene expression data into a single metric which can be used for comparison purposes.
from a biological perspective pif analysis is useful as it gives more meaningful functional enrichments than ranking on de, including better characterising muscle fibre composition shifts in muscle transcriptome data  <cit> .
also, pif is a foundation of the rif differential network analysis  where the target de molecules are first defined by extreme pif, and then the regulators differential co-expression to those targets is integrated into a differential connectivity metric.
we ranked on pif and assessed functional enrichment at the extremes of the list using gorilla  <cit> .
vmus3d muscle model
to visualise spatial locations of proteins encoded by de genes derived from the hfe versus lfe birds we utilised the virtual muscle 3d .
vmus3d, described in detail here  <cit> , is an annotated 3-dimensional representation of key contractile protein complexes, including thick, thin, z-disc and costameric complexes, and their relative protein spatial locations.
vmus3d is a database driven visualisation tool that enables mapping of gene expression changes onto their encoded proteins, represented as an extensible 3d object, via the document object model.
here, we mapped m-values onto the vmus3d and coloured change as follows: blue , green , yellow , orange  and red  .fig.
<dig> 3d muscle model visualisation illustrating the downregulation of slow subunits in hfe birds consistent with a whiter, sprint muscle phenotype.
a. longitudinal view of the contractile apparatus with coordinated de changes of slow unit proteins highlighted.
b. z-disk proteins differentially expressed.
c. protein locations across the contractile apparatus and their corresponding values.
colours correspond to the level of de  as follows: blue , green , yellow , orange  and red
permut matrix clustering on extreme pif
the top  <dig> up and top  <dig> downregulated genes  were identified as described above.
a matrix was formed comprising as many rows as genes  and as many columns as birds  with each cell containing the normalised gene expression.
this matrix was imported into permutmatrix software  <cit>  for hierarchical clustering .fig.
<dig> a hierarchically clustered heatmap showing the expression patterns of the  <dig> most extreme pif mrna.
this panel of  <dig> genes which includes mrna encoding muscle structural proteins and mitochondrial proteins can correctly discriminate the two treatment groups
regulatory impact factor  analysis
the purpose of rif is to identify regulatory molecules whose activity can change independent of any change in gene expression level.
this includes transcription factors whose activity may be influenced by ligand binding, co-factor binding, cellular localisation and other post-translational processes.
the method works by establishing different patterns of network connectivity in the two states .
rif <dig> and rif <dig> are two versions of essentially the same analysis.
in both cases the abundance and differential expression of the ‘target genes’ is exploited in conjunction with the differential co-expression of the ‘regulators’ to those ‘targets.’ the output of both versions has been plotted and reported here.
we computed rif <dig> and rif <dig> as previously described .
this procedure exploits global patterns of differential co-expression  to infer those regulatory molecules whose behaviour is systematically different in a contrast of interest, in this case the hfe versus lfe birds.
herein, the experimental contrast was hfe vs. lfe and the rif metrics for the r-th regulator  were computed using the following formulae: rif1r=1nde∑j=1j=ndexj×dj×dcrj <dig>
and rif2r=1nde∑j=1j=ndexjhfe×rrjhfe2−xjlfe×rrjlfe <dig> where nde represented the number of de genes; x
j was the average expression of the j-th de gene across all time points; d
j was the de of the j-th gene in the hfe vs. lfe contrast; dcrj was the differential co-expression between the r-th regulator and the j-th de gene, and computed from the difference between r
rjhfe and r
rjlfe, the correlation co-expression between the r-th regulator and the j-th de gene in the hfe and lfe samples, respectively; finally, x
jhfe and x
jlfe represented the average expression of the j-th de or ts gene in the hfe and lfe samples, respectively.
the rif <dig> and rif <dig> output were simultaneously plotted .fig.
<dig> the extremely differentially connected tf as illustrated by rif <dig> and rif <dig> scores.
progesterone receptor and  <dig> other tf involved in progesterone signalling are awarded extreme scores by both metrics
quantitative pcr technical validation
to independently validate our rnaseq normalisation method and method to detect de, we selected a subset of nine representative genes for qpcr amplification.
these were mrna encoding slow twitch subunits , slow muscle associated proteins  and muscle fat content .
we used pcr cycling conditions as previously described  <cit> .
primers are detailed in additional file  <dig>
in brief, one microgram of total rna was obtained from  <dig> muscle samples each for high and low fe pedigree broilers for general validation of the rnaseq results and for specific confirmation.
rna was converted into cdna with qscripttm cdna supermix  following the manufacturer’s instructions.
the cdna samples were diluted in a 1: <dig> ratio and a portion  of the cdna was subjected to a real-time quantitative pcr  reaction using an abi prism 7500ht system  with poweruptm sybr® green master mix .
the specific oligonucleotide primers were designed by the primer <dig> program .
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  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  <dig> times and the values of average cycle threshold  were determined for each sample, and 2 − ΔΔct values for the comparison of hfe and lfe muscles were used for relative quantification by fold-change and statistical significance.
quail muscle  cell culture and treatments
qm <dig> cells were grown in m <dig> medium  with 10% fetal bovine serum , 10% tryptose phosphate , and 1% penicillin-streptomycin  at 37 °c under a humidified atmosphere of 5% co <dig> and 95% air.
at 80–90% confluence, cells were synchronized overnight in serum-free medium.
progesterone receptor immunofluorescence
immunofluorescence was performed as previously described  <cit>  with modifications.
briefly, cells were grown to 50–60% confluence in chamber slides  followed by methanol fixation for 10 min at 20 °c.
serum free blocking buffer , and the cells incubated with rabbit anti-progesterone receptor antibody overnight at 4 °c.
cells were treated with mitotracker deep red  and visualized with alexa fluor 488-conjugated secondary antibody .
after dapi counterstaining, a cover slip was placed on slides with vectashield .
images were obtained using the zeiss imager m <dig> with a 20x plan-apochromat 20x/ <dig>  objective and a 100x ec planneofluor 100x/ <dig>  oil objective.
the alexa fluor  <dig> fluorophore was observed through filter set  <dig> 1031- <dig> with an excitation of bp 470/ <dig>  beamsplitter of ft  <dig>  and emission spectrum of bp 525/ <dig>  differential interference contrast images were collected using dic m <dig> condensers.
the alexa fluor  <dig> fluorophore was excited for 500 ms prior to capturing each image using an axio cam mr <dig> camera.
all analysis was performed using axiovision se <dig>  <dig> .
<dig> sp <dig> software .
statistical analysis
given the two treatment groups are so similar and that multiple testing penalties can be very severe in this kind of experimental design we elected not to identify significantly de genes on a gene-by-gene basis.
rather, we used a range of systems-wide analyses to interpret our data.
this approach has the advantage of not being reliant on the performance of any given gene, but rather draws on information present in entire functional groups of genes.
for basic functional enrichment analysis we used gorilla  <cit> .
here, we examined functional enrichment in the extreme  <dig> up and downregulated genes  genes ranked by pif, compared to a background list consisting of all genome-wide data; both p and fdr q-values reported.
this nominal 5% de genes also represent the target de genes for the rif analysis.
for the mitoproteome analysis we used binomial statistics to quantify the skew for all mrna encoding mitochondrial proteins.
for the rif analysis we report on extreme outlier regulators only.
RESULTS
twelve rnaseq libraries were constructed using breast muscle rna samples from  <dig> hfe and  <dig> lfe broilers.
in total,  <dig> million 100 bp sequences were obtained with an average of  <dig> million reads per sample and 80% of the reads were mapped to the g. gallus reference genome assembly.
after mapping, counting and normalising we returned data for  <dig>  genes.
we screened for those mrna with no missing values across the  <dig> samples.
this gave a total  <dig>  genes which were used for all subsequent analysis.
both rpm and rpkm values were computed and compared by correlation analysis .
to achieve this we calculated the pif for each gene  and correlated the pif scores by the two methods.
given the high correlation  and the concern that the additional normalisation used in generating rpkm values can theoretically introduce a variety of problems  <cit> , we elected to proceed with the normalised rpm values.
principal component analysis
the raw read count data could not discriminate the birds into treatment group of origin.
log <dig> stabilising the data allowed separation of individuals into the two groups with one exception.
this separation by treatment group was reinforced by a subsequent quantile normalisation which standardises the distributions of the reads for each sample .
this separation provided a quality check that our normalisation approach had maintained the treatment differences without introducing any systematic bias.
further validation that the normalisation approach did not introduce bias was provided by determination that the ma plot was centred on  <dig>
differential expression between hfe and lfe birds
we used a modified de, named pif for: 1) functional enrichment analysis in gorilla, 2) to generate a panel of genes that can discriminate the birds into treatment group  and, 3) as the prelude step identifying the targets for the rif analysis.
pif is computed by multiplying de by abundance i.e.
the two sources of variation in expression data.
the top enrichment among the  <dig> upregulated genes  in hfe by pif is “cholesterol biosynthesis”  based on the presence of dhcr <dig>  insig <dig>  fdps, sc5d and mvd.
the top enrichment among the  <dig> downregulated genes in hfe by pif is “extracellular matrix organisation”  based on the presence of fmod, itga <dig>  actn <dig>  smoc <dig>  fn <dig>  itga <dig>  mfap <dig>  cma <dig>  tnc, gsn, tll <dig>  nrxn <dig>  col1a <dig>  col3a <dig>  ccdc <dig>  dcn, ltbp <dig>  cyp1b <dig>  itgb <dig>  ncam <dig>  mmp <dig>  lox, ctsk, col6a <dig>  tgfb <dig>  col6a <dig>  col22a <dig>  cd <dig>  ctss, col5a <dig>  lcp <dig>  col6a <dig>  loxl <dig>  comp, anxa <dig>  thbs <dig>  abi3bp, col12a <dig>  ecm <dig>  col11a <dig>  vcan, postn and nid <dig>
the top  <dig> upregulated and downregulated genes by pif are tabulated with brief functional descriptions and absolute fold changes in de .table  <dig> the  <dig> most extreme genes by phenotypic impact factor in high fe versus low fe birds.
these  <dig> genes have been reordered manually into functional groupings.
the absolute de expressed as a fold change has also been reported
mb
myh15
tpm3
myoz2
tnni1
pln
col3a1
col12a1
fn1
bmp5
rsfr
ly86
mgat5b
nes
apoa1
mgp
igj
sepp1
ca3
ctss
h2afy
hig2a
tesk1
c11orf89
*
ckmt1
cntfr
enkd1
ca4
insig1
ttpa
lingo1
dcun1d5
sncg1
c16orf45
lipg
bfsp1
myl2_ <dig>
mybph
* we further explored the possible function of the gene encoding an unannotated open reading frame  c11orf <dig>  given it was in the top  <dig> upregulated genes in the hfe birds.
the entry in genecards  predicts the orf to translate into a large protein containing regions of low complexity and comparative analysis suggests that the protein was present in the ancestor of the chordates.
it is located in the g. gallus genome adjacent to tnnt <dig> and tnni <dig>  and cellular compartment analysis  suggests there is moderate evidence it is a cytoskeletal protein
we also cross-referenced the genes and direction of change in table  <dig> to the previous microarray study on the same samples  <cit> .
of the  <dig> genes for which we could find exact matches,  <dig> showed the same direction of change by the two different mrna quantitation technologies.
mrna encoding the mitoproteome
of the  <dig> genes encoding proteins localised to the mitochondrion,  <dig> were upregulated in hfe and  <dig> were downregulated in hfe.
this is a very significant skew  implying the hfe birds have an elevated mitochondrial content.
a similar skew was also observed in upregulated proteins in hfe breast muscle in a recent proteomics study  <cit> .
this indicates that in general terms the upregulated genes tended to translate into upregulated mitochondrial proteins.
muscle model visualisation
overlaying de genes between hfe and lfe birds  on the vmus3d, identified a coordinated down-regulation of thick , thin  and z-disc  proteins encoded by slow twitch isoform genes in hfe birds .
alterations of genes encoding myosin heavy and light chain proteins  together with the tropomyosin : troponin  are components of the key protein complex in determining biomechanical efficiency of striated muscles  <cit>  .
furthermore, the z-disk and costamere have been documented as key complexes involved in sensing mechanical stress and converting this to chemical signals  <cit> .
although no costameric proteins were identified as de using vmus3d, the z-disk protein csrp <dig> has been associated with both structural and signal transduction roles in cardiac muscle  <cit>  and myoz <dig> has been linked to calcineurin signalling, a key protein phosphatase involved in fibre type regulation  <cit> .
together these results indicate a difference in contractile protein composition and key contractile proteins involved in the regulation of fibre type.
it is worth noting that in the previous microarray study that used pooled rna samples from the hfe and lfe muscle tissue, csrp <dig>  myoz <dig>  and tnni <dig> were the number  <dig>   <dig>  and  <dig> upregulated genes in the lfe phenotype.
thus, there appears to be good agreement between the microarray and rnaseq analysis data.
permut matrix clustering on extreme pif
after normalising on rows, clustering on columns correctly ascribed each bird into the correct group of origin  .
interestingly, the dendrogram suggests the hfe birds have lower dispersion than the lfe birds.
functional groups contained in the  <dig> downregulated genes in hfe include slow twitch contractile subunits  and associated slow twitch metabolic machinery ; components of the extracellular matrix ; negative regulation of igf <dig> signalling ; and cholesterol metabolism .
the observed downregulation of bmp <dig>  a negative regulator of igf <dig> signalling, implies igf <dig> signalling is increased in hfe birds.
functional groups contained in the  <dig> upregulated genes include mitochondrial metabolism ; fast twitch muscle ; and acid-base balance and response to ph .
regulatory impact factor analysis
the rif <dig> and rif <dig> scores were plotted and manually explored for outlier transcription factors  .
most of the data is centred close to  <dig>  implying that the majority of tfs are not differentially connected and therefore perform the same function in the same manner in the two groups.
however, there are a small number of tf that are highly differentially connected.
the top  <dig> encode progesterone receptor  and two tfs that are components of progesterone signalling pathways .
quantitative pcr technical validation
the comparison between the rnaseq and qpcr data is outlined in table  <dig>  the direction of expression change was confirmed in all nine cases.
in seven cases the qpcr detected a larger de than the rnaseq.table  <dig> the correspondence between rnaseq and qpcr data for mrna encoding slow twitch subunits, proteins associated with slow fibres and intramuscular fat content
myh15
tpm3
myoz2
tnni1
mybpc1
pln
fabp4
mb
ca3
progesterone receptor immunofluorescence
in a preliminary attempt to explore the localisation of the progesterone receptor in birds, we undertook an immunofluorescence approach in a quail muscle cell line.
this qualitative approach showed a near complete overlap between the presence of mitochondria and the progesterone receptor.
this indicates that most, if not all, progesterone receptor proteins are localised to the mitochondria in this particular in vitro avian muscle cell line .fig.
<dig> the strong co-localisation of the progesterone receptor in the mitochondria of a quail cell line
discussion
muscle structure
at a gross histological level, broiler pectoralis muscle is unique among the vertebrates, being almost exclusively built of fast, glycolytic type iib fibres  <cit> .
this muscle tissue would be expected to have a sparse capillary network, and the individual fibres to possess a low mitochondrial content and few lipid droplets.
this expectation was essentially borne out in the 1970’s by  <cit> , where a mitochondrial content of just 4% was determined by transmission electron microscopy and morphometric analysis in chicken pectoralis muscle.
this is approximately 10-fold less than the 35% mitochondrial content of very athletic hummingbird pectoralis muscle  <cit> .
collectively, these observations indicate that chickens have an unusually homogeneous breast muscle structure bearing a very low aerobic capacity.
broiler chickens are the most feed efficient vertebrates known with a dry weight in to wet weight out ratio of ~ <dig>   <cit> .
it has been previously argued that the very low aerobic capacity of their muscle is a likely driver of this phenotype through avoiding paying for ‘expensive’ spare aerobic capacity and considerations of economic physiological design  <cit> .
indeed, it has been observed in many independent circumstances that the most feed efficient breeds, selected on either fe or muscle mass, transition from a redder, slow, aerobic ancestral muscle to a faster, whiter, less aerobic glycolytic one in the derived state.
examples include mstn mutant piedmontese  <cit> , belgium blue  <cit>  and blond d’aquitaine  <cit>  cattle, callipyge sheep  and both large white  <cit>  and german landrace pigs  <cit> .
the transition from high myoglobin content red to low myoglobin content white fibres in pigs has been so dramatic that they are no longer considered a rich source of dietary iron and are viewed as the second ‘white’ meat after chicken, rather than the third red meat after beef and lamb.
we know that these broiler muscle transcriptomes are derived  almost exclusively from white, fast type iib fibres  <cit> .
in line with anatomical and physiological expectation the two most highly expressed genes in the muscle are the fast myosin heavy chain isoform myh1e and the regulator of anaerobic glycolysis gapdh.
however, we also detected gene expression of a number of slow contractile subunits  in addition to the expected fast ones.
this somewhat surprising finding implies that a subset of remnant slow subunits make up a small but measurable proportion of the sarcomeric structure of the type iib fibres, perhaps such that they are a mosaic at the molecular level.
the subtle but coordinated reduction of these slow twitch isoforms and associated slow machinery  in the higher fe birds is in line with the broad expectation of muscle ‘whitening’ being associated with high fe described above.
it also concurs with our previous observations based on a cdna microarray  <cit> .
the divergent pattern we observe in ca <dig> and ca <dig> may also be a simple consequence of fast versus slow fibre subunit abundance  <cit> .
intriguingly, while the majority of slow subunits were downregulated in hfe including myl2b, a separate myl <dig> isoform  was upregulated in hfe birds.
we do not have an explanation for this apparent inconsistency.
one would expect the expression of mitochondrial proteins to track the fibre composition data.
however, our observations on the gene expression of mitochondrial proteins are not congruent with the reduction in the slow muscle fibre subunit and myoglobin/phospholamban data in hfe birds.
mitochondrial metabolism
in terms of mitochondrial content we found that the mrna encoding the mitoproteome is modestly but highly significantly skewed towards the hfe birds.
we propose that this indicates a slightly higher mitochondrial content and therefore a somewhat higher aerobic capacity.
this is prima facie contrary to the reduction in slow fibre subunits and reduced myoglobin/phospholamban outlined above.
it is also at odds with theoretical expectation  <cit>  and mrna/protein abundance data from numerous other species, particularly recent genome-wide screens performed on pigs divergent for rfi  <cit> .
it also appears contradictory to previous work showing in chickens divergent for mass, the more efficient larger birds have a lower, not higher, mass-specific metabolic rate  <cit> .
having said this, we believe these specific samples have been correctly characterised.
these mrna data are in agreement with a mass spectrometry-based proteomic analysis carried out on a subset of the same breast muscle tissues , as the protein expression of the mitochondrial proteome was also somewhat elevated in the hfe phenotype  <cit> .
therefore, while the ‘big picture’ finding in production animals appears to be that selection on mass and fe reduces mitochondrial content and aerobic capacity, the data from this particular fe contrast between highly selected birds tends to point in the opposite direction.
the underlying reason for this notable disparity are the source of continued investigations.
one possibility is that the commercial birds have such a low mitochondrial content that they are on the verge of expressing metabolically driven muscle pathologies.
this would be consistent with the recent emergence of concerning phenotypes such as ‘white striping’ and ‘wooden breast.’ if a low mitochondrial content drives these pathologies, one might speculate that subtly increasing mitochondrial content may increase fe through impact on general health.
ultimately, the mitochondrial content finding has been approximated using gene expression.
ideally, this finding would be validated with an independent technology such as qpcr of mtdna expressed relative to ndna, but unfortunately we no longer have adequate sample material.
the general upregulation of mrna encoding the mitoproteome in hfe birds is modest in absolute terms and there are many genes that are downregulated as well as upregulated.
extreme among the upregulated molecules in hfe are tdrkh, mrpl <dig>  ldhd, ckmt1a encoding a gene silencer, a mitoribosome subunit, lactate dehydrogenase d and mitochondrial creatine kinase 1a.
the biological implications of tdrkh for skeletal muscle are unknown, but mrpl <dig>  ldhd and ckmt1a indicate an increase in mitochondrial protein synthesis , a change in redox status and increasing the transfer of high energy phosphate from the mitochondria to the cytosolic carrier, creatine.
the downregulation of pdk <dig>  acacb and ldhb  indicate the hfe birds have a preference for burning carbohydrate instead of fat, a reduction in fatty acid uptake and combustion by mitochondria and a reduced reliance on anaerobic glycolysis.
the first two point to a change in fuel combustion pattern that is consistent with a subtle fibre subunit transition towards a whiter muscle phenotype–after all, glycolytic fibres are geared to burn glycogen not fat during exercise.
the reduction in the anaerobic ldhb in hfe is perhaps consistent with an elevated mitochondrial content, as the presence of higher aerobic metabolism may reduce the need to resort to anaerobic metabolism for atp production.
intramuscular fat and extracellular matrix content
although arguably not statistically significant sensu stricto, we wish to highlight the expression profiles of several genes previously implicated in intramuscular fat  content given the important role imf plays in both meat quality and energetics.
in the hfe birds there is a trend towards a small reduction in expression of two genes  previously associated with amount and activity of intramuscular fat in cattle , sheep  <cit>  and pigs  <cit> .
in absolute terms plin <dig> and fabp <dig> are  <dig>  fold and  <dig>  fold downregulated in hfe respectively.
plin <dig> is ranked in position  <dig>  when the  <dig>   <dig> genes are ranked by pif from negative  to positive  values.
together, the expression of these two genes implies the hfe birds possess a slightly leaner breast muscle.
in cattle muscle, these two genes tend to be tightly co-expressed with thrsp, s100g, cidea and cidec  <cit> , but we did not detect expression for these additional genes in the broiler data.
the tendency towards higher muscle tissue leanness in the high fe birds is unsurprising, relating to the selection for overall animal mass in the calculation of residual feed intake.
using basic animal mass as the selection phenotype will inevitably select against light, anhydrous but energy-dense adipose tissue and promote heavy, watery lean muscle tissue.
further, very lean meat is texturally dry and relatively flavourless–a universal problem in production animals or breeds with relatively high fe such as callipyge sheep  <cit> , myostatin mutant sheep  <cit>  and large white pigs  <cit> .
additionally, a very strongly coordinated reduction in expression of  <dig> genes encoding extracellular matrix proteins in hfe birds is consistent with less ‘white striping’, a meat defect that has recently appeared in poultry, particularly large birds grown too quickly  <cit> .
white striping comprises seams of collagenous material containing some embedded adipocytes.
the coordinate down regulation in hfe birds of genes encoding ecm components in addition to those relating to lipid content can also be seen as consistent with a reduction in this pathology.
progesterone signalling and fe
we have characterised the anatomical and physiological differences between high and low fe broiler muscle.
the hfe bird muscle exhibits a reduction in slow muscle subunits and associated machinery and a subtly elevated mitochondrial content.
but what is the ultimate cause of this muscle phenotype shift?
the differential network rif analysis clearly points to alterations in progesterone signalling as being the major driver.
the  <dig> most differentially connected transcription factors are cited <dig>  pgr and gzf <dig>  that is, the progesterone receptor itself, and two molecules that communicate progesterone signalling or are influenced by it  <cit> .
however, the rif analysis does not imply activation or repression of the pathway, the direction of change needs to be inferred by different means.
in a concurrent proteomics study being conducted on the pedm broiler line, progesterone was predicted to be activated using a commercial software program , based on expression of downstream proteins  <cit> .
thus, two independent analytical approaches, one based on global gene expression and the other based on global protein expression both point towards progesterone as a driver of the fe phenotype.
progesterone is synthesised from cholesterol, so our additional identification of “cholesterol biosynthesis”  as enriched in the upregulated genes in hfe birds is noteworthy.
in mammals, progesterone synthesis is known to occur in reproductive organs, adrenal glands, nervous tissue and adipose tissue.
to the best of our knowledge endogenous synthesis has not been detected in mammalian skeletal muscle.
the situation in birds is unknown.
in fact, we could not find any published information on progesterone production nor mode of action in avian muscle.
however, in cattle exogenous administration of progesterone is often used in hormone growth promotant mixes  that raise fe by ~20%  <cit> .
indeed, the synthetic gestagen , melengestrol acetate, has been found to stimulate muscle growth in yearling heifers  <cit> .
furthermore, in pregnant humans progesterone alters muscle phenotype from a slow oxidative to fast glycolytic state  <cit> .
this slow to fast contractile skeletal isoform transition is certainly consistent with what we observe in our hfe broiler data .
progesterone also promotes protein synthesis in human skeletal muscle and rat cardiac muscle  <cit> .
finally, a truncated isoform of the progesterone receptor localises to the mitochondrion in mammals  <cit> .
here, we present a first line of evidence that the progesterone receptor also localises to the mitochondria in birds, based on immunofluorescence in a quail muscle cell line.
consequently, the posited role of progesterone signalling via the progesterone receptor in driving variation in fe is a further line of evidence consistent with the mitochondrion being a central player.
previous work on muscle gene expression and poultry feed efficiency
we have previously documented gene expression patterns in these same samples using microarray technology  <cit> , finding a similar reduction in the hfe birds of genes encoding a set of slow twitch contractile subunits  plus slow machinery pln and atp2a <dig>  another group  <cit>  examined  <dig> breast muscle samples from a different population of commercial broilers also using rna sequencing technology.
surprisingly, some of their basic structural findings are exactly opposite to ours, in that they found the slow contractile subunits  and intramuscular fat content  as upregulated in the hfe birds.
the reason for this very fundamental contradiction is a mystery at this stage, but may be attributable to different genetics or some other factor.
finally, we wish to emphasise that the progesterone signalling prediction made by rif on our mrna data has also been independently borne out by ingenuity pathway analysis ‘upstream analysis’ on proteomics data derived from the same muscle samples  <cit> .
thus, analyses performed at two levels of biological organisation  and two independent inference strategies  have come to the same conclusion.
CONCLUSIONS
selection on fe at the whole bird level culminates in measurable changes in breast muscle phenotype.
despite the subtlety of the changes, we have identified a panel of  <dig> genes whose expression level can discriminate the two groups and whose functions include slow muscle structure  and mitochondrial activity.
hierarchical clustering shows the hfe birds show less dispersal than the lfe, reflecting a reduction in gene expression variation.
this concept of reduced variation  has been observed before as a real signal of pathway constraint  <cit> , and therefore can be taken as informative in the same way as measures of differential expression and differential connectivity.
the gene expression changes in the hfe birds suggest a whiter, leaner tissue which presumably reflects the well established relationship between elevating fe on the one hand but increasing the likelihood of deficits in meat quality on the other.
this exact same issue has been observed in many other highly selected production animals such as large white pigs and callipyge sheep.
the differential network analysis clearly points to modifications in progesterone signalling as a key driver of the difference between the two groups.
the next challenge is to harness this finding, perhaps through diet or genetics, to help drive broiler fe even higher without exogenous administration of progesterone.
any attempt to increase fe will need to be cogniscent of possible negative effects on associated meat quality traits or the health of the birds.
additional files
additional file 1: pca analysis of  the raw rpm read count data,  the log <dig> transformed data and  the quantile normalised log <dig> data.
green denote hfe, red denote lfe and blue are barred rock outgroup animals.
the barred rock birds were used for collective normalisation but not formally analysed in this study.
the outlier red sample is lfe bird  <dig>
additional file 2: the qpcr primers used to validate the de of a subset of genes prioritised by rnaseq.
additional file 3: the very strong relationship between rpm and rpkm in our data, using pif as the metric for comparison.
the pif based on rpm has a correlation of  <dig>  to the pif based on rpkm.
additional file 4: the rif <dig> and rif <dig> output listed in full for the  <dig> transcription factors.
abbreviations
fefeed efficiency
hfehigh feed efficiency
lfelow feed efficiency
pedmpedigree broiler males
drs.
ron okimoto , tony vuocolo  and gene wijffels  provided insightful reviews of this manuscript prior to submission.
funding
this research was based upon work that is supported, in part, by the national institute of food and agriculture, u.s. department of agriculture, under award number 2013- <dig>  the commonwealth science and industrial research organisation mcmaster fund supported wb on a visiting fellowship to australia, during which time njh and wb conceived this project.
the funders had no say in the design of the experiment or its biological interpretation.
availability of data and materials
the data sets analysed during the current study are available from the corresponding author upon reasonable request.
authors’ contributions
wb and njh designed the study.
njh drafted the manuscript.
bwk performed the wet lab work.
njh and bwk mapped reads, counted reads and normalised the data.
aw generated the 3d muscle model and interpreted the output.
ar performed the rif analysis.
njh, bwk and wb biologically interpreted the data.
kl performed the quail muscle immunofluorescence study.
wb managed the project resources and provided oversight.
all authors read, contributed suggestions and approved the submitted version.
competing interests
the authors declare that they have no competing interests.
consent for publication
not applicable.
ethics approval
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 : protocol # <dig>  the animal handling and euthanasia for this project were approved by the institutional animal care and use committee .

 in a resource-constrained world supporting a rapidly growing human population there is great interest in enhancing the production efficiency of our major animal and plant food industries. feed is the single largest commercial cost in animal production. any inefficiency not only affects the bottom line but also negatively impacts resource usage and waste production. consequently, there is considerable economic and environmental interest in increasing animal feed conversion efficiency. small increases in feed efficiency across numerous animals can have a large industrial and environmental impact. broiler chickens are the most feed efficient of all the vertebrates, with some companies boasting a conversion of just ~ units dry weight intake per unit wet weight gain. this makes broilers particularly valuable as biological models for understanding the mechanistic basis of feed efficiency. along with other avian meat producing species, cf. turkeys, pheasants, partridge, grouse and quails, chickens are members of the phasianidae taxonomic clade, the largest branch of the galliformes. the members of this group tend to be sedentary, resident ground-dwelling birds that use short, burst flights to escape predators. this behaviour explains the functionally unusual breast muscle metabolism of the ancestral birds such as the red junglefowl progenitor of modern domestic chickens, which is dominated by explosive, fast twitch contractile isoforms and a relatively low oxidative capacity metabolism. the birds we chose for this analysis are pedigree broiler males from a population of highly inbred, commercial animals subject to very strong historical selection for high muscle mass and elevated feed efficiency. the significant difference in fe between the two treatment groups is in the order of fold in absolute terms. some of the variation in the fe phenotype presumably derives from biological modifications in the non-muscle tissues. this logic, coupled with the prior application of microarray technology on the same individuals, gave us the expectation that any molecular differences we may detect in the breast muscle between the two groups would be fairly subtle. nevertheless, it has been established for some time that commercial broilers possess a muscle structure that is, histologically speaking, composed almost exclusively of low mitochondrial content type iib glycolytic fibres. this fibre homogeneity is advantageous as it provides a very stable transcriptional background against which to perform the various molecular analyses. here, we have used illumina rna sequencing technology to screen the breast muscle transcriptomes of the high fe and low fe groups. we aimed to both test existing hypotheses as well as generate new ones. we have previously argued for a role of both isolated mitochondrial physiology and overall tissue mitochondrial content in driving variation in production animal feed efficiency. therefore, as part of our analyses we have specifically harvested the data for genes encoding mitochondrial proteins, finding expression was skewed towards the hfe birds. furthermore, we have also used a combination of basic differential expression and a recently developed differential connectivity algorithm called regulatory impact factor analysis to screen the data in an unbiased fashion. the rif algorithm detects regulatory molecules that are transcriptionally rewired between the two states, irrespective of whether that molecule is de itself. this is a potent approach for the identification of transcription factors that tend to be lowly and stably expressed at the mrna level, but whose functional activity is controlled at the protein level through changes in cellular localisation, ligand binding or post-translational modifications. this analysis points to progesterone signalling through the progesterone receptor as a key molecular driver. methods animal resources and phenotypes we used breast skeletal muscle samples from male broilers divergent in feed efficiency as previously described. male broilers were raised under standard industry practices and euthanized using carbon dioxide asphyxiation. cobb-vantress inc. adheres to internal industrial standard guidelines for animal welfare and has standard operating procedures for euthanasia, dissection and other experimental techniques. in brief, the broilers were derived from a single genetic line of a commercial cobb-vantress inc. population. both the hfe and lfe birds were sub-sampled from the extremes of an original group of birds. in turn, those birds were selected as the most efficient from a larger population of birds. the high fe group had greater body weight gain from to weeks but consumed the same amount of feed as the lfe group. for the purposes of this manuscript we have defined fe as the gain in bird mass relative to feed intake, such that a higher number equates to a more efficient bird. in our case, the hfe birds yielded efficiencies of  ±  compared to  ±  for the lfe birds. this animal resource is noteworthy in the sense that the two groups are not divergently selected. this is an important distinction compared to other research populations in this area. for example, the observed de in divergent selection experiments could be a consequence of functionally irrelevant founder effects that have accumulated over subsequent matings. in our case the birds are from the same line and generation, but separated on the basis of 1 week of fcr testing in a controlled experimental environment. mapping, counting and normalising the mrna reads after phenol chloroform rna extraction we submitted our samples to the research support facility at michigan state university for illumina hiseq base pair paired end read sequencing. in brief, we used clc genomics workbench to map the reads to gallus gallus genome assembly version the clc software follows the analytical pipeline recommended in. we log transformed the rpm data to stabilise the variance and then performed a further quantile normalisation. the impact of the two levels of normalisation on whether the individual samples could discriminate into the treatment groups of origin was confirmed by principal component analysis. differential expression to compare the treatment groups we plotted the ma plot comparing hfe birds to lfe birds. the x axis, a, represents the average abundance in the groups, and the y axis, m, is hfe average transcript abundance minus lfe average transcript abundance. a subset of functionally relevant outlier de genes was chosen for annotation on the plot.fig. a the slow contractile subunits highlighted are myh tpm myoz tnni myl tnnt and csrp consistent with a shift in fibre composition towards the fast glycolytic and away from the slow aerobic fibres in hfe birds. b the mrna encoding the mitoproteome is significantly skewed towards the hfe consistent with a slight increase in mitochondrial content. c the extreme 5% pif used as input for the rif analysis have been highlighted. these de genes are distributed evenly across the ma plot, closely tracking the overall shape of the plot mrna encoding the mitoproteome we downloaded the complete human mitoproteome from http://mitominer.mrc-mbu.cam.ac.uk/release- /begin.do generating a list of nuclear and mitochondrial genes encoding mitochondrial proteins. we converted the protein names to gene names, and found matches in our chicken rnaseq data. the mitoproteome was then plotted in ma format, and all those with higher values in the hfe group were colour coded blue and those with lower values in the hfe colour coded red. the skew in distribution away from the null expectation of 50: was quantified by binomial statistics. this approach formalises the extent to which the mitoproteome data is biased to one or the other of the two groups. phenotypic impact factor we next computed a modified de called pif which we have found to have a number of appealing characteristics, both numerical and biological. from a numerical perspective it can be used to establish extreme de genes in a way that accounts for the structure of the distribution of the data. it de-emphasises lowly abundant genes that are inherently noisy as they approach the detection limit of the technology. it also draws together the major sources of variation in gene expression data into a single metric which can be used for comparison purposes. from a biological perspective pif analysis is useful as it gives more meaningful functional enrichments than ranking on de, including better characterising muscle fibre composition shifts in muscle transcriptome data. also, pif is a foundation of the rif differential network analysis where the target de molecules are first defined by extreme pif, and then the regulators differential co-expression to those targets is integrated into a differential connectivity metric. we ranked on pif and assessed functional enrichment at the extremes of the list using gorilla. vmus3d muscle model to visualise spatial locations of proteins encoded by de genes derived from the hfe versus lfe birds we utilised the virtual muscle 3d. vmus3d, described in detail here, is an annotated 3-dimensional representation of key contractile protein complexes, including thick, thin, z-disc and costameric complexes, and their relative protein spatial locations. vmus3d is a database driven visualisation tool that enables mapping of gene expression changes onto their encoded proteins, represented as an extensible 3d object, via the document object model. here, we mapped m-values onto the vmus3d and coloured change as follows: blue, green, yellow, orange and red.fig. 3d muscle model visualisation illustrating the downregulation of slow subunits in hfe birds consistent with a whiter, sprint muscle phenotype. a. longitudinal view of the contractile apparatus with coordinated de changes of slow unit proteins highlighted. b. z-disk proteins differentially expressed. c. protein locations across the contractile apparatus and their corresponding values. colours correspond to the level of de as follows: blue, green, yellow, orange and red permut matrix clustering on extreme pif the top up and top downregulated genes were identified as described above. a matrix was formed comprising as many rows as genes and as many columns as birds with each cell containing the normalised gene expression. this matrix was imported into permutmatrix software for hierarchical clustering.fig. a hierarchically clustered heatmap showing the expression patterns of the most extreme pif mrna. this panel of genes which includes mrna encoding muscle structural proteins and mitochondrial proteins can correctly discriminate the two treatment groups regulatory impact factor analysis the purpose of rif is to identify regulatory molecules whose activity can change independent of any change in gene expression level. this includes transcription factors whose activity may be influenced by ligand binding, co-factor binding, cellular localisation and other post-translational processes. the method works by establishing different patterns of network connectivity in the two states. rif and rif are two versions of essentially the same analysis. in both cases the abundance and differential expression of the ‘target genes’ is exploited in conjunction with the differential co-expression of the ‘regulators’ to those ‘targets.’ the output of both versions has been plotted and reported here. we computed rif and rif as previously described. this procedure exploits global patterns of differential co-expression to infer those regulatory molecules whose behaviour is systematically different in a contrast of interest, in this case the hfe versus lfe birds. herein, the experimental contrast was hfe vs. lfe and the rif metrics for the r-th regulator were computed using the following formulae: rif1r=1nde∑j=1j=ndexj×dj×dcrj and rif2r=1nde∑j=1j=ndexjhfe×rrjhfe2−xjlfe×rrjlfe where nde represented the number of de genes; x j was the average expression of the j-th de gene across all time points; d j was the de of the j-th gene in the hfe vs. lfe contrast; dcrj was the differential co-expression between the r-th regulator and the j-th de gene, and computed from the difference between r rjhfe and r rjlfe, the correlation co-expression between the r-th regulator and the j-th de gene in the hfe and lfe samples, respectively; finally, x jhfe and x jlfe represented the average expression of the j-th de or ts gene in the hfe and lfe samples, respectively. the rif and rif output were simultaneously plotted.fig. the extremely differentially connected tf as illustrated by rif and rif scores. progesterone receptor and other tf involved in progesterone signalling are awarded extreme scores by both metrics quantitative pcr technical validation to independently validate our rnaseq normalisation method and method to detect de, we selected a subset of nine representative genes for qpcr amplification. these were mrna encoding slow twitch subunits, slow muscle associated proteins and muscle fat content. we used pcr cycling conditions as previously described. primers are detailed in additional file in brief, one microgram of total rna was obtained from muscle samples each for high and low fe pedigree broilers for general validation of the rnaseq results and for specific confirmation. rna was converted into cdna with qscripttm cdna supermix following the manufacturer’s instructions. the cdna samples were diluted in a 1: ratio and a portion of the cdna was subjected to a real-time quantitative pcr reaction using an abi prism 7500ht system with poweruptm sybr® green master mix. the specific oligonucleotide primers were designed by the primer program. 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 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 times and the values of average cycle threshold were determined for each sample, and 2 − ΔΔct values for the comparison of hfe and lfe muscles were used for relative quantification by fold-change and statistical significance. quail muscle cell culture and treatments qm cells were grown in m medium with 10% fetal bovine serum, 10% tryptose phosphate, and 1% penicillin-streptomycin at 37 °c under a humidified atmosphere of 5% co and 95% air. at 80–90% confluence, cells were synchronized overnight in serum-free medium. progesterone receptor immunofluorescence immunofluorescence was performed as previously described with modifications. briefly, cells were grown to 50–60% confluence in chamber slides followed by methanol fixation for 10 min at 20 °c. serum free blocking buffer, and the cells incubated with rabbit anti-progesterone receptor antibody overnight at 4 °c. cells were treated with mitotracker deep red and visualized with alexa fluor 488-conjugated secondary antibody. after dapi counterstaining, a cover slip was placed on slides with vectashield. images were obtained using the zeiss imager m with a 20x plan-apochromat 20x/ objective and a 100x ec planneofluor 100x/ oil objective. the alexa fluor fluorophore was observed through filter set 1031- with an excitation of bp 470/ beamsplitter of ft and emission spectrum of bp 525/ differential interference contrast images were collected using dic m condensers. the alexa fluor fluorophore was excited for 500 ms prior to capturing each image using an axio cam mr camera. all analysis was performed using axiovision se. sp software. statistical analysis given the two treatment groups are so similar and that multiple testing penalties can be very severe in this kind of experimental design we elected not to identify significantly de genes on a gene-by-gene basis. rather, we used a range of systems-wide analyses to interpret our data. this approach has the advantage of not being reliant on the performance of any given gene, but rather draws on information present in entire functional groups of genes. for basic functional enrichment analysis we used gorilla. here, we examined functional enrichment in the extreme up and downregulated genes genes ranked by pif, compared to a background list consisting of all genome-wide data; both p and fdr q-values reported. this nominal 5% de genes also represent the target de genes for the rif analysis. for the mitoproteome analysis we used binomial statistics to quantify the skew for all mrna encoding mitochondrial proteins. for the rif analysis we report on extreme outlier regulators only.  twelve rnaseq libraries were constructed using breast muscle rna samples from hfe and lfe broilers. in total, million 100 bp sequences were obtained with an average of million reads per sample and 80% of the reads were mapped to the g. gallus reference genome assembly. after mapping, counting and normalising we returned data for genes. we screened for those mrna with no missing values across the samples. this gave a total genes which were used for all subsequent analysis. both rpm and rpkm values were computed and compared by correlation analysis. to achieve this we calculated the pif for each gene and correlated the pif scores by the two methods. given the high correlation and the concern that the additional normalisation used in generating rpkm values can theoretically introduce a variety of problems, we elected to proceed with the normalised rpm values. principal component analysis the raw read count data could not discriminate the birds into treatment group of origin. log stabilising the data allowed separation of individuals into the two groups with one exception. this separation by treatment group was reinforced by a subsequent quantile normalisation which standardises the distributions of the reads for each sample. this separation provided a quality check that our normalisation approach had maintained the treatment differences without introducing any systematic bias. further validation that the normalisation approach did not introduce bias was provided by determination that the ma plot was centred on differential expression between hfe and lfe birds we used a modified de, named pif for: 1) functional enrichment analysis in gorilla, 2) to generate a panel of genes that can discriminate the birds into treatment group and, 3) as the prelude step identifying the targets for the rif analysis. pif is computed by multiplying de by abundance i.e. the two sources of variation in expression data. the top enrichment among the upregulated genes in hfe by pif is “cholesterol biosynthesis” based on the presence of dhcr insig fdps, sc5d and mvd. the top enrichment among the downregulated genes in hfe by pif is “extracellular matrix organisation” based on the presence of fmod, itga actn smoc fn itga mfap cma tnc, gsn, tll nrxn col1a col3a ccdc dcn, ltbp cyp1b itgb ncam mmp lox, ctsk, col6a tgfb col6a col22a cd ctss, col5a lcp col6a loxl comp, anxa thbs abi3bp, col12a ecm col11a vcan, postn and nid the top upregulated and downregulated genes by pif are tabulated with brief functional descriptions and absolute fold changes in de.table the most extreme genes by phenotypic impact factor in high fe versus low fe birds. these genes have been reordered manually into functional groupings. the absolute de expressed as a fold change has also been reported mb myh15 tpm3 myoz2 tnni1 pln col3a1 col12a1 fn1 bmp5 rsfr ly86 mgat5b nes apoa1 mgp igj sepp1 ca3 ctss h2afy hig2a tesk1 c11orf89 * ckmt1 cntfr enkd1 ca4 insig1 ttpa lingo1 dcun1d5 sncg1 c16orf45 lipg bfsp1 myl2_ mybph * we further explored the possible function of the gene encoding an unannotated open reading frame c11orf given it was in the top upregulated genes in the hfe birds. the entry in genecards predicts the orf to translate into a large protein containing regions of low complexity and comparative analysis suggests that the protein was present in the ancestor of the chordates. it is located in the g. gallus genome adjacent to tnnt and tnni and cellular compartment analysis suggests there is moderate evidence it is a cytoskeletal protein we also cross-referenced the genes and direction of change in table  to the previous microarray study on the same samples. of the genes for which we could find exact matches, showed the same direction of change by the two different mrna quantitation technologies. mrna encoding the mitoproteome of the genes encoding proteins localised to the mitochondrion, were upregulated in hfe and were downregulated in hfe. this is a very significant skew implying the hfe birds have an elevated mitochondrial content. a similar skew was also observed in upregulated proteins in hfe breast muscle in a recent proteomics study. this indicates that in general terms the upregulated genes tended to translate into upregulated mitochondrial proteins. muscle model visualisation overlaying de genes between hfe and lfe birds on the vmus3d, identified a coordinated down-regulation of thick, thin and z-disc proteins encoded by slow twitch isoform genes in hfe birds. alterations of genes encoding myosin heavy and light chain proteins together with the tropomyosin: troponin are components of the key protein complex in determining biomechanical efficiency of striated muscles. furthermore, the z-disk and costamere have been documented as key complexes involved in sensing mechanical stress and converting this to chemical signals. although no costameric proteins were identified as de using vmus3d, the z-disk protein csrp has been associated with both structural and signal transduction roles in cardiac muscle and myoz has been linked to calcineurin signalling, a key protein phosphatase involved in fibre type regulation. together these results indicate a difference in contractile protein composition and key contractile proteins involved in the regulation of fibre type. it is worth noting that in the previous microarray study that used pooled rna samples from the hfe and lfe muscle tissue, csrp myoz and tnni were the number and upregulated genes in the lfe phenotype. thus, there appears to be good agreement between the microarray and rnaseq analysis data. permut matrix clustering on extreme pif after normalising on rows, clustering on columns correctly ascribed each bird into the correct group of origin. interestingly, the dendrogram suggests the hfe birds have lower dispersion than the lfe birds. functional groups contained in the downregulated genes in hfe include slow twitch contractile subunits and associated slow twitch metabolic machinery; components of the extracellular matrix; negative regulation of igf signalling; and cholesterol metabolism. the observed downregulation of bmp a negative regulator of igf signalling, implies igf signalling is increased in hfe birds. functional groups contained in the upregulated genes include mitochondrial metabolism; fast twitch muscle; and acid-base balance and response to ph. regulatory impact factor analysis the rif and rif scores were plotted and manually explored for outlier transcription factors. most of the data is centred close to implying that the majority of tfs are not differentially connected and therefore perform the same function in the same manner in the two groups. however, there are a small number of tf that are highly differentially connected. the top encode progesterone receptor and two tfs that are components of progesterone signalling pathways. quantitative pcr technical validation the comparison between the rnaseq and qpcr data is outlined in table  the direction of expression change was confirmed in all nine cases. in seven cases the qpcr detected a larger de than the rnaseq.table the correspondence between rnaseq and qpcr data for mrna encoding slow twitch subunits, proteins associated with slow fibres and intramuscular fat content myh15 tpm3 myoz2 tnni1 mybpc1 pln fabp4 mb ca3 progesterone receptor immunofluorescence in a preliminary attempt to explore the localisation of the progesterone receptor in birds, we undertook an immunofluorescence approach in a quail muscle cell line. this qualitative approach showed a near complete overlap between the presence of mitochondria and the progesterone receptor. this indicates that most, if not all, progesterone receptor proteins are localised to the mitochondria in this particular in vitro avian muscle cell line.fig. the strong co-localisation of the progesterone receptor in the mitochondria of a quail cell line discussion muscle structure at a gross histological level, broiler pectoralis muscle is unique among the vertebrates, being almost exclusively built of fast, glycolytic type iib fibres. this muscle tissue would be expected to have a sparse capillary network, and the individual fibres to possess a low mitochondrial content and few lipid droplets. this expectation was essentially borne out in the 1970’s by, where a mitochondrial content of just 4% was determined by transmission electron microscopy and morphometric analysis in chicken pectoralis muscle. this is approximately 10-fold less than the 35% mitochondrial content of very athletic hummingbird pectoralis muscle. collectively, these observations indicate that chickens have an unusually homogeneous breast muscle structure bearing a very low aerobic capacity. broiler chickens are the most feed efficient vertebrates known with a dry weight in to wet weight out ratio of ~. it has been previously argued that the very low aerobic capacity of their muscle is a likely driver of this phenotype through avoiding paying for ‘expensive’ spare aerobic capacity and considerations of economic physiological design. indeed, it has been observed in many independent circumstances that the most feed efficient breeds, selected on either fe or muscle mass, transition from a redder, slow, aerobic ancestral muscle to a faster, whiter, less aerobic glycolytic one in the derived state. examples include mstn mutant piedmontese, belgium blue and blond d’aquitaine cattle, callipyge sheep and both large white and german landrace pigs. the transition from high myoglobin content red to low myoglobin content white fibres in pigs has been so dramatic that they are no longer considered a rich source of dietary iron and are viewed as the second ‘white’ meat after chicken, rather than the third red meat after beef and lamb. we know that these broiler muscle transcriptomes are derived almost exclusively from white, fast type iib fibres. in line with anatomical and physiological expectation the two most highly expressed genes in the muscle are the fast myosin heavy chain isoform myh1e and the regulator of anaerobic glycolysis gapdh. however, we also detected gene expression of a number of slow contractile subunits in addition to the expected fast ones. this somewhat surprising finding implies that a subset of remnant slow subunits make up a small but measurable proportion of the sarcomeric structure of the type iib fibres, perhaps such that they are a mosaic at the molecular level. the subtle but coordinated reduction of these slow twitch isoforms and associated slow machinery in the higher fe birds is in line with the broad expectation of muscle ‘whitening’ being associated with high fe described above. it also concurs with our previous observations based on a cdna microarray. the divergent pattern we observe in ca and ca may also be a simple consequence of fast versus slow fibre subunit abundance. intriguingly, while the majority of slow subunits were downregulated in hfe including myl2b, a separate myl isoform was upregulated in hfe birds. we do not have an explanation for this apparent inconsistency. one would expect the expression of mitochondrial proteins to track the fibre composition data. however, our observations on the gene expression of mitochondrial proteins are not congruent with the reduction in the slow muscle fibre subunit and myoglobin/phospholamban data in hfe birds. mitochondrial metabolism in terms of mitochondrial content we found that the mrna encoding the mitoproteome is modestly but highly significantly skewed towards the hfe birds. we propose that this indicates a slightly higher mitochondrial content and therefore a somewhat higher aerobic capacity. this is prima facie contrary to the reduction in slow fibre subunits and reduced myoglobin/phospholamban outlined above. it is also at odds with theoretical expectation and mrna/protein abundance data from numerous other species, particularly recent genome-wide screens performed on pigs divergent for rfi. it also appears contradictory to previous work showing in chickens divergent for mass, the more efficient larger birds have a lower, not higher, mass-specific metabolic rate. having said this, we believe these specific samples have been correctly characterised. these mrna data are in agreement with a mass spectrometry-based proteomic analysis carried out on a subset of the same breast muscle tissues, as the protein expression of the mitochondrial proteome was also somewhat elevated in the hfe phenotype. therefore, while the ‘big picture’ finding in production animals appears to be that selection on mass and fe reduces mitochondrial content and aerobic capacity, the data from this particular fe contrast between highly selected birds tends to point in the opposite direction. the underlying reason for this notable disparity are the source of continued investigations. one possibility is that the commercial birds have such a low mitochondrial content that they are on the verge of expressing metabolically driven muscle pathologies. this would be consistent with the recent emergence of concerning phenotypes such as ‘white striping’ and ‘wooden breast.’ if a low mitochondrial content drives these pathologies, one might speculate that subtly increasing mitochondrial content may increase fe through impact on general health. ultimately, the mitochondrial content finding has been approximated using gene expression. ideally, this finding would be validated with an independent technology such as qpcr of mtdna expressed relative to ndna, but unfortunately we no longer have adequate sample material. the general upregulation of mrna encoding the mitoproteome in hfe birds is modest in absolute terms and there are many genes that are downregulated as well as upregulated. extreme among the upregulated molecules in hfe are tdrkh, mrpl ldhd, ckmt1a encoding a gene silencer, a mitoribosome subunit, lactate dehydrogenase d and mitochondrial creatine kinase 1a. the biological implications of tdrkh for skeletal muscle are unknown, but mrpl ldhd and ckmt1a indicate an increase in mitochondrial protein synthesis, a change in redox status and increasing the transfer of high energy phosphate from the mitochondria to the cytosolic carrier, creatine. the downregulation of pdk acacb and ldhb indicate the hfe birds have a preference for burning carbohydrate instead of fat, a reduction in fatty acid uptake and combustion by mitochondria and a reduced reliance on anaerobic glycolysis. the first two point to a change in fuel combustion pattern that is consistent with a subtle fibre subunit transition towards a whiter muscle phenotype–after all, glycolytic fibres are geared to burn glycogen not fat during exercise. the reduction in the anaerobic ldhb in hfe is perhaps consistent with an elevated mitochondrial content, as the presence of higher aerobic metabolism may reduce the need to resort to anaerobic metabolism for atp production. intramuscular fat and extracellular matrix content although arguably not statistically significant sensu stricto, we wish to highlight the expression profiles of several genes previously implicated in intramuscular fat content given the important role imf plays in both meat quality and energetics. in the hfe birds there is a trend towards a small reduction in expression of two genes previously associated with amount and activity of intramuscular fat in cattle, sheep and pigs. in absolute terms plin and fabp are fold and fold downregulated in hfe respectively. plin is ranked in position when the genes are ranked by pif from negative to positive values. together, the expression of these two genes implies the hfe birds possess a slightly leaner breast muscle. in cattle muscle, these two genes tend to be tightly co-expressed with thrsp, s100g, cidea and cidec, but we did not detect expression for these additional genes in the broiler data. the tendency towards higher muscle tissue leanness in the high fe birds is unsurprising, relating to the selection for overall animal mass in the calculation of residual feed intake. using basic animal mass as the selection phenotype will inevitably select against light, anhydrous but energy-dense adipose tissue and promote heavy, watery lean muscle tissue. further, very lean meat is texturally dry and relatively flavourless–a universal problem in production animals or breeds with relatively high fe such as callipyge sheep, myostatin mutant sheep and large white pigs. additionally, a very strongly coordinated reduction in expression of genes encoding extracellular matrix proteins in hfe birds is consistent with less ‘white striping’, a meat defect that has recently appeared in poultry, particularly large birds grown too quickly. white striping comprises seams of collagenous material containing some embedded adipocytes. the coordinate down regulation in hfe birds of genes encoding ecm components in addition to those relating to lipid content can also be seen as consistent with a reduction in this pathology. progesterone signalling and fe we have characterised the anatomical and physiological differences between high and low fe broiler muscle. the hfe bird muscle exhibits a reduction in slow muscle subunits and associated machinery and a subtly elevated mitochondrial content. but what is the ultimate cause of this muscle phenotype shift? the differential network rif analysis clearly points to alterations in progesterone signalling as being the major driver. the most differentially connected transcription factors are cited pgr and gzf that is, the progesterone receptor itself, and two molecules that communicate progesterone signalling or are influenced by it. however, the rif analysis does not imply activation or repression of the pathway, the direction of change needs to be inferred by different means. in a concurrent proteomics study being conducted on the pedm broiler line, progesterone was predicted to be activated using a commercial software program, based on expression of downstream proteins. thus, two independent analytical approaches, one based on global gene expression and the other based on global protein expression both point towards progesterone as a driver of the fe phenotype. progesterone is synthesised from cholesterol, so our additional identification of “cholesterol biosynthesis” as enriched in the upregulated genes in hfe birds is noteworthy. in mammals, progesterone synthesis is known to occur in reproductive organs, adrenal glands, nervous tissue and adipose tissue. to the best of our knowledge endogenous synthesis has not been detected in mammalian skeletal muscle. the situation in birds is unknown. in fact, we could not find any published information on progesterone production nor mode of action in avian muscle. however, in cattle exogenous administration of progesterone is often used in hormone growth promotant mixes that raise fe by ~20%. indeed, the synthetic gestagen, melengestrol acetate, has been found to stimulate muscle growth in yearling heifers. furthermore, in pregnant humans progesterone alters muscle phenotype from a slow oxidative to fast glycolytic state. this slow to fast contractile skeletal isoform transition is certainly consistent with what we observe in our hfe broiler data. progesterone also promotes protein synthesis in human skeletal muscle and rat cardiac muscle. finally, a truncated isoform of the progesterone receptor localises to the mitochondrion in mammals. here, we present a first line of evidence that the progesterone receptor also localises to the mitochondria in birds, based on immunofluorescence in a quail muscle cell line. consequently, the posited role of progesterone signalling via the progesterone receptor in driving variation in fe is a further line of evidence consistent with the mitochondrion being a central player. previous work on muscle gene expression and poultry feed efficiency we have previously documented gene expression patterns in these same samples using microarray technology, finding a similar reduction in the hfe birds of genes encoding a set of slow twitch contractile subunits plus slow machinery pln and atp2a another group examined breast muscle samples from a different population of commercial broilers also using rna sequencing technology. surprisingly, some of their basic structural findings are exactly opposite to ours, in that they found the slow contractile subunits and intramuscular fat content as upregulated in the hfe birds. the reason for this very fundamental contradiction is a mystery at this stage, but may be attributable to different genetics or some other factor. finally, we wish to emphasise that the progesterone signalling prediction made by rif on our mrna data has also been independently borne out by ingenuity pathway analysis ‘upstream analysis’ on proteomics data derived from the same muscle samples. thus, analyses performed at two levels of biological organisation and two independent inference strategies have come to the same conclusion.  selection on fe at the whole bird level culminates in measurable changes in breast muscle phenotype. despite the subtlety of the changes, we have identified a panel of genes whose expression level can discriminate the two groups and whose functions include slow muscle structure and mitochondrial activity. hierarchical clustering shows the hfe birds show less dispersal than the lfe, reflecting a reduction in gene expression variation. this concept of reduced variation has been observed before as a real signal of pathway constraint, and therefore can be taken as informative in the same way as measures of differential expression and differential connectivity. the gene expression changes in the hfe birds suggest a whiter, leaner tissue which presumably reflects the well established relationship between elevating fe on the one hand but increasing the likelihood of deficits in meat quality on the other. this exact same issue has been observed in many other highly selected production animals such as large white pigs and callipyge sheep. the differential network analysis clearly points to modifications in progesterone signalling as a key driver of the difference between the two groups. the next challenge is to harness this finding, perhaps through diet or genetics, to help drive broiler fe even higher without exogenous administration of progesterone. any attempt to increase fe will need to be cogniscent of possible negative effects on associated meat quality traits or the health of the birds. additional files additional file 1: pca analysis of the raw rpm read count data, the log transformed data and the quantile normalised log data. green denote hfe, red denote lfe and blue are barred rock outgroup animals. the barred rock birds were used for collective normalisation but not formally analysed in this study. the outlier red sample is lfe bird additional file 2: the qpcr primers used to validate the de of a subset of genes prioritised by rnaseq. additional file 3: the very strong relationship between rpm and rpkm in our data, using pif as the metric for comparison. the pif based on rpm has a correlation of to the pif based on rpkm. additional file 4: the rif and rif output listed in full for the transcription factors. abbreviations fefeed efficiency hfehigh feed efficiency lfelow feed efficiency pedmpedigree broiler males drs. ron okimoto, tony vuocolo and gene wijffels provided insightful reviews of this manuscript prior to submission. funding this research was based upon work that is supported, in part, by the national institute of food and agriculture, u.s. department of agriculture, under award number 2013- the commonwealth science and industrial research organisation mcmaster fund supported wb on a visiting fellowship to australia, during which time njh and wb conceived this project. the funders had no say in the design of the experiment or its biological interpretation. availability of data and materials the data sets analysed during the current study are available from the corresponding author upon reasonable request. authors’ contributions wb and njh designed the study. njh drafted the manuscript. bwk performed the wet lab work. njh and bwk mapped reads, counted reads and normalised the data. aw generated the 3d muscle model and interpreted the output. ar performed the rif analysis. njh, bwk and wb biologically interpreted the data. kl performed the quail muscle immunofluorescence study. wb managed the project resources and provided oversight. all authors read, contributed suggestions and approved the submitted version. competing interests the authors declare that they have no competing interests. consent for publication not applicable. ethics approval 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: protocol # the animal handling and euthanasia for this project were approved by the institutional animal care and use committee.
functional groups contained in this de gene list include metabolic proteins, mitochondrial proteins, oxidative status and cholesterol homeostasis. this analysis clearly points to alterations in progesterone signalling as the major driver. application of a virtual muscle model to gene expression data pointed to a coordinated reduction in slow twitch muscle isoforms of the contractile apparatus, consistent with diminishment in associated slow machinery in the high fe animals. we contrast the pectoralis muscle transcriptomes of broilers selected from within a single genetic line expressing divergent feed efficiency in an effort to improve our understanding of the mechanistic basis of fe. surprisingly, we found that the expression of genes encoding the broiler mitoproteome is modestly–but significantly–biased towards the high fe group, suggesting a slightly elevated mitochondrial content. we applied a differential network method whose aim is to use patterns of differential co-expression to detect regulatory molecules transcriptionally rewired between the groups. we show the progesterone receptor localises to the mitochondria in a quail muscle cell line. because the progesterone receptor can localise to avian mitochondria, our data continue to point to muscle mitochondrial metabolism as an important component of the phenotypic expression of variation in broiler fe.
 we contrast the pectoralis muscle transcriptomes of broilers selected from within a single genetic line expressing divergent feed efficiency in an effort to improve our understanding of the mechanistic basis of fe.  application of a virtual muscle model to gene expression data pointed to a coordinated reduction in slow twitch muscle isoforms of the contractile apparatus, consistent with diminishment in associated slow machinery in the high fe animals. these data are in line with the repeated transition from red slow to white fast muscle fibres observed in agricultural species selected on mass and fe. surprisingly, we found that the expression of genes encoding the broiler mitoproteome is modestly–but significantly–biased towards the high fe group, suggesting a slightly elevated mitochondrial content. this is contrary to expectation based on the slow muscle isoform data and theoretical physiological capacity arguments. reassuringly, the extreme most de genes can successfully cluster the individuals into the appropriate fe treatment group. functional groups contained in this de gene list include metabolic proteins, mitochondrial proteins, oxidative status and cholesterol homeostasis. we applied a differential network method whose aim is to use patterns of differential co-expression to detect regulatory molecules transcriptionally rewired between the groups. this analysis clearly points to alterations in progesterone signalling as the major driver. we show the progesterone receptor localises to the mitochondria in a quail muscle cell line.  progesterone is sometimes used in the cattle industry in exogenous hormone mixes that lead to a ~20% increase in fe. because the progesterone receptor can localise to avian mitochondria, our data continue to point to muscle mitochondrial metabolism as an important component of the phenotypic expression of variation in broiler fe. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords feed efficiencymitochondriaprogesteronehttp://dx.doi.org/ /100000199u.s. department of agriculture2013-01953bottje walter issue-copyright-statement© the author 2017 