Kong B and Hudson and Feed Efficiency

Research article
Open Access
Published: 24 February 2017

Progesterone signalling in broiler skeletal muscle is associated with divergent feed efficiency

Byung-Whi Kong¹,
Antonio Reverter²,
Ashley J. Waardenberg^2,3,
Kentu Lassiter¹ &
…
Nicholas J. Hudson⁴

BMC Systems Biology volume 11, Article number:29 (2017) Cite this article

2777 Accesses
21 Citations
Metrics details

Abstract

Background

We contrast the pectoralis muscle transcriptomes of broilers selected from within a single genetic line expressing divergent feed efficiency (FE) in an effort to improve our understanding of the mechanistic basis of FE.

Results

Application of a virtual muscle model to gene expression data pointed to a coordinated reduction in slow twitch muscle isoforms of the contractile apparatus (MYH15, TPM3, MYOZ2, TNNI1, MYL2, MYOM3, CSRP3, TNNT2), consistent with diminishment in associated slow machinery (myoglobin and phospholamban) 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 699 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 40 most DE genes can successfully cluster the 12 individuals into the appropriate FE treatment group. Functional groups contained in this DE gene list include metabolic proteins (including opposing patterns of CA3 and CA4), mitochondrial proteins (CKMT1A), oxidative status (SEPP1, HIG2A) and cholesterol homeostasis (APOA1, INSIG1). We applied a differential network method (Regulatory Impact Factors) 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 (via the receptor PGR) as the major driver. We show the progesterone receptor localises to the mitochondria in a quail muscle cell line.

Conclusions

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.

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 [1]. Any inefficiency not only affects the bottom line but also negatively impacts resource usage (e.g. water and energy) and waste production (e.g. urine). 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 ~1.5 units dry weight intake per unit wet weight gain [2]. 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 [3, 4]. This behaviour explains the functionally unusual breast muscle metabolism of the ancestral birds such as the Red Junglefowl (Gallus gallus) 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 (PedM) from a population of highly inbred, commercial animals subject to very strong historical selection for high muscle mass and elevated feed efficiency (FE) [5]. The significant difference in FE (g gain/g feed) between the two treatment groups is in the order of 1.4 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 [5], 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 [3]. 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 (HFE) and low FE (LFE) 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 [5–20] and overall tissue mitochondrial content [21] 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 (DE) and a recently developed differential connectivity (DC) algorithm called Regulatory Impact Factor (RIF) analysis to screen the data in an unbiased fashion [22–24]. 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 [25], 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 12 male broilers divergent in feed efficiency as previously described [5]. 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 12 broilers were derived from a single genetic line of a commercial Cobb-Vantress Inc. population. Both the HFE (n = 6) and LFE (n = 6) birds were sub-sampled from the extremes of an original group of 100 birds. In turn, those 100 birds were selected as the most efficient from a larger population of 300 birds. The high FE group had greater body weight gain from 6 to 7 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 0.65 ± 0.01 compared to 0.46 ± 0.01 (g gain to g feed intake) 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 (East Lansing, MI) for Illumina HiSeq 100 base pair paired end read sequencing. In brief, we used CLC Genomics Workbench 8 to map the reads to Gallus gallus genome assembly version 4. The CLC software follows the analytical pipeline recommended in [26]. We log2 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 12 individual samples could discriminate into the treatment groups of origin was confirmed by Principal Component Analysis (Additional file 1).

Differential expression

To compare the 2 treatment groups we plotted the MA (i.e. Minus Average) plot comparing HFE birds to LFE birds (Fig. 1 panel a). The x axis, A, represents the average abundance in the 2 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.

mRNA encoding the mitoproteome

We downloaded the complete human mitoproteome from http://mitominer.mrc-mbu.cam.ac.uk/release-3.1/begin.do generating a list of 1046 nuclear and mitochondrial genes encoding mitochondrial proteins. We converted the protein names to gene names, and found 699 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 (Fig. 1 panel b). The skew in distribution away from the null expectation of 50: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 (PIF)

We next computed a modified DE called PIF (DE multiplied by abundance) 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 (Fig. 1 panel c). It de-emphasises lowly abundant genes that are inherently noisy as they approach the detection limit of the technology. It also draws together the 2 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 [23]. Also, PIF is a foundation of the RIF differential network analysis (below) 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 [27].

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). VMus3D, described in detail here [28, 29], 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 (≥2), green (< 2 and ≥ 1), yellow (> −1 and < 1), orange (≤ −1 and > −2) and red (≤ −2) (Fig. 2).

Permut Matrix clustering on extreme PIF

The top 20 up and top 20 downregulated genes (by PIF) were identified as described above. A matrix was formed comprising as many rows as genes (40) and as many columns as birds (12) with each cell containing the normalised gene expression. This matrix was imported into PermutMatrix software [30] for hierarchical clustering (Fig. 3).

Regulatory Impact Factor (RIF) 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 (here HFE versus LFE). RIF1 and RIF2 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 RIF1 and RIF2 as previously described [22–24]. This procedure exploits global patterns of differential co-expression (or 'differential wiring') 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 (r = 1, 2, …898) were computed using the following formulae:

$$ \mathrm{R}\mathrm{I}\mathrm{F}{1}_r=\frac{1}{{\mathrm{n}}_{\mathrm{DE}}}{\displaystyle \sum_{j=1}^{j={\mathrm{n}}_{\mathrm{DE}}}{x}_j\times {d}_j\times {\mathrm{DC}}_{rj}^2} $$

and

$$ \mathrm{R}\mathrm{I}\mathrm{F}{2}_r=\frac{1}{{\mathrm{n}}_{\mathrm{DE}}}{\displaystyle \sum_{j=1}^{j={\mathrm{n}}_{\mathrm{DE}}}\left[{\left({x}_j^{\mathrm{HFE}}\times {r}_{rj}^{\mathrm{HFE}}\right)}^2-{\left({x}_j^{\mathrm{LFE}}\times {r}_{rj}^{\mathrm{LFE}}\right)}^2\right]} $$

where n_DE 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; DC_rj was the differential co-expression between the r-th regulator and the j-th DE gene, and computed from the difference between r ^HFE _rj and r ^LFE _rj , the correlation co-expression between the r-th regulator and the j-th DE gene in the HFE and LFE samples, respectively; finally, x ^HFE _j and x ^LFE _j represented the average expression of the j-th DE or TS gene in the HFE and LFE samples, respectively. The RIF1 and RIF2 output were simultaneously plotted (Fig. 4).

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 (MYH15, TPM3, MYOZ2, TNNI1, MYBPC1), slow muscle associated proteins (MB, CA3) and muscle fat content (PLN and FABP4). We used PCR cycling conditions as previously described [5, 11]. Primers are detailed in Additional file 2.

In brief, one microgram of total RNA was obtained from 6 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 (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 PowerUpTM SYBR® Green Master Mix (ThermoFisher Scientific, Waltham, MA). The specific oligonucleotide primers were designed by the PRIMER3 program (http://frodo.wi.mit.edu). 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 HFE and LFE muscles were used for relative quantification by fold-change and statistical significance.

Quail muscle (QM7) cell culture and treatments

QM7 cells were grown in M199 medium (Life Technologies, Grand Island, NY) with 10% fetal bovine serum (Life Technologies), 10% tryptose phosphate (Sigma-Aldrich, St. Louis, MO), and 1% penicillin-streptomycin (Biobasic, Amherst, NY) at 37 °C under a humidified atmosphere of 5% CO2 and 95% air. At 80–90% confluence, cells were synchronized overnight in serum-free medium.

Progesterone receptor immunofluorescence

Immunofluorescence was performed as previously described [31] with modifications. Briefly, cells were grown to 50–60% confluence in chamber slides (Lab-Tek, Hatfield, PA) followed by methanol fixation for 10 min at 20 °C. Serum free blocking buffer (Dako, Carpinteria, CA), and the cells incubated with rabbit anti-progesterone receptor antibody overnight at 4 °C. Cells were treated with MitoTracker Deep Red (Molecular Probes, Inc., Eugene OR) and visualized with Alexa Fluor 488-conjugated secondary antibody (Molecular Probes, Life Technologies). After DAPI counterstaining, a cover slip was placed on slides with Vectashield (Vector Laboratories, Burlingame, CA). Images were obtained using the Zeiss Imager M2 with a 20X Plan-APOCHROMAT 20X/0.8 objective and a 100X EC PLANNEOFLUOR 100X/1.3 oil objective. The Alexa Fluor 488 fluorophore was observed through filter set 38 1031-346 with an excitation of BP 470/40, beamsplitter of FT 495, and emission spectrum of BP 525/50. Differential interference contrast images were collected using DIC M27 condensers. The Alexa Fluor 488 fluorophore was excited for 500 ms prior to capturing each image using an Axio Cam MR3 camera. All analysis was performed using AxioVision SE64 4.9.1 SP1 software (Carl Zeiss Microscopy 2006–2013).

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 [27]. Here, we examined functional enrichment in the extreme 260 up and downregulated genes (nominal 5% or 520 out of the total 10, 416) 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 6 HFE and 6 LFE broilers. In total, 804 million 100 bp sequences were obtained with an average of 67 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 12,814 genes. We screened for those mRNA with no missing values across the 12 samples. This gave a total 10,412 genes which were used for all subsequent analysis. Both RPM and RPKM values were computed and compared by correlation analysis (Additional file 3). To achieve this we calculated the PIF for each gene (abundance multiplied by DE) and correlated the PIF scores by the two methods. Given the high correlation (0.93) and the concern that the additional normalisation used in generating RPKM values can theoretically introduce a variety of problems [32], 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. Log2 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 (Additional file 1). 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 0.

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 (below) 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 260 upregulated genes (top 2.5% out of 10,412) in HFE by PIF is "cholesterol biosynthesis" (P = 0.00027; FDR q value = 0.1) based on the presence of DHCR24, INSIG1, FDPS, SC5D and MVD. The top enrichment among the 260 downregulated genes in HFE by PIF is "extracellular matrix organisation" (P = 1.03 E-26; FDR q value = 1.24 E-22) based on the presence of FMOD, ITGA4, ACTN1, SMOC2, FN1, ITGA8, MFAP5, CMA1, TNC, GSN, TLL1, NRXN1, COL1A2, COL3A1, CCDC80, DCN, LTBP3, CYP1B1, ITGB2, NCAM1, MMP2, LOX, CTSK, COL6A2, TGFB1, COL6A3, COL22A1, CD44, CTSS, COL5A2, LCP1, COL6A1, LOXL2, COMP, ANXA2, THBS1, ABI3BP, COL12A1, ECM2, COL11A1, VCAN, POSTN and NID1.

The top 20 upregulated and downregulated genes by PIF are tabulated with brief functional descriptions and absolute fold changes in DE (Table 1).

Table 1 The 40 most extreme genes by Phenotypic Impact Factor in High FE versus Low FE birds. These 40 genes have been reordered manually into functional groupings. The absolute DE expressed as a fold change has also been reported

Full size table

We also cross-referenced the genes and direction of change in Table 1 to the previous microarray study on the same samples [5]. Of the 14 genes for which we could find exact matches, 12 showed the same direction of change by the two different mRNA quantitation technologies.