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

,


INTRODUCTION
The management of muscle mass and/or of meat-eating qualities requires a good understanding of the molecular drivers 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 variations affect sensorial and nutritional meat-eating qualities (41).This has motivated the age-dependent high-throughput molecular 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 consensus on the molecular signature of muscle fiber growth relative to connective and adipose tissues (15,18,21,63).Metaanalysis (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.

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, hypermuscled, hypermuscularity, double-muscled, 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.

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 the from 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

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 metaanalysis 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.cutreefunction 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 pheatmap 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 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.
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.

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).

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 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.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).

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 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, ATRAID, 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 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 meateating 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.62DoM.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 (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 mitogenactivated 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.6fold 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 metaanalysis 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.revised manuscript; J.B., F.J., P.D., M.R., C.K., D.L., and M.B. approved final version of manuscript.

Fig. 1 .
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 Table1.n, Number of data sets/publications; GEO, Gene Expression Omnibus; WoS, Web of Science.

Fig. 3 .
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.

Table 1 .
Descriptive data of GSE used