# Molecular signatures of muscle growth and composition deciphered by the meta-analysis of age-related public transcriptomics data

Jeanne Bazile, Florence Jaffrezic, Patrice Dehais, Matthieu

Matthieu.Reichstadt@inrae.Fr Reichstadt, Christophe C. Kropp, Denis Laloë, Muriel Bonnet

► To cite this version:

Jeanne Bazile, Florence Jaffrezic, Patrice Dehais, Matthieu Matthieu.Reichstadt@inrae.Fr Reichstadt, Christophe C. Kropp, et al.. Molecular signatures of muscle growth and composition deciphered by the meta-analysis of age-related public transcriptomics data. Physiological Genomics, 2020, 52 (8), pp.322-332. 10.1152/physiolgenomics.00020.2020. hal-02937790

HAL Id: hal-02937790

https://hal.inrae.fr/hal-02937790

Submitted on 14 Sep 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 | *-Omic Approaches to Understanding Muscle Biology*

Molecular signatures of muscle growth and composition deciphered by the
meta-analysis of age-related public transcriptomics data

Jeanne Bazile,¹ Florence Jaffrezic,² Patrice Dehais,³,⁴ Matthieu Reichstadt,¹ Christophe Klopp,³,⁴
Denis Laloe,² and Muriel Bonnet¹

¹INRAE, UMR Herbivores, Université Clermont Auvergne, VetAgro Sup, Saint-Genès-Champanelle, France; ²INRAE, UMR1313 Génétique Animale et Biologie Intégrative, Jouy-en-Josas, France; ³Plate-forme bio-informatique Genotoul, Mathématiques et Informatique Appliquées de Toulouse, INRAE, Castanet Tolosan, France; and ⁴SIGENAE, GenPhySE, Université de Toulouse, INRAE, ENVT, Castanet Tolosan, France

Submitted 27 February 2020; accepted in final form 7 July 2020

Bazile J, Jaffrezic F, Dehais P, Reichstadt M, Klopp C, Laloe D, Bonnet M. Molecular signatures of muscle growth and composition deciphered by the meta-analysis of age-related public transcriptomics data. *Physiol Genomics* 52: 322–332, 2020. First published July 13, 2020; doi:10.1152/physiolgenomics.00020.2020.—The lean-to-fat ratio is a major issue in the beef meat industry from both carcass and meat production perspectives. This industrial perspective has motivated meat physiologists to use transcriptomics technologies to decipher mechanisms behind fat deposition within muscle during the time course of muscle growth. However, synthetic biological information from this volume of data remains to be produced to identify mechanisms found in various breeds and rearing practices. We conducted a meta-analysis on 10 transcriptomic data sets stored in public databases, from the longissimus thoracis of five different bovine breeds divergent by age. We updated gene identifiers on the last version of the bovine genome (UCD1.2), and the 715 genes common to the 10 studies were subjected to the meta-analysis. Of the 238 genes differentially expressed (DEG), we identified a transcriptional signature of the dynamic regulation of glycolytic and oxidative metabolisms that agrees with a known shift between those two pathways from the animal puberty. We proposed some master genes of the myogenesis, namely MYOG and MAPK14, as probable regulators of the glycolytic and oxidative metabolisms. We also identified overexpressed genes related to lipid metabolism (APOE, LDLR, MXRA8, and HSP90AA1) that may contribute to the expected enhanced marbling as age increases. Lastly, we proposed a transcriptional signature related to the induction (YBX1) or repression (MAPK14, YWAH, ERBB2) of the commitment of myogenic progenitors into the adipogenic lineage. The relationships between the abundance of the identified mRNA and marbling values remain to be analyzed in a marbling biomarkers discovery perspectives.

bovine; meta-analysis; muscle; public data; transcriptomic

INTRODUCTION

The management of muscle mass and/or of meat-eating
qualities requires a good understanding of the molecular driv-
ers of muscle growth and composition. Indeed, in cattle as in
other vertebrates, muscle is composed of several tissues, such
as muscle fibers as well as connective and adipose tissues (41).
Genetic (breed, sex, etc.) and nongenetic (age, nutrition, etc.)
factors strongly influence the growth of each cell type within
muscle (9). For example, aging, the most influential nongenetic

factor, decreases muscular fiber and increases intramuscular
adipose and connective tissues proportions (9, 26, 41, 48).
Consequently, these age-related muscular composition varia-
tions affect sensorial and nutritional meat-eating qualities (41).
This has motivated the age-dependent high-throughput molec-
ular characterization of bovine muscle (4, 14, 30, 33, 43, 50,
55, 57, 72). The aim of these studies was to further knowledge
important in managing an optimal trade-off between muscle
growth and composition and to enhance beef production and
consumer satisfaction. However, to date, there is little consen-
sus on the molecular signature of muscle fiber growth relative
to connective and adipose tissues (15, 18, 21, 63). Meta-
analysis (52), as well as the emerging concept of data reuse for
transcriptomic [microarray or RNA-Seq (56)] or proteomic
results (71), offers the opportunity to produce knowledge from
available data with several advantages: reduced cost, reduced
use of animals, and reliable and robust results that have been
recorded in several independent studies. The first and rare
attempts of bovine public “omics” data meta-analyses have
successively provided signatures of the lactation process in
dairy cattle (23), predictors of oocyte quality in the bovine
(35), and conserved gene expression patterns between species
for normal and pathological muscles (6). We thus hypothesized
that the meta-analysis of available muscular transcriptomic
data would allow us to identify robust molecular signatures of
muscle growth and composition, especially marbling. We use
variation of postnatal ages across studies as a proxy for
physiological variation in muscular mass, composition, and
intramuscular lipid content. We use, reuse, reannotate and
analyze 10 muscular transcriptomic data sets available from
meat-producing bovines to describe transcriptomic signatures
of muscle according to muscle growth and breed. To our
knowledge, this is the first attempt to reuse public bovine
muscular transcriptomic data to identify conserved patterns
between breeds and rearing practices and to propose potential
indicators of muscular growth and composition in bovine.

MATERIAL AND METHODS

The work presented in this article was conducted in accordance
with guidelines previously proposed for meta-analysis of microarrays
data sets (52, 70). All data analysis was performed using R statistical
software (version 3.5.2) (51).
## Data Retrieval

Data sets and their corresponding raw expression data were collected from the public repositories Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) (7) or PubMed (https://pubmed.ncbi.nlm.nih.gov) of the National Center for Biotechnology Information, as well as from the Web of Science (WoS; https://www.webofknowledge.com/). These publications (43, 50) and data sets (from GSE5659, GSE48136, GSE21782) were identified using keywords related to muscle (muscle, myo*, satellite cells; where * allows one to broaden the search to every word starting with "myo" such as myofiber, myogenesis, etc.), growth (growth, hypertrophy, atrophy, hypermuscle, hypermuscularity, double-muscle, marbling), bovine species (bovine, calf, steers, heifers, cows, beefs, bulls), methods (transcript*, microarray, chip, gene array), and exclusion criteria (fetal, disease, cardio*, heart). Of the hundreds of data sets retrieved, we selected 25 transcriptomic data sets (18 from GEO, 7 from the WOS, see the PRISMA flowchart in Fig. 1), obtained in the longissimus thoracis muscle. These 25 transcriptomic data sets reported comparisons between bovines differing within the study only by the age (similar breed and rearing practice and same sex). We removed duplicate data from the GSM132170 and GSM132194 (38) (a GSM stands for GEO samples and corresponds to 1 data set, a GSE

(GEO Series Experiment) may comprise several GSM) that have been reanalyzed in the most recent GSE5659. We then excluded a data set that has been obtained by RT-PCR, as there would be too few genes to include in the meta-analysis (37), an RNA-Seq data set (29), and a microarray data set due to the lack of gene identifiers (GSE60844). Six data sets comparing a fetal and a postnatal age were excluded because cellular and biological pathways involved in pre- and postnatal muscle growth are very different (47). It was not possible to update gene identifiers of three data sets on the current version of the bovine genome due to the lack of sequences for the microarrays. We excluded these three data sets. Thus, 10 data sets were finally subjected to data preprocessing, gene identifier update, and meta-analysis as described thereafter (Fig. 1). These 10 data sets were stored in GEO as three GSE, namely GSE21782 (1 data set), GSE48136 (3 data sets), and GSE5659 (6 data sets) and reported mRNA abundance assayed by microarrays in muscles from five bovine pure- or cross-breeds aged from 31 to 732 days after birth (Table 1). The GSE5659 reported data from 147 samples differing by the breed, the nutrition or the age. Thus, six subdata sets were created from the GSE5659 (with 2-4 technical replicates per animal) to perform the comparison between two ages. These subdata sets were named E2, E3, E4, E6, E8, and E9 and were produced from a noncommercial microarray (39). We have created one subdata set from GSE21782, namely the study E5, with data from three muscles produced from the Affymetrix Bovine Genome Array (50). Three subdata sets were created from GSE48136, named E1, E7, and E10 with data from 13 or 14 muscles.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Prisma flowchart of the data selection. In total, 10 studies comparing muscular transcripts abundance between 2 ages of bovine were included in the meta-analysis. The characteristics of the included studies are listed in Table 1. n, Number of data sets/publications; GEO, Gene Expression Omnibus; WoS, Web of Science.</figcaption></figure>

## Data Renormalization

Most of the time, raw data from GEO correspond to already preprocessed and/or normalized data. However, to remove artefacts and to ensure that each study follows the same procedure we normalized the preprocessed microarray data. To achieve this goal, each selected microarray data set was verified manually to take into account if they were labeled according to a dye swap or a dye switch protocol, if the data were paired, normalized, and log-transformed. According to these indications, the values from the GSE21782 and GSE48136 were normalized with the anapuce package normalization function (73), converted into log base 2, and the sign was reversed in case of a dye swap or dye switch. The abundances of replicated probes or of technical replicates were averaged. Because data from GSE5659 (6 data sets, most of half of the available data) were expressed as ratio, data from the remaining four data sets were also expressed as mRNA ratio calculated as the log abundance in the muscle of the oldest bovine divided by the log abundance in the youngest animal. Based on the similarity of variance abundances between the biological replicates, 59 biological replicates from 10 data sets were included in the meta-analysis.

## Gene Identifier Update to Match Probes with the Last Version of the Bovine Genome

As the data sets were obtained over the last decade and were performed with probes that targeted genes putatively different from the current version of bovine genome, the gene IDs of the probes were updated. Thus, probe IDs were matched with their corresponding official gene symbol updated on the current bovine genome (UCD1.2) using the pipeline sigReannot as previously described by Casel et al. (12). In brief, a BLAST (blastall 2.2.26) (2) was performed using the FASTA files provided within the raw data GPL file of each microarray. The criteria applied for successful BLAST results were continuous stretch of at least 15 base pairs for 50 or 20 bp for 70 mers and a global identity percentage of at least 85%. The sequences of probes were blasted against various banks of the current bovine genome, namely verified bovine transcripts, bovine ncRNA, against the whole genome, and against annex databases namely RefSeq *Bos taurus*, RefSeq *Sus scrofa*, RefSeq *Mus musculus*, RefSeq *Homo sapiens*, and ab initio bovine cDNA (all versions available in June 2018). The
<table><caption>Table 1 Descriptive data of GSE used</caption><thead><tr><th>Study No.</th><th>GSE No.</th><th>Age of the Older Animals, days</th><th>Age of the Younger Animals, days</th><th>Sex</th><th>Breed</th><th>Animals Used for Each Age, n</th><th>Live Body Weight for Older Animals, kg</th><th>Live Body Weight for Younger Animals, kg</th><th>Theoretical Adult Weight, kg</th><th>Physiological Age for Older Animals</th><th>Physiological Age for Younger Animals</th><th>Group</th><th>Maturity</th><th>Genes Analyzed after ID Update, n</th></tr></thead><tbody><tr><td>1</td><td>48136</td><td>215</td><td>155</td><td>C</td><td>Angus</td><td>7</td><td>265</td><td>200</td><td>900</td><td>0.29</td><td>0.22</td><td>1</td><td>early</td><td>7054</td></tr><tr><td>2</td><td>5659</td><td>366</td><td>244</td><td>C</td><td>Belmont Red</td><td>9</td><td>325</td><td>250</td><td>950</td><td>0.34</td><td>0.26</td><td>1</td><td>late</td><td>999</td></tr><tr><td>3</td><td>5659</td><td>214</td><td>92</td><td>F</td><td>Piedmontese × Hereford</td><td>9</td><td>186.2</td><td>93.6</td><td>533</td><td>0.35</td><td>0.18</td><td>1</td><td>inter</td><td>999</td></tr><tr><td>4</td><td>5659</td><td>214</td><td>92</td><td>F</td><td>Wagyu × Hereford</td><td>9</td><td>185.9</td><td>92.1</td><td>570</td><td>0.36</td><td>0.18</td><td>1</td><td>early</td><td>999</td></tr><tr><td>5</td><td>21782</td><td>732</td><td>8</td><td>M</td><td>Chinese Red Steppe</td><td>8</td><td>NS</td><td>NS</td><td>1000</td><td>0.4</td><td>0.1</td><td>2</td><td>early</td><td>10638</td></tr><tr><td>6</td><td>5659</td><td>458</td><td>31</td><td>C</td><td>Belmont Red</td><td>9</td><td>400</td><td>325</td><td>950</td><td>0.42</td><td>0.34</td><td>2</td><td>late</td><td>999</td></tr><tr><td>7</td><td>48136</td><td>275</td><td>155</td><td>C</td><td>Angus</td><td>7</td><td>387</td><td>200</td><td>900</td><td>0.43</td><td>0.22</td><td>2</td><td>early</td><td>6993</td></tr><tr><td>8</td><td>5659</td><td>366</td><td>214</td><td>F</td><td>Piedmontese × Hereford</td><td>9</td><td>269.3</td><td>186.2</td><td>533</td><td>0.51</td><td>0.35</td><td>2</td><td>inter</td><td>999</td></tr><tr><td>9</td><td>5659</td><td>366</td><td>214</td><td>F</td><td>Wagyu × Hereford</td><td>9</td><td>267.5</td><td>185.9</td><td>570</td><td>0.52</td><td>0.36</td><td>2</td><td>early</td><td>995</td></tr><tr><td>10</td><td>48136</td><td>375</td><td>155</td><td>C</td><td>Angus</td><td>7</td><td>560</td><td>200</td><td>900</td><td>0.62</td><td>0.22</td><td>3</td><td>early</td><td>6996</td></tr></tbody></table>

The comparison between 2 ages from a unique GSE were numbered; comparison concerned data from the older bovine divided by data from the younger one. In the sex column: C, castrated; F, female; M, male. The physiological age is a value between 0 and 1, calculated as the live body weight at slaughter divided by the assumed adult weight. This physiological age is an indicator of animal maturity. The type of bovine is described by their rearing purpose as meat production (beef), milk production, or dual (meat and milk production) and by their maturity: either early-maturing breed (early), late-maturing breed (late) or intermediate (inter). The number of remaining genes after ID update shows those that were subjected to meta-analysis.

success rate of gene name update ranged from 23.7% (E6) to 64.7% (E5) of probes IDs updated with a gene name.

## Meta-analysis

Since the number of genes were different between studies (from 995 to 10,638 gene names), the 715 common genes across the 10 data sets, were gathered and loaded in the meta-analysis. The meta-analysis was performed using the mRNA ratio abundances based on *P* value combination by the inverse normal method. Meta-analysis calculates for each study, the *P* value of differential ratio for paired data with moderated *t* test method (Limma) using the package metaMA (version 3.1.2) (42). *P* values were corrected for multiple testing with a Benjamini-Hochberg correction and were considered as significant when *P* < 0.05. Whenever a gene was found several times in a same sample, the mean of the abundance ratios for this gene was calculated and used for the meta-analysis. By the end, meta-analysis was performed for 715 unique genes that had expression values across all the 10 studies. Based on the Benjamini-Hochberg adjusted *P* values (< 0.05) a list of differentially abundant mRNA between two postnatal ages was realized.

## Cluster Analysis and Annotation according to Gene Ontology

Heatmap and clustering methods were applied to group mRNA together based on the similarity of their abundance pattern. Coexpressed genes were identified by a hierarchical clustering procedure (hclust) implemented in the R package stats (version 3.6.0) on which the best cutree function was used (https://github.com/larmarange/JLutils/blob/master/man/best cutree.Rd) to determine the best partition to cut the dendrogram based on the highest relative loss of inertia criterion. Best cutree was applied to the average ratio values of the differentially abundant mRNA between two postnatal ages. Then a graphical representation was performed with the p Heatmap package (version 1.0.12 (36)). Gene names within each cluster were analyzed for Gene Ontology (GO) enrichment using ProteINSIDE [https://www.proteinside.org/ (32)]. To take advantage of the better annotations within human, the GO enrichment analysis was performed in the human species, and GO terms with a Benjamini-Hochberg adjusted *P* value < 0.001 were considered.

## Functional Analyses

The topological analysis of networks constructed between the gene products of the DEG and the human proteins [because protein-protein interactions (PPI) were much more described in human than bovine] declared in Psicquic (3) was realized with ProteINSIDE (32). We recorded the number of interactions, as well as centralities that were proven as efficient to reveal proteins that play important roles in a network, namely betweenness and closeness centralities. In brief, betweenness centrality quantifies how frequently a node is on the shortest path between every pair of nodes to detect bottlenecks in a network. Closeness centrality quantifies how short minimal paths are from a given node to all others; a large closeness indicates that a node is close to the topological center of the network. Centralities were calculated on a PPI network that linked the gene products identified as differentially abundant (DEG) in the meta-analysis and interacting proteins declared in IntAct, BioGrid, and UniProt databases. Values of betweenness and closeness centralities were downloaded to identify proteins with the highest values of centralities. Only PPI that were reported by experiments were considered; thus 6,903 PPI were identified between 3,249 proteins; of these 214 were DEG from the meta-analysis (version of ProteINSIDE Database was 1.2.11 and of Psicquic 1.4.11).

## Analysis of Genes Mapped in Quantitative Trait Locus Regions

The genes that were found as differentially expressed between two ages were mapped against quantitative trait locus (QTL) regions using
the ProteINSIDE webservice that uses resources from the Cattle QTL database (https://www.animalgenome.org/cgi-bin/QTldb/BT/index). We mapped genes in regions where QTL related to muscle growth (lean meat yield, longissimus muscle area), intramuscular fat content (intramuscular fat, marbling score), and carcass adiposity (fat thickness at the 12th rib, subcutaneous fat, rib fat) were reported.

## Data Analysis according to the Physiological Age of Bovine

Across lifespan, the rate of hypertrophy and growth of longissimus thoracis muscle, as well as the percent of intramuscular fat are very different among bovine breeds (26) because of differences in adult weight. Regarding muscle composition, the amount of intramuscular lipids strongly increases from 8 to 14 mo of age in early- and late-maturing breeds (26) corresponding to a degree of maturity (DoM) of 0.30 ± 0.07. Thus, to ensure relevant data-mining across the five pure- or cross-breeds under the meta-analysis we calculated the DoM (Table 1) of the animals, which is the ratio between slaughter body weight (kg) and adult body weight (kg). The 10 studies used for the meta-analysis were then categorized according to the early (group 1), mid (group 2), and more than the mid (group 3) animal growth, keeping in mind that marbling deposition begins from 0.30 DoM and is high after 0.50 DoM (Table 1). The first group is composed from studies E1, E2, E3, and E4, which compared muscles from bovines from 0.18 up to 0.36 of adult weight or DoM. The second group is composed from studies E5, E6, E7, E8, and E9, which compared bovines from 0.1 to 0.52 of DoM. The third group with the study E10 reported results from bovine muscles from 0.22 to 0.62 of DoM. An average ratio of mRNA abundance per group of DoM was calculated as the mean of the individual abundance ratio between old and young bovines per studies and by group.

# RESULTS

## Identification of DEG

We considered 125 bovines from Angus, Belmont Red, Chinese Red Steppe pure breeds, as well as from Piedmontese × Hereford and Wagyu × Hereford crossbreds over the time course of longissimus thoracis muscle growth, i.e., from 31 to 732 days of postnatal life (Table 1). A total of 10 microarray data sets have been used, which provided 715 genes that were quantified for their mRNA abundance in all the data sets. Of these, 238 genes were identified as differentially expressed (DEG, *P* value adjusted by the Benjamini-Hochberg correction < 0.05) between two ages and whatever the breed. Of these, 214 mRNA were overabundant and 24 were less abundant when older bovines were compared with the youngest animals. The top 20 mRNA with the highest significant (*P* < 0.05) fold-change were FMO5, CYCS, TCAP, SMG7, RPL29, CTSS, PGAM2, STRA6, CYC1, STXBP1, ASPH, MAPK14, GPRC5C, ARPC5, DUSP1, RPS12, ACSL3, TIMP1, and HSPA5, all of them having a ratio >2.5 (log 2 fold change >1.35 when old bovines were compared with the youngest ones).

Among the 238 DEGs, 97 mRNA were identified as differentially abundant between two ages in at least two studies and the meta-analysis. Of these, 15 DEGs were identified in more than half of the studies and the meta-analysis: namely, MB

was DEG in six studies and the meta-analysis; ASB12, FN1, LRRC20, PKM in five studies and the meta-analysis; and ACTN3, ATF4, CLNS1A, CLU, COPG1, CYC1, GPI, PFKFB3, PGAM2, PGK1 in four studies and the meta-analysis

Among the 238 DEGs between two ages, 17 were exclusively identified in the meta-analysis, namely APOE, CAMTA2, CTSH, HMBS, HSP90AA1, LDLR, MGC148714, MXRA8, MYL12A, NNT, RPL32, RPLP0, RPS13, SGCE, SOD2, STAC3, and UQCR11.

## Cluster Analysis and GO Annotation

To group DEG with similar abundance variation, we applied a hierarchical clustering to the 238 mRNA differentially abundant according to the age. Three clusters were identified with two global profiles consisting of up- (clusters 1 and 3) and downregulated (cluster 2) gene expression when older animals were compared with younger ones (Figs. 2 and 3).

Cluster 1 contains 78 genes with the highest abundances in studies 10, 5, and 3 comparing Angus at 375 versus 155 days, Chinese Red Steppe at 732 versus 31 days, and Angus at 215 versus 155 days. The lowest abundances were recorded in studies 8 and 9 comparing Piedmontese × Hereford and Wagyu × Hereford at 366 versus 214 days (Fig. 3). The top biological processes related to 78 genes that are annotated by enriched (*P* < 1.13E-03) GO terms are reported in Fig. 2. Of the 78 annotated genes, ACTN3, TNNT1, TNNC2, MYBPC1, TPM2, MYL3, STAC3, MB, and MYH1 are annotated by five terms related to muscle organization and contraction (muscle filament sliding, skeletal muscle contraction, slow-twitch skeletal muscle fiber contraction, muscle contraction, transition between fast and slow fiber, sarcomere organization, and skeletal muscle cell differentiation). PGK1, LDHA, OGDH, PKM, GPI, MDH2, ATF4, GOT2, PGK1, PDHA1, BSG, and PLN are annotated by six GO terms related to glucose metabolism and signaling (glycolytic process, gluconeogenesis, canonical glycolysis, pyruvate metabolic process, response to insulin). DHRS3, NNT, MDH2, PDHA1, OGDH, CYC1, SOD2, LDHA, ATP5F1B, COX7B, ATP2A1, GOT2 are annotated by seven terms related to oxidative metabolism (oxidation-reduction process, tricarboxylic acid cycle, proton transmembrane transport, oxaloacetate metabolic process, 2-oxoglutarate metabolic process, mitochondrial ATP synthesis coupled proton transport, and ATP biosynthetic process). Less numerous DEG are annotated by terms related to cell volume homeostasis (CLNS1A, AQP1), skeletal muscle (FOS, ASB2), or fat cell (MB, FABP4) differentiation, as well as negative (ANGPTL4, NNT, AQP1, CD74) and positive (LDHA, TXNIP, ATF4, CTSD) regulation of apoptotic process.

Cluster 2 contains 24 genes with the highest abundance ratios in studies 4 and 9 comparing Wagyu × Hereford at 214 versus 92 days and at 366 versus 214 days. The lowest abundance ratios were recorded in studies 1 and 10 comparing Angus at 215 versus 155 and at 375 versus 155 days (Fig. 3). The top biological processes related to the 24 genes that are

Fig. 2. Heatmap and clustering of 238 differentially expressed genes (DEGs) from the meta-analysis, have highlighted 3 clusters. Gene names by cluster were subjected to a Gene Ontology (GO) annotation (ProteINSIDE). Results of enrichment are expressed as -log10 (*P* value) to visually plot them on graphs, which means that -log10 (*P* value) of 3, 2, and 1.3 correspond to a *P* values of 0.001, 0.01, and 0.05, respectively. The genes or the number of genes that have been annotated by a GO term are written next to the bar.
<figure><img src="image_3.png" /></figure>

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

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

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

<figure><img src="image_7.png" /></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 3. Radar plot of mean log ratio profiles of the mRNA abundance by cluster and by study. A log ratio of -1 or 1 indicates a twofold variation in abundance. Shown near the name of the studies (E1-E10) are the breed, and the physiological ages are compared. DoM, degree of maturity.</figcaption></figure>

annotated by enriched (*P* < 3.69E-03) GO terms are reported in Fig. 2. Of these 24 genes, RPS9, RPL18, RPL32, RPLP0, RPS13, IGFBP3, ATG5, RBBP7 are annotated by terms related to translation (translational initiation, translation, cytoplasmic translation, posttranslational protein modification). RBBP7, HYAL2, ENO1, H3F3A are annotated by GO terms related to cell cycle (negative regulation of cell growth, negative regulation of chromosome condensation, telomeric heterochromatin assembly, negative regulation of G0 to G1 transition). IGFBP3, CYP19A1, RPLP0 are annotated by terms related to hormone or cytokine signaling (positive regulation of insulin-like growth factor receptor signaling pathway, positive regulation of estradiol secretion, interleukin-12-mediated signaling pathway).

*Cluster 3* contains 136 DEG with the highest positive abundance ratios in *studies 1, 7, and 10* comparing Angus at 215 versus 155 days, at 155 versus 275 days, and at 375 versus 155 days of life. The lowest negative abundance ratios were recorded in *studies 8* and 2 comparing Piedmontese × Hereford at 366 versus 214 days and Belmont Red at 366 versus 244 days of life (Fig. 3). The top biological processes related to the 132 genes that are annotated by enriched (*P* < 1.61E-6) GO terms are reported in Fig. 2. MAPRE2, MYOG, DDIT3, GNAI1, USP16, PIK3C3, ENSA, TERF1, NR3C1, ANAPC11, DUSP1, CSF2, CTSH, CXCL12, FN1, CAPNS1, NPM1, TIMP1, ERBB2, HAS2, SFRP4, PHB, ASPH, TIMP2 are annotated by GO terms related to cell cycle, and positive or

negative regulation of cell population proliferation. FSTL1, TNF, CSRP3, PMCH, SFRP4, MGP, BIN1, CREG1, MYOG, MAPK14, ELN, CSRP3, SGCE, UTRN are annotated by GO terms related to muscle development (multicellular organism development, skeletal muscle tissue development, muscle organ development). ADRB2, PCBP1, MYOG, TNF, NR3C1, YBX1, CAMTA2, NFIC, SERPINE1, MAPK14, NFKB2, DDIT3, NPM1, USP16, CSRP3, SMG7, TERF1, HSP90AA1, APOE, PHB, ASPH, YWHAH, NCOR1, YBX3, CREG1, CTDSP2 and ERBB2 are annotated by GO terms related to the regulation of transcription (positive regulation of transcription by RNA polymerase II, telomere maintenance via telomerase, positive regulation of transcription, DNA-templated, negative regulation of transcription by RNA polymerase II, regulation of transcription by RNA polymerase II). GSPT1, RPS2, MRPL2, RPL19, EIF4G2, RPL29, RPL27A, MRPL35, RPL15, RPS12, RPL24, RPS12, HSPA14, CST3, APOE, FSTL1, TIMP1, FN1, MXRA8 are annotated by GO terms related to the regulation of translation (translation, translational initiation, cytoplasmic translation, posttranslational protein modification). Other biological processes relevant to muscle structure and metabolism were enriched within the DEG from *cluster 3*: oxidation-reduction process (NDUFB9, IFI30, CYCS, SDHA, FMO5, ACOX2, UQCR11, ACADS, ASPH, IDH3B), lipid metabolic process (ACADS, FDFT1, ACOX2, ACSL3, LDLR, LIPE, PLA2G2A, APOE, PLA2G12A), muscle contraction (TPM1, FXYD1, MYL12A, ASPH, UTRN), and extra-
cellular matrix disassembly (TIMP2, CTSK, CAPNS1, TIMP1, CTSS).

## Network Analysis and Linkage with Muscle Growth

**Identification of genes with the highest regulatory potential through network analysis.** To identify the genes with the highest regulatory potential we evaluated the closeness and betweenness centralities of the DEG product within a PPI network. Of the 50 top genes with the highest closeness and betweenness centralities, 29 were found in the two lists of centralities. Of these YBX1, HSPA5, STXBP1, YWHAH, SDHA, HSP90AA1, and ERBB2 are among the genes with the highest fold change. Six additional DEG (RPL29, CYCS, MAPK14, FDFT1, MYOG, AQP1) are among the genes with the highest fold change, and closeness or betweenness centralities (Fig. 4).

**Identification of QTL related to muscle growth and marbling around DEG identified by meta-analysis.** Of the 238 DEG identified by the meta-analysis, OGDH, EIF4A2, ANGPTL4, RBM38, ROMO1, CYC1, AQP1, SLA, DES, IGFBP3 are near QTL for marbling fat or marbling score. Moreover, CLNS1A, PKM, ACAT1, COQ8A, ATG101, CSDE1, CLU, GOT2, TXNIP, LDHA, CAMK2D, ANGPTL4, HMBS, CYC1, IRF8, SBDS, MYH1, PDLIM1, SLA, LRRC20, ATP2A1, VEGFC are near QTL for lean meat yield. STAC3, CKM, OGDH, ATG101, GPI, TPPP3, EIF4A2, NNT, AQP1, RTN2, DES, ATG5, IGFBP3 were shown to be near QTL related to longissimus muscle area. Of these, GOT2, TXNIP (longissimus muscle area), CYC1 (marbling and lean meat yield QTL), and AQP1 (marbling and longissimus muscle area) are among the top 50 DEG with the highest fold changes.

**Identification of the most differentially abundant mRNA regarding the phases of muscular growth.** Genes with ratio of abundance between two ages (not log2 ratio) < 0.5 and > 2 were mined according to the growth phase of the bovine. We

listed mRNA with the top fold change during the early (up to 0.32 DoM, group 1), mid (up to 0.52 DoM, group 2), and more than the mid (up to 0.62 DoM, group 3) animal growth. Among the 238 DEG of the meta-analysis, the number of mRNA that are twofold higher in older than younger muscle are 42, 20, and 204 for groups 1, 2, and 3. Of these and when a maximum of 30 mRNA was considered by group, TCAP and SMG7 were identified in groups 1 and 2, ASPH, STXBP1, RPS12, CYCS, MYOG in groups 2 and 3, and PFKFB3, ATRAUD, SDHA, ARPC5, FMO5, PHB, ELN, STRA6, TIMP1 in groups 1 and 3 (data not shown). The number of genes downregulated more than twofold are 10, 3, and 10 for groups 1, 2, and 3, respectively. Of these, TMED9 was identified in groups 2 and 3, while ENO1 and MATR3 were identified among the most regulated genes in groups 1 and 3 (data not shown).

## DISCUSSION

Postnatal muscle growth requires the concerted regulation of numerous genes to ensure hypertrophy, but also metabolic and contractile differentiations. Moreover, muscle is not just composed of muscular contractile cells; satellite cells, adipocytes, fibrocytes, endothelial cells, pericytes, nonsatellite stem cells within the muscle also contribute to muscle growth by hypertrophy (increase in cell volume) and hyperplasia [increase in cell number (9)]. The rates and waves of growth for the different muscular cell types vary according to age and at similar age according to the breed (26). Thus, the proportion of muscular, fibrous, and adipose (called marbling) cells determines the muscle composition with consequences on meat-eating qualities and thus on the economic value of bovine meat. Using publicly available microarray data from muscles of both high muscling (Belmont Red, Piedmontese × Hereford), high marbling (Angus, Wagyu × Hereford), and intermediate (Chinese Red Steppe) bovine breeds, we carried out a meta-analysis to identify DEGs between two postnatal ages distributed from 10 to 62% of mature weight, also cited as 0.1–0.62 DoM. As main results, we identified 238 DEG between two ages whatever the bovine breeds. Of these, 97 were identified as DEG in at least two independent studies and in the meta-analysis, when we carried out a one-to-one data set analysis, while 17 were never described before as DEG between two ages. Among the DEG, we emphasized mRNA that could be transcriptional signatures of muscular lipid deposition as the result of regulated metabolisms or of the differentiation of news adipocytes.

### Transcriptional Signature of the Dynamic Regulation of Glycolytic and Oxidative Metabolisms

Besides protein synthesis, glycolytic and oxidative metabolisms are the pathways that sustain postnatal muscle growth. Muscles use both carbohydrates and fatty acids as energy sources, according to a ratio that differs between glycolytic and oxidative fibers of a muscle and according to the age. An increased oxidative metabolism is generally concomitant to an increase in marbling (8, 11). Accordingly, we detected a muscular twofold mean overexpression for genes related to glycolysis (GPI, PDHA1, PGAM2, PGK1, PKM) and tricarboxylic acid cycle (IDH3B, MDH2, OGDH, SDHA, SUCLG1) when bovines from 0.32 and 0.18 DoM were compared (group 1), and for the five breeds. However, the abundances of these 10 mRNA were similar between bovines at 0.52 or 0.10 DoM

<figure><img src="image_9.png" /><figcaption>Fig. 4. Venn diagram (5) highlighting genes among the 238 differentially expressed genes (DEGs) of the meta-analysis that were found with the highest betweenness centrality (BC), closeness centrality (CC) and fold change (FC).</figcaption></figure>
(group 2), while they were 2.5 (glycolytic)- to 5.5 (oxidative)- fold higher in muscle from bovine at 0.62 than 0.22 DoM (group 3). These results may reflect known variation in the metabolic properties of muscle at puberty that generally occurs in cattle at a DoM of ~0.5. The DoM near puberty represents an inflection point, as has been reported in the quadratic and cubic relationships of ICDH (tricarboxylic acid cycle enzyme) and LDH (glycolytic enzyme) activities regressed with the DoM of the bovines in several breeds (58). The dynamic nature of the abundance of mRNA related to glycolytic and oxidative metabolisms with DoM may reflect the use of glycolytic metabolism to supply energy during periods of rapid growth before puberty and the use of glycolytic and more importantly oxidative metabolism during muscle aging probably in relationship to an increased proportion of lipid deposit (9). We observed a similar dynamic nature of mRNA abundance related to fast glycolytic or slow oxidative fibers. Indeed, MYH1 mRNA (coding for MyHC-2x isoform), related to fast glycolytic fiber, was overexpressed in *groups 1* and *3*, while the slow-type mRNA such as TNNT1 related to oxidative fibers was overexpressed only in *group 3*.

From the current meta-analysis, we proposed that dynamic regulation of the glycolytic and oxidative metabolisms may result from the signaling of two muscular master genes MYOG and MAPK14. The myogenin (MYOG) and the mitogen-activated protein kinase p38 alpha (MAPK14) mRNA are among the top fold change according to the age and with a high regulatory potential suggested by high centralities within the network. MYOG and MAPK14 are key transcription factors of muscle-specific genes, and they are recognized markers of differentiated fibers and of differentiating satellite cells (27, 44, 54, 61). However, increased MYOG and MAPK14 mRNA abundance during muscle growth and after puberty (top fold change in *group 3* with an increase of 9- and 10-fold) is much less documented than an overexpression during myogenesis, which raises the question of the significance of MYOG's and MAPK14's persistent expression in adult bovine muscle deciphered by the current study. One possible explanation for MYOG overexpression in bovine muscle, especially after puberty, is to shift muscle enzymatic activity from glycolytic toward oxidative metabolism, as has been reported in several mice models (24).

In differentiated contractile muscular fibers, MAPK14 was shown to be a key regulatory gene of glucose and insulin signaling pathways in muscle from monogastric species (66, 69). Additionally, TXNIP may contribute to or be a sign of enhanced oxidative metabolism. TXNIP was found as a mRNA with a high fold-change and as a gene close to a chromosomal region with a QTL of meat yield. The gene TXNIP encodes the thioredoxin-interacting protein that has been shown to regulate both insulin-dependent and insulin-independent pathways of glucose uptake in human skeletal muscle (45). The current meta-analysis quantified TXNIP mRNA overexpression of 2.6- and 5.6-fold in *groups 1* and *3* of DoM, which parallels the overexpression of glycolytic mRNA however with a higher fold-change. Thus, MOYG, MAPK14, and TXNIP may contribute to the dynamic regulation of glycolytic and oxidative metabolisms within bovine muscle from 0.1 to 0.62 of DoM.

Transcriptional Signature Related to Lipid Metabolism

Of the 17 DEG between two ages that were identified exclusively by the meta-analysis, APOE, LDLR, MXRA8, and HSP90AA1 upregulation may be related to lipid metabolism and enhanced marbling. Indeed, circulating triacylglycerols such as VLDL and chylomicrons contribute to ATP production within the contractile muscular cells and to energy storage within adipocytes in bovine muscle through lipoprotein lipase (LPL) activity (10). APOE and LDLR are among the regulatory proteins involved in lipid delivery to muscle. Apolipoprotein E (APOE) is synthesized primarily by the liver and the adipose tissue but also by the muscle, and it regulates blood lipid and lipoprotein levels in multiple ways. For example, APOE acts as a ligand for low-density lipoprotein receptor (LDLR) and affects the activities of LPL (34). Furthermore, accumulated data report that APOE contributes to adipogenesis by inducing the assembly of triglycerides in adipocytes (40). LDLR mRNA was shown to be highly abundant in the heart but also present in skeletal muscle and adipose tissues according to an age-related pattern in mice (68). Studies using LDLR-deficient and transgenic mice have provided evidence for the contribution of LDLR in VLDL-triglyceride metabolism in peripheral tissue including muscle, and for triglyceride storage in adipocytes (64). Interestingly, another major physiological regulator of triacylglycerol metabolism through LPL activity in muscle is angiopoietin-like protein 4 (ANGPTL4) protein. The current study reports a 2 (APOE and LDLR)- or 4.5 (ANGPTL4)-fold overexpression of mRNA abundance in 0.62 compared with 0.22 DoM bovines (group 3). These overexpressions may be a sign of regulation of triacylglycerol breakdown and storage within muscular adipose cells, if we assume that in bovines as in monogastric species these mRNA are much more abundant in adipose than muscular cells. In agreement with this hypothesis of the transcriptional signature of adipose cells, LIPE and FABP4 mRNA abundances followed the same pattern as APOE and LDLR mRNA and were among the 238 DEG identified by the current meta-analysis. LIPE is an adipose enzyme involved in the lipolysis of triacylglycerol stored within adipocytes, and a balance between triacylglycerol synthesis and lipolysis was repeatedly observed during lipid deposition (65). COPG1 is a coatomer GTPase involved in lipolysis signaling within skeletal muscle (17) and was also overexpressed with age in the present study. FABP4 is a well-known adipose cell marker within the muscle, and its mRNA abundance was shown to be overexpressed with increased marbling of bovine muscle (31, 36).

MXRA8 and HSP90AA1 mRNA abundance was also identified among the overexpressed DEG between two ages (2- and 2.5-fold for MXRA8 and HSP90AA1 between bovines of 0.62 and 0.22 or 0.52 and 0.1 DoM). MXRA8 is a transmembrane protein that can modulate activity of various signaling pathways (28). However, to date a link between MXRA8 and lipid metabolism has been documented by only one publication. MXRA8 was indeed shown as an upregulated gene when bovines with high carcass adiposity and marbling were compared with slightly adipose bovines (16). The heat shock protein HSP 90-alpha (HSP90AA1) is related to muscular physiology by regulating transforming growth factor β (TGF-β), by folding the myosin motor domain and organizing the filament structure of myofibers (22). Inhibition of
HSP90AA1 function was reported to block TGF-β-induced signaling and transcriptional responses through enhanced ubiquitination and degradation of Smad3 in various cell lines (62). In bovine reduced myostatin abundance, a member of the TGF-β family was associated with higher muscle mass and lower intramuscular fat content (1). Thus, the upregulation of HSP90AA1 (one of the mRNA with the highest fold change and centralities) in the muscle of aging bovine may contribute to the upregulation of TGF-β with a putative positive impact on adipogenesis or lipid metabolism. In agreement with these various signaling mechanisms related to an increased lipid deposition with age, we reported an overexpression of TNF-βI mRNA in 0.62 compared with 0.22 DoM bovines. In bovine muscle, decreased abundance of TNF-βI mRNA was associated with increased longissimus muscle area and decreased marbling (13).

### Transcriptional Signature of a Balance between Adipogenic and Myogenic Lineages

Muscle grows by hypertrophy of muscular contractile cells and by hyperplasia and hypertrophy of other myogenic cell types, among them the most important being satellite cells. Indeed, the number of fibers being fixed at birth in bovine (47), satellite cells are a well-known muscle resident myogenic progenitor recruited for muscle growth. Once recruited, satellite cells express muscle-specific transcriptional factors including MYOG and MAPK14, found upregulated in *groups 2* and mainly *3* for MYOG and in *group 3* for MAPK14. However, an age-related decrease in satellite cell activity is thought to exist in bovines as in monogastric species and may partly result from MAPK14 upregulation. Indeed, inhibition of MAPK14 signaling was shown to maintain bovine satellite cells stemness (20) and to decrease the adipogenic transdifferentiation of C2C12 myoblasts (49). Beside satellite cells, fibrocyte/adipocyte progenitors are a well-known muscle resident cell population recruited for muscle growth under the control of several proteins, among these, bone morphogenetic protein 7 (BMP 7) (9, 59). Y-box-binding protein (YBX1) depletion was shown to lead to a loss of Bmp7 expression and decreased adipogenesis in mice (46). The protein 14-3-3 eta encoded by the gene YWHAH is involved in a large spectrum of both general and specialized signaling pathways. siRNA-mediated silencing of YWHAH leads to an increased number of adipose cells differentiated from mesenchymal skeletal stem cells (19). The gene ERBB2 encodes the receptor tyrosine-protein kinase erbB-2, which has been shown to regulate the expression of acetylcholine receptors and glucose transport within normal skeletal muscle, as well as to favor the activation of quiescent satellite cells in mice (25) and protein synthesis in culture bovine satellite cells treated with testosterone analog (67). Moreover, ErbB2 signaling downregulation was proposed to enhance the differentiation of adipogenic progenitors through upregulation of the Akt pathway and downregulation of Erk-1 (60). MAPK14, YBX1, YWAH, and ERBB2 were among the mRNA with the highest fold change and with the highest centralities and are thus genes with high regulatory potential through various pathways. Strong upregulation of those genes, especially in muscle from 0.62 compared with 0.22 DoM bovines, may be favorable (YBX1) or repressive (MAPK14, YWAH, ERBB2), depending on the presence of

their downstream targets, the differentiation of myogenic progenitors into adipogenic cells. Collectively, these results highlight a probable fine-tuning of gene expressions for muscle growth probably involving several cell populations and cross talk between cell types.

### Conclusion

The present study was designed to identify DEG between two different ages in growing bovine muscle, by using meta-analysis of multiple microarray data sets throughout the animal lifespan in five breeds or crossbreeds. The difference between two ages was used as a proxy of intramuscular lipid content variation, in the absence of enough muscular microarray data sets carried out in bovines divergent only by marbling. In total, we identified 238 DEG between two ages; of these, many are related to the known dynamic regulation of contractile and metabolic properties of muscle. Network and fold-change ranking analysis highlights some potential “hub” genes that were proposed to contribute to the regulation of dynamic contractile and metabolic pathways. We also highlight transcriptional signatures of lipid metabolism and of the contribution of muscle-resident progenitor cells to the myogenic or the adipogenic growth with possible impact on the lean-to-fat ratio within muscle and thus on marbling. These findings have practical implications as they pave the way through the second step of the marbling biomarkers discovery pipeline; that is the qualification step (53). Indeed, we provide a list of DEG valid for different breeds that are regulated in muscle when bovines, after and before puberty, are compared as a proxy of marbling. The qualification step requires the demonstration that the differential abundance of a candidate biomarker observed in the current meta-analysis driven discovery phase is seen by an alternative method. In the case of mRNA, quantitative PCR should be a valuable method to confirm current DEG as differentially abundant between groups of bovine divergent for marbling and also to investigate the relationship between the mRNA abundance and the intramuscular fat content. Thus, the present meta-analysis achieved the expected goals: the verification and validation of independent results and the production of conclusive findings related to transcriptional signatures of lipid deposition. Numerous genes were not discussed in this paper and are made available to meat physiologists interested in the regulation of protein deposition or in the regulation of growth of the connective tissue.

### ACKNOWLEDGMENTS

We thank Isabelle Hue, Isabelle Cassar-Malek and Brigitte Picard for helpful discussions.

### GRANTS

We thank Region Auvergne-Rhône-Alpes and INRAe-Phase Division for a J. Bazile grant. This study was realized in the framework of the French government IDEX-ISITE initiative 16-IDEX-0001 (CAP 20-25).

### DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the authors.

### AUTHOR CONTRIBUTIONS

J.B. and M.B. conceived and designed research; J.B. and F.J. performed experiments; J.B., F.J., P.D., M.R., C.K., D.L., and M.B. analyzed data; J.B. and M.B. interpreted results of experiments; J.B. prepared figures; J.B. and M.B. drafted manuscript; J.B., F.J., P.D., C.K., D.L., and M.B. edited and
revised manuscript; J.B., F.J., P.D., M.R., C.K., D.L., and M.B. approved final version of manuscript.

REFERENCES

1. Aiello D, Patel K, Lasagna E. The myostatin gene: an overview of mechanisms of action and its relevance to livestock animals. *Anim Genet* 49: 505–519, 2018. doi:10.1111/age.12696.

2. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. *J Mol Biol* 215: 403–410, 1990. doi:10.1016/S0022-2836(05)80360-2.

3. Aranda B, Blankenburg H, Kerrien S, Brinkman FSL, Ceol A, Chautard E, Dana JM, De Las Rivas J, Dumousseau M, Galeota E, Gaulton A, Goll J, Hancock REW, Isserlin R, Jimenez RC, Kersse- makers J, Khadake J, Lynn DJ, Michaut M, O'Kelly G, Ono K, Orchard S, Prieto C, Razick S, Rigina O, Salwinski L, Simonovic M, Velankar S, Winter A, Wu G, Bader GD, Cesareni G, Donaldson IM, Eisenberg D, Kleywegt GJ, Overington J, Ricard-Blum S, Tyers M, Albrecht M, Hermjakob H. PSICQUIC and PSISCORE: accessing and scoring molecular interactions. *Nat Methods* 8: 528–529, 2011. doi:10.1038/nmeth.1637.

4. Baik M, Kang HJ, Park SJ, Na SW, Piao M, Kim SY, Fassah DM, Moon YS. TRIENNIAL GROWTH AND DEVELOPMENT SYMPOSIUM: Molecular mechanisms related to bovine intramuscular fat deposition in the longissimus muscle. *J Anim Sci* 95: 2284–2303, 2017. doi:10.2527/jas2016.1160.

5. Bardou P, Mariette J, Escudié F, Djemiel C, Klopp C. jvenn: an interactive Venn diagram viewer. *BMC Bioinformatics* 15: 293, 2014. doi:10.1186/1471-2105-15-293.

6. Baron D, Dubois E, Bihouée A, Teusan R, Steenman M, Jourdon P, Magot A, Péréon Y, Veitia R, Savagner F, Ramstein G, Houlgatte R. Meta-analysis of muscle transcriptome data using the MADMuscle database reveals biologically relevant gene patterns. *BMC Genomics* 12: 113, 2011. doi:10.1186/1471-2164-12-113.

7. Barrett T, Wilhite SE, Ledoux P, Evangelista C, Kim IF, Tomashevsky M, Marshall KA, Phillippy KH, Sherman PM, Holko M, Yefanov A, Lee H, Zhang N, Robertson CL, Serova N, Davis S, Soboleva A. NCBI GEO: archive for functional genomics data sets–update. *Nucleic Acids Res* 41, D1: D991–D995, 2013. doi:10.1093/nar/gks1193.

8. Bazile J, Picard B, Chambon C, Valais A, Bonnet M. Pathways and biomarkers of marbling and carcass fat deposition in bovine revealed by a combination of gel-based and gel-free proteomic analyses. *Meat Sci* 156: 146–155, 2019. doi:10.1016/j.meatsci.2019.05.018.

9. Bonnet M, Cassar-Malek I, Chilliard Y, Picard B. Ontogenesis of muscle and adipose tissues and their interactions in ruminants and other species. *Animal* 4: 1093–1109, 2010. doi:10.1017/S1751731110000601.

10. Bonnet M, Faulconnier Y, Hocquette J-F, Bocquier F, Leroux C, Martin P, Chilliard Y. Nutritional status induces divergent variations of GLUT4 protein content, but not lipoprotein lipase activity, between adipose tissues and muscles in adult cattle. *Br J Nutr* 92: 617–625, 2004. doi:10.1079/BJN20041240.

11. Bonnet M, Faulconnier Y, Leroux C, Jurie C, Cassar-Malek I, Bauchart D, Boulesteix P, Pethick D, Hocquette JF, Chilliard Y. Glucose-6-phosphate dehydrogenase and leptin are related to marbling differences among Limousin and Angus or Japanese Black x Angus steers. *J Anim Sci* 85: 2882–2894, 2007. doi:10.2527/jas.2007-0062.

12. Casel P, Moreews F, Lagarrigue S, Klopp C. sigReannot: an oligo-set re-annotation pipeline based on similarities with the Ensembl transcripts and Unigene clusters. *BMC Proc* 3, Suppl 4: S3, 2009. doi:10.1186/1753-6561-3-s4-s3.

13. Casey TM, Walker JF, Bhide K, Thimmapuram J, Schoonmaker JP. Global transcriptional differences in myokine and inflammatory genes in muscle of mature steer progeny are related to maternal lactation diet and muscle composition. *Physiol Genomics* 50: 884–892, 2018. doi:10.1152/physiolgenomics.00060.2018.

14. Chaze T, Meunier B, Chambon C, Jurie C, Picard B. In vivo proteome dynamics during early bovine myogenesis. *Proteomics* 8: 4236–4248, 2008. doi:10.1002/pmic.200701101.

15. Chaze T, Meunier B, Chambon C, Jurie C, Picard B. Proteome dynamics during contractile and metabolic differentiation of bovine foetal muscle. *Animal* 3: 980–1000, 2009. doi:10.1017/S17517311109004315.

16. Clark DL, Boler DD, Kutzler LW, Jones KA, McKeith FK, Killefer J, Carr TR, Dilger AC. Muscle gene expression associated with increased

marbling in beef cattle. *Anim Biotechnol* 22: 51–63, 2011. doi:10.1080/10495398.2011.552031.

17. Covington JD, Galgani JE, Moro C, LaGrange JM, Zhang Z, Rustan AC, Ravussin E, Bajpeyi S. Skeletal muscle perilipin 3 and coatomer proteins are increased following exercise and are associated with fat oxidation. *PLoS One* 9: e91675, 2014. doi:10.1371/journal.pone.0091675.

18. Dayanidhi S, Lieber RL. Skeletal muscle satellite cells: mediators of muscle growth during development and implications for developmental disorders. *Muscle Nerve* 50: 723–732, 2014. doi:10.1002/mus.24441.

19. De S, Kline D. Evidence for the requirement of 14-3-3eta (YWHAH) in meiotic spindle assembly during mouse oocyte maturation. *BMC Dev Biol* 13: 10, 2013. doi:10.1186/1471-213X-13-10.

20. Ding S, Swennen GNM, Messmer T, Gagliardi M, Molin DGM, Li C, Zhou G, Post MJ. Maintaining bovine satellite cells stemness through p38 pathway. *Sci Rep* 8: 10808, 2018. doi:10.1038/s41598-018-28746-7.

21. Dodson MV, Allen RE, Du M, Bergen WG, Velleman SG, Poulos SP, Fernyhough-Culver M, Wheeler MB, Duckett SK, Young MRI, Voy BH, Jiang Z, Hausman GJ. Evolution of meat animal growth research during the past 50 years. *J Anim Sci* 93: 457–481, 2015. doi:10.2527/jas2014-8221.

22. Etard C, Armant O, Roostalu U, Gourain V, Ferg M, Strähle U. Loss of function of myosin chaperones triggers Hsf1-mediated transcriptional response in skeletal muscle cells. *Genome Biol* 16: 267, 2015. doi:10.1186/s13059-015-0825-8.

23. Farhadian M, Rafat SA, Hasanpur K, Ebrahimie E. Transcriptome signature of the lactation process, identified by meta-analysis of microarray and RNA-Seq data. *BioTechnologia* 99: 153–163, 2018. doi:10.5114/bta.2018.75659.

24. Flynn JM, Meadows E, Fiorotto M, Klein WH. Myogenin regulates exercise capacity and skeletal muscle metabolism in the adult mouse. *PLoS One* 5: e13535, 2010. doi:10.1371/journal.pone.0013535.

25. Golding JP, Calderbank E, Partridge TA, Beauchamp JR. Skeletal muscle stem cells express anti-apoptotic ErbB receptors during activation from quiescence. *Exp Cell Res* 313: 341–356, 2007. doi:10.1016/j.yexcr.2006.10.019.

26. Gotoh T, Albrecht E, Teuscher F, Kawabata K, Sakashita K, Iwamoto H, Wegner J. Differences in muscle and fat accretion in Japanese Black and European cattle. *Meat Sci* 82: 300–308, 2009. doi:10.1016/j.meatsci.2009.01.026.

27. Guo L, Cui H, Zhao G, Liu R, Li Q, Zheng M, Guo Y, Wen J. Intramuscular preadipocytes impede differentiation and promote lipid deposition of muscle satellite cells in chickens. *BMC Genomics* 19: 838, 2018. doi:10.1186/s12864-018-5209-5.

28. Han SW, Jung YK, Lee EJ, Park HR, Kim GW, Jeong JH, Han MS, Choi JY. DICAM inhibits angiogenesis via suppression of AKT and p38 MAP kinase signalling. *Cardiovasc Res* 98: 73–82, 2013. doi:10.1093/cvr/cvt019.

29. He H, Liu X. Characterization of transcriptional complexity during longissimus muscle development in bovines using high-throughput sequencing. *PLoS One* 8: e64356, 2013. doi:10.1371/journal.pone.0064356.

30. Jeong J, Bong J, Kim GD, Joo ST, Lee HJ, Baik M. Transcriptome changes favoring intramuscular fat deposition in the longissimus muscle following castration of bulls. *J Anim Sci* 91: 4692–4704, 2013. doi:10.2527/jas.2012-6089.

31. Jurie C, Picard B, Hocquette JF, Dransfield E, Micol D, Listrat A. Muscle and meat quality characteristics of Holstein and Salers cull cows. *Meat Sci* 77: 459–466, 2007. doi:10.1016/j.meatsci.2007.04.014.

32. Kaspric N, Picard B, Reichstadt M, Tournayre J, Bonnet M. ProteIN-SIDE to Easily Investigate Proteomics Data from Ruminants: Application to Mine Proteome of Adipose and Muscle Tissues in Bovine Foetuses. *PLoS One* 10: e0128086, 2015. doi:10.1371/journal.pone.0128086.

33. Keady SM, Kenny DA, Ohlendieck K, Doyle S, Keane MG, Waters SM. Proteomic profiling of bovine M. longissimus lumborum from Crossbred Aberdeen Angus and Belgian Blue sired steers varying in genetic merit for carcass weight. *J Anim Sci* 91: 654–665, 2013. doi:10.2527/jas.2012-5850.

34. Kersten S. Physiological regulation of lipoprotein lipase. *Biochim Biophys Acta* 1841: 919–933, 2014. doi:10.1016/j.bbalip.2014.03.013.

35. Khan DR, Landry DA, Fournier É, Vigneault C, Blondin P, Sirard M-A. Transcriptome meta-analysis of three follicular compartments and its correlation with ovarian follicle maturity and oocyte developmental competence in cows. *Physiol Genomics* 48: 633–643, 2016. doi:10.1152/physiolgenomics.00050.2016.
36. Komolka K, Ponsuksili S, Albrecht E, Kühn C, Wimmers K, Maak S. Gene expression profile of Musculus longissimus dorsi in bulls of a Charolais × Holstein F2-cross with divergent intramuscular fat content. *Genom Data* 7: 131–133, 2015. doi:10.1016/j.gdata.2015.12.023.

37. Lee SH, Park EW, Cho YM, Kim SK, Lee JH, Jeon JT, Lee CS, Im SK, Oh SJ, Thompson JM, Yoon D. Identification of differentially expressed genes related to intramuscular fat development in the early and late fattening stages of hanwoo steers. *J Biochem Mol Biol* 40: 757–764, 2007. doi:10.5483/BMBRep.2007.40.5.757.

38. Lehnert SA, Reverter A, Byrne KA, Wang Y, Nattrass GS, Hudson NJ, Greenwood PL. Gene expression studies of developing bovine longissimus muscle from two different beef cattle breeds. *BMC Dev Biol* 7: 95, 2007. doi:10.1186/1471-213X-7-95.

39. Lehnert SA, Wang YH, Byrne KA. Development and application of a bovine cDNA microarray for expression profiling of muscle and adipose tissue. *Aust J Exp Agric* 44: 1127–1133, 2004. doi:10.1071/EA03238.

40. Li YH, Liu L. Apolipoprotein E synthesized by adipocyte and apolipoprotein E carried on lipoproteins modulate adipocyte triglyceride content. *Lipids Health Dis* 13: 136, 2014. doi:10.1186/1476-511X-13-136.

41. Listrat A, Lebret B, Louveau I, Astruc T, Bonnet M, Lefaucheur L, Bugeon J. How muscle structure and composition determine meat quality. *ScientificWorldJournal* 2016: 3182746, 2016. doi:10.1155/2016/3182746.

42. Marot G, Foulley J-L, Mayer C-D, Jaffréric F. Moderated effect size and P-value combinations for microarray meta-analyses. *Bioinformatics* 25: 2692–2699, 2009. doi:10.1093/bioinformatics/btp444.

43. Moisá SJ, Shike DW, Graugnard DE, Rodriguez-Zas SL, Everts RE, Lewin HA, Faulkner DB, Berger LL, Loor JJ. Bioinformatics analysis of transcriptome dynamics during growth in angus cattle longissimus muscle. *Bioinform Biol Insights* 7: 253–270, 2013. doi:10.4137/BBI.S12328.

44. Pallafacchina G, Blaauw B, Schiaffino S. Role of satellite cells in muscle growth and maintenance of muscle mass. *Nutr Metab Cardiovasc Dis* 23, Suppl 1: S12–S18, 2013. doi:10.1016/j.numecd.2012.02.002.

45. Parikh H, Carlsson E, Chutkow WA, Johansson LE, Storgaard H, Poulsen P, Saxena R, Ladd C, Schulze PC, Mazzini MJ, Jensen CB, Krook A, Björnholm M, Tornqvist H, Zierath JR, Ridderstråle M, Altshuler D, Lee RT, Vaag A, Groop LC, Mootha VK. TXNIP regulates peripheral glucose metabolism in humans. *PLoS Med* 4: e158, 2007. doi:10.1371/journal.pmed.0040158.

46. Park JH, Kang HJ, Kang SI, Lee JE, Hur J, Ge K, Mueller E, Li H, Lee BC, Lee SB. A multifunctional protein, EWS, is essential for early brown fat lineage determination. *Dev Cell* 26: 393–404, 2013. doi:10.1016/j.devcel.2013.07.002.

47. Picard B, Lefaucheur L, Berri C, Duclos MJ. Muscle fibre ontogenesis in farm animal species. *Reprod Nutr Dev* 42: 415–431, 2002. doi:10.1051/rnd:2002035.

48. Purslow PP, Archile-Contreras AC, Cha MC. Meat Science and Muscle Biology Symposium: manipulating meat tenderness by increasing the turnover of intramuscular connective tissue. *J Anim Sci* 90: 950–959, 2012. doi:10.2527/jas.2011.4448.

49. Qi R, Liu H, Wang Q, Wang J, Yang F, Long D, Huang J. Expressions and Regulatory Effects of P38/ERK/JNK Maps in the Adipogenic Trans-Differentiation of C2C12 Myoblasts. *Cell Physiol Biochem* 44: 2467–2475, 2017. doi:10.1159/000486169.

50. Qin L, Zhang G, Cao Y, Zhang J, Zhao Y, Zhao Z. Microarray Analysis on the Differences of Gene Expression in Longissimus Dorsi Muscle Tissue Between 1 and 24 Months Chinese Red Steppes. *J Anim Vet Adv* 10: 428–436, 2011. doi:10.3923javava.2011.428.436.

51. R Core Team. R: A Language and Environment for Statistical Computing [Online]. R Foundation for Statistical Computing. 2018. https://www.r-project.org/ [9 Apr. 2019].

52. Ramasamy A, Mondry A, Holmes CC, Altman DG. Key issues in conducting a meta-analysis of gene expression microarray datasets. *PLoS Med* 5: e184, 2008. doi:10.1371/journal.pmed.0050184.

53. Rifai N, Gillette MA, Carras SA. Protein biomarker discovery and validation: the long and uncertain path to clinical utility. *Nat Biotechnol* 24: 971–983, 2006. doi:10.1038/nbt1235.

54. Roberts MD, Haun CT, Mobley CB, Mumford PW, Romero MA, Roberson PA, Vann CG, McCarthy JJ. Physiological differences between low versus high skeletal muscle hypertrophic responders to resistance exercise training: Current perspectives and future research directions. *Front Physiol* 9: 834, 2018. doi:10.3389/fphys.2018.00834.

55. Romao JM, Jin W, He M, McAllister T, Guan LL. Elucidation of molecular mechanisms of physiological variations between bovine subcutaneous and visceral fat depots under different nutritional regimes. *PLoS One* 8: e83211, 2013. doi:10.1371/journal.pone.0083211.

56. Rung J, Brazma A. Reuse of public genome-wide gene expression data. *Nat Rev Genet* 14: 89–99, 2013. doi:10.1038/nrg3394.

57. Sadkowski T, Jank M, Oprzadek J, Motyl T. Age-dependent changes in bovine skeletal muscle transcriptomic profile. *J Physiol Pharmacol* 57, Suppl 7: 95–110, 2006.

58. Schreurs NM, Garcia F, Jurie C, Agabriel J, Micol D, Bauchart D, Listrat A, Picard B. Meta-analysis of the effect of animal maturity on muscle characteristics in different muscles, breeds, and sexes of cattle. *J Anim Sci* 86: 2872–2887, 2008. doi:10.2527/jas.2008-0882.

59. Schulz TJ, Huang TL, Tran TT, Zhang H, Townsend KL, Shadrach JL, Cerletti M, McDougall LE, Giorgadze N, Tchkonia T, Schrier D, Falb D, Kirkland JL, Wagers AJ, Tseng YH. Identification of inducible brown adipocyte progenitors residing in skeletal muscle and white fat. *Proc Natl Acad Sci USA* 108: 143–148, 2011. doi:10.1073/pnas.1010929108.

60. Scioli MG, Bielli A, Gentile P, Mazzaglia D, Cervelli V, Orlandi A. The biomolecular basis of adipogenic differentiation of adipose-derived stem cells. *Int J Mol Sci* 15: 6517–6526, 2014. doi:10.3390/ijms15046517.

61. Segalés J, Perdiguero E, Muñoz-Cánoves P. Regulation of muscle stem cell functions: A focus on the p38 MAPK signaling pathway. *Front Cell Dev Biol* 4: 91, 2016. doi:10.3389/fcel.2016.00091.

62. Shang Y, Xu X, Duan X, Guo J, Wang Y, Ren F, He D, Chang Z. Hsp70 and Hsp90 oppositely regulate TGF-β signaling through CHIP/Spb1. *Biochem Biophys Res Commun* 446: 387–392, 2014. doi:10.1016/j.bbrc.2014.02.124.

63. Song Y. Function of membrane-associated proteoglycans in the regulation of satellite cell growth. In: *Growth Factors and Cytokines in Skeletal Muscle Development, Growth, Regeneration and Disease*, edited by White J, Smythe G. Springer International Publishing, p. 61–95, 2016.

64. Tacken PJ, Hofker MH, Havekes LM, van Dijk KW. Living up to a name: the role of the VLDL receptor in lipid metabolism. *Curr Opin Lipidol* 12: 275–279, 2001. doi:10.1097/00041433-200106000-00006.

65. Taga H, Bonnet M, Picard B, Zingaretti MC, Cassar-Malek I, Cinti S, Chilliard Y. Adipocyte metabolism and cellularity are related to differences in adipose tissue maturity between Holstein and Charolais or Blond d'Aquitaine fetuses. *J Anim Sci* 89: 711–721, 2011. doi:10.2527/jas.2010-3234.

66. Tang X, Deng L, Xiong H, Li G, Lin J, Liu S, Xie J, Liu J, Kong F, Tu G, Peng H, Liang S. Expression profile of mitogen-activated protein kinase (MAPK) signaling genes in the skeletal muscle & liver of rat with type 2 diabetes: role in disease pathology [Online]. *Indian J Med Res* 140: 744–755, 2014.

67. Thornton KJ, Kamanga-Sollo E, White ME, Dayton WR. Active G protein-coupled receptors (GPCR), matrix metalloproteinases 2/9 (MMP2/9), heparin-binding epidermal growth factor (hbEGF), epidermal growth factor receptor (EGFR), erbB2, and insulin-like growth factor 1 receptor (IGF-1R) are necessary for trenbolone acetate-induced alterations in protein turnover rate of fused bovine satellite cell cultures. *J Anim Sci* 9: 2332–2343, 2016. doi:10.2527/jas.2015-0178.

68. Tiebel O, Oka K, Robinson K, Sullivan M, Martinez J, Nakamura M, Ishimura-Oka K, Chan L. Mouse very low-density lipoprotein receptor (VLDLR): gene structure, tissue-specific expression and dietary and developmental regulation. *Atherosclerosis* 145: 239–251, 1999. doi:10.1016/S0021-9150(99)00068-4.

69. Tremblay F, Dubois MJ, Marette A. Regulation of GLUT4 traffic and function by insulin and contraction in skeletal muscle. *Front Biosci* 8: d1072–d1084, 2003. doi:10.2741/1137.

70. Tseng GC, Ghosh D, Feingold E. Comprehensive literature review and statistical considerations for microarray meta-analysis. *Nucleic Acids Res* 40: 3785–3799, 2012. doi:10.1093/nar/gkr1265.

71. Vaudel M, Verheggen K, Csordas A, Raeder H, Berven FS, Martens L, Vizcaíno JA, Barsnes H. Exploring the potential of public proteomics data. *Proteomics* 16: 214–225, 2016. doi:10.1002/pmic.201500295.

72. Wang YH, Bower NI, Reverter A, Tan SH, De Jager N, Wang R, McWilliam SM, Cafe LM, Greenwood PL, Lehnert SA. Gene expression patterns during intramuscular fat development in cattle. *J Anim Sci* 87: 119–130, 2009. doi:10.2527/jas.2008-1082.

73. Yang YH, Dudoit S, Luu P, Lin DM, Peng V, Ngai J, Speed TP. Normalization for cDNA microarray data: a robust composite method addressing single and multiple slide systematic variation. *Nucleic Acids Res* 30: e15, 2002. doi:10.1093/nar/30.4.e15.