Transcriptional Profiles of Skeletal Muscle Associated With Increasing Severity of White Striping in Commercial Broilers

Development of the white striping (WS) abnormality adversely impacts overall quality of broiler breast meat. Its etiology remains unclear. This study aimed at exploring transcriptional profiles of broiler skeletal muscles exhibiting different WS severity to elucidate molecular mechanisms underlying the development and progression of WS. Total RNA was isolated from pectoralis major of male 7-week-old Ross 308 broilers. The samples were classified as mild (n = 6), moderate (n = 6), or severe (n = 4), based on number and thickness of the white striations on the meat surface. The transcriptome was profiled using a chicken gene expression microarray with one-color hybridization technique. Gene expression patterns of each WS severity level were compared against each other; hence, there were three comparisons: moderate vs. mild (C1), severe vs. moderate (C2), and severe vs. mild (C3). Differentially expressed genes (DEGs) were identified using the combined criteria of false discovery rate ≤ 0.05 and absolute fold change ≥1.2. Differential expression of 91, 136, and 294 transcripts were identified in C1, C2, and C3, respectively. There were no DEGs in common among the three comparisons. Based on pathway analysis, the enriched pathways of C1 were related with impaired homeostasis of macronutrients and small biochemical molecules with disrupted Ca2+-related pathways. Decreased abundance of the period circadian regulator suggested the shifted circadian phase when moderate WS developed. The enriched pathways uniquely obtained in C2 were RNA degradation, Ras signaling, cellular senescence, axon guidance, and salivary secretion. The DEGs identified in those pathways might play crucial roles in regulating cellular ion balances and cell-cycle arrest. In C3, the pathways responsible for phosphatidylinositol 3-kinase-Akt signaling, p53 activation, apoptosis, and hypoxia-induced processes were modified. Additionally, pathways associated with a variety of diseases with the DEGs involved in regulation of [Ca2+], collagen formation, microtubule-based motor, and immune response were identified. Eight pathways were common to all three comparisons (i.e., calcium signaling, Ras-associated protein 1 signaling, ubiquitin-mediated proteolysis, vascular smooth muscle contraction, oxytocin signaling, and pathway in cancer). The current findings support the role of intracellular ion imbalance, particularly Ca2+, oxidative stress, and impaired programmed cell death on WS progression.


INTRODUCTION
White striping (WS) abnormality, as depicted by the appearance of white lines running parallel to muscle fibers on the surface of chicken breast meat, has become of great concern in the poultry industries. Based on their aberrant visual characteristics, the breasts with severe WS condition have reduced appeal to consumers and are downgraded for use as processed meat products, thus resulting in significant economic losses for the poultry industry (Kuttappan et al., 2012(Kuttappan et al., , 2016Petracci et al., 2014). The WS breasts consistently exhibit technological challenges, particularly inferior water holding capacity (Kuttappan et al., 2013;Petracci et al., 2013). The nutritional profile of the WS meat is shifted toward higher fat and lower protein proportions compared to non-WS chicken breasts (Petracci et al., 2014;Zambonelli et al., 2016;Malila et al., 2018).
The specific etiology of the WS abnormality is still under investigation. In order to elucidate the molecular mechanisms underlying development of WS myopathy, Zambonelli et al. (2016) compared the gene expression patterns of pectoralis major muscle from Ross 708 broiler hybrids exhibiting severe WS together with wooden breast (WB) against those without any macroscopic lesions. The differentially expressed genes (DEGs) identified were associated with muscle development, inflammation, polysaccharide metabolism, and calcium signaling. It is worth noting that, in the study of Zambonelli et al. (2016), the abnormal group was affected with both WS and WB conditions. Therefore, the biological pathways directly related with WS development alone could be diluted by the strong effect of WB myopathy and might not be clearly shown in the results. Transcriptome profiles of severe WS (white line thickness >1 mm) in pectoralis major collected from 42-day-old Cobb 500 broilers showed significant increases in transcripts involved in the activation of hypoxia, immune responses, and angiogenesis compared with that of nondefective samples (Marchesi et al., 2019). Dysregulated metabolic processes were also identified at the metabolomic level (Boerboom et al., 2018) with opposite directionalities of some intermediates of tricarboxylic acid cycle as well as incomplete breakdown products of protein catabolism. Recently, using droplet digital polymerase chain reaction (ddPCR), an increased absolute transcript abundance of hypoxia-inducible factor 1 alpha (HIF1A) along with differential expression of HIF1A-related genes were identified in pectoralis major of the WS-and WS+WB-affected Ross 308 broilers (Malila et al., 2019). So far, the growing evidence is pointing toward oxidative stress as a primary contributor to development of WS muscle.
While the previous studies focused on exploring different transcriptomic or metabolomic profiles of severe WS in respect to nondefective muscles, the molecular aspects underlying the progression in severity of WS lesion still remain unclear. The objective of this study was to examine differences in transcriptome profiles and define molecular mechanisms associated with elevating WS severity in pectoralis major collected from 7-week-old male commercial Ross 308 broilers. Our findings regarding the progression of the WS myopathy at molecular level provides complementary information to existing studies concerning the development of WS in commercial broilers. The current findings could lay foundation in the development of an intervention to prevent or delay the progression of WS abnormality.

Samples and Sample Collection
All breast samples were collected within 20 min postmortem from 83 broiler carcasses (7-week-old male Ross 308), raised at a facility of a local poultry meat producer and processed in the industrial processing plant under the routine Halal standard practice. One side of each breast was immediately dissected, snapfrozen, and kept in liquid nitrogen during transportation back to the laboratory (Pathum Thani, Thailand). The muscle samples were stored at -80 • C until total RNA isolation. The other side of the breast was dissected, placed in a plastic bag, kept on ice during transportation, and stored at 4 • C until use.
In this study, all samples were purchased in the form of whole carcasses from the commercial processing plant. Neither experimental treatments or scientific procedures were subjected to the living animals. Thereby, according to BIOTEC Institutional Animal Care and Use Committee, the ethical approval was not required.

White Striping Defect Classification
The WS abnormality and its severity were inspected using the criteria previously described (Malila et al., 2018). Three WS severity levels (Figure 1), including "mild WS" (1-40 white lines with the thickness of ≤0.5 mm), "moderate WS" [more than 40 white lines or 1-5 line(s) with the thickness of 1-1.9 mm], and "severe WS" (more than 5 lines with the thickness of 1-1.9 mm thickness or at least 1 white line with thickness >2 mm), were observed among the samples. The samples showing other abnormalities, e.g., WB, were excluded from this experiment in order to avoid any confounding interpretations from such myopathies. Accordingly, 29 mild, 38 moderate, and 4 severe WS samples were available for this study. Six biological replicates were randomly selected from "mild" and "moderate" groups. For FIGURE 1 | Classification of white striping (WS) abnormality by degree of severity. Mild, 1-40 white lines with the thickness of ≤0.5 mm; moderate, more than 40 white lines or 1-5 line(s) with the thickness of 1-1.9 mm; and severe, more than 5 lines with the thickness of 1-1.9 mm thickness or at least 1 white line with thickness >2 mm. Black arrows, black triangles, and red triangles indicate white lines with thickness <1.0 mm, 1.0-1.9 mm, and ≥2.0 mm, respectively.
"severe" samples, all four samples were used. Although our study was initially designed to include the analysis of non-WS samples, there were no nondefective samples (without any macroscopic trait associated with the development of WS) found among our total of 83 samples of the 7-week-old broilers collected for this trail. This phenomenon emphasized the high prevalence of WS abnormality in 7-week-old commercial broilers.

Total RNA Isolation
Total RNA was isolated using TRIzol TM Reagent (Invitrogen), subsequently treated with DNase I (Thermo Scientific, Inc.), and repurified using GeneJET RNA Cleanup and Concentration Micro Kit (Thermo Scientific, Inc.). The samples with an RNA integrity number exceeding 7.0 were analyzed by microarray hybridization at Molecular Genomics Pte. Ltd. (Singapore) and quantitative real-time polymerase chain reaction (qPCR).

Microarray Platform and Experimental Design
The Agilent SurePrint G3 Custom GE 8 × 60K chicken gene expression microarray (Agilent Technologies, Inc.) used in this study was designed based on the National Center for Biotechnology Information (NCBI) Gallus gallus Annotation Release 103. The detailed information of the array is available at NCBI Gene Expression Omnibus (GEO) with the platform accession GPL24307.
To determine differential gene expression associated with WS severity, expression patterns obtained among the WS levels were compared against one another (Figure 2A). Overall, 16 arrays were hybridized in order to accomplish three comparisons (C1, moderate against mild; C2, severe against moderate; C3, severe against mild).

Microarray Hybridizations
Total RNA was labeled using one-color low-input Quick Amp labeling kit (Agilent) following the company's instruction. In brief, 100 ng total RNA was reverse transcribed into double-stranded complementary DNA (cDNA) by priming with oligo(dT) primer. The synthesized cDNA was subjected to an in vitro transcription using T7 RNA polymerases to produce cyanine 3-CTP-labeled complementary RNA (cRNA). After purification, 600 ng of the purified cRNA was hybridized onto the Agilent SurePrint array at 65 • C for 17 h. The array was washed, gently blotted dry, and scanned using an Agilent High-Resolution Microarray Scanner (C Model, Agilent). The TIFF image was saved and analyzed using Agilent Feature Extraction Software version V10.7.1.1 (Agilent).

Microarray Data Analysis and Gene Annotation
Raw microarray signal values were generated by Agilent Feature Extraction Software version V10.7.1.1 (Agilent). Feature extraction and probe quality control were processed using GeneSpring software. The data were further analyzed in R version 3.4.0. The signal values of probes deemed to be suspicious or faulty were quarantined using flags and expression values. Flags were categorical indicators including "detected, " "compromised, " and "not detected" from the scanner. Probes that were flagged as "not detected" and "compromised" were assigned to not available values (NA). Raw signal intensities of probes below 20 for that microarray were assigned as NA values. Probes with missing intensity data were discarded. The filtered raw data were normalized in R using quantile normalization (Bolstad, 2017) and using combat normalization afterward (Müller et al., 2016;Leek et al., 2017). Log 2 transformation was applied to normalized data for statistical analysis. Statistically significant difference of each transcript between two treatments was identified by analysis of independent two-group t test. The DEGs were identified using combined criteria of false discovery rate (FDR) ≤ 0.05 and absolute fold change (|FC|) ≥ 1.2. Positive and negative FC values represent increased or decreased expression of a particular gene in the samples exhibiting a greater severity level relative to its counterparts. All microarray data have been submitted to the NCBI GEO repository with GEO accession number GSE107362.

Confirmation of Differential Gene Expression
Twenty-one genes, resulting in a total of 63 comparison counterparts, were chosen for confirming the microarray data. Of the 63 comparison counterparts, 11 and 18 showed decreased and increased expressions, respectively, with |FC| ≥ 1.2, whereas the rest of the genes exhibited |FC| < 1.2. Primers (Supplementary Table S1) were designed using Primer-BLAST 1 .
Changes in expression patterns of the selected genes observed in the current microarray analysis were confirmed using qPCR (Malila et al., 2015). Threshold cycle (Ct) was analyzed using BioRad CFX Manager 2.1 software (BioRad). Hypoxanthine-guanine phosphoribosyltransferase (HPRT) was chosen as the reference gene due to its unchanged expression across WS groups. Relative messenger RNA (mRNA) abundance of each gene in the sample group relative to its designated counterparts was calculated based on 2 − Ct method (Livak and Schmittgen, 2001). To determine statistical difference between two treatments, Ct was used in calculation using Student's t test. The significance was set at p value of 0.05. All relative abundances obtained from qPCR were plotted against log 2 FC analyzed from microarrays.

Biological Function and Pathway Analysis of Differentially Expressed Genes
All DEGs of each comparison were grouped into known networks and pathways using Ingenuity Pathway Analysis R (IPA) (trial version, Qiagen) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (Kanehisa et al., 2012). Mapping takes the identifiers from the DEGs and compares them to the database. In IPA, the selected DEGs were mapped to the Ingenuity Knowledge Base database and were then ranked by score, which is calculated based on the hypergeometric distribution and Fisher's exact test. The significance was set at a p value of 0.05. The enrichment of genes based on KEGG pathways was the number that has highest number of genes mapped on the pathway.

Differentially Expressed Genes Associated With Increasing WS Severity
The current microarray analysis revealed 91, 136, and 294 differentially expressed genes in C1, C2, and C3 comparisons, respectively (Figures 2B,C). The complete DEGs lists for each comparison are available in Supplementary Table S2. No unique genes were common to all three comparisons. However, C1 shared 8 and 10 DEGs with C2 and C3, respectively; for C2 and C3, there were 29 common DEGs (Supplementary Table S3). Approximately 75% of the DEGs in all comparisons showed |FC| ranging between 1.2 and 1.5 ( Figure 2D), except for the C2 of which all 33 downregulated DEGs fell within such |FC| range. Numbers of the up-and downregulated transcripts were comparable within each comparison. The most up-and downregulated genes in C1 were chromosome 3 open reading frame (C1ORF95, FC = 2.2) and nucleolin-like (LOC107057290, FC = -3.3), respectively. Mitochondrial translational factor 2 (MTIF2, FC = 2.3) and poly(A)-specific ribonuclease regulatory subunit (PAN3, FC = -9.6) were the most changed DEGs in C2. For C3, general transcription factor IIF subunit 1-like (LOC107051030, FC = 2.1) and neuronal pentraxin 2 (NPTX2, FC = -2.5) were the DEGs showing the largest changes.
The changes in gene expression patterns from the present microarray data were further validated using qPCR (Figure 3). Of the 63 comparison counterparts, the changes in expression of 45 counterparts accounted for 71% of the qPCR data, exhibiting similar trends as that acquired from the microarray platform. As expected, the |FC| values obtained from qPCR of some genes were greater than those reported in microarray analysis. This is due mainly to better sensitivity of qPCR over microarray (Camarillo et al., 2011).

Functional and Pathway Analysis
In order to review the enriched biological functions, the DEGs were analyzed using a trial version of IPA. Based on the IPA knowledge base, 50 (55%), 89 (65%), and 191 (65%) DEGs of C1, C2, and C3, respectively, were recognized and mapped into the respective canonical pathways. Considering the top altered physiological functions (Table 1), the increment of WS severity from mild to moderate (C1) and from moderate to severe (C2) similarly impacted organ morphology and development as well as nervous system function. On the other hand, in C3, the large abundance of DEGs involved cardiovascular system and connective tissue. For cellular functions (Table 2), cell-to-cell FIGURE 3 | Confirmation of differential gene expression. Scatter plot illustrates the fold changes of 21 selected genes determined by microarray and quantitative real-time PCR (qPCR). Of the 63 comparisons, 45 counterparts (accounted for 71.4%) show comparable results between microarray and qPCR.  signaling and interaction was the overlapping category among the three comparisons. The increasing WS severity from mild to moderate (C1) mainly affected amino acid, carbohydrate, and lipid metabolisms. Although the altered cellular functions in C2 were somewhat similar with those in C1, the DEGs involved in cellular morphology and cellular assembly and organization were also enriched. As for C3, the large enrichment of DEGs was related with cellular compromise and the changes in cellular growth, proliferation, and cell cycle. All DEGs in the top altered cellular functions are listed in Supplementary Table S3.
To visualize the functional interactions among the DEGs, all DEGs of each comparison were subsequently submitted to KEGG pathway database. The top 30 enriched pathways annotated by KEGG are depicted in

DISCUSSION
The severity of the WS abnormality varies among the meat depending on the age, the gender, and the size of the chickens as well as depending on their management, e.g., feeding regime. In our previous investigation (Malila et al., 2018), the incidence of the WS abnormality was as high as 98% of the samples with broad spectrum of severity. The majority of the abnormal meat samples were categorized as mild (55.7% of 183 collected broilers) which showed little to no effect on both nutritional and technological properties of the meat. In contrast, severe WS was detected at a small prevalence of 3.3% but significantly impact the quality of the resulting meat (Malila et al., 2018). Once downgraded and used for manufacturing processed meat products, the severe WS breasts, even in small proportions, might alter the quality and consistency of an entire batch of processed products (Petracci and Cavini, 2012). In this study, we profiled transcriptomes within the pectoralis major of 7week-old broilers exhibiting three different WS levels (i.e., mild, moderate, and severe) with the aim of understanding how changes in the transcriptome are associated with increasing WS severity. Since we did not obtain any nondefective breast samples from our sampling sites, the experiment was modified to compare each pair of defective based on their WS severity (Figure 2A). Table S3), the eight shared DEG responses between C1 and C2 exhibited comparable fold changes but in opposite directions. The results indicated that differential gene expression patterns were associated with the progression of the myopathic conditions from mild to moderate, followed by a return to non-significant levels similar to those of the mild samples when the WS severity increased from moderate to a severe degree. These expression patterns suggest the crucial roles of those DEGs in development of moderate WS as well as the potent metabolic switch when the myopathy progressed from moderate to severe level. Two out of the eight shared DEGs, the transcripts of dimethyladenosine transferase 1 (DIMT1) and zinc finger protein 251 (LOC101750329), are involved in RNA processing. The other two transcripts, the plateletspecific isoform of phosphofructokinase (PFKP) and phytanoyl-CoA dioxygenase domain-containing protein 1 (LOC100858951), play essential roles in carbohydrate and lipid metabolisms, respectively. The gene CALM encodes the Ca 2+ -binding protein calmodulin that mediates a wide range of cellular response to the changes in intracellular [Ca 2+ ] (Hamilton et al., 2000). The relevant roles in maintenance of cellular structure of mitochondrial ribosome-associated GTPase 2 (mtg2) and spastin (SPAST) have been documented (Datta et al., 2005;Roll-Mecak and Vale, 2008). For eva-1 homolog B (EVA1B), although the specific function of this member of EVA1 family is still scarce, the roles of EVA1 in mediating muscle cell extensions were previously demonstrated in the nematode Caenorhabditis elegans (Chan et al., 2014).

Regarding the shared DEGs between each comparison pair (Supplementary
Based on the fold changes of the 10 DEGs shared between C1 and C3, the expressions of those DEGs were altered when the WS severity progressed from mild to the other severity degrees. Among those 10 shared DEGs, increased expressions of prolyl 4-hydroxylase, transmembrane (P4HTM), and sperm autoantigenic protein 17 (SPA17) were of interest. The P4HTM encodes an endoplasmic reticulum-anchored prolyl 4-hydroxylase possessing a transmembrane domain. Unlike the other prolyl 4-hydroxylases that act on collagens, this specific isoform catalyzes hydroxylation of the critical prolines in HIF-1α but not in type I procollagen, suggesting the crucial participation of this isoform in adaptive response under reduced oxygen availability (Koivunen et al., 2007). Transcription of this gene was induced in human cell lines exposed to hypoxic condition (Koivunen et al., 2007), supporting the greater extent of hypoxia within the moderate and severe WS breasts in comparison to mild WS. As for SPA17, the potential pathogenic role of this transcript, but not the protein, in highly proliferating and neoplastic cells has been addressed (De Jong et al., 2002). The increased SPA17 mRNA abundance in moderate and severe WS might reflect increased cell proliferation within those samples. Considering the DEGs shared between C2 and C3, similar fold changes between those two comparisons indicated the upregulation or, in the case of neuronal pentraxin 2 (NPTX2), downregulation only in severe WS samples. Among those 29 shared DEGs between C2 and C3, six transcripts were annotated as noncoding RNAs (ncRNA), whereas the other two were annotated as miscellaneous RNAs (miscRNA). The regulatory roles of ncRNAs in oxidative stress response and development of cardiac hypertrophy and apoptosis, the biological pathways frequently identified to be associated with myopathies in broilers, were also addressed (Dong et al., 2019;Turunen et al., 2019). In addition, Zhang et al. (2017) identified a set of ncRNAs involved in regulation of intramuscular preadipocyte differentiation in breast muscle of Chinese indigenous chickens, which was similar to the deposition of fat on connective tissue of the WS breasts. Still, specific biological functions of those ncRNAs associated with WS development remain to be elucidated. Overall, the list of the shared DEGs between the comparisons suggests the key genes partly responsible for each WS severity degree.

Comparison Between Mild and Moderate WS Degrees (C1)
The small numbers of DEGs in C1 relative to the other comparisons suggest slight differences at the cellular level between mild and moderate WS degree. Considering the unique enriched biological pathways altered specifically for C1 (Table 3), the majority of the pathways are linked with phosphatidylinositol phospholipase C (PLCB) and CALM. The PLCB-encoding enzyme catalyzes the hydrolysis of phosphatidylinositol 4,5-bisphosphate into two secondary messengers, inositol triphosphate (IP 3 ) and diacylglycerol (DAG), which regulate various transduction cascades, particularly homeostasis of macronutrients and intracellular [Ca 2+ ] (Kadamar and Ross, 2013). The CALM, the second most differentially expressed gene in C1, plays roles in a vast diversity of cellular functions (Berchtold and Villalobo, 2014). In this study, an increased PLCB expression coupled with reduced CALM abundance identified in moderate WS compared with those of mild WS appeared to alter metabolisms of macronutrients as well as the biological processes, including circadian clock, dopaminergic synapse, long-term potentiation, gonadotrophin-releasing hormone (GnRH) secretion, melanogenesis, and glucagon signaling ( Table 3). Alteration of lipid and amino acid metabolisms could be partly supported by an increase in hepatic aminoadipate aminotransferase (AADAT), currently highlighted in the pathway of antibiotic biosynthesis. Upregulation of AADAT, reported in liver of fattening pigs, was speculated to promote triglyceride biosynthesis from degraded lysine (Kojima et al., 2018).
The association between the dysregulated intracellular [Ca 2+ ] and development of moderate WS was emphasized by the 1.4fold decreased chemokine (C motif) ligand (XCL1), encoding lymphotactin as shown in the chemokine signaling pathway ( Table 3). Similar to other chemokines, lymphotactin recruits T and natural killer (NK) cells and induces intracellular Ca 2+ mobilization by binding a specific G-protein-coupled receptor (Volkman et al., 2009).
Among the pathways relevant to dopaminergic synapse, sphingolipid signaling, and Chagas disease, transcript abundance of the beta isoform of regulatory subunit 55 of protein phosphatase 2A (PPP2R2B) was decreased 1.2-fold in moderate WS relative to mild samples. Decreased expression of PPP2R2B has been linked with the promotion of deregulated oncogene signaling pathways required for survival and proliferation of cancer cells (Tan et al., 2010), suggesting the different cell proliferation rates between mild and moderate WS.
Regarding the circadian entrainment, the changes in the clock gene networks not only associated with the disrupted immunity and metabolic disorders (Turek et al., 2005) but also the structure within the muscle (Harfmann et al., 2015). In this study, apart from PLCB and CALM, an increased  (Aggarwal et al., 2017). Although the expression of PER3 was upregulated in the C1 comparison in the current study, differential expression of this gene likely affects adipogenesis in these birds and may be partly responsible for the increased fat content within the moderate WS breasts. As the light is the dominant signal for the biological clock-controlled system, the findings imply that the lighting management might partly influence the entanglement of biological cascades, including circadian clock and metabolic functions during the progression of WS from mild to moderate level. Further investigation is required to elucidate such association. Overall, the current transcriptome profiling consistently suggested the association between the progression of mild to moderate WS and the dysregulated metabolic processes of macronutrients along with the shifts of Ca 2+ homeostasis. The impacts on amino acids, lipids, and carbohydrate were in accordance with the changes in glycolytic enzymes (Kuttappan et al., 2017) and metabolites from tricarboxylic acid cycle (Boerboom et al., 2018) previously reported in pectoralis major of WS + WB samples, and the differential expression of various genes related to carbohydrate and fatty acid metabolisms identified in WS-affected Cobb 500 broilers (Marchesi et al., 2019).

Comparison Between Moderate and Severe WS Degrees (C2)
Similar to C1, differential gene expression identified in C2 in which the WS abnormality progressed from moderate to severe levels were related to metabolisms and physiological activities within the abnormal broilers through the signals of phosphatidylinositol, cyclic guanosine monophosphate (cGMP)dependent protein kinase, and apelin (Figure 4). Consistent with C1, an altered intracellular ion homeostasis, particularly Ca 2+ , was observed in C2. An increased CALM abundance in severe WS relative to moderate WS might exert relevant roles in several unique enriched biological pathways, including Ras signaling, cellular senescence, and axon guidance ( Table 4). In addition to changes in [Ca 2+ ], a difference in [K + ] between severe and moderate WS might be hypothesized as a consequence of the 1.5-fold increase in transcriptional level of the calciumactivated potassium channel subunit alpha-1 (KCNMA1). In salivary secretion pathway, the product of KCNMA1 functions as a component of the voltage-gated potassium channel, which is activated by either an increased intracellular [Ca 2+ ] or changes in membrane electrical potential (Miller, 2000).
Unlike C1, the modules containing sequential activities regulating cell proliferation, differentiation, and death, particularly through the pathways of mitogen-activated protein kinase (MAPK) necroptosis, RNA degradation, and cellular senescence (Figure 4), were impacted in C2. In addition, RNA degradation and cellular senescence were among the five unique enriched pathways identified for C2 (Table 4). Considering the complete list of the unique enriched pathways identified in C2 ( Table 4), most of them transmit signals in governing cell proliferation and survival. In Ras signaling, Ras functions as a molecular switch for signaling pathways through a cycle of guanidine triphosphates (GTP)-and guanidine diphosphate (GDP)-bound states (Milburn et al., 1990). One of the mechanisms by which GTP-Ras operates includes the recruitment of Raf-1 proteins, encoded by RAF1, to form Raf-1 heterodimers, which in turn stimulate various cellular signaling cascades that control cell proliferation, differentiation, and migration (Zenonos and Kyprianou, 2013). As for the axon guidance, its signal cues non-canonically regulate cell proliferation, cell migration, and gene expression (Russell and Bashaw, 2018). It is also interesting that several DEGs in C2 exhibited consistent fold change with the expression required for proliferation and survival of various cancer cell types as addressed in previous studies. For instance, suppressed Ras guanyl-releasing protein 3 (RASGRP3) was previously identified in transgenic hypertrophic mouse hearts of which malignant mutations of cardiac myosin regulatory light chain (Huang et al., 2016). In the study of Ji et al. (2011), the abundance of ephrin type-B receptor 3 (EPHB3) was positively correlated with severe pathological characteristics, particularly tumor size, in lung cancer cells. Overexpression of KCNMA1 promoted tumor progression in mouse model (Khaitan and Ningaraj, 2019). Taken together, the underlined molecular shifts when the WS abnormality elevating from moderate to severe lesion was more likely related with dysregulated cell proliferation.
In addition to the potential abnormality in cell proliferation, the transcriptional change specific for C2 exerted more deleterious impact than those observed in C1. The speculation was supported by the two unique enriched pathways, including RNA degradation and cellular senescence. The RNA degradation is required for monitoring RNA species in living cells. Transcriptional changes of a regulatory subunit of poly(A)nuclease (PAN3, decreased by 9.6-fold) and scavenging decapping enzyme (DCPS, increased by 1.2-fold) suggest the fluctuations in RNA degradation pathway to a greater extent in severe WS in compared to moderately affected cases. Even small fluctuations are adequate to initiate cell-to-cell variation and interfere with the precision of cell division cycle (Baudrimont et al., 2019). Knockdown of PFKP stimulates survival of breast cancer cells (Kim et al., 2017). In addition to PFKP, a delayed poly(A) tail degradation associated with reduced PAN3 protein level was observed in HeLa cells subjected to oxidative stress (Yamagishi et al., 2014). The findings underline the more pronounced deviation of cell cycle and accelerated cell proliferation with the increment of WS from moderate to severe level. As for cellular senescence, this pathway is an irreversible cell-cycle arrest induced by a number of endogenous and exogeneous stresses arising from oncogene, DNA damage, and oxidative stress (Zhao et al., 2016). Accumulation of the senescent cells interferes with the functions of stem cells, attenuates muscle regeneration, and promotes an imbalance between bone formation and bone resorption (Baar et al., 2018). Sun et al. (2017) reported an increase in senescent biomarkers in association with isoproterenol-induced pathological cardiac hypertrophy, which also exhibited fibrosis similar to those found in severe WS muscle (Kuttappan et al., 2013;Malila et al., 2019).
The findings in the current microarray study may reflect the uncontrolled proliferation of muscle cells and the disruption of cell cycle within the hypertrophic breast affected with severe WS in comparison to moderate level.

Comparison Between Mild and Severe WS Degrees (C3)
Microarray data of C3 revealed the greatest numbers of DEGs due potentially to the profound distinctive phenotypic characteristics between the severe and mild groups. Approximately twothirds of the DEGs in C3 exhibited positive fold changes. Differences, at gene expression level, in cardiovascular system, connective tissue (Table 1), cell cycle, and cellular compromise ( Table 2) are highlighted in between severe and mild cases. The enriched biological pathways specific for C3 (Table 5) could reflect the transcriptional alterations of signaling hubs, including phosphatidylinositol 3-kinase (PI3K)/Akt pathway and p53 cascade, and the downstream outputs of such stress-induced adaptive mechanisms.
Based on the DEGs enlisted in the PI3K/Akt signaling, it is interesting that several DEGs, including angiopoietin 1-receptor tyrosine kinase (TEK), integrin subunit alpha 3 (ITGA3), membrane associated guanylate kinase with inverted orientation (MAGI1), and nitric oxide synthase 3 (NOS3), are found mainly in endothelial cells (Laura et al., 2002;Zhen et al., 2008). The TEK-encoding enzyme regulates vascular morphogenesis and homeostasis (Augustin et al., 2009). Targeted deletion of TEK contributed to postnatal degeneration of newly formed veins and development of hemangioma-like vascular tufts in mouse retina (Chu et al., 2016). Upregulated ITGA3 has been linked with a progression of renal failure (Steenhard et al., 2012) and pathological features of liver cancer (Huang et al., 2018). Depletion of MAGI1 induced upregulation of endoplasmic reticulum (ER) stress-related genes, which in turn signal apoptosis in the endothelial cells (Abe et al., 2019). As demonstrated in vitro and in vivo, NOS3 was highly expressed upon the exposure to reactive oxygen species (Zhen et al., 2008). An increase in synthesis of nitric oxide to increase blood flow due to the limited vascularization within WS/WB affected muscle has been speculated . As shown in the study of Christov et al. (2007), endothelial cells stimulate myogenic cell (MPCs) growth, whereas the MPCs showed proangiogenic effect upon differentiation. The altered expression of those DEGs supported the hypothesis regarding the imbalanced interplay between the endothelial cells and MPCs within the hypertrophic myofibers in breast muscle of broilers exhibiting severe WS. Still, further studies are necessary to obtain a better understanding underlying specific functions of those DEGs in broiler skeletal muscle. One of the downstream processes from PI3K-Akt or p53 pathways is the onset of apoptosis with the aim to remove any damaged cells via caspase activity (Chang et al., 2003). Transcript abundance of caspase-8 (CASP8), the principal upstream apoptotic caspase effector of deathinducing receptor signaling (Kruidering and Evan, 2000), was decreased by 1.2-fold in severe relative to mild WS samples. A reduced mRNA level of CASP8 was recently reported in breast tumors collected from the subjects diagnosed with breast cancer (Aghababazadeh et al., 2017). Apart from CASP8, the decreases of caspase 8 and Fas-associated death domain-like apoptosis regulator (CFLAR) along with the increased cathepsin D (CTSD) abundance were identified in severe WS compared to mild samples. The reduced CFLAR expression accelerates inflammation, ER stress, and cell death in ischemia-induced mouse brain tissues , whereas overexpression of CTSD was linked with various diseases, e.g., type 2 diabetes (Liu et al., 2017), acute coronary syndromes (Gonçalves et al., 2016), and malignant epithelial tumor . Additionally, the extensive alteration of apoptosis in severe WS is emphasized by the increased abundance of dynactin subunit 1 (DCTN1) coupled with decreased axonemal heavy chain 3 of dynein (DNAH3), highlighted in Huntington disease pathway ( Table 5). Differential expression of those two genes inhibits apoptotic membrane trafficking during apoptosis (Lane et al., 2001). Overall, the modulated programmed cell death could be anticipated in severe WS.
The importance of altered intracellular [Ca 2+ ] in the development of WS was again underlined by the DEGs identified in C3. The Ca 2+ -/calmodulin-dependent protein kinase type II delta chain (CAMK2D) is part of a complex that regulates ER-mediated Ca 2+ influx in either skeletal or cardiac muscles. The Ca 2+ voltage-gated channel subunit alpha 1D (CACNA1D) is specifically part of the channel responsible for conferring low-voltage activation and slowly inactivating Ca 2+ currents. The Ca 2+ release-activated calcium channel protein 1 (ORAI1), a Ca 2+ -selective channel particularly in T cells, mediates Ca 2+ influx following the depletion of intracellular Ca 2+ stores (Gwozdz et al., 2012). Overexpression of CAMK2D and CACNA1D induces cardiac enlargement, with mild fibrosis resembling the severe WS lesion, in transgenic mice (Zhang et al., 2003). In addition, a decreased expression of the aminobutyrate aminotransferase (ABAT), highlighted in GABAergic signaling, is associated with increased intracellular [Ca 2+ ], leading to an increased Ca 2+ -sensing nuclear translocation associated with aggressive breast cancer . The changes in transcription of those DEGs likely lead to the excessive Ca 2+ influx in severe WS samples (Kuttappan et al., 2013;Zambonelli et al., 2016;Marchesi et al., 2019).
Increased abundances of LOC101747382, LOC107055590, and LOC107055658, encoding type I collagen, coupled with upregulated NOS3 and ORAI1, were highlighted under the pathway of platelet activation. The pathway is responsible primarily for coagulation cascades upon the interaction between platelet membrane glycoproteins and exposed collagens and secondarily for intercellular communication that mediates inflammatory and immunomodulatory activities, particularly with the triggers from increasing intracellular [Ca 2+ ] and extracellular matrix (Yun et al., 2016). Generally, the ratio between types 1 and 3 determines mechanical properties of collagen matrix, hence the tissues. An increase in proportion of collagen type 1 has been demonstrated in healing wounds, hypertensive arteries, and fibrotic tissues (Friedman et al., 1993), which was in accordance with fibrotic characteristics of severe WS previously described (Trocino et al., 2015).
Regarding the arrhythmogenic right ventricular cardiomyopathy, the pathological characteristics of this heart muscle disease includes progressive myocyte loss and fibrofatty replacement, which were resembled to the histological lesions of severe WS. Under this enriched pathway, increased transcript abundance of desmin (DES) was significant in severe WS in comparison to that of mild WS. This observation was consistent with the previous report on the induced DES at the protein levels in WS/WB breast muscle samples compared with nondefective ones (Zambonelli et al., 2016;Kuttappan et al., 2017;Soglia et al., 2020). Based on the functions of desmin in structural organization of myofibers during myogenesis, the greater DES in severe WS suggested the occurrence of extensive muscle regeneration within the severe WS in comparison with those of found in mild WS samples.
The majority of the DEGs enlisted in the unique pathways of C3 were also associated with growth and proliferation of cancer cells as observed in C1 and C2. Apart from those DEGs mentioned above, transcript abundance of plekstrin homology domain leucine-rich repeat-containing protein phosphatase 2 (PHLPP2), a negative regulator of PI3K/Akt pathway, was decreased by 1.2-fold. A loss of PHLPP2 is associated with progression of pancreatic (Smith et al., 2016) and prostatic cancer (Nowak et al., 2019). Differential expressions of cyclin G (CCNG) and increased RING finger and WD repeat domain protein 2 (RFWD2), the two proteins involved in negative autoregulatory feedback loops of the p53 response observed in this study, were concomitant with the modified transcriptional direction addressed in gastric cancer tissues and in ovarian and breast cancer tissues, respectively (Dornan et al., 2004;Gao et al., 2018). Decreased expression of the gene encoding B-Raf (BRAF) by 1.2fold in severe compared to mild WS might imply the impaired proliferative arrest in the severe WS.
HIF-1 signaling is also highlighted in C3. The pathway, consistently identified in the WS/WB myopathies (Malila et al., 2019;Marchesi et al., 2019;Pampouille et al., 2019), stimulates an adaptive mechanism via an activation of transcription factor HIF-1 when cells sense low oxygen availability. Despite no difference in expression of HIF-1, the shifts of TEK and NOS3 expression suggested modified angiogenesis and antioxidant activities under hypoxia within the severe WS compared with those in mild WS samples. The current results correspond well with our previous report (Malila et al., 2019) in which an increased HIF1A was only detected when absolute transcription level of WS and WS + WB muscle were compared with that of nondefective ones; however, no change was observed among the defective samples. This might explain the absence of HIF1A in the DEGs obtained in this report.
Collectively, the changes in gene expression patterns combined with the unique biological pathways in C3 were in the same direction related to development and proliferation of tumor and cancer in several tissue types. The molecular resemblance underlined the dysregulation of damaged cells as well as the aberrant cell proliferation potentially during muscle regeneration within the muscle exhibiting severe WS myopathy.

Altered Biological Pathways in Common Among All Comparisons
Most of the overlapping enriched pathways comprised the same set of the DEGs discussed in each comparison, suggesting the importance of those DEGs in the progression of WS condition. The shared altered pathways emphasize the impact of dysregulated intracellular ion homeostasis and modulated cellular cascades that remove damaged cells in relevance with WS abnormality. In this study, transcript abundances of PLCB, CACNA1D, CAMK2, and ORAI1 gradually increased during WS progression. The findings were consistent with previous studies of Zambonelli et al. (2016) and Marchesi et al. (2019) identifying upregulation of those genes in WS-related muscle compared with those of nondefective ones. Increases in expression of those genes have been related with elevated intracellular [Ca 2+ ]. As Ca 2+ is one of the crucial secondary messengers associated with phosphoinositol system, [Ca 2+ ] disturbance may affect both the molecular functions and the viability of the cell. Elevated intracellular [Ca 2+ ] modulates apoptosis and protein degradation, compromises mitochondrial function, and interferes with cell membrane integrity (Florea et al., 2005;Mutryn et al., 2015;Petracci et al., 2015). This could explain the extensive muscle protein degradation observed in WS samples collected from Cobb 500 broilers in the previous report of Vignale et al. (2017). Likewise, NOS3 was increased as WS severity increased. This result supported the speculation of Petracci et al. (2019), suggesting the elevated hypoxic condition as WS severity elevated. In addition, the gene encoding solute carrier family 25 (SLC25A4) slightly decreased as the myopathy progressed. Expression of this gene promoted oxidative stress and apoptosis (Lena et al., 2010). A progressive decrease in PHKB was detected, which agreed well with the previous findings of Marchesi et al. (2019), suggesting an unusual muscle contraction process within the affected muscle.
Consistent with the observation in C2 and C3, the numbers of the DEGs obtained in the overlapped pathways were changed in the direction that facilitated proliferation and progression of cancer cells and tumors (Figure 6). Apart from the DEGs discussed earlier, the DNA mismatch repair protein, encoded by MSH6, gradually decreased while the WS severity increased. Previously, suppressed MSH6 was detected under an experimental hypoxic condition (Koshiji et al., 2005). The shifts of glutathione (GST)-encoding gene among the WS severity levels, corresponding with our recent ddPCR results (Malila et al., 2019), implied responsive mechanism of the muscle cells against increased oxidative stress. In addition to impaired apoptosis, ubiquitination-mediated proteolysis (Supplementary Figure S1) was disrupted. An accumulation of the unrepaired oxidative damaged biomolecules could induce a cascade of molecular events that ultimately contribute to promotion of cancer or tumor growth (Conde-Pérezprina et al., 2012). The current microarray results strengthen the case for oxidative stress-induced damage within the WS breast that induce a cascade of molecular events contributing to aberrant cell proliferation.
Although the actual etiology regarding development of WS abnormality remains under investigation, the current transcriptome analyses ascertained an association between disturbances of calcium homeostasis during progression of WS severity. In complement to the previous report, our results underlined the major shifts of biological processes while the abnormality progressed. When WS abnormality advanced from mild to moderate WS, the microarray data revealed transcriptional changes of the key upstream mediators of Ca 2+ . The impacts of imbalanced intracellular [Ca 2+ ] coupled with alterations of the other secondary messengers may exert an overall metabolic shift and inflammation. As the WS myopathy progresses to the severe degree, a more profound perturbation of intracellular [Ca 2+ ] along with oxidative stress could be anticipated. One of the major origins of the oxidative stress might potentially be the limited capillarization. The differences in metabolisms of macronutrients and small biochemical molecules might be shadowed by the extensive results of impaired cellcycle arrest and cell proliferation. With the unrepaired DNA mismatches along with modulated programmed cell death, accumulation of abnormal genomic materials and intermediate metabolites could trigger dysregulation of physiological functions and excessive cell proliferation. The aberrant muscle cells underwent necrosis as supported by the previous histological reports ( Kuttappan et al., 2013;Malila et al., 2018;Marchesi et al., 2019) and ultimately developed severe WS characteristics.
In conclusion, the current transcriptome analysis strengthens the association between development of WS abnormality in commercial broilers and aberrant molecular mechanisms responsible for controlling ion homeostasis and muscle repair process. The transcriptional alteration, when WS progressed from mild to moderate, also influenced metabolisms of macronutrients and small biomolecules. In the development of severe WS lesion, besides the potential leak of the ions, the changes in expression of several genes resemble those found in highly proliferating cells. In addition to the link with an intensive selection for fast-growing and enlarged muscle mass, this study has shown an alteration of the circadian clock system at transcriptional level. The current findings offer evidence that both the intensive, selective breeding and overall management of the modern broiler may play roles in development and severe progression of this myopathy. Optimizing rather than maximizing growth performance of broilers should be prioritized to ensure security of this inexpensive protein source and the sustainability of poultry industry.

ETHICS STATEMENT
All samples were purchased in the form of whole carcasses from the commercial processing plant. Neither experimental treatments nor scientific procedures were subjected to the living animals. Thereby, according to BIOTEC Institutional Animal Care and Use Committee, the ethical approval was not exempted.