# Molecular regulation of high muscle mass in developing Blonde d'Aquitaine cattle fetuses

Isabelle Cassar-Malek, Céline Boby, Brigitte Picard, Antonio Reverter-Gomez, Nicholas Hudson

► To cite this version:

Isabelle Cassar-Malek, Céline Boby, Brigitte Picard, Antonio Reverter-Gomez, Nicholas Hudson. Molecular regulation of high muscle mass in developing Blonde d'Aquitaine cattle fetuses. Biology Open, 2017, 6 (10), pp.bio.024950. 10.1242/bio.024950. hal-01607200

HAL Id: hal-01607200

https://hal.science/hal-01607200

Submitted on 26 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

# Molecular regulation of high muscle mass in developing Blonde d'Aquitaine cattle foetuses

Isabelle Cassar-Malek<fnref n="1" />, Céline Boby<fnref n="1" />, Brigitte Picard<fnref n="1" />, Antonio Reverter<fnref n="2" /> and Nicholas J. Hudson³,*

## ABSTRACT

The Blonde d'Aquitaine (BA) is a French cattle breed with enhanced muscularity, partly attributable to a *MSTN* mutation. The BA *m. Semitendinosus* has a faster muscle fibre isoform phenotype comprising a higher proportion of fast type IIX fibres compared to age-matched Charolais (CH). To better understand the molecular network of modifications in BA compared to CH muscle, we assayed the transcriptomes of the *m. Semitendinosus* at 110, 180, 210 and 260 days postconception (dpc). We used a combination of differential expression (DE) and regulatory impact factors (RIF) to compare and contrast muscle gene expression between the breeds. Prominently developmentally regulated genes in both breeds reflected the replacement of embryonic myosin isoforms (*MYL4*, *MYH3*) with adult isoforms (*MYH1*) and the upregulation of mitochondrial metabolism (*CKMT2*, *AGXT2L1*) in preparation for birth. However, the transition to a fast, glycolytic muscle phenotype in the MSTN mutant BA is detectable through downregulation of various slow twitch subunits (*TNNC1*, *MYH7*, *TPM3*, *CSRP3*) beyond 210 dpc, and a small but consistent genome-wide reduction in mRNA encoding the mitoproteome. Across the breeds, *NRIP2* is the regulatory gene possessing a network change most similar to that of *MSTN*.

**KEY WORDS:** Skeletal muscle, Myostatin, Mitochondria, Genomics, Myogenesis

## BACKGROUND

Feed is the single largest cost in many intensive livestock farming enterprises including beef feedlotting (Robinson and Oddy, 2004). Consequently, the food production industry is interested in enhancing the feed conversion efficiency of its major production species. Accelerating improvements in these areas beyond what can be achieved by traditional artificial selection on phenotype first requires a base level of understanding. What physiological and developmental processes are taking place, and what are their functional limits? Fortunately, domesticated species provide opportunity to uncover the molecular mechanisms underpinning phenotypes of interest (Andersson and Georges, 2004). Agricultural

scientists have numerous breeds at their disposal that, through generations of selective breeding, possess strong phenotypes expressed against a relatively stable background genome. The mechanistic basis of feed efficiency is not known, but variation in feed intake, digestion, metabolism, activity and thermogenesis have all been implicated (Herd and Arthur, 2009). Focussing on skeletal muscle tissue, mammals possess a range of fibre types with different functional properties. Interestingly, we already know that highly feed-efficient breeds of all the major production species, such as Piedmontese cattle (Lehnert et al., 2007), Belgian Blue cattle (Bouley et al., 2005; Deveaux et al., 2001; Hanset et al., 1987), Large White pigs (Lefaucheur and Ecolan, 2005; Lefaucheur et al., 2004, 2011), German Landrace pigs (Rehfeldt et al., 2008) and Callipyge sheep (Jackson et al., 1997a,b,c), collectively possess exactly the same shift towards a whiter type IIB/IIX phenotype compared to their wild counterparts and less efficient breeds. Broiler chicken muscle can be seen as the most extreme example of this phenotype (Kiessling, 1977).

A bioenergetic feature of type IIB/X fibres is a lower mitochondrial content than type I and IIA. This implies that the glycolytic fibre shift documented above will result in muscle tissue mitochondrial content decreasing. Direct measures of mitochondrial content have not been routinely taken. However, Vincent et al. (2015) recently discovered that when divergently selecting on residual feed intake, the more efficient pigs possessed a lower muscle oxidative metabolism and mitochondrial content. Similarly, in Yorkshire boars, low muscle mitochondrial content has recently been associated with higher feed efficiency (Fu et al., 2017). In addition, double-muscled cattle or cattle selected for high growth potential and feed efficiency also possess less oxidative muscles (Sudre et al., 2005; Picard et al., 2006), as do Charolais (CH) divergent for muscle mass (Sudre et al., 2005). A reduction in oxidative capacity has previously been argued to improve metabolic efficiency on the grounds that it is inefficient to possess too much fuel burning capacity (Hudson, 2009). The engineering argument runs as follows: it is not economical to pay for the construction (biogenesis), maintenance (trans-membrane proton gradient) and load (space that could be occupied by other cellular structures) of energetically expensive mitochondrial capacity that is not needed. Further, Bottje and Carstens (2009) have reviewed the impact of isolated mitochondrial performance and its impact on poultry and livestock feed, efficiency documenting associations in several species and circumstances. These lines of evidence support a role for mitochondrial function in contributing to variation in feed efficiency.

However, the genetic bases of the ‘whitening’ of livestock muscle and other biological aspects of efficiency remain a mystery. Other than some progress in identifying causal mutations such as the various myostatin (*MSTN*) mutants, the detailed rewiring of gene regulatory networks that manifest these anatomical and physiological changes through muscle development are largely unknown. Here, we attempt to shed more light on this question by

<figure><img src="image_2.png" /></figure>

<footnote n="1">¹UMR1213 Herbivores, Institut National de la Recherche Agronomique, VetAgro Sup, 63122 Saint Genès Champanelle, Clermont-Ferrand F-63122, France.</footnote>

<footnote n="2">²Agriculture, Commonwealth Science and Industrial Research Organisation, Queensland Bioscience Precinct, St. Lucia, Brisbane 4075, Australia. ³School of Agriculture and Food Sciences, University of Queensland, Brisbane 4075, Australia.</footnote>

<footnote>*Author for correspondence (n.hudson@uq.edu.au)</footnote>

<footnote>N.J.H., 0000-0002-3549-9396</footnote>

<footnote>This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.</footnote>
contrasting the genome-wide muscle transcriptomes of two breeds of cattle. A *MSTN* mutation has been identified in the Blonde d'Aquitaine (BA) breed (Bouyer et al., 2014) that culminates in greater muscle mass and feed efficiency than the comparison CH, but not all individuals carry the mutation. In line with the other examples given above, BA also has a more glycolytic, less oxidative muscle metabolism and a higher proportion of fast glycolytic IIX fibres (Listrat et al., 2001; Picard and Cassar-Malek, 2009). We analysed *m. Semitendinosus* samples at three prenatal time points coincident with secondary myogenesis (110 days postconception (dpc), onset of functional differentiation (180 dpc), completion of functional differentiation (210 dpc), in addition to a later time point prior to birth when mature muscle has been established (260 dpc). Through the application of several methods including a differential network algorithm we have identified molecules behaving differently between the two breeds during pre-natal development, when adult muscle potential is defined.

## RESULTS

### Histology

Representative histology sections of the *m. Semitendinosus* for the two breeds of cattle were analysed at each time point (data not shown). By 260 dpc, comparison of fibre frequencies in the muscle of BA versus CH at 260 dpc (Fig. S1) showed that BA muscle contained a lower proportion of slow (I) fibres than CH (8.2% versus 15.8%), of which 77% still expressed the foetal myosin heavy chain (MyHC). Moreover, BA muscle contained a higher proportion of fast (II) fibres than CH (92.8% versus 76.4%), of which 47% still expressed the foetal MyHC. Lastly, the proportion of IICF fibres (transition fibres containing I, IIa and F MyHC) was lower in the muscle of BA versus CH. This clearly indicated a predominant fast glycolytic phenotype in the BA. Accordingly, at 260 dpc, the activity of oxidative enzymes was lower in the BA than in the CH. The higher content of foetal MyHC also suggested that the differentiation of BA fibres originating from both the first and secondary generation was delayed compared to that of CH fibres.

### Developmental patterns of differential expression (DE) in both breeds

We determined normalized mean expression values for 45,220 probes mapped to 17,646 bovine protein coding genes. After quality checking for a minimal level of expression in seven of the eight

treatments, we were left with a refined list containing 19,100 probes mapping to 10,138 genes. When comparing the breeds for patterns of high developmental regulation (including but not limited to stepwise increases and decreases across development), a number of genes were identified. This was achieved by computing the standard deviation of expression values across the four time points on a breed-specific basis. High standard deviations reflect high developmental regulation. The replacement of embryonic myosin heavy chain isoforms (*MYL4*, *MYH3*) with adult isoforms (*MYH1*) was observed. Further, we detected an upregulation of creatine kinase (*CKMT2*), a marker of mitochondrial metabolism, in preparation for birth (Table 1). We also detected strong developmental regulation of an unannotated Open Reading Frame (ORF), *C13H20ORF194*, *AGXT2L1* also has high developmental regulation, but missing values at some developmental stages in both breeds meant it could not be reported in Table 1. It can be seen that the gene expression pattern change in these genes is dominated by upregulation as development progresses. However, we wish to emphasize that this is not a result forced by the analysis. Further exploration of the data finds that *SNAP91* ranked in 11th position has a pattern of decline in gene expression over development. This gene encodes a protein enriched in nervous tissue so its expression pattern in muscle samples presumably reflects some aspect of the neuromuscular junction and associated motor innervation.

A subset of these genes has also been highlighted in the MA plots represented in Fig. 1. In these plots, in line with all our analysis, the DE is computed by BA minus CH. Patterns of developmental expression change (as opposed to breed DE) can be monitored by following the change in position on the x axis (average abundance) of the plots as one proceeds through development, 110 days to 260 days.

### Breed contrast using DE

As part of our strategy to detect patterns of DE between the breeds at each time point individually we plotted the MA plots as four separate panels (Fig. 1). We only plotted those data minimally differentially expressed (>0.5 log2 DE) to simplify presentation. This 0.5 log2 DE is an arbitrary decision and is not meant to reflect a judgment of statistical significance. In all the plots we have highlighted functional groups of genes we have previously found to be tightly co-expressed in bovine muscle and where the encoded proteins have fundamental roles in embryonic muscle development (*MYL4* and *MYH3*), adult

<table><caption>Table 1. Top 10 highly developmentally regulated genes in both breeds</caption><thead><tr><th>Gene</th><th>Probe</th><th>BA110</th><th>BA180</th><th>BA210</th><th>BA260</th><th>CH110</th><th>CH180</th><th>CH210</th><th>CH260</th><th>BA developmental regulation (estimated by standard deviation)</th><th>CH developmental regulation (estimated by standard deviation)</th></tr></thead><tbody><tr><td><strong>MYH1</strong></td><td>A_73_P332346</td><td>0.60</td><td>2.41</td><td>7.38</td><td>10.84</td><td>0.16</td><td>3.73</td><td>8.71</td><td>10.64</td><td>4.67</td><td>4.75</td></tr><tr><td><strong>MYBPC2</strong></td><td>A_73_P068061</td><td>0.56</td><td>5.47</td><td>8.44</td><td>10.39</td><td>-</td><td>6.34</td><td>8.16</td><td>10.51</td><td>4.28</td><td>*</td></tr><tr><td><strong>C13H20ORF194</strong></td><td>A_73_107199</td><td>1.40</td><td>4.40</td><td>6.84</td><td>8.88</td><td>1.23</td><td>5.07</td><td>6.80</td><td>8.67</td><td>3.22</td><td>3.16</td></tr><tr><td><strong>ACTN3</strong></td><td>A_73_P416836</td><td>2.83</td><td>6.56</td><td>8.63</td><td>10.03</td><td>3.04</td><td>6.99</td><td>8.60</td><td>10.10</td><td>3.12</td><td>3.03</td></tr><tr><td><strong>RPL3L</strong></td><td>A_73_118617</td><td>4.01</td><td>6.90</td><td>9.17</td><td>10.59</td><td>4.02</td><td>7.38</td><td>9.05</td><td>10.34</td><td>2.87</td><td>2.73</td></tr><tr><td><strong>TMOD4</strong></td><td>A_73_P350016</td><td>3.68</td><td>7.17</td><td>8.71</td><td>9.81</td><td>4.09</td><td>7.65</td><td>8.54</td><td>10.05</td><td>2.67</td><td>2.52</td></tr><tr><td><strong>ART3</strong></td><td>A_73_P262561</td><td>2.35</td><td>5.18</td><td>6.45</td><td>8.68</td><td>2.59</td><td>5.64</td><td>6.94</td><td>8.76</td><td>2.63</td><td>2.60</td></tr><tr><td><strong>CKMT2</strong></td><td>A_73_P251541</td><td>0.50</td><td>1.57</td><td>4.02</td><td>6.38</td><td>0.76</td><td>2.45</td><td>4.26</td><td>6.79</td><td>2.62</td><td>2.58</td></tr><tr><td><strong>BTC</strong></td><td>A_73_P041286</td><td>1.40</td><td>5.49</td><td>6.72</td><td>7.05</td><td>2.11</td><td>5.87</td><td>6.14</td><td>6.77</td><td>2.60</td><td>2.11</td></tr><tr><td><strong>SLC16A3</strong></td><td>A_73_112950</td><td>2.17</td><td>6.62</td><td>7.69</td><td>7.46</td><td>2.30</td><td>6.67</td><td>7.90</td><td>6.48</td><td>2.58</td><td>2.44</td></tr></tbody></table>

Gene expression values are expressed as log2. Developmental regulation was estimated by computing the standard deviation of the mean values across the four time points for each breed separately.

<footnote>*Not computable because of one missing value.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** MA plots at all four time points, highlighting a subset of developmentally regulated and breed-affected genes. These include muscle structural isoforms, mitochondrial metabolism and the causal mutation *MSTN*. The slow myosin isoforms being downregulated in the BA is consistent with a shift to a faster muscle phenotype in the mutant BA by 260 days. The strong developmental regulation in both breeds of *C13H20ORF194* and *MYH1* is also apparent. *MSTN* is consistently upregulated in the more muscular BA. We only included data from those genes >0.5 log DE in an effort to simplify the plots. The 0.5 threshold is an entirely arbitrary decision and is not meant to represent significance.</figcaption></figure>

muscle fibre composition (*TNNC1*, *MYH7* and *MYL2*), mitochondrial
metabolism (*CKMT2* and *AGXT2L1*) and triacylglyceride synthesis
(*FABP4*, *PLIN1* and *DGAT2*). The coordination in expression of
these independently measured genes lends more credibility to the
biological implications. Further, we also annotated on each plot an
adult muscle myosin (*MYH1*) and an ORF (*C13H20ORF194*) with
strong developmental regulation. The latter was previously identified
in an independent bovine muscle data set. Finally, we further
annotated genes with extreme DE values in at least one of the plots.
This includes *XIST*, a non-coding RNA involved in X chromosome
silencing.

In order to get a broader perspective on DE between the breeds, we
established which genes showed cumulatively the most DE between
breeds across all time points combined. We summed the absolute
mean log2 DE at each of the four developmental stages and ranked
the output (see Table S1 for the full data; Table 2 contains the extreme
20). The entire ranked list (10,138 unique genes) was imported into
GOrilla (Eden et al., 2009) for functional enrichment analysis
assessed by hypergeometric statistics. Gorilla used 49% of the genes
in the analysis based on those gene names it recognized and for which
it had accompanying functional information. GOrilla found ‘type I
interferon signalling pathway’ to be highly enriched based on genes
<table><caption>Table 2. Top 20 most differentially expressed genes summed across all four developmental time points when contrasting the two breeds</caption><thead><tr><th>Gene</th><th>Probe</th><th>110DE</th><th>180DE</th><th>210DE</th><th>260DE</th><th>Cumulative DE</th><th>Gene Ontology</th></tr></thead><tbody><tr><td>IFI27</td><td>A_73_P036551</td><td>3.78</td><td>2.34</td><td>2.68</td><td>0.36</td><td>9.18</td><td>Apoptotic signalling; mitochondrial</td></tr><tr><td>OAS2</td><td>A_73_P035131</td><td>4.13</td><td>1.91</td><td>2.07</td><td>-0.13</td><td>8.26</td><td>Innate immune response</td></tr><tr><td><strong>MSTN</strong></td><td>A_73_120254</td><td>1.87</td><td>1.82</td><td>2.18</td><td>2.20</td><td>8.09</td><td>Causal mutation in BA, negative regulator of muscle mass</td></tr><tr><td>XAF1</td><td>A_73_112973</td><td>3.40</td><td>2.71</td><td>1.37</td><td>0.43</td><td>7.93</td><td>Apoptotic signalling; mitochondrial</td></tr><tr><td>ISG15</td><td>A_73_P122741</td><td>3.81</td><td>1.77</td><td>1.76</td><td>0.15</td><td>7.51</td><td>Cytokine-mediated pathway</td></tr><tr><td>PLAC8</td><td>A_73_115944</td><td>1.52</td><td>2.84</td><td>1.66</td><td>1.29</td><td>7.32</td><td>Apoptotic signalling</td></tr><tr><td>IFI44L</td><td>A_73_106244</td><td>3.43</td><td>1.76</td><td>1.60</td><td>0.50</td><td>7.30</td><td>Immune response</td></tr><tr><td>LOC616942</td><td>A_73_P098701</td><td>1.75</td><td>2.13</td><td>1.83</td><td>1.28</td><td>7.01</td><td>Unknown function</td></tr><tr><td>BOLA-A</td><td>A_73_P032576</td><td>1.70</td><td>2.11</td><td>1.58</td><td>1.28</td><td>6.68</td><td>Major histocompatibility</td></tr><tr><td>MX1</td><td>A_73_103867</td><td>3.15</td><td>1.39</td><td>1.67</td><td>0.32</td><td>6.55</td><td>Apoptosis</td></tr><tr><td>BOLA</td><td>A_73_P062321</td><td>1.79</td><td>2.19</td><td>1.36</td><td>1.12</td><td>6.47</td><td>Major histocompatibility</td></tr><tr><td>PSMB9</td><td>A_73_118729</td><td>2.18</td><td>2.04</td><td>1.03</td><td>1.07</td><td>6.34</td><td>Immunoproteosome, processing class I MHC</td></tr><tr><td>PMEL</td><td>A_73_107037</td><td>-1.23</td><td>-1.91</td><td>-1.82</td><td>-1.30</td><td>6.27</td><td>Melanin biosynthesis</td></tr><tr><td>CCL5</td><td>A_73_100958</td><td>1.24</td><td>2.23</td><td>1.38</td><td>1.32</td><td>6.19</td><td>Inflammation</td></tr><tr><td>CTDSP2</td><td>A_73_106483</td><td>1.31</td><td>2.01</td><td>1.15</td><td>1.62</td><td>6.10</td><td>Cell cycle regulation</td></tr><tr><td>LEAP2</td><td>A_73_P051936</td><td>1.80</td><td>0.53</td><td>0.38</td><td>-2.88</td><td>5.61</td><td>Antimicrobial peptide</td></tr><tr><td>KYNU</td><td>A_73_105245</td><td>1.39</td><td>1.43</td><td>1.68</td><td>0.71</td><td>5.22</td><td>Tryptophan metabolism</td></tr><tr><td>MIC1</td><td>A_73_P062286</td><td>1.40</td><td>1.74</td><td>1.18</td><td>0.84</td><td>5.17</td><td>MHC class I</td></tr><tr><td>UXS1</td><td>A_73_P164082</td><td>1.27</td><td>1.43</td><td>1.34</td><td>1.07</td><td>5.13</td><td>UDP-glucuronate decarboxylase activity</td></tr><tr><td>CACNA1A</td><td>A_73_P101386</td><td>-3.07</td><td>0.17</td><td>-1.05</td><td>0.75</td><td>5.05</td><td>Voltage-dependent calcium channel</td></tr></tbody></table>

present near the top of the list (FDR Q value<1.25 e⁻¹¹) in the ‘Process’ ontology. This enrichment was derived from the following high ranking genes: XAF1, IFITM1, RSAD2, ISG15, USP18, MX1, IFI27, PSMB8, OAS2 and IFI6.

Overall, *MSTN* was the third most differentially expressed gene using this approach. *MSTN* was expressed at relatively higher levels in the BA animals at all time points. The breed difference in expression level of the *MSTN* downstream canonical signalling pathway could not be determined because the *MSTN* receptor *ACTVIIB* and major downstream signalling molecule *SMAD7* did not pass the quality checking criteria in at least one of the time points. The interferon signalling molecules are prominent among the differentially expressed genes. It is possible that the patterns of DE across development are quite different in these extreme performers identified by this metric. For example, using this approach, one very divergent time point could give the same ranking as more modest differences at each time point. Consequently, all extreme values need assessment on a case by case basis to better infer the biological implications. The highlighting of the extreme 20 in Table 2 is an arbitrary decision but falls well within a nominal *P*-value of 0.01, which could be calculated by prioritising the extreme 1% from the full data set.

The downregulation of mRNA encoding slow twitch subunits *TNNC1* and *MYH7* is apparent in the mutant BA at 210 days (Table S2) and 260 days (Table 3). While the DE of these mRNA may appear modest in absolute terms, it should be emphasized that (1) they do rank as extreme by a modified DE that accounts for data distribution [see Phenotypic Impact Factor (PIF) below]; (2) they are coordinately expressed so the collective observation is arguably more convincing than a single mRNA taken in isolation; (3) the same gene expression change is observed in other *MSTN* mutant models; (4) we have phenotyped these animals at the protein level, and the mRNA and protein data are mutually reinforcing. Finally, the highly divergent expression of the female X chromosome silencing non-coding RNA *XIST* at many time points reflects a gender imbalance in the experimental design.

## Breed contrast using PIF

We have previously found that multiplying the differentially expressed genes by their average abundance in the two treatments

(PIF) is a robust preliminary step for functional enrichment analysis, as well as being an analytical step in the regulatory impact factor (RIF) analysis differential connectivity (DC) analysis (Hudson et al., 2009a, 2012) described later. PIF emphasizes those genes that are substantially differentially expressed as well as abundant and de-emphasizes those molecules that may be highly differentially expressed but so lowly expressed that they approach the detection limit of the technology where the data are noisier. First of all, we computed PIF on a time-point specific basis for all genes (Table S3). The functional enrichments by GOrilla for the extreme 1% from the entire gene lists ranked on PIF are listed in Table S2.

We next focussed more specifically on 260 dpc, when mature muscle has been developed and the postnatal growth potential for the animal has been set up. *MSTN* receives the fifth highest score for genes upregulated in BA (Table 3). Of the four genes awarded higher scores than *MSTN*, *XIST* generates a non-coding RNA expressed in female tissue to silence expression from one of the two X chromosomes. *JSP.1* encodes an MHC class I protein. *UPK2* encodes a membrane protein that is present in endoplasmic reticulum and plasma membrane. *CTDSP2* encodes a protein that co-regulates androgen receptor activity. The fibre shift to a faster glycolytic phenotype in the *MSTN* mutant BA is apparent through relative downregulation of *TNNC1* and *MYH7* (Table 3).

## Comparison of muscle physiology to MSTN mutant Piedmontese

We attempted to find commonalities between the observed muscle structural DE (driven by the *TNNC1*, *MYH7*, *TPM3* and *CSRP3* observations described above) with that made in a different *MSTN* mutant system from an entirely independent animal trial. Exploration of the birth time point in an equivalent comparison (Hudson et al., 2009a) between *MSTN* mutant Piedmontese cross bred animals versus *MSTN* wild-type Wagyu cross bred animals found the same slow twitch subunits downregulated in the mutant breed (*TPM3*, *CSRP3*, *MYH7B*, *TNNC1*). The exception is *MYL2*, which is not differentially expressed in the BA versus CH (Table S4).
<table><caption>Table 3. Top 20 abundant up- and downregulated genes in BA minus CH at 260 dpc</caption><thead><tr><th>Gene</th><th>Probe</th><th>Abundance (log2)</th><th>DE (log2)</th><th>PIF</th></tr></thead><tbody><tr><td>XIST</td><td>A_73_119243</td><td>6.80</td><td>5.53</td><td>37.6</td></tr><tr><td>JSP.1</td><td>A_73_P296616</td><td>8.81</td><td>3.47</td><td>30.6</td></tr><tr><td>UPK2</td><td>A_73_P329141</td><td>5.90</td><td>2.32</td><td>13.7</td></tr><tr><td>CTDSP2</td><td>A_73_106483</td><td>8.08</td><td>1.62</td><td>13.1</td></tr><tr><td>MSTN</td><td>A_73_P036291</td><td>5.30</td><td>2.09</td><td>11.1</td></tr><tr><td>BOLA</td><td>A_73_P062321</td><td>8.16</td><td>1.12</td><td>9.2</td></tr><tr><td>SCAMP5</td><td>A_73_112884</td><td>5.34</td><td>1.68</td><td>9.0</td></tr><tr><td>ANKRD1</td><td>A_73_103476</td><td>9.60</td><td>0.89</td><td>8.5</td></tr><tr><td>LOC616942</td><td>A_73_P108026</td><td>8.14</td><td>1.00</td><td>8.1</td></tr><tr><td>OTUD1</td><td>A_73_P143321</td><td>6.94</td><td>1.11</td><td>7.7</td></tr><tr><td>PWWP2B</td><td>A_73_119461</td><td>8.20</td><td>0.91</td><td>7.5</td></tr><tr><td>GPR37L1</td><td>A_73_117230</td><td>11.66</td><td>0.64</td><td>7.4</td></tr><tr><td>U2AF2</td><td>A_73_P379036</td><td>9.98</td><td>0.73</td><td>7.3</td></tr><tr><td>NA</td><td>A_73_P359701</td><td>10.75</td><td>0.66</td><td>7.1</td></tr><tr><td>ZMAT5</td><td>A_73_106955</td><td>9.91</td><td>0.71</td><td>7.1</td></tr><tr><td>TAF1C</td><td>A_73_P037996</td><td>10.39</td><td>0.67</td><td>7.1</td></tr><tr><td>DNAJC15</td><td>A_73_P042026</td><td>6.70</td><td>1.04</td><td>7.0</td></tr><tr><td>CS</td><td>A_73_P107236</td><td>6.38</td><td>1.10</td><td>7.0</td></tr><tr><td>GATAD2A</td><td>A_73_P244505</td><td>10.48</td><td>0.66</td><td>6.9</td></tr><tr><td>MS4A15</td><td>A_73_P061221</td><td>10.72</td><td>0.64</td><td>6.9</td></tr><tr><td>PYCR1</td><td>A_73_113219</td><td>7.46</td><td>-1.45</td><td>-10.8</td></tr><tr><td>LEAP2</td><td>A_73_105018</td><td>3.35</td><td>-2.64</td><td>-8.8</td></tr><tr><td>NA</td><td>A_73_P190802</td><td>3.72</td><td>-2.36</td><td>-8.8</td></tr><tr><td>GARS</td><td>A_73_P034521</td><td>7.97</td><td>-1.07</td><td>-8.5</td></tr><tr><td>CEBPG</td><td>A_73_P035406</td><td>6.85</td><td>-1.24</td><td>-8.5</td></tr><tr><td>AARS</td><td>A_73_108830</td><td>9.67</td><td>-0.85</td><td>-8.3</td></tr><tr><td>C14H8ORF22</td><td>A_73_P031046</td><td>5.41</td><td>-1.46</td><td>-7.9</td></tr><tr><td>ATF4</td><td>A_73_P033631</td><td>11.47</td><td>-0.65</td><td>-7.5</td></tr><tr><td>MYH7</td><td>A_73_P344101</td><td>10.36</td><td>-0.71</td><td>-7.3</td></tr><tr><td>SARS</td><td>A_73_P042826</td><td>8.41</td><td>-0.86</td><td>-7.3</td></tr><tr><td>SLC1A5</td><td>A_73_111842</td><td>5.47</td><td>-1.31</td><td>-7.1</td></tr><tr><td>PSPH</td><td>A_73_P268996</td><td>4.45</td><td>-1.58</td><td>-7.0</td></tr><tr><td>CES1</td><td>A_73_107919</td><td>4.43</td><td>-1.53</td><td>-6.8</td></tr><tr><td>SLC1A4</td><td>A_73_103367</td><td>5.82</td><td>-1.16</td><td>-6.8</td></tr><tr><td>TNNC1</td><td>A_73_P475968</td><td>12.06</td><td>-0.56</td><td>-6.8</td></tr><tr><td>MRPL37</td><td>A_73_108914</td><td>6.86</td><td>-0.91</td><td>-6.3</td></tr><tr><td>C9H6ORF225</td><td>A_73_P081161</td><td>7.06</td><td>-0.88</td><td>-6.2</td></tr><tr><td>GPT2</td><td>A_73_108464</td><td>8.67</td><td>-0.71</td><td>-6.2</td></tr><tr><td>PRKG1</td><td>A_73_102930</td><td>5.96</td><td>-1.03</td><td>-6.2</td></tr><tr><td>IARS</td><td>A_73_105874</td><td>7.60</td><td>-0.81</td><td>-6.1</td></tr></tbody></table>

We have reported average abundance in BA and CH (abundance), DE (BA minus CH) and abundance multiplied by DE (PIF).

## Mitochondrial physiology

Citrate synthase (CS) is a nuclear-encoded mitochondrial gene, indicative of aerobic capacity and mitochondrial content. The presence of CS mRNA as upregulated in BA at 260 dpc (Table 3) is contrary to our expectation by fibre composition (>fast glycolytic) and also contrary to previously documented biochemical assay data for the oxidative enzymes citrate synthase and isocitrate dehydrogenase. This conflicting observation prompted a deeper exploration of the mitochondrial proteome to determine whether this was a general trend or specific to this particular gene.

In order to determine whether there are substantial differences in muscle mitochondrial physiology between BA and CH animals by birth, we quantified the changes in the mitochondrial proteome at the mRNA level. A mitoproteome list including both mitochondrial and nuclear encoded transcripts contained 1289 proteins. We identified 514 matches including duplicates in our data (20,354 probes). The underestimate reflects mismatches due to gene aliases and/or gene names versus protein names. We attempted to get a global perspective on mitochondrial expression between the two breeds (Fig. 2). Comparing all the data with respect to the 0 baseline,

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** mRNA encoding the mitoproteome. Overlaying the mRNA encoding the mitoproteome onto the 260 days MA plot illustrates a significant skew towards the CH, consistent with the expectation of lower oxidative capacity in the fast, glycolytic musculature of the BA.</figcaption></figure>

169 mRNA were higher in BA but more than double that were higher in CH (345), a highly significant deviation from the null expectation of 50:50 equilibrium (P=3.21 e⁻¹⁵; binomial test statistic).

The transcript abundance of the mRNA encoding the mitoproteome was therefore skewed towards the slow, oxidative CH. While it is not clear exactly how patterns of mitochondrial gene expression reflects mitochondrial content (as opposed to mitochondrial activity), this finding is certainly consistent with a reduction in mitochondrial content in the MSTN mutant fast glycolytic BA muscle tissue. This result meets our expectation given the CS biochemical data (an established proxy of mitochondrial content), and as indicated by the fibre histology data being skewed away from the oxidative mitochondrial rich fibres. Nevertheless, a substantial number of mRNA-encoding mitochondrial proteins are more highly expressed in the BA despite the lower mitochondrial content. These include *GMPPB*, *ACYP2*, *TFAM*, *ACSF2*, *UXS1* and *MRS2*. They imply some changes to isolated mitochondrial performance/metabolism in the BA versus CH animals that cannot be simply explained by reduced mitochondrial content. Most prominently upregulated in the CH include *GSTK1*, *ACACB*, *PXMP2* and *ADCK3*.

## DC

The two complementary DC scores (RIF scores) were computed in a manner taking into account correlations computed across the four developmental time points for 870 regulatory molecules identified in the data (Table S5). They were plotted and manually explored to identify regulatory molecules at the extremes. The RIF data become more sparse at the very extremes, indicating that most regulatory molecules are not behaving differently between the BA and CH but a smaller subset are substantially rewired. We focussed on those molecules that were deemed extreme by both versions of RIF.

MSTN, which harbours a causal double-muscling mutation in BA (Bouyer et al., 2014), possessed strong negative RIF1 and RIF2 values, so when graphed was positioned at one of the extremes in the bottom left quadrant in Fig. 3. Given that there is a MSTN mutation in BA, this observation can be used to ground-truth the utility of our approach. The regulatory molecule for which behaviour was most like MSTN was *NRIP2* (Fig. 3). Other molecules prominently displayed by the RIF plot include *TAF12*, *MBD1*, *NMI*, *PROX1* and *OVGP1*. The muscle regulatory factors (*MYOD1*, *MYOG*, *MYF5* and
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** RIF1 and RIF2 plotted for all transcriptional regulators. Molecules at the extremes of the distribution are highly differentially connected in developing skeletal muscle when comparing the two breeds. They may or may not be differentially expressed. Regulators of mitochondrial content and *MSTN* have been highlighted, along with a variety of extreme performers.</figcaption></figure>

*MEF2C*) showed only modest evidence for rewiring, implying the MSTN mutation's impact acts relatively independent of them. Finally, we overlaid known master regulators (including both drivers and suppressors) of mitochondrial content harvested from a number of literature searches. Of the 29 for which we have RIF data, the most prominently rewired were *TFB1M* (Falkenberg et al., 2002) and *DGUOK* (Mousson de Camaret et al., 2007). Interestingly, 23 of the 29 mitochondrial content regulators were awarded positive (rather than negative) RIF2 scores, a significant enrichment (binomial P=0.0013) for which we have no explanation.

## Hierarchical clustering

Finally, we represented a number of outlier genes, both highly developmentally regulated, highly differentially expressed and highly differential connected as a heat map hierarchically clustered on rows (Fig. 4). The clusters highlight patterns of tight co-expression across the eight treatments. Notable highly co-expressed genes include *TNNC1*, *TPM3* and *MYH7*. These three genes encode contractile proteins all of which are highly expressed in slow twitch fibres. The clustering of mRNA that encode clearly functionally related proteins is an indication of data robustness.

We also wished to explore whether the individual normalized biological replicates are representative of the group means we have been using. Reassuringly, a small panel of five genes, the expression of which was identified as highly developmentally regulated in both breeds (*MYL4*, *MYH3*, *MYH1*, *CKMT2* and *C13H20ORF194*) were collectively able to cluster the samples into the developmental treatment group of origin i.e. the individual samples for d110 clustering with other samples from d110 and so on through each developmental group (Fig. S2). Moreover, the top 20 most upregulated and top 20 most downregulated genes in BA versus CH animals at d260 were able to cluster the animals by breed (Fig. 5). Both matrices were also clustered on rows to identify genes co-expressed across treatments of interest. In addition to the

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** A hierarchically clustered heat map of normalized mean expression values with treatments as columns and genes as rows. The genes selected for hierarchical clustering have been emphasized in this manuscript because of high developmental regulation (*MYL4*, *MYH3*, *C13H20ORF194*, *CKMT2* and *MYH1*), high DE between the breeds (*TPM3*, *MYH7*, *TNNC1*, *CSRP3*, *MSTN*) or high DC between the breeds (*MSTN*, *NRIP2* and *TFB1M*). Red represents high expression and green represents low expression. A number of tight co-expression relationships are apparent, including that of *TNNC1*, *TPM3* and *MYH7*, all of which encode contractile subunits predominantly expressed in slow twitch fibres.</figcaption></figure>

clustering of mRNA encoding slow twitch proteins as described above, we also determined the very tight co-expression of mRNA (*AARS*, *SARS* and *IARS* with *GARS* only slightly removed) encoding different tRNA synthetases. The biologically meaningful nature of the clusters on both columns and rows is an indication that the individual samples perform in a similar manner to the treatment means and therefore can be considered reliable in terms of signal to noise.

## DISCUSSION

The BA is a feed efficient, hypermuscular breed whose muscularity is not present at birth but rather appears during early postnatal development. However, like most mammals the muscle potential is set by birth because of fibre number determination (Bonnet et al., 2010). This means hypermuscularity is achieved by hypertrophy of a set of existing fibres, not by the addition of new fibres (B.P., unpublished data). Compared to the CH, it is possible that the BA is still undergoing hyperplasia at 180 dpc. If true, this extension of hyperplasia in BA may increase the total number of fibres available for postnatal hypertrophy, with implications for enhancing muscle mass. We have characterized the molecular differences between the two breeds during prenatal development – when the muscle program is being set up – using a set of genome-wide transcriptome analyses. Strictly speaking, our approach on ranked genome-wide lists only allowed us to interpret biological processes on the 16 foetuses we studied. Whether these data can be extrapolated to a larger population of BA and CH foetuses remains to be evaluated, although many of our molecular findings do fit with our biological expectations of the two breeds.
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** The hierarchically-clustered heatmap for the individual normalized mean expression values for the 40 most differentially expressed genes between BA and CH at d260. This panel of 40 DE genes is able to correctly cluster the individual samples into the appropriate breed treatment group at this time point. This suggests that the treatment means we used as the basis for much of the analysis are reasonable representations of the full data set.</figcaption></figure>

Before describing the contrasts between the breeds, we will briefly comment on developmental programs they have in common, and which have previously been observed in other species of cattle.

## Developmental profiles common to both breeds

In both breeds the most developmentally regulated genes include embryonic muscle subunit isoforms (e.g. *MYL4*, *MYH3*), the expression of which drops dramatically prior to birth when they are replaced by adult isoforms (*MYH1*); metabolic enzymes (e.g. the mitochondrially localized *CKMT2* and *AGXT2L1*) that prepare the foetus for functional independence postpartum, as seen at the protein level (Chaze et al., 2008); and a gene of unknown function (*C13H20ORF194*). Overall, these muscle subunit and metabolic changes are consistent with previous results on MyHC abundances (Picard et al., 2002) and also on *m. Semitendinosus* proteomic analyses (Chaze et al., 2008; Chaze et al., 2009) observed during the foetal stages in CH. Similarly, these basic changes were also observed in an independent mRNA developmental analysis of

Wagyu by Hereford breed cross *Longissimus* muscle (Hudson et al., 2013). The dramatic increase in *C13H20ORF194* has also been described in the prenatal developing muscle of Wagyu by Hereford animals (Hudson et al., 2013). To better understand the significance of *C13H20ORF194*, we performed some basic bioinformatic analysis using The Human Protein Atlas (http://www.proteinatlas.org) and the TIGER tissue-specific expression database (Liu et al., 2008). It appears that *C13H20ORF194* encodes a large protein that is highly conserved across the vertebrates. There is some evidence for cytoplasmic localization and it possesses a strong coiled-coil domain at the C-terminus. According to the Mammalian Phenotype Browser, the encoded protein plays a role in metabolism (http://www.informatics.jax.org/marker/MGI:1923029) and there is evidence for expression in skeletal muscle based on both NGS and SAGE (http://www.genecards.org/cgi-bin/carddisp.pl?gene=C20orf194). In light of these prior findings and our own expression data, we hypothesize that there would be value in functionally characterising this protein as it looks to play a role in developing production animal muscle.

## BA versus CH

The two major structural and functional changes associated postnatally with the MSTN mutation in addition to an increase in muscle mass are (1) a shift to a faster muscle isoforms and (2) a decrease in mitochondrial content. The latter is because the faster muscle fibres are predominantly geared towards glycolytic rather than oxidative metabolism. In relaxed or gently exercising muscle the 'glycolytic' fibres would still be expected to convert energy aerobically, given aerobic metabolism is the only sustainable source of ATP, and can provide an adequate ATP supply under conditions of low demand. We fibre typed the animals using immunohistochemistry, finding substantially less oxidative type I fibres and substantially more type II fibres in the BA. These basic physiology results are similar to data obtained in double-muscled cattle compared to non-double-muscled cattle and to cattle selected on their muscle growth potential (Bouley et al., 2005; Deveaux et al., 2001; Gagnière et al., 1997; Picard et al., 2006). Both these functional changes are evident at the mRNA level. First, we find a decrease in expression of a number of slow twitch subunits in BA cattle in one or both of the later developmental time points (*TPM3*, *CSRP3*, *MYH7* and *TNNC1*), some of which (e.g. *TNNC1* and *MYH7*) are highly positively co-expressed. This suggests not only that the BA possesses a shift towards fast, glycolytic muscle but that this is already apparent at the 210 dpc time point prior to birth. Interestingly, all of these subunits are also relatively reduced at similar time points in another MSTN mutant, the Piedmontese (Hudson et al., 2009a). One subtle difference is the expression of *MYL2* which surprisingly is not relatively downregulated in BA. This may reflect the different genetic backgrounds against which the MSTN mutation is expressed.

Second, we observe that the majority of mRNA encoding the mitoproteome for which we have data (345 of 514 proteins or 68% of the total) are downregulated in BA. This is consistent with a reduction in mitochondrial content in BA muscle and aligns with our expectation based on both CS biochemical activity and fibre histology. The consistent upregulation of interferon signalling apparatus in the BA compared to the CH at multiple time points is an unresolved mystery at this stage. However, one possible explanation relates to the known ability of interferon to promote muscle healing through formation of new muscle fibres (Cheng et al., 2008). Alternatively, the result may reflect the immune systems utility in a monitoring and signalling capacity as previously been observed for toll receptors and hematopoietic homeostasis (Cannova et al.,
2015). A third more prosaic possibility is that these data simply reflect breed-specific MHC haplotypes being manifested in a muscle context.

Interestingly, at d260, we found coordinate downregulation of four transcripts (AARS, SARS, IARS and GARS) encoding different tRNA synthetases, perhaps indicating reduced protein synthesis in the BA. At 260 dpc MSTN is the fifth most upregulated gene in BA. Of the four genes more upregulated than *MSTN* (*XIST*, *JSP.1*, *UPK2* and *CTDSP2*), *XIST* (X chromosome inactivation) and *CTDSP2* (androgen receptor regulation) may be a reflection of gender imbalance in the experimental design. A subset of the extreme DE genes show high individual variability, and their potential role in the process under investigation can be considered less likely.

Next, we were interested in the regulatory machinery that may be
driving the changes in fibre composition and mitochondrial content.
We made use of our DC analysis (RIF), which was developed on a
cattle muscle developmental system similar to this one (Hudson
et al., 2009a), but subsequently applied in a diverse range of
circumstances to unravel the molecular basis of cancer (Reverter
et al., 2010), cell determination pathways (Reverter et al., 2010), the
impact of hormone growth promotants (De Jager et al., 2011) and
several human neurodegenerative diseases (Rhinn et al., 2013;
Rhinn et al., 2012). The basis of the analysis is to ask the question
'which transcriptional regulator changes its network topology the
most between the two breeds?' Regulators with markedly different
network topologies may or may not be differentially expressed
themselves, so the approach can be considered complementary to
conventional DE-based analyses.

We find *MSTN* to be extremely differentially connected between the breeds as expected. Further, we find that another regulator (*NRIP2*) has a very similar change in behaviour to *MSTN*, implying a hitherto unrecognized role in driving the muscle mass, fibre composition and mitochondrial content phenotype change in *MSTN* mutant animals. *NRIP2* is a poorly characterized molecule, the protein of which is thought to interact with regulators of metabolism such as NR1F2, RARA and THRB. We hypothesize that *NRIP2* deserves more scrutiny in a production muscle context. In this experiment, the anticipated change in *MSTN* activity is also detectable by DE as the mRNA is prominently differentially expressed between the breeds at all time points. However, the direction of change (higher expression in the more muscular breed) is opposite to that predicted by its functional outcome (as it is a repressor of muscle mass).

MSTN higher expression in the BA may represent an attempt to restore tissue mass homeostasis, while its DC likely reflects the loss of growth repression by the dysfunctional protein. We wish to point out that these BA foetuses have not been genotyped for the mutation. Typically, 80% of BA carry the mutated transcript, and it works in correspondence with the wild-type transcript to drive muscle mass (I.C.-M., unpublished data). In a similar experiment comparing *MSTN* mutant Piedmontese and Wagyu (Hudson et al., 2009a), we made slightly different observations: while we were able to identify the differential activity/behaviour of *MSTN* through DC, this was in the absence of any DE. Finally, in the specific context of reduced mitochondrial content in BA versus CH we found evidence for substantial rewiring of the master regulator of mitochondrial transcription and biogenesis, *TFB1M*.

CONCLUSIONS

We have used global gene expression analyses to better understand the
molecular basis of high muscling in the *MSTN* mutant BA. A number
of regulatory genes have been identified as behaving differently

between the two breeds, some through conventional DE and some
through DC. Interestingly, the influence of the causal *MSTN* mutation
was highlighted through both DE and DC. A number of molecular
observations made here in mutant BA are concordant with a previously
published comparison between the *Longissimus* muscle of *MSTN*
mutant Piedmontese by Hereford crosses versus Wagyu by Hereford
crosses (Hudson et al., 2009a). This includes the gene expression
evidence for a fibre composition shift towards the low mitochondrial
content glycolytic fibres. Overall, these comparisons imply some
common muscle biology across these different hyper-muscular breeds
that transcend the particular muscle selected for analysis.

MATERIALS AND METHODS

Animal resources

CH and BA fetuses were generated by artificial insemination of CH and BA primiparous cows with CH and BA bull sperm, respectively, at the INRA experimental farm, located at Theix near Clermont Ferrand, France. Fetuses were collected from cows slaughtered at the INRA experimental abattoir at 110 (n=3/breed), 180 (n=3/breed), 210 (n=3/breed) and 260 (n=3/breed) dpc. The collection of fetuses at these stages was completed to target the proliferation of foetal myoblasts (110 dpc), the end of proliferation (180 dpc), the contractile and metabolic differentiation of fibres, and the perinatal period, respectively. The *m. Semitendinosus* muscle was sampled for histochemistry, biochemistry and transcriptomics according to (Cassar-Malek et al., 2010).

All experimental procedures were approved by the Animal Care
Committee in accordance with the Use of Vertebrates for Scientific
Purposes Act 1985. The foetuses were collected respecting broadly
recognized ethical standards, although a formal ethical evaluation was not
mandatory at the time of this experiment.

Histology

Serial sections (10 µm) were cut perpendicular to the *m. Semitendinosus* muscle fibres on a cryostat at -25°C. Contractile type was determined using anti-MyHC antibodies. Briefly, serial transverse sections were incubated with specific monoclonal antibodies (Picard et al., 1998) raised either against the slow (5B9) or the fast (F113 15F4) or a combination of slow and fast IIx (F5 8H2) or of the foetal *F* (4C10) MyHC categories. Immune complexes were revealed using a fluorescein isothiocyanate (FITC)-conjugated secondary antibody (Interchim, Montiluçon, France). Immunofluorescent staining was analysed in two zones under a microscope (Labophot-2, Nikon, Tokyo, Japan) with a fluorescent light (excitation filter 450 to 490 nm and stop filter 520 nm). The determination of fibre types and their frequencies expressed as a percentage was based on 200 fibres per cross section. Muscle areas from exactly the same position were analysed.

RNA extraction and microarray hybridization

Total RNA was extracted from snap-frozen *m. Semitendinosus* tissue stored at -80°C using Trizol following the manufacturer's instructions. A Purelink RNA Mini Kit (Invitrogen) was used followed by a Dnase treatment (Qiagen). Then, 50 ng total RNA was amplified and fluorescently labelled to produce Cy3 cRNA using the Low Input Quick Amp Labeling One Color Kit (Agilent, Santa Clara, CA, USA). After purification (RNeasy Mini Kit, Qiagen), 1650 ng labelled cRNA was hybridized onto G2519F 023647 Bovine v2 GE 4x44 K microarray (Agilent) with hybridization buffer (Gene Expression Hybridization Kit, Agilent) according to the manufacturer's instructions. After 17 h of hybridization, slides were washed with the Gene Expression Wash Pack (Agilent). Microarrays were scanned with an Agilent G2505 scanner. Based on considerations of quality control we generated microarray data for the following subset of samples: 110d (BA *n*=3, CH *n*=2), 180 days (BA *n*=2, CH=3), 210 days (BA *n*=3, CH *n*=3), 260 days (BA *n*=2, CH *n*=3).

Data normalization

Within and between array normalization was performed as previously
described (Hocquette et al., 2012) Data were extracted and normalized with
Feature Extraction 10.7 software (Agilent) using background and multiplicative detrend (no background subtraction). The mean intensity between the technical replicates was then calculated for each probe. Normalized log2 mean intensity (Cy3) were used. Data were submitted to GEO and published under accession number GSE60844.

## Calculation of DE

DE was calculated using the LIMMA method (Smyth et al., 2005) with a Benjamini-Hochberg multiple testing correction. Very few significant DE genes were identified using this method. Consequently, we proceeded with a less stringent approach, relying on ranked genome-wide lists to biologically interpret the data. These approaches are less dependent on the significance of any one gene, but rather explore patterns of functional enrichment detected in large lists. All the downstream analyses can then focus on either the extremes of ranked lists using nominal P-values (top 1% or 5% of the genes in the list) or exploiting the entire ranked list for functional enrichment analysis using GOrilla (Eden et al., 2009) in which case no nominal threshold needs to be set. We focussed on the 'Process' level ontology for reporting in this manuscript.

To compute the MA plots between the breeds the average abundance and the DE (BA-CH) was computed for each gene at each time point. To ensure a more manageable data set for plotting we filtered to include all those genes minimally DE (a log2 difference of 0.5 equating to a fold change of 1.41). This yielded 2911, 1599, 966 and 1031 genes for 110, 180, 210 and 260 days, respectively, reflecting a reduction in breed difference as the muscle developmental program unfolds. We also computed a cumulative 'overall' DE by summing the absolute DE at each time point and then ranking on the sum. We generated a ranked list of developmentally regulated genes by ranking on standard deviation across the four time points for the BA.

## Comparison to MSTN mutant Piedmontese versus wild-type Wagyu

We attempted to find commonalities between the BA versus CH muscle structural DE with that made in an independent MSTN mutant system. We made use of the published data comparing MSTN mutant Piedmontese by Hereford crosses versus MSTN wild-type Wagyu by Hereford crosses (Hudson et al., 2009a), examining DE (Piedmontese minus Wagyu) of the slow muscle subunits at the 280 days peri-natal period.

## Quantifying the mitoproteome in BA versus CH

We filtered the 260 days data with a comprehensive mitoproteome list derived from http://mitominer.mrc-mbu.cam.ac.uk/release-3.1/begin.do. This list includes mitochondrially encoded proteins as well as the much more numerous nuclear-encoded proteins. We generated a gene list with no duplicates and based our gene list on the human data (as opposed to the mouse data) in the mitominer database.

## Calculation of PIF

Next we computed PIFs for each gene at each time point by multiplying average abundance by breed DE as previously described (Hudson et al., 2009a). We ranked the lists and explored for functional enrichment by hypergeometric statistics using GOrilla analysis. We used the two list import option, with a background list of 19,100 genes and target lists derived from the extreme 191 (nominal 1%) up- and downregulated.

## Calculation of DC

We calculated RIF1 and RIF2 DC scores using the methods outlined in Hudson et al. (2009a) and Reverter et al. (2010). In brief, to identify transcription factors and other regulatory molecules that have been rewired in the BA we applied our recently developed RIF analysis. We identified the regulatory molecules using Genomatix software as detailed in Hudson et al. (2009b). As we previously described (Hudson et al., 2015), the RIF procedure exploits global patterns of differential co-expression (or 'differential wiring') to infer those regulatory molecules for which behaviour is systematically different in a contrast of interest, in this case BA versus CH *m. Semitendinosus* muscle. RIF metrics for the *r*-th regulator

were computed using the following formulae:

$$RIF1_r = \frac{1}{n_{DE}} \sum_{j=1}^{j=n_{DE}} x_j \times d_j \times DC_{rj}^2$$

and

$$RIF2_r = \frac{1}{n_{DE}} \sum_{j=1}^{j=n_{DE}} [(x_j^{BA} \times r_{rj}^{BA})^2 - (x_j^{CH} \times r_{rj}^{CH})^2]$$

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 BA versus CH 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_{rj}^{BA}$ and $r_{rj}^{CH}$, the correlation co-expression between the *r*-th regulator and the *j*-th DE gene in the BA and CH samples, respectively; and $x_j^{BA}$ and $x_j^{CH}$ represented the average expression of the *j*-th DE or TS gene in the BA and CH samples, respectively.

## Hierarchical clustering of both developmentally and breed regulated genes

We created a matrix of normalized mean expression values with as many columns as treatments (8) and as many rows as genes (12). The 12 genes were selected for deeper analysis on the basis of extreme developmental regulation (*MYH1*, *MYL4*, *MYH3*, *CKMT2*, *C13H20ORF194*) or extreme breed DE (*MSTN*, *CSRP3*, *TNNC1*, *MYH7* and *TPM3*) or extreme DC (*MSTN*, *NRIP2*, *TFB1M*). The matrix was imported into Permut Matrix software (Caraux and Pinloche, 2005), normalized on rows, and then hierarchically clustered on rows to highlight genes co-expressed across the eight treatments. Clustered genes were explored for functional relatedness of the encoded proteins.

In a complementary analysis, we also explored the numerical behaviour of the individual biological replicates to see if they were representative of the treatment means. Two matrices of individual normalized mean expression were generated to be imported into Permut Matrix for hierarchical cluster analysis. The first matrix was to test the performance of those genes previously highlighted as developmentally regulated. Are those individual normalized expression values correctly able to discriminate the samples into developmental group of origin? This matrix had 21 columns (all individual values across the eight treatments) and five rows (representing the developmentally regulated genes *MYL4*, *MYH3*, *MYH1*, *CKMT2* and *C13H20ORF194*).

The second was to test the performance of those genes previously highlighted as having substantial DE between the breeds at d260 when the muscle has reached its mature form. Are those individual normalized expression values correctly able to discriminate the samples into the breed of origin? This matrix had five columns (individual values across d260 samples only) and 40 rows (the 20 most up- and downregulated genes at d260 in BA versus CH as determined by PIF). In both cases the matrices were imported into Permut Matrix software, normalized on rows, then hierarchically clustered on both rows and columns.

## Acknowledgements

We thank Anne de La Foye, who planned the hybridization experiment, completed the normalization of microarray data and performed LIMMA differential analysis. We also thank G. Gentés for assistance with RNA extraction and microarray hybridization, C. Barboiron for assistance with foetus dissection and immunohistochemistry, and D. Chadeyron for the measurement of citrate synthase activity.

## Competing interests

The authors declare no competing or financial interests.

## Author contributions

Conceptualization: I.C.-M., C.B., B.P., A.R., N.J.H.; Methodology: I.C.-M, C.B., B.P., A.R., N.J.H.; Software: B.P.; Formal analysis: C.B., A.R., N.J.H.; Investigation: B.P., N.J.H.; Resources: I.C.-M; Writing - original draft: N.J.H.; Writing - review & editing: I.C.-M, C.B., B.P., A.R.; Supervision: I.C.-M; Project administration: I.C.-M; Funding acquisition: I.C.-M
Funding

This work was supported by the Commonwealth Science and Industrial Research Organisation (JM Rendel Visiting Fellowship to N.J.H.). The funder had no role in study design nor any other aspect of the research.

Data availability

The gene expression data have been deposited into the Gene Expression Omnibus database under GEO GSE60844.

Supplementary information

Supplementary information available online at
http://bio.biologists.org/lookup/doi/10.1242/bio.024950.supplemental

References

Andersson, L. and Georges, M. (2004). Domestic-animal genomics: deciphering the genetics of complex traits. *Nat. Rev. Genet.* **5**, 202-212.

Bonnet, M., Cassar-Malek, I., Chilliard, Y. and Picard, B. (2010). Ontogenesis of muscle and adipose tissues and their interactions in ruminants and other species. *Animal* **4**, 1093-1109.

Bottje, W. G. and Carstens, G. E. (2009). Association of mitochondrial function and feed efficiency in poultry and livestock species. *J. Anim. Sci.* **87**, E48-E63.

Bouley, J., Meunier, B., Chambon, C., De Smet, S., Hocquette, J. F. and Picard, B. (2005). Proteomic analysis of bovine skeletal muscle hypertrophy. *Proteomics* **5**, 490-500.

Bouyer, C., Forestier, L., Renand, G. and Oulmouden, A. (2014). Deep intronic mutation and pseudo exon activation as a novel muscular hypertrophy modifier in cattle. *PLoS ONE* **9**, e97399.

Cannova, J., Breslin, S. J. P. and Zhang, J. (2015). Toll-like receptor signaling in hematopoietic homeostasis and the pathogenesis of hematologic diseases. *Front Med.* **9**, 288-303.

Caraux, G. and Pinloche, S. (2005). PermutMatrix: a graphical environment to arrange gene expression profiles in optimal linear order. *Bioinformatics* **21**, 1280-1281.

Cassar-Malek, I., Picard, B., Jurie, C., Listrat, A., Guillomot, M., Chavatte-Palmer, P. and Heyman, Y. (2010). Myogenesis is delayed in bovine fetal clones. *Cell. Reprogram.* **12**, 191-201.

Chaze, T., Meunier, B., Chambon, C., Jurie, C. and Picard, B. (2008). In vivo proteome dynamics during early bovine myogenesis. *Proteomics* **8**, 4236-4248.

Chaze, T., Meunier, B., Chambon, C., Jurie, C. and Picard, B. (2009). Proteome dynamics during contractile and metabolic differentiation of bovine foetal muscle. *Animal* **3**, 980-1000.

Cheng, M., Nguyen, M.-H., Fantuzzi, G. and Koh, T. J. (2008). Endogenous interferon-gamma is required for efficient skeletal muscle regeneration. *Am. J. Physiol. Cell Physiol.* **294**, C1183-C1191.

De Jager, N., Hudson, N. J., Reverter, A., Wang, Y. H., Nagaraj, S. H., Cafe, L. M., Greenwood, P. L., Barnard, R. T., Kongsuwan, K. P. and Dalrymple, B. P. (2011). Chronic exposure to anabolic steroids induces the muscle expression of oxytocin and a more than fiftyfold increase in circulating oxytocin in cattle. *Physiol. Genomics.* **43**, 467-478.

Deveaux, V., Cassar-Malek, I. and Picard, B. (2001). Comparison of contractile characteristics of muscle from Holstein and double-muscled Belgian Blue foetuses. *Comp. Biochem. Physiol. A Mol. Integr. Physiol.* **131**, 21-29.

Eden, E., Navon, R., Steinfeld, I., Lipson, D. and Yakhini, Z. (2009). GOrilla: a tool for discovery and visualization of enriched GO terms in ranked gene lists. *BMC Bioinformatics* **10**, 48.

Falkenberg, M., Gaspari, M., Rantanen, A., Trifunovic, A., Larsson, N.-G. and Gustafsson, C. M. (2002). Mitochondrial transcription factors B1 and B2 activate transcription of human mtDNA. *Nat. Genet.* **31**, 289-294.

Fu, L., Xu, Y., Hou, Y., Qi, X., Zhou, L., Liu, H., Luan, Y., Jing, L., Miao, Y., Zhao, S. et al. (2017). Proteomic analysis indicates that mitochondrial energy metabolism in skeletal muscle tissue is negatively correlated with feed efficiency in pigs. *Sci. Rep.* **7**, 45291.

Gagnière, H., Picard, B., Jurie, C. and Geay, Y. (1997). Comparative study of metabolic differentiation of foetal muscle in normal and double-muscled cattle. *Meat Sci.* **45**, 145-152.

Hanset, R., Michaux, C. and Stasse, A. (1987). Relationships between growth rate, carcass composition, feed intake, feed conversion ratio and income in four biological types of cattle. *Genet. Sel. Evol.* **19**, 225-248.

Herd, R. M. and Arthur, P. F. (2009). Physiological basis for residual feed intake. *J. Anim. Sci.* **87**, E64-E71.

Hocquette, J. F., Bernard-Capel, C., Vidal, V., Jesson, B., Leveziel, H., Renand, G. and Cassar-Malek, I. (2012). The GENOTEND chip: a new tool to analyse gene expression in muscles of beef cattle for beef quality prediction. *BMC Vet. Res.* **8**, 135.

Hudson, N. J. (2009). Symmorphosis and livestock bioenergetics: production animal muscle has low mitochondrial volume fractions. *J. Anim. Physiol. Anim. Nutr. (Berl.)* **93**, 1-6.

Hudson, N. J., Reverter, A. and Dalrymple, B. P. (2009a). A differential wiring analysis of expression data correctly identifies the gene containing the causal mutation. *PLoS Comput. Biol.* **5**, e1000382.

Hudson, N. J., Reverter, A., Wang, Y. H., Greenwood, P. L. and Dalrymple, B. P. (2009b). Inferring the transcriptional landscape of bovine skeletal muscle by integrating co-expression networks. *PLoS ONE* **4**, e7249.

Hudson, N. J., Dalrymple, B. P. and Reverter, A. (2012). Beyond differential expression: the quest for causal mutations and effector molecules. *BMC Genomics* **13**, 356.

Hudson, N. J., Lyons, R. E., Reverter, A., Greenwood, P. L. and Dalrymple, B. P. (2013). Inferring the in vivo cellular program of developing bovine skeletal muscle from expression data. *Gene. Expr. Patterns* **13**, 109-125.

Hudson, N. J., Reverter, A., Greenwood, P. L., Guo, B., Cafe, L. M. and Dalrymple, B. P. (2015). Longitudinal muscle gene expression patterns associated with differential intramuscular fat in cattle. *Animal* **9**, 650-659.

Jackson, S. P., Green, R. D. and Miller, M. F. (1997a). Phenotypic characterization of rambouillet sheep expressing the callipyge gene: I. Inheritance of the condition and production characteristics. *J. Anim. Sci.* **75**, 14-18.

Jackson, S. P., Miller, M. F. and Green, R. D. (1997b). Phenotypic characterization of Rambouillet sheep expressing the callipyge gene: II. Carcass characteristics and retail yield. *J. Anim. Sci.* **75**, 125-132.

Jackson, S. P., Miller, M. F. and Green, R. D. (1997c). Phenotypic characterization of rambouillet sheep expression the callipyge gene: III. Muscle weights and muscle weight distribution. *J. Anim. Sci.* **75**, 133-138.

Kiessling, K.-H. (1977). Muscle structure and function in the goose, quail, pheasant, guinea hen, and chicken. *Comp. Biochem. Physiol. B* **57**, 287-292.

Lefaucheur, L. and Ecolan, P. (2005). Pattern of muscle fiber formation in Large White and Meishan pigs. *Arch. Anim. Breed.* **48**, 117-122.

Lefaucheur, L., Milan, D., Ecolan, P. and Le Callennec, C. (2004). Myosin heavy chain composition of different skeletal muscles in Large White and Meishan pigs. *J. Anim. Sci.* **82**, 1931-1941.

Lefaucheur, L., Lebret, B., Ecolan, P., Louveau, I., Damon, M., Prunier, A., Billon, Y., Sellier, P. and Gilbert, H. (2011). Muscle characteristics and meat quality traits are affected by divergent selection on residual feed intake in pigs. *J. Anim. Sci.* **89**, 996-1010.

Lehnert, S. A., Reverter, A., Byrne, K. A., Wang, Y., Nattrass, G. S., Hudson, N. J. and Greenwood, P. L. (2007). Gene expression studies of developing bovine longissimus muscle from two different beef cattle breeds. *BMC Dev. Biol.* **7**, 95.

Listrat, A., Picard, B., Jailler, R., Collignon, H., Peccatte, J. M., Micol, D., Geay, Y. and Dozias, D. (2001). Grass valorisation and muscular characteristics of blonde d'Aquitaine steers. *Anim. Res.* **50**, 105-118.

Liu, X., Yu, X., Zack, D. J., Zhu, H. and Qian, J. (2008). TIGER: a database for tissue-specific gene expression and regulation. *BMC Bioinform.* **9**, 271.

Mousson De Camaret, B., Taanman, J. W., Padet, S., Chassagne, M., Mayencon, M., Clerc-Renaud, P., Mandon, G., Zabot, M. T., Lachaux, A. and Bozon, D. (2007). Kinetic properties of mutant deoxyguanosine kinase in a case of reversible hepatic mtDNA depletion. *Biochem. J.* **402**, 377-385.

Picard, B. and Cassar-Malek, I. (2009). Evidence for expression of Ilb myosin heavy chain isoform in some skeletal muscles of Blonde d'Aquitaine bulls. *Meat Sci.* **82**, 30-36.

Picard, B., Duris, M. P. and Jurie, C. (1998). Classification of bovine muscle fibres by different histochemical techniques. *Histochem. J.* **30**, 473-479.

Picard, B., Lefaucheur, L., Berri, C. and Duclos, J.-M. (2002). Muscle fibre ontogenesis in farm animal species. *Reprod. Nutr. Dev.* **42**, 415-431.

Picard, B., Jurie, C., Duris, M. P. and Renand, G. (2006). Consequences of selection for higher growth rate on muscle fibre development in cattle. *Livestock Sci.* **102**, 107-120.

Rehfeldt, C., Henning, M. and Fiedler, I. (2008). Consequences of pig domestication for skeletal muscle growth and cellularity. *Livestock Sci.* **116**, 30-41.

Reverter, A., Hudson, N. J., Nagaraj, S. H., Pérez-Enciso, M. and Dalrymple, B. P. (2010). Regulatory impact factors: unraveling the transcriptional regulation of complex traits from expression data. *Bioinformatics* **26**, 896-904.

Rhinn, H., Qiang, L., Yamashita, T., Rhee, D., Zolin, A., Vanti, W. and Abeliovich, A. (2012). Alternative alpha-synuclein transcript usage as a convergent mechanism in Parkinson's disease pathology. *Nat. Commun.* **3**, 1084.

Rhinn, H., Fujita, R., Qiang, L., Cheng, R., Lee, J. H. and Abeliovich, A. (2013). Integrative genomics identifies APOE epsilon4 effectors in Alzheimer's disease. *Nature* **500**, 45-50.

Robinson, D. L. and Oddy, V. H. (2004). Genetic parameters for feed efficiency, fatness, muscle area and feeding behaviour of feedlot finished beef cattle. *Livestock Prod. Sci.* **90**, 255-270.

Smyth, G. K., Michaud, J. and Scott, H. S. (2005). Use of within-array replicate spots for assessing differential expression in microarray experiments. *Bioinformatics* **21**, 2067-2075.

Sudre, K., Cassar-Malek, I., Listrat, A., Ueda, Y., Leroux, C., Jurie, C., Auffray, C., Renand, G., Martin, P. and Hocquette, J.-F. (2005). Biochemical and transcriptomic analyses of two bovine skeletal muscles in Charolais bulls divergently selected for muscle growth. *Meat Sci.* **70**, 267-277.

Vincent, A., Louveau, I., Gondret, F., Tréfeu, C., Gilbert, H. and Lefaucheur, L. (2015). Divergent selection for residual feed intake affects the transcriptomic and proteomic profiles of pig skeletal muscle. *J. Anim. Sci.* **93**, 2745-2758.