Differential expression and co-expression gene network analyses reveal molecular mechanisms and candidate biomarkers involved in breast muscle myopathies in chicken

The broiler industry is facing an increasing prevalence of breast myopathies, such as white striping (WS) and wooden breast (WB), and the precise aetiology of these occurrences remains poorly understood. To progress our understanding of the structural changes and molecular pathways involved in these myopathies, a transcriptomic analysis was performed using an 8 × 60 K Agilent chicken microarray and histological study. The study used pectoralis major muscles from three groups: slow-growing animals (n = 8), fast-growing animals visually free from defects (n = 8), or severely affected by both WS and WB (n = 8). In addition, a weighted correlation network analysis was performed to investigate the relationship between modules of co-expressed genes and histological traits. Functional analysis suggested that selection for fast growing and breast meat yield has progressively led to conditions favouring metabolic shifts towards alternative catabolic pathways to produce energy, leading to an adaptive response to oxidative stress and the first signs of inflammatory, regeneration and fibrosis processes. All these processes are intensified in muscles affected by severe myopathies, in which new mechanisms related to cellular defences and remodelling seem also activated. Furthermore, our study opens new perspectives for myopathy diagnosis by highlighting fine histological phenotypes and genes whose expression was strongly correlated with defects.

Fine and quantitative histological traits. To quantify muscular defects, quantitative histological traits were measured in pectoralis major muscle cross sections obtained from the three groups (SG, FG-C and FG-WSWB) ( Table 1). A quantification of fibrosis and adiposis was performed using a collagen VI-bodipy co-labelling technique (Fig. 1A-C). The percentage of the area labelled with collagen VI (representative of fibrosis) on the microscopic field was 4.2 and 6.9 times higher in FG-C and FG-WSWB muscles, respectively, than in SG pectoralis major muscles (p-values ≤ 0.0001). Thus, compared to SG muscles, both FG-C and FG-WSWB muscles showed extended endomysial and perimysial connective tissues, but the rate was 1.6 lower in the FG muscles macroscopically unaffected by WS and WB than in affected muscles. It is interesting to note that the percentage of the area labelled with collagen VI in one sample of the FG-WSWB group and one sample of the FG-C group was very different than the average value obtained in their respective groups, i.e., it was lower for the FG-WSWB  Coefficient of variation of fibre size, % 25.5 ± 2.0 a 37.8 ± 7.9 a 50.5 ± 10.7 b 3.14 × 10 −3 Table 1. Means and standard deviations for histological traits measured in pectoralis major muscle crosssections. * Relative area (%) occupied by collagen VI or adipose tissue. ** Calculated from the Feret's minimum diameter. a,b,c Analysis of variance performed with the group factor for each trait. Different superscripts in each line indicate significant differences between groups, as determined by pairwise comparison (p-values ≤ 0.05).
sample (8.1%) and higher for the FG-C sample (13.6%). This suggests that these samples were either misclassified macroscopically or phenotypically intermediate between the FG-C and the FG-WSWB classes. However, the percentage of the area occupied by bodipy 493 staining (representative of adiposis) was approximately three times larger in FG-WSWB than in FG-C pectoralis major muscles (p-value ≤ 0.01) and was almost non-existent in SG pectoralis major muscles. These results indicate that muscle histological architecture is modified according to the selection on growth rate and that fibrosis and adiposis establishment are exacerbated in cases of severe defects. Capillary density and fibre size were measured using laminin-fibronectin co-immunolabelling ( Fig. 1D-F). The number of capillaries in the microscopic field was 1.5 and 2 times lower in severely affected muscles (FG-WSWB) than in macroscopically unaffected muscles (FG-C) (p-value ≤ 0.05) and normal muscles (SG) (p-value ≤ 0.0001), respectively. In addition, the structure and location of the capillaries appeared to be modified in FG-WSWB muscles as we noticed the presence of lumens in their capillaries, likely because of vasodilatation, and of connective tissue extensions that altered the connection between capillaries and muscle fibres (Fig. 1D-F).
As expected, the number of fibres in a microscopic field was significantly lower in FG than in SG muscles (p-values ≤ 0.0001) due to the difference in fibre size (measured as the minimal Feret's diameter), which was significantly lower in SG than in FG muscles (p-values ≤ 0.0001) ( Table 1). The number of fibres per capillary was approximately 2 to 3 for FG muscles and approximately 6 for SG muscles. However, if we consider the fibre surface supplied by one capillary, it was almost two times larger in FG (20169 µm² for FG-C and 22150 µm² for FG-WSWB) than in SG muscles (12961 µm²) (p-values ≤ 0.05) due to the difference in fibre size and conjunctive tissue extension. Interestingly, neither the number of fibres per capillary nor the averaged fibre size was different between FG-C and FG-WSWB muscles. However, there was significantly more fibre size variation in FG-WSWB muscles than in SG and FG-C muscles (p-values ≤ 0.01), probably due to the strong regeneration process indicated by an increase in frequency of regenerative fibre foci that contain more basophilic fibres with a great proportion of central nuclei. lying pectoralis major muscle susceptibility to WS and WB defects, a custom chicken 8 × 60 K Gene Expression Agilent Microarray was used to perform the transcriptomic analysis of muscles obtained from the SG, FG-C and FG-WSWB groups. Among the 61,657 probes spotted on the array, 43,688 (71%) were expressed, and 36,473 (59%) were annotated (corresponding to 12,272 genes) based on the NCBI EntrezGene database and Galgal5 chicken genome assembly. Based on quantitative histological characterization and a multidimensional scaling plot of gene expression, three samples (i.e., the two samples previously identified in the histological analysis as misclassified into groups according to the macroscopic examination of defects and another sample from the FG-C group) were considered outliers and were not kept in the differential analysis. The analysis of 21 samples was performed using a linear model and pairwise comparisons and revealed that 10,482 of the annotated genes were differentially expressed (DE) genes with a fold-change |FC| ≥ 1.2 in at least one of the three comparisons (adjusted p-value < 0.05). Among these, 6,725 genes were DE between SG and FG-C, 7,887 were DE between FG-C and FG-WSWB, and 9,608 were DE between SG and FG-WSWB. The log2-transformed FC of the DE genes ranged from −3.48 to 6.64 for FG-C vs. SG, from −6.81 to 9.19 for FG-WSWB vs. SG and from −3.87 to 4.22 for FG-WSWB vs. FG-C, indicating the presence of substantial gene expression differences between these groups.
The top 10 DE genes in the three comparisons based on adjusted p-value and log2 FC are presented in Table 2  and Supplementary Table S2, respectively. These genes were mainly involved in the extracellular matrix, collagen fibrils and cytoskeleton organization, muscle contraction, the inflammation response, calcium homeostasis, oxidoreductase activity and angiogenesis.
A hierarchical clustering was applied to the list of 10,482 DE genes with an |FC| ≥ 1.2 to group them into clusters according to their expression profile. A modular height cut-off value in the hierarchical tree allowed the identification of six clusters of DE genes that showed similar profiles of expression within clusters and a clear discrimination among the three groups (Fig. 2a). The number of genes and the expression profile of each cluster is presented in Fig. 2b. Clusters c1 and c4, which represented 12% of the DE genes, discriminated SG from FG muscles (regardless of the presence of defects). In all, 79% of DE genes were found in clusters c2 and c5, which were associated with both growth rate variations and the severity of defects. The expression of the 3,117 DE genes in cluster c5 gradually increase, while that of the 5,168 DE genes belonging to cluster c2 decreased from SG to FG-C and FG-C to FG-WSWB. Among all DE genes, only 9% appeared to be regulated between FG muscles free or affected by WS and WB regardless of the effect of growth. They were expressed at either higher levels (cluster c6) or lower levels (cluster c3) in the FG-WSWB group than in the FG-C and SG groups.

Functional analysis of DE genes.
All expressed genes on the array were annotated by Gene Ontology for Biological Process (BP) based on NCBI EntrezGene chicken IDs (taxonomic id: 9031) and orthologues. Among the 12,272 expressed genes on the array, 6,232 were annotated with at least one BP GO term (51%). The enrichment analysis was performed independently for each DE genes cluster with all expressed genes used as the background (Supplementary Data S1). Two hundred and sixty-seven GO terms were enriched in at least one cluster (p-value < 0.01), and 264 were unique, indicating that clusters of genes were very specific and different from each other ( Table 3). The number of enriched GO terms in each cluster varied between 8 (cluster c3) and 128 (cluster c5).
Co-expression network analysis and module-trait relationships. The histological study identified fine phenotypes to accurately quantify the degree of muscle damage and significantly improve diagnosis, which is currently based on visual notation or muscle palpation after slaughter. Thus, studying the correlations between gene expression levels and histological quantitative phenotypes using gene network construction served as a complementary approach to the differential analysis. Moreover, this approach made it possible to make no a priori conclusions about the groups since it was based on continuous phenotypic variables and thus kept all the samples (n = 24) and took into account the variability in cell disorders highlighted at the histological level in the correlation analysis. From the normalized gene expression matrix, a weighted correlation network analysis (WGCNA) was performed to investigate the relationship between modules of co-expressed genes and histological traits and identify candidate biomarkers of the defects. Only the histological phenotypes allowed us to discriminate muscles severely affected by WS and WB from the two other groups based on, i.e., collagen VI and adipose tissue contents, the number of capillaries per field and the coefficient of variation of fibre size (Table 1) were kept for this analysis. Fifteen co-expression modules containing 50 to 6516 genes were constructed. In Supplementary  Fig. S1, correlations between eigengenes and quantitative phenotypes show that several modules are related to one or more histological traits. We chose to focus on modules that had a mean gene significance (i.e., mean of correlations between gene expression and phenotype, GS) greater than 0.5 ( Supplementary Fig. S2) and absolute value of the module-histological trait correlation higher than 0.7 ( Supplementary Fig. S1). Hence, the light green (55 genes), royal blue (140 genes), brown (6516 genes), black (420 genes), and blue (3393 genes) modules were selected for further analysis. The light green module was strongly positively related to collagen VI content (0.89) and to variation in fibre size (0.85) but negatively related to the number of capillaries (−0.71). Brown and royal blue modules were negatively related to collagen VI content (−0.87 and −0.76, respectively) but strongly positively related to the number of capillaries (0.73 and 0.79, respectively). The blue module was negatively related to collagen VI content (−0.79) and variation in fibre size (−0.74), while the black module was positively correlated with variation in fibre size (0.84).
Genes for which expression was both highly correlated with a histological phenotype (GS ≥ 0.7) and significant modules eigengene (Module Membership, MM ≥ 0.7) were considered hub genes. Hence, 4,894 hub genes were associated with collagen VI labelling, 1,456 with the number of capillaries and 915 with variation in fibre size. www.nature.com/scientificreports www.nature.com/scientificreports/ Identification of candidate biomarkers for severe myopathies. Genes that were DE between FG-C and FG-WSWB muscles (adjusted p-value < 0.05) and hub genes identified in the weighted correlation network analysis were combined in a Venn diagram to identify the most confident biomarkers related to muscle disorders involved in severe WS and WB pathologies (Fig. 3).
Among the 6,418 DE genes between FG-C and FG-WSWB muscles, more than one-third (36%) were not found to be strongly correlated with any of the three histological phenotypes studied. The vast majority of hub genes were DE between FG-WSWB and FG-C muscles. Only 17, 4, and 8% of hub genes that were related to collagen, number of capillaries, and variation of the fibre size, respectively, were not DE between FG-WSWB and FG-C muscles. Among the 7,679 genes included in the Venn diagram, 17 were common among all the lists, 1,611 (21%) between 3 lists, and 2,731 (36%) between two lists. Except for the 17 genes that were common to all the lists, no gene was common between the lists of genes correlated with the number of capillaries and variation in fibre size or between the lists of genes correlated with each of the three histological phenotypes. A total of 548 genes were common between the list of genes related to collagen VI and variation in fibre size, of which 527 were DE between FG-WSWB and FG-C muscles. Finally, 1,336 genes were highly correlated with both collagen VI and the number of capillaries, of which 1,084 were DE between FG-WSWB and FG-C muscles.
Validation using RT-qPCR for a subset of genes. Among the genes found to be DE among the groups, a subset of 12 genes was chosen for validation by RT-qPCR. Genes were selected according to several criteria, including their functional role, their belonging to clusters, and, for some of them, their location within QTL regions recently identified for WS defect in broiler chickens 18 . The expression level of DE genes obtained by RT-qPCR was compared with expression obtained by microarrays. For each tested gene, the DE between groups observed on the microarray was confirmed with the exception of one gene, LRSAM1, the expression of which was only slightly correlated between RT-qPCR and microarray results (0.34). Except for this gene, expression www.nature.com/scientificreports www.nature.com/scientificreports/ measured by RT-qPCR and microarray was highly correlated, as shown by Pearson's correlation coefficients between 0.71 and 0.99 (Fig. 4). Moreover, the gene expression values between groups obtained by qPCR were consistent with the categorization of these genes to different clusters. www.nature.com/scientificreports www.nature.com/scientificreports/ Among the identified genes, the expression of the myosin heavy chain 15 (MYH15), which is involved in the contraction, development, and regeneration of avian skeletal muscles 24 , was dramatically higher (×771) in FG-WSWB than in SG muscles, while the difference between FG-C and SG muscles was much less dramatic (×60). As for MYH15, fibronectin 1 (FN1) expression was the highest in FG-WSWB and the lowest in SG, but the difference between the FG groups was much lower than that for MYH15. Platelet-derived growth factor receptor alpha (PDGFRA), collagen type VI alpha 3-chain (COL6A3) and leucine-rich repeat and sterile alpha motif-containing 1 (LRSAM1) were expressed at significantly highest levels in FG-WSWB than in FG-C and SG muscles. Conversely, the expression levels of adult (MYH1E) and neonatal (MYH1F) myosin heavy chain isoforms were lower in the FG-WSWB group than in the SG and FG-C groups. Other genes varied between the slow-growing and fast-growing groups, regardless of the presence of WS and WB defects. This was the case for caveolin 3 (CAV3) and embryonic myosin heavy chain isoform (MYH1B), whose expression was higher in FG muscles, and for myocardin (MYOCD), calpain 3 (CAPN3) and protein phosphatase-1 regulatory subunit 3 A (PPP1R3A), whose expression was lower in FG muscles than in SG muscles.

Discussion
In recent years, the broiler industry has faced the emergence of serious meat quality defects, such as WS or WB, which affect both the competitiveness and the image of poultry production. Although research has intensified on this subject in recent years, the aetiology of WS and WB defects is still not understood except that it is now well-established that WS and WB conditions relate to an abnormal extent of connective and adipose tissues, which is associated with muscle fibre necrosis and regeneration processes 2,4 . One hypothesis is that these processes are the result of a long-term effect of genetic selection and intensification of production methods that have led to huge progress in terms of broiler growth and breast muscle development, which have also likely reached physiological limits. The comparison of muscles obtained from slow-growing chickens with muscles obtained from fast-growing broilers either affected or not affected by WS and WB defects revealed more than 10,000 DE genes. Part of the differential gene expression is due to change in the relative cell types present in the muscle, e.g. accumulation of adipocytes, fibroblasts and immune cells, and not strictly due to changes within myocytes. Our study revealed however many changes related to muscle fibre metabolism and contractile activity. Altogether, our observations allowed specifying structural and metabolic changes, which have progressively led to the onset of myopathies in fast-growing chicken (Supplementary Data S1).
Long-term effects of selection on muscle structure and metabolism. Clusters c1 and c4 include genes whose expression in muscle has been altered either positively or negatively in response to long-term selection for growth (SG ≠ FG-C) but were not associated with the occurrence of WB and WS syndromes (FG-C = FG-WSWB). Their functions relate to non-specific biological processes and general functions, such as metabolic processes, signal transduction and the immune system (Table 3). Among these genes, those participating in the "positive regulation of glycogen metabolic process" (GO:0070875, Supplementary Data S1) function were expressed at lower levels in FG than in SG muscles, suggesting a decrease in the use of glycogen in fast-growing muscles independent of the presence of WS and WB defects. This finding is consistent with previous studies that clearly demonstrated the negative impact of selection for growth and/or breast meat yield on muscle glycogen content 25,26 . Among the genes likely to be involved in such changes, we identified PPP1R3A, whose expression was much lower in FG than in SG muscles (Fig. 4), as an important actor. PPP1R3A encode a muscle-specific regulatory subunit of protein phosphatase 1, which binds glycogen with high affinity, activates glycogen synthase, and inhibits glycogen phosphorylase. It has been reported that the overexpression of PPP1R3A leads to increased glycogen content 27 , while its KO results in a drastic reduction (90%) in the glycogen content in muscle 28 . Moreover, PPP1R3A was recently identified as the most interesting candidate gene for a major QTL region for ultimate pH (pHu) 29 that reflects muscle glycogen reserves 30 . The depletion of glycogen reserves in the muscle of fast-growing lines could lead to a predisposing metabolic environment favourable to the development of myopathic syndromes. Indeed, many studies have reported that the occurrence of WS and WB  www.nature.com/scientificreports www.nature.com/scientificreports/ defects is often associated with high meat ultimate pH 3,5,17,31-33 or low muscle glycogen content 19,34 . Consistent with these findings, the occurrence and severity of WS were higher in a broiler line selected for high ultimate pH (pHu + ) than in its counterpart selected for low ultimate pH (pHu-) 17 . Moreover, in a recent genetic study, several QTL-controlling WS were identified in the pHu + line but not in the pHu-one 18 , suggesting a possible interaction between the genetic control of WS and the metabolic status of muscle.
In addition to a decrease in glycogen stores, it has been shown that the haem pigment content of muscles also decreases with selection for growth and breast meat yields 25 , probably indicating that there is less muscle vasculature in fast than in slow-growing lines. In the present study, we observed that the expression of MYOCD was 3-fold lower in FG than in SG muscles. MYOCD is involved in extracellular matrix modelling and expressed in smooth muscle cell-containing tissues, in which it may play a crucial role in the differentiation of the smooth muscle cell lineage and vascularization. Interestingly, MYOCD was also recently identified in two different studies of QTL regions that control breast muscle yield 35 and WS occurrence 18 in slow-and fast-growing chicken lines, respectively. In parallel with the decrease in the expression of MYOCD observed in FG muscles (affected or not affected by defects), we report that compared to SG muscles, in FG muscles, there was a strong increase in the fibre surface supplied by one capillary (Table 1). Although this parameter was similar between FG-C and FG-WSWB chickens, we found that the extension of endomysial and perimysial connective tissues observed in the latter moved the capillaries away from the fibre membrane ( Fig. 1D-F). This could be an aggravating factor that could further lower oxygen and nutrient supply to fibres in muscles affected by WS and WB defects.
Altogether, our observations formally demonstrate a recurrent assumption presented in the literature 2,19,26,36 , namely the suggestion that increasing breast muscle weight and yield via selection through hypertrophy of existing fibres has led to a lower capillary density. This may have led to an inadequate supply of oxygen and nutrients to the muscle fibres, which likely promotes the excessive use of glycogen and probably increases oxidative stress in fast-growing muscles.
In the present study, we also observed that CAPN3 is expressed at 3-fold lower levels in FG than in SG muscles (Fig. 4). CAPN3 is an intracellular protease specific to muscle that is involved in the calcium-dependent proteolytic system. A deficiency in this enzyme results in sarcomere disorganization and muscular limb-girdle muscular dystrophy type R1 (LGMDR1 calpain3-related) in humans 37 . In mice, its expression decreases following muscle denervation or denervation-devascularization injury, which induces muscle degeneration and regeneration 38 . Its lower expression in FG muscles is more consistent with the lower amount of vasculature observed in this line, as atrophy induced by denervation was not observed in such muscles. Although the embryonic myosin heavy chain isoform MYH1B is normally no longer expressed in adult muscle, its expression was increased in FG muscles (Fig. 4) and was one of the top 10 DE genes between the FG-C and SG groups (Supplementary Table S2). The expression of CAV3, which plays a key role in the fusion of myoblasts into myotubes during the maturation process of muscle fibres, was also higher in FG muscles than in SG muscles. These observations support the notion that selection for growth and breast muscle yield contributes to the occurrence of degeneration-regeneration processes in FG muscles, independent of the development of WS and WB myopathies.
How has progressive metabolic modification due to selection resulted in myopathic disorders in broiler breast muscle tissues? If not directly related to myopathies, the metabolic and structural changes induced by long-term selection on growth performance and breast muscle yield and evidenced by functional analysis of clusters c1 and c4 have likely created a favourable physiological context for their onset. Clusters c2 and c5 indicate what types of biological processes are involved in the progression of myopathic defects in relation to growth variations. These two clusters accounted for almost 80% of the DE genes (Fig. 2). Several genes belonging to cluster c5, in which gene expression progressively increased from SG to FG-C and from FG-C to FG-WSWB, are involved in the generation of superoxide anion (GO:0042554, Supplementary Data S1), a reactive oxygen species (ROS) issued from the first electronic reduction of oxygen during oxidative phosphorylation 39 . Among these genes, some contribute to the oxidative stress response (ARNT, HIF1A, GPX8, ALDOC, and CRYAB), and these include heat-shock proteins (HSBP1, HSPA5, HSF2, HSPB7, and HSPH1); these genes are known to be expressed at substantially higher levels during oxidative stress 40 . These results strongly support the notion that selection for growth progressively resulted in increased oxidative stress and that this process was exacerbated in fast-growing muscle affected by WS and WB. It is worth noting that elevated ROS levels contribute to human muscular dystrophies, such as Duchenne muscular dystrophy 41 , and that excessive ROS production alters calcium homeostasis via the release of calcium from the sarcoplasmic reticulum 19,36,42 and leads to skeletal muscle necrosis and atrophy 43 . Indeed, several genes belonging to cluster c5 (such as P2RX4, ATP2B1, ANXA6, ATP2A3, SLC8A1, ANXA2, ATP2A2, and SLC24A2) participate in calcium ion transport, indicating the possible release of calcium at the origin of muscle calcium-dependent proteolysis and tissue degeneration 41 , which has been reported in WS and WB muscles.
Cluster c5 includes a great number of genes involved in the production of inflammatory molecules, such as cytokines, chemokines, cell adhesion molecules and tissue-degrading enzymes. They are grouped into several GO terms, including "inflammatory response" (GO:0006954), "positive regulation of cytokine secretion" (GO:0050715), "regulation of cytokine-mediated signalling pathway" (GO:0001959), "regulation of interleukin-1 secretion" (GO:0050704), "positive regulation of interleukin-12 production" (GO:0032735), and "positive regulation of tumour necrosis factor (TNF) production" (GO:0032760) (Supplementary Data S1). TNF is a multifunctional proinflammatory cytokine, and one of its isoforms (i.e., TNFAIP6) was in the top 10 DE genes between the FG-WSWB and FG-C groups ( Table 2, Supplementary Table S2). It induces the expression of metalloproteinases, which are known to degrade several components of the extracellular matrix-cytoskeleton linkage in skeletal muscle and to thereby lead to extracellular alterations and myopathy 43 . Among the 14 matrix metalloproteinases included on the microarray, seven belong to cluster c5 (MMP2, MMP10, MMP14, MMP16, MMP17, MMP24, MMP27, and MMP23B). Consistent with the increase in the inflammatory response found in cluster c5, we observed the activation of the immune system cascade. This included 20% of the genes in cluster c5, with one of the three most significant functions the "innate immune response" (GO:0045087, Supplementary Data S1). We identified several chemokines, cell adhesion molecules and extracellular proteasomes that can favour the infiltration of macrophages, lymphocytes and neutrophils. Such infiltrations, observed in cases of muscle injury, are also frequently observed in muscles affected by WS and WB 1,2,4,14,44 . There are also several cytokines that are produced during the inflammatory response, such as TNF-α or IL-6, that play an important role in muscle regeneration. www.nature.com/scientificreports www.nature.com/scientificreports/ TNF-α attracts muscle satellite cells to the damaged site, promotes satellite cell proliferation by activating NF-κB signalling, and stimulates differentiation by activating the p38 signalling pathway. IL-6 stimulates the migration, proliferation and differentiation of myoblasts 45 .
After injury, muscle regeneration mediated by muscle stem cells is switched on. Muscle stem cells (or satellite cells) are the major contributors to muscle regeneration and remain quiescent under normal conditions 45 . After muscle trauma, they are activated, giving rise to a proliferative myoblast population that differentiates to restore damaged myofibres. The "regulation of skeletal muscle satellite cell proliferation" (GO:0014842, Supplementary Data S1) was significantly enriched in cluster c5. It includes genes involved in the activation of muscle satellite cell proliferation (FGF2, SDC4, GPC1, PAX7, and STAT3), indicating the activation of the muscle repair process. This is consistent with the decrease observed in the expression of the adult and neonatal myosin heavy chain isoforms, MYH1E and MYH1F, which belong to cluster c2, and the increase observed in the expression of MYH15, which is involved in the regeneration of avian skeletal muscles 24 and was one of the top 10 DE genes in the three comparisons (Supplementary Table S2). Moreover, cluster c5 also included several genes that contribute to collagen formation (GO:0032964, GO:0032967, GO:0030199), the GO term "collagen fibril organization" being one of the three most enriched terms of this cluster (Supplementary Data S1). All these observations were consistent with our histological study that showed that there was an increase in both the relative surface occupied by collagen VI across SG, FG-C and FG-WSWB muscles and the coefficient of variation of fibre size for FG-WSWB muscles compared to the other groups (Table 1). Collectively, our molecular and histological observations support the hypothesis that regeneration processes are more activated in fast than in slow-growing muscles, with this phenomenon exacerbated in muscles affected by WS and WB (Fig. 4).
Functional analysis of clusters c2 and c5 also revealed several enriched GO terms related to energetic pathways, especially glycogen metabolism, glucose use, cellular respiration and ATP production (Supplementary Data S1). This finding reveals that there was a progressive decrease (SG > FG-C> FG-WSWB) in the expression of genes encoding glycolytic enzymes (PFKM, ALDOB, TPI1, GAPDH, PGK1, PGAM1, ENO3 and PKM), suggesting a lower level of anaerobic glycolytic activity, which is already observed in muscles affected by WS and WB 11,19,34,46 . This decrease in anaerobic glycolytic activity is concomitant with an alteration in mitochondrial cellular respiration. Indeed, several genes belonging to cluster c2 (SG> FG-C> FG-WSWB) were related to "cellular respiration" (GO:0045333), "aerobic respiration" (GO:0009060), "ATP synthesis coupled electron transport" (GO:0042773), and "mitochondrial respiratory chain complex I and IV assembly" (GO:0032981-GO:0033617) metabolic processes (Supplementary Data S1). These genes are primarily involved in the production of electron transport chain components in mitochondria (CYCS, COX16, SCO1, SURF1, TIMM21, BCS1L, FOXRED2…) and the Krebs cycle (ACO2, PDHB, MDH1, MDH2, PDHA1). Moreover, several genes belonging to cluster c5 (SG <FG-C <FG-WSWB) are involved in carbohydrate metabolic processes, such as gluconeogenesis (FBP1, PCK2, GALE, and ALDOC), carbohydrate derivative transport, nucleotide-sugar biosynthesis, the pentose phosphate pathway (RBKS, RPE, TKT, ALDOC, DERA, FBP1, PGD, and PGM2) and glycogenolysis (GCG, INPP5K). All these results strongly suggested that selection by itself has markedly modified the metabolic pathways used to produce energy in muscles and that these metabolic modifications are even more pronounced in muscles affected by WS and WB. This metabolic reorientation may be related to the progressive decrease in carbohydrate (especially glycogen) availability that has accompanied selection for growth and breast meat yield 25,26 and that is also observed in muscles affected by WS or WB 5,21,34 . It is likely that these metabolic modifications have contributed to the appearance of muscle disorders in fast-growing birds. Indeed, it was recently shown that muscles exhibiting a heightened ability to store glycogen used carbohydrates as the main source of energy, whereas those lacking glycogen used energy produced from amino acid catabolism and lipid oxidation, leading to an adaptive response to oxidative stress and muscle regeneration processes 21,47 . Mechanisms of cellular defence and remodelling are specifically involved in case of severe myopathies. Functional analysis of clusters c3 and c6 highlighted the biological processes specifically activated in fast-growing muscles that are severely affected by WS and WB. Genes belonging to these clusters are regulated between FG-C and FG-WSWB regardless of the effect of growth. They account for 9% of the total DE genes and can be distinguished from that DE genes of clusters 2 and 5 (79% of DE genes) whose expression gradually changes in relation to growth variation and severity of defects. Cluster c6 contained genes that are overexpressed in cases of severe WS and WB. Its most enriched GO term is related to oxygen transport (GO:0015671, Supplementary Data S1) and encompasses genes encoding for oxygen carriers (HBE, HBBA, HBM, HBA1, and  HBBR), likely indicating insufficient oxygen supply in muscles. Another group of functions relates to glycosaminoglycan and proteoglycan synthesis ("protein O-linked glycosylation via threonine" (GO:0018243) and "proteoglycan metabolic process" (GO:0006029)), which are both components of fibrous tissue (Supplementary Data S1). This is fully in agreement with histological observations that showed significant infiltration and extension of connective tissues within muscles (Fig. 1). At the molecular level, we noticed a strong increase in FG-WSWB muscles of that expressed PDGFRA (Fig. 4), a marker of fibro-adipogenic precursors (FAPs), which were recently characterized as playing a role in muscle regeneration and repair processes. Indeed, FAPs differentiate into fibroblasts and adipocytes under pathological conditions, such as muscular dystrophy, leading to connective and fat tissue infiltration 48 . Similarly, the expression of COL6A3 and FN1, which respectively contribute to the formation of collagen VI and are found in most connective tissues, fibres and extracellular matrix (ECM) during tissue repair after degeneration, were expressed at substantially higher levels in FG-WSWB muscles (Fig. 4), in accordance with the high proportion of ECM that characterizes these muscles. FN1 was also one of the top 10 DE genes between FG-WSWB and FG-C (Table 2).
Interestingly, cluster c6 was associated with several GO terms related to neurogenesis, such as "positive regulation of neurogenesis" (GO:0050769), "axonogenesis" (GO:0007409), "neuron differentiation" (GO:0030182), www.nature.com/scientificreports www.nature.com/scientificreports/ "positive regulation of synapse assembly" (GO:0051965), "regulation of neuron differentiation" (GO:0045664), and "neuron projection extension" (GO:1990138) (Supplementary Data S1). Genes belonging to these functions are mostly involved in axon guidance and the formation or maintenance of neuromuscular junctions (such as  AGRN, SLIT2, SEMA3A, EPHB3 and NTRK2). This is likely more a sign of the involvement of neuromuscular junction maintenance and repair than of motor neuron injury, which would result in foci containing atrophic fibres 49 , which were not observed in FG-WSWB muscles. Indeed, because neuromuscular junction structure alterations occur in affected muscles, axons would need to be attracted and redirected to neuromuscular junctions. Many factors that promote myofibre regeneration and repair, such as semaphorin 3A (SEMA3A), which is expressed by satellite cells in injured muscle, also serve as promoters of axonal sprouting and guidance and neuromuscular junction repair as nerves grow towards newly regenerated myofibres 50 . If cluster c6 is associated with the activation of several mechanisms, such as oxygen transport, axon guidance, and neuromuscular junction repair, then in the case of severe myopathy, cluster c3 could be associated with limitations in the development of cellular defence systems. Indeed, it contains genes whose expression is downregulated in muscles affected by WS and WB, with some of these genes participating in apoptosis and cytokinesis, two processes involved in cell defence or survival. This could be a sign that these cells fail to develop effective mechanisms of regeneration, and this defect is probably responsible for the substantial alteration of the structure of muscle and muscle functions observed in muscles severely affected by WS and WB.

Identification of quantitative histological phenotypes and molecular markers for the diagnosis of muscle disorders in chicken.
Our study, like other recent studies 19,20,22,35 , helps to advance our understanding of the biological processes that are involved in the development of myopathies in the breast muscle tissues of modern broiler lines. It also provides several finely defined histological phenotypes, which were used to precisely quantify the degree of muscle damage and significantly improve diagnosis, which is currently based on visual or palpation scoring of muscle after slaughter. The study of these fine histological phenotypes helps researchers to better describe the major structural changes that occur in muscles affected by WS and WB, and we also used them to determine, through a combination of DE analysis and weighted correlation network analysis, the most confident molecular markers related to muscle disorders involved in severe WS and WB pathologies (Fig. 3). Contrary to connective tissue content, fibre size variation and capillary density, no WGCNA module was found to be highly correlated with adipose tissue content. This is probably due to the great variability observed between samples from a same category (see standard deviations in Table 1). It is also possible that changes in fat metabolism could occur earlier during development, as recently shown in Papah et al. 23 , and that its later extension could interfere with that of connective tissue.
Our integrative analysis led to a set of 17 relevant biomarkers that were DE between the FG-C and FG-WSWB groups and highly correlated with histological phenotypes related to myopathy ( Table 4). Three of these genes have not yet been characterized. Among the others, most are involved in muscle regeneration phenomena (PDLIM1), fibrosis (COL6A2), hypoxia (ARNT2, HIF1AN, and PTGS1) and axon guidance (CADM1, SLIT3 and NAV3). These functions are the main biological processes that are specifically activated in fast-growing muscles severely affected by WS and WB defects (see the functional annotation of clusters c3 and c6). We also propose to enrich this set of biomarkers by 5 additional genes (MYH15, FN1, COL6A3, MYH1E, and MYH1F) that were strongly DE between the FG-C and FG-WSWB groups, were highly correlated to two of the three histological phenotypes studied, and were located in the QTL region controlling the WS defect in the breast muscle tissues of chicken 18 (Table 4). Among these five genes, MYH15 was the only one that was highly correlated with the number of capillaries. This slow-type myosin, also called ventricular myosin in chicken because it is mainly expressed in heart, is transiently expressed in embryonic chick skeletal muscles and then re-expressed during muscle regeneration 24 . Interestingly, it is associated with bovine pulmonary hypertension related to high altitude, suggesting that MYH15 could be involved in the response to hypoxic conditions 51 . The ability of this set of 22 genes to discriminate muscles with severe myopathies from non-affected muscles should now be tested in larger populations obtained from different genetic backgrounds. These biomarkers based on gene expression provide fine diagnostic tools to quantify WS and WB defects. Although they require invasive muscle sampling, they can serve as intermediate phenotypes for the search of non-invasive predictive measurements, based on blood test or imaging on live animal for instance.
In conclusion, the present study provides new information on the biological processes that lead to the occurrence of severe myopathy in the breast muscles of modern broiler lines. By comparing the muscle transcriptome of slow-growing to fast-growing broilers that were affected or not affected by myopathies, we attempt to specify molecular pathways that were altered by the long-term selection of growth performance to those directly involved in the development of muscle disorders. Our data suggest that selection for fast growing and breast meat yield has progressively led to inadequate supplies of oxygen and the nutrients that promote carbohydrate depletion in muscle. Such changes would have provided favourable conditions for metabolic modifications that shift cells toward use of alternative catabolic pathways to produce energy, leading to an adaptive response to oxidative stress and the first signs of inflammatory, regenerative and fibrotic processes. All these processes appeared to be dramatically intensified in muscles affected by myopathies, as are new mechanisms related to cellular defence and remodelling to allow the muscles to function (oxygen transport, neurogenesis). Finally, our study highlights the fine histological phenotypes and genes that have expression profile that are highly correlated with muscle disorders involved in the development of myopathies, thus opening new perspectives for a more fine diagnosis and a better screening tool for broilers affected these such defects.

Methods
Animals, sample collection, muscle and meat quality measurements. The present study was carried out on pectoralis major muscles obtained from a fast-and a slow-growing chicken genotype. The experimental slow-growing (SG) chicken line was created by INRA in the seventies' from three original strains: www.nature.com/scientificreports www.nature.com/scientificreports/ Bresse-Blanche, New-Hampshire and White-American breeds. Although this line has been selected for improved body weight, its growth rate is close to that of chicken strains currently used for the French "Label Rouge" (freerange) production 52 . The fast-growing (FG) line is a grand-parental commercial modern line highly selected by Hubbard for growth and breast muscle yield. Birds of the two genotypes were reared following the same standard practices. They were raised separately but received the same diet. All animals were slaughtered at 6 weeks of age at PEAT INRA Poultry Experimental Facility (2018, https://doi.org/10.15454/1.5572326250887292E12). The slaughter and carcass processing conditions of the FG broilers are fully described in Alnahhas et al. 53 . Samples of the pectoralis major muscle were collected fifteen minutes after slaughter and either directly snap-frozen in liquid nitrogen for molecular analysis or stored in isopentane cooled with liquid nitrogen for histological purposes. All muscle samples were then stored at −80 °C until further analysis. One day after slaughter, three trained people visually graded pectoralis major muscles for WS and WB. Both WS and WB were categorized separately on a three-point scale: 0 = absence of defects, 1 = moderate severity of defects, and 2 = high severity of defects. Defects were scored according to the widths of white stripes based on the scale of Kuttappan et al. 6 for WS in addition to visual characteristics and hardness on palpation for WB.
Body composition was characterized for all SG and FG birds through the measurement of pectoralis major yield (PMY) expressed in relation to body weight (BW). Breast meat yield (BMY), abdominal fat yield (AFP), and pectoralis minor yield (PmY) were also measured in FG birds as were several breast meat quality traits that were evaluated in pectoralis major muscle: ultimate pH (pHu) and the colour parameters L*, a*, b*, which were measured 24 hours post-mortem, drip loss (DL) during a 4-day storage period, cooking loss (CL), Warner Bratzler shear force (SF) of cooked meat, curing-cooking yield (CCY) and the thiobarbituric acid-reactive substance index (TBA-RS), used as markers of lipid peroxidation and measured after 8 days of storage. All these measurements were obtained as described in Alnahhas et al. 17 .
Three groups of eight animals were selected for further molecular and histological analyses: eight animals obtained from the SG line and two groups of eight animals selected from among the 200 FG birds according to their extreme visual score for WS and WB defects and DL, CL and CCY phenotypes. The 8 FG-C birds were visually free of any defects. They were selected within their group based on their high water retention capacity (with below-average DL and CL and above-average CCY). The 8 FG-WSWB birds were affected by both the WS and WB defects (average score of 2 for both WS and WB) and have a CCY below the average of their group.
Histochemical traits. Immunohistochemistry was performed on 10 µm-thick pectoralis major muscle cross-sections obtained from the FG-C, FG-WSWB and SG groups. All primary antibodies were obtained from the Developmental Studies Hybridoma Bank, created by the NICHD of the NIH, maintained at The University of Iowa, Department of Biology, Iowa City, IA, and developed by Douglas M. Fambrough, as follows: 39 (anti-collagen VI), B3/D6 (anti-fibronectin), and 31/31-2 (anti-laminin). Secondary antibodies and conjugated  www.nature.com/scientificreports www.nature.com/scientificreports/ streptavidin were obtained from Southern Biotech, Birmingham, AL, USA, as follows: goat anti-mouse IgG1 conjugated to a TR antibody, goat anti-mouse IgG H+L conjugated to a biotin antibody, and streptavidin conjugated to Cy2. The washing procedure consisted of incubating cross-sections for 3 × 5 minutes in phosphate-buffered saline solution. Slides were mounted in Mowiol medium (Sigma, Lyon, France) and stored at 4 °C in a dark chamber before they were observed. Whatever the labelling procedure, six images per sample were selected randomly (captured at 10× and/or 20× magnification using a Leica MC170 camera; Leica Microsystems, Nanterre, France) on a Leica DMRB microscope.
Fibrosis and adiposis quantification. The Collagen VI-Bodipy co-labelling procedure was used to quantify the extension of connective and adipose tissues within muscle. Muscle cross-sections were fixed with 4% paraformaldehyde solution, incubated in 10% goat serum (Sigma, Lyon, France) for 30 minutes, and then incubated with primary Collagen VI antibody (1/200) for 2 hours. After the tissues were washed, muscle cross-sections were incubated with a secondary goat anti-mouse IgG1 conjugated to a TR antibody (1/500) for 1 hour, washed again, and incubated with Bodipy 493 (4,4-Difluoro-1,3,5,7,8-Pentamethyl-4-Bora-3a,4a-Diaza-s-Indacene, 1/3000, Sigma, Lyon, France) for 15 minutes. The relative surface areas labelled with collagen VI and Bodipy were quantified using the Area Fraction Parameter function in ImageJ software after transformation into black and white images (8-bit) and threshold adjustment.
Capillary density and fibre size measurement. A fibronectin-laminin co-labelling procedure was used to label capillary membranes (fibronectin) and basal lamina surrounding muscle fibres (laminin). Muscle cross-sections were incubated in 10% goat serum for 30 minutes and with primary fibronectin antibody (1/100) for 1 hour. After they were washed, cross-sections were incubated with secondary goat anti-mouse IgG H+L conjugated to biotin antibody (1/1000) for 45 minutes, washed, and incubated with streptavidin conjugated to Cy2 (1/500) for 45 minutes. After they were washed, the cross sections were incubated with primary laminin antibody (1/50) for 1 hour, washed, and finally incubated with secondary goat anti-mouse IgG1 conjugated to a TR antibody (1/1000) for 45 minutes. Capillaries were counted manually, and fibre size was determined using the minimum Feret's diameter parameter function of the ImageJ software after transformation into black and white images (8-bit) and threshold adjustment.
Microarray expression profiling. RNA extraction. Total RNA was extracted from pectoralis major muscle samples using the RNeasy ® Mini kit (Qiagen, Valencia, CA, USA). Residual genomic DNA was removed by DNase I treatment (Qiagen, Valencia, CA, USA). RNA concentration was measured using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), and RNA integrity was checked according to an RNA Integrity Number (RIN > 8) using RNA 6000 Nano chips run on a Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA).
Labelling and hybridization. Muscle transcriptional profiling was performed using an 8 × 60 K Agilent Gallus gallus custom-commercial array (Agilent Technologies, Santa Clara, CA, USA) containing 61,657 probes (Platform GPL20588 in the U.S. National Center for Biotechnology Information Gene Expression Omnibus (GEO) microarray database). RNA labelling and microarray processing were performed by the @BRIDGe platform (INRA, UMR GABI, Jouy-en-Josas) using the procedure described in Jacquier et al. 54 . The microarray datasets were submitted to the GEO microarray database with the accession number GSE127806.
Re-annotation of the microarray chip on the current Galgal5 chicken genome. To compare the DE genes identified in this study and the genetic markers of WS provided by our previous study 18 , we re-annotated the chip.
The transcriptome alignment on the Galgal5 chicken assembly (EntrezGene database) of all the probes deposited on the chip was carried out using the blastn algorithm (hybridization tolerance of two mismatches and a single probe-associated gene as parameters), which is available in the BLAST + suite (ncbi-blast-2.6.0+) 55 , on the cluster obtained from the Genotoul bioinfo platform (http://bioinfo.genotoul.fr/). Probes were annotated on a Galgal5 assembly as reported in Supplementary Data S2.
Probes with a hybridization tolerance of two mismatches and a single probe-associated gene were considered correctly annotated and were kept for further analysis.
Differential expression analysis. First-step filtering was applied to keep only those probes that were expressed in this study (43,688 probes). The fluorescence signal for each probe was log2-transformed and normalized by median centring for each array. The differences in expression among the three groups (SG, FG-C, and FG-WSWB) were computed using the log2 transformation of the fold-change (log2FC). The differential analysis was performed using a linear model to identify DE probes among all pairwise comparisons of the three groups implemented in R package limma version 3.32.7 (Linear Model for Microarray Data) 56 . The p-values were adjusted for multiple testing by the Benjamini-Hochberg method 57 to control the False Discovery Rate (FDR). Probes were considered DE when the adjusted p-value was ≤ 0.05 and |FC| ≥ 1.2.
Hierarchical clustering. To identify clusters of DE genes, hierarchical clustering was performed on the gene expression matrix. The correspondence between expressed probes and annotated genes was determined using the re-annotation of the chip. Genes with at least one Galgal5 annotated DE probe, |FC| ≥ 1.2 and the highest mean expression probe level in one of groups were kept. Normalized expression (after log2-transformation) was scaled for each gene. Hierarchical clustering on the scaled gene expression matrix was based on the Pearson correlation distance and an average link aggregation distance. A modular height cut-off value in the hierarchical tree obtained using the cutreeDynamic function of dynamicTreeCut package 58 (with the following parameters: Module detection. A hierarchical clustering of the genes based on the TOM dissimilarity measure followed by a modular height cut-off value of branches in the hierarchical tree using the cutreeDynamic function (deep-Split = 4, minClusterSize = 30) was performed to detect modules of co-expressed genes. The module eigengene, which was the first principal component of each module and represented the expression value of each module, was calculated. Modules with expression profiles that were very similar (correlation = 0.90) were merged since there was a high probability that genes belonging to these modules are highly co-expressed.
Module-trait relationships. Only the histological phenotypes that allowed discrimination of muscles severely affected by WS and WB from other muscles (i.e., collagen VI labelling, adipose tissue, number of capillaries and coefficient of variation of fibre size) were kept for the analysis. The module eigengene was used to detect biologically relevant modules. Indeed, module-trait relationships were estimated using Pearson's correlation between the module eigengene and the trait of interest.
Genes for each module with high Gene Significance (GS ≥ 0.7) corresponding to the absolute value of the correlation between gene expression and the histological trait and high Module Membership (MM ≥ 0.7), which was defined as the absolute value of the correlation of the module eigengene and the gene expression profile, were defined as hub genes.

Validation by RT-qPCR.
To validate the results obtained on the microarray, RT-qPCR was performed to examine the relative expression of 12 genes selected according to several criteria, including their functional role, their belonging to each of the identified clusters, and, for some of them, their location within QTL regions recently identified as controlling the WS defect in broiler chickens 18 . Total RNA was extracted from pectoralis major muscle samples using RNA NOW (Ozyme, St Quentin en Yvelines, France). Ten micrograms of RNA from each sample were reverse-transcribed using RNase H − MMLV reverse transcriptase (Superscript II, Invitrogen, Illkirch, France) and random primers (Promega, Charbonnières les Bains, France). Primers targeting the studied genes were designed with Primer3 version 4.0.0 62,63 . The list of primer sequences is available in Supplementary  Table S3. The products of amplification were analysed by electrophoresis and further sequenced.
Gene expression was quantified by qPCR using a Roche LightCycler ® 480 II (Roche Applied Science) and Takyon ® mix (Eurogentec, Seraing, Belgium) according to the manufacturer's recommendations. Quantitative PCR conditions were set at 95 °C for 5 min, followed by forty-five cycles of 10 s at 95 °C, 20 s at 60 °C and 10 s at 72 °C. 18S ribosomal RNA was used as a housekeeping gene to normalize the CT values as its expression was invariant between groups. A mix of chicken pectoralis major muscle cDNA served as a reference sample (control cDNA). The calculation of absolute mRNA levels was based on the PCR efficiency and the threshold cycle (CT) deviation of an unknown cDNA versus the control cDNA according to the equation proposed by Pfaffl 64 . Ethics approval. Experimentation was conducted in PEAT INRA Poultry Experimental Facility (2018, https://doi.org/10.15454/1.5572326250887292E12) and authorized for animal experimentation by the French veterinary service, which is the competent authority (D-37-175-1). The scientist in charge of the experimentation received training and personal authorization . In agreement with the Ethical Committee for Animal Experimentation in Val de Loire (C2EA-19), the experiments reported here did not require approval from a specific ethical committee since only classical rearing practices for broilers were used. No permissions were required to collect the chicken lines owned by INRA.

Data availability
The microarray data were submitted to the Gene Expression Omnibus (GEO) microarray database (accession number GSE127806).