Transcriptome Response of Differentiating Muscle Satellite Cells to Thermal Challenge in Commercial Turkey

Early muscle development involves the proliferation and differentiation of stem cells (satellite cells, SCs) in the mesoderm to form multinucleated myotubes that mature into muscle fibers and fiber bundles. Proliferation of SCs increases the number of cells available for muscle formation while simultaneously maintaining a population of cells for future response. Differentiation dramatically changes properties of the SCs and environmental stressors can have long lasting effects on muscle growth and physiology. This study was designed to characterize transcriptional changes induced in turkey SCs undergoing differentiation under thermal challenge. Satellite cells from the pectoralis major (p. major) muscle of 1-wk old commercial fast-growing birds (Nicholas turkey, NCT) and from a slower-growing research line (Randombred Control Line 2, RBC2) were proliferated for 72 h at 38 °C and then differentiated for 48 h at 33 °C (cold), 43 °C (hot) or 38 °C (control). Gene expression among thermal treatments and between turkey lines was examined by RNAseq to detect significant differentially expressed genes (DEGs). Cold treatment resulted in significant gene expression changes in the SCs from both turkey lines, with the primary effect being down regulation of the DEGs with overrepresentation of genes involved in regulation of skeletal muscle tissue regeneration and sarcomere organization. Heat stress increased expression of genes reported to regulate myoblast differentiation and survival and to promote cell adhesion particularly in the NCT line. Results suggest that growth selection in turkeys has altered the developmental potential of SCs in commercial birds to increase hypertrophic potential of the p. major muscle and sarcomere assembly. The biology of SCs may account for the distinctly different outcomes in response to thermal challenge on breast muscle growth, development, and structure of the turkey.


Introduction
Early muscle development is a multifaceted process that involves the proliferation of undifferentiated stem cells in the mesoderm (myoblasts) followed by differentiation and fusion to form multinucleated myotubes. Both early and subsequent muscle growth and repair is dependent upon a population of muscle stem cells (satellite cells, SCs). In the mammalian neonate, SCs comprise~30% of muscle micronuclei but this number declines to only 1-5% in adult skeletal muscle [1,2]. Based on study of chickens, a similar progression of SCs is thought to exist in birds [3,4]. Satellite cells are also a heterogeneous population, with activation, proliferation and differentiation potentials that change with age [5]. For example, in turkeys maximum SC activity occurs immediately after hatch; proliferation declines after wk 4 and differentiation significantly decreases as early as wk 4 post hatch. Cell culture studies suggest that SCs are also multi-potential as they can be induced to follow osteogenic or adipogenic cellular pathways in addition to myogenesis [6][7][8].
Satellite cells are modulated by the cellular microenvironment and their active state is defined by antagonistic Notch and Wnt signaling, with Notch maintaining expression was removed, cells were rinsed with PBS, collected into RNAzol RT (Sigma Aldrich) and held at −80 • C until RNA isolation.

RNA Isolation and Sequencing
Total RNA extraction, sample QC, library preparation were as described in Reed et al. [26]. Briefly, Total RNA was isolated by RNAzol RT extraction, DNase-treated (Turbo DNA-freeTM Kit, Ambion, Inc., Austin, TX, USA), and stored at −80 • C. Concentration and quality of RNAs were assessed by spectrophotometry (Nanodrop 1000). Each sample was then quantified by RiboGreen Assay (Invitrogen Corp., Waltham, MA, USA) and RNA integrity confirmed on the 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). For each treatment group, replicate samples were sequenced (n = 2 replicates per treatment group). Library preparation and sequencing was conducted at the University of Minnesota Genomics Center on the NovaSeq SP platform using v1.5 chemistry (Illumina, Inc., San Diego, CA, USA) to produce 51-bp paired-end reads.

RNAseq Data Analyses
Sequence processing and mapping were as described in Reed et al. [26] using the turkey genome (UMD 5.1, ENSEMBL Annotation 104). Principal component analysis (PCA) was performed in CLC Genomics Workbench (CLCGWB v. 10.1, Qiagen, Redwood City, CA, USA) on normalized read counts. Gene IDs were determined as described in Reed et al. [26]. Differential gene expression and ANOVA (Bonferroni and False Discovery Rate [FDR] corrected) was performed in CLCGWB using standard work flow and with the RBC2 cells set as control for between-line comparisons and the 38 • C treatment groups set as control for temperature/treatment comparisons (Wald test). Pair-wise comparisons between treatment groups were made in the Bioconductor (3.2) R package DESeq2 [29]. For all comparisons, p-values less than 0.05 were considered significant. Affected genes and pathways were investigated with Ingenuity Pathway Analysis [30], Gene Ontology (GO analysis) and PANTHER Overrepresentation Tests (GO Ontology database doi:10.5281/zenodo.6799722 Released 1 July 2022, [31]) http://geneontology.org/ (accessed 8 October 2022).

Results
Total RNA was used to construct individual barcoded libraries for the SC cultures (2 lines × 3 treatments × 2 replicates = 12 total). Sequencing produced over 273M pairedend reads (Table 1)   For each library the total number of raw paired-end reads, median read qualities (R1 and R2), the number of observed genes (mapped reads > 1) by library and treatment group and the percentage of expressed genes (group read count > 1.0) are given.

Gene Expression
Evidence for expression (at least one mapped read in at least one treatment group) was observed for 15,771 annotated genes representing 87.8% of the ENSEMBL turkey gene set (~90% of protein-coding genes) (Supplementary Table S1). The number of observed genes per library ranged from 13,537 to 14,000 (average 13,752.0), and the number of observed genes per treatment ranged from 13,653 to 13,949 (~80% of gene set) ( Table 1). Distribution of expressed protein-coding genes (shared and unique) in treatment group comparisons is summarized in Table 2. The total number of shared genes (expressed in all treatments) was 12,487 with the number of unique genes ranging from 367 to 651 in within-group comparisons. Uniquely expressed genes were higher in the RBC2 cells at 33 • C and 38 • C, but greater in the NCT SCs at 43 • C. For each comparison of the treatment groups (Temperature: 33 • C, 38 • C or 43 • C, Line: RBC2 = R or NCT = N), the total number of expressed and uniquely expressed genes (genes expressed uniquely in each group comparison), the number of genes with significant False Discovery Rate (FDR) p-value, and the numbers of significant genes also with |Log 2 FC (fold change)| > 1.0 and > 2.0 are given. Only those genes with treatment group mean normalized read counts > 1.0 are counted as expressed.
Normalized read counts were used to conduct a principal component analysis (PCA) to partition variation among treatment groups ( Figure 1). Within the first two principal components, treatment groups clustered distinctly irrespective of line and replicate treatment pairs occurred as nearest neighbors within the PCA space.
Normalized read counts were used to conduct a principal component analysis (PCA) to partition variation among treatment groups ( Figure 1). Within the first two principal components, treatment groups clustered distinctly irrespective of line and replicate treatment pairs occurred as nearest neighbors within the PCA space. Several genes that serve as markers of the muscle differentiation cascade were included in the RNAseq data. Contrasting expression in differentiating cells with that of proliferating cells [26] found significant up regulation for MYOD1, PAX7, MYOG (myogenin), DLL1 (delta like canonical Notch ligand 1), and MSTN (myostatin) ( Table 3). Significantly down-regulated genes included MYF6 (myogenic factor 6), PPARγ (peroxisome proliferator activated receptor γ), and LIF (LIF interleukin 6 family cytokine). Table 3. Differential expression of skeletal muscle marker genes in SCs following 48 h of differentiation at 38 °C compared to 72 h of proliferation [26]. For each gene, the directional change in gene expression between differentiating and proliferating cells, the fold change (Log2FC) and False Discovery Rate (FDR) p-value are given for the comparisons within each turkey line (RBC2 and NCT). Directional expression of genes is indicated by red (up) or green (down).

Differential Expression
Seven two-way gene expression contrasts were made based on temperature treatment (cold and hot) and turkey line (RBC2 and NCT). Temperature effects were examined in four pairwise, within-line comparisons: cold-versus control (33 • C vs. 38 • C) and hotversus control (43 • C vs. 38 • C) (Supplementary Table S2). On average, more genes were significantly affected (FDR p-value < 0.05) by heat (43 • C) treatment than by cold (33 • C) although a greater number of DEGs (|Log 2 FC| > 2.0) were seen in the NCT cells with cold treatment and a greater number in the RBC2 cells with heat. A large portion of DEGs identified in the treatment comparisons (36-56%) were unique to treatment groups (temperature/line) (Supplementary Figure S1).

Effect of Cold Treatment
A total of 575 DEGs were identified in the 33 •  significantly affected (FDR p-value < 0.05) by heat (43 °C) treatment than by cold (33 °C) although a greater number of DEGs (|Log2FC| > 2.0) were seen in the NCT cells with cold treatment and a greater number in the RBC2 cells with heat. A large portion of DEGs identified in the treatment comparisons (36-56%) were unique to treatment groups (temperature/line) (Supplementary Figure S1).

Effect of Cold Treatment
A total of 575 DEGs were identified in the 33 °C vs. 38   In the RBC2 SCs 244 DEGs were identified (|Log2FC| > 2.0) with 181 (74.2%) being down regulated by cold treatment. The 40 DEGs with the greatest expression change in each comparison are presented in Figure 3. DEGs with the greatest down regulation included ENSMGAG00000018824 (a novel protein coding gene with BLAST homology to Proline-Rich Protein 33, PRR33), ENSMGAG00000001884 (a myosin heavy chain homolog), ADCY8 (adenylate cyclase 8), and TRABD2B (TraB domain containing 2B). In humans, PRR33 is predicted to act within the response to wounding, ADCY8 catalyzes the formation of cyclic AMP from ATP [32], and TRABD2B is involved in several processes, including negative regulation of Wnt signaling pathway [33]. DEGs with the greatest up regulation included IGSF21 (immunoglobin superfamily member 21), SOX10 (SRY-box transcription factor 10), and LOC100545586, (microsomal triglyceride transfer protein large subunit-like). As a transcription factor SOX10 is thought to be involved in the regulation of embryonic development and in the determination of the cell fate [34]. Microsomal triglyceride transfer proteins are involved in lipoprotein assembly [35]. In the RBC2 SCs 244 DEGs were identified (|Log 2 FC| > 2.0) with 181 (74.2%) being down regulated by cold treatment. The 40 DEGs with the greatest expression change in each comparison are presented in Figure 3. DEGs with the greatest down regulation included ENSMGAG00000018824 (a novel protein coding gene with BLAST homology to Proline-Rich Protein 33, PRR33), ENSMGAG00000001884 (a myosin heavy chain homolog), ADCY8 (adenylate cyclase 8), and TRABD2B (TraB domain containing 2B). In humans, PRR33 is predicted to act within the response to wounding, ADCY8 catalyzes the formation of cyclic AMP from ATP [32], and TRABD2B is involved in several processes, including negative regulation of Wnt signaling pathway [33]. DEGs with the greatest up regulation included IGSF21 (immunoglobin superfamily member 21), SOX10 (SRY-box transcription factor 10), and LOC100545586, (microsomal triglyceride transfer protein large subunit-like). As a transcription factor SOX10 is thought to be involved in the regulation of embryonic development and in the determination of the cell fate [34]. Microsomal triglyceride transfer proteins are involved in lipoprotein assembly [35].
In the NCT SCs, 456 DEGs were identified (|Log 2 FC| > 2.0) with 313 (68.6%) being down regulated by cold treatment (Figure 2). DEGs with the greatest down regulation included ENSMGAG00000018824 and ENSMGAG00000001884 shared with the RBC2 SCs ( Figure 3, Supplementary Table S2). Down regulated DEGs included myosin heavy chain genes LOC100544198, LOC100544198, and MYH1E (average Log 2 FC = −6.8) and the myosin light chain gene MYL3 (Log 2 FC = −6.64). MYL3 was also significantly down regulated in the RBC2 SCs as was stathmin 4 (STMN4), a gene thought to enable tubulin binding activity. Highly up-regulated genes include two (ADAMTS3 and ENSMGAG00000021208) with predicted function in collagen formation ( Figure 3). ADAMTS3 is a protease that may play a role in the processing of type II fibrillar collagen [36]. ENSMGAG00000021208, annotated as a novel protein coding has BLAST similarity to translocation associated membrane protein 2 (TRAM2). TRAM2 is predicted to couple the activity of the ER Ca 2+ pump SERCA2B to increase local Ca 2+ concentration at the site of collagen synthesis [37].
In the comparison of cold treated cells, 125 DEGs were shared between the RBC2 and NCT SCs. Of these all responded similarly with 113 being down regulated and 12 up regulated by heat treatment in both lines (Supplementary Table S2). Comparison analysis of the cold-treated RBC2 and NCT DEGs in IPA identified Oxidative phosphorylation as the canonical pathway with the greatest average |z|-score (−7.21 and −5.10 in NCT and RBC2 cells, respectively, Figure 4). Directional activation was generally consistent between lines, except for Cholesterol biosynthesis, where positive z-scores occurred in NCT and negative scores in RBC2.
genes LOC100544198, LOC100544198, and MYH1E (average Log2FC = −6.8) and the myosin light chain gene MYL3 (Log2FC = −6.64). MYL3 was also significantly down regulated in the RBC2 SCs as was stathmin 4 (STMN4), a gene thought to enable tubulin binding activity. Highly up-regulated genes include two (ADAMTS3 and ENSMGAG00000021208) with predicted function in collagen formation ( Figure 3). ADAMTS3 is a protease that may play a role in the processing of type II fibrillar collagen [36]. ENSMGAG00000021208, annotated as a novel protein coding has BLAST similarity to translocation associated membrane protein 2 (TRAM2). TRAM2 is predicted to couple the activity of the ER Ca 2+ pump SERCA2B to increase local Ca 2+ concentration at the site of collagen synthesis [37].
In the comparison of cold treated cells, 125 DEGs were shared between the RBC2 and NCT SCs. Of these all responded similarly with 113 being down regulated and 12 up regulated by heat treatment in both lines (Supplementary Table S2). Comparison analysis of the cold-treated RBC2 and NCT DEGs in IPA identified Oxidative phosphorylation as the canonical pathway with the greatest average |z|-score (−7.21 and −5.10 in NCT and RBC2 cells, respectively, Figure 4). Directional activation was generally consistent between lines, except for Cholesterol biosynthesis, where positive z-scores occurred in NCT and negative scores in RBC2.  common to the two genetic lines. In contrast to the cold treated cells, the number of DEGs up regulated by heat was greater in both lines than those down regulated.
In the RBC2 SCs, 525 DEGs were identified (|Log 2 FC| > 2.0) with 328 (62.5%) being up regulated. The 40 significant DE genes with the greatest expression change in each comparison are presented in Figure 5. DEGs with the greatest down regulation in the RBC2 SCs included ND4L (NADH dehydrogenase subunit 4 L), and G-coupled receptors ADRA1D (adrenoceptor α 1D), NTSR1 (neurotensin receptor 1) and ADGRL4. These receptors can activate mitogenic responses. For example, ADRA1D has been implicated in growth-promoting pathways [39].  In the RBC2 SCs, 525 DEGs were identified (|Log2FC| > 2.0) with 328 (62.5%) being up regulated. The 40 significant DE genes with the greatest expression change in each comparison are presented in Figure 5. DEGs with the greatest down regulation in the RBC2 SCs included ND4L (NADH dehydrogenase subunit 4 L), and G-coupled receptors ADRA1D (adrenoceptor α 1D), NTSR1 (neurotensin receptor 1) and ADGRL4. These receptors can activate mitogenic responses. For example, ADRA1D has been implicated in growth-promoting pathways [39].  In the NCT SCs, 395 DEGs were identified (|Log 2 FC| > 2.0) with 248 (62.8%) being up regulated by heat treatment (Figure 2, Supplementary Table S2). DEGs with the greatest down regulation included TMED6 (transmembrane p24 trafficking protein 6), AKR1D1 (aldo-keto reductase family 1 member D1) and ENSMGAG00000020839 (novel protein with BLAST similarity to collagen α-2(I) chain-like, LOC104913712). Two of the 20 most down-regulated DEGs (ENSMGAG00000019963 and ENSMGAG00000022949) were shared with the DEGs found in the RBC2 SCs. Both of these are annotated in ENSEMBL as novel protein coding genes. ENSMGAG00000019963 has BLAST similarity to uncharacterized protein C9orf152 and ENSMGAG00000022949 is likely a ncRNA. Highly up-regulated genes included LOC100546217 (protein TENP, also known as BPIFB2 BPI fold containing family B member 2), TACR1 (tachykinin receptor 1), and ENSMGAG00000020551 (annotated novel protein with BLAST similarity to CLASP2, CLIP-associating protein 2). BPIFB2 is a lipid transfer/lipopolysaccharide binding protein [40]. The tachykinin receptor 1 (TACR1) is characterized by interactions with G proteins [41], and this DEG was also significantly up regulated in the heat-treated RBC2 SCs. CLASP2 is a tracking protein that promotes the stabilization of dynamic microtubules [42].
Analysis of the DEGs with |Log 2 FC| > 2.0 (205 IDs mapped to G. gallus gene set) in PANTHER found greatest enrichment in the GO Biological Process category Cell adhesion (4.00-fold, p = 1.00 × 10 −05 ). Based on IPA analysis, the canonical pathway of Oxidative phosphorylation was significantly affected in the NCT SCs as it was for the RBC2 cells ( Figure 6). However, in contrast to the RBC2 cells, the activation z-score was negative (z-score = −6.6, −log(pval) = 10.6). Significant positive z scores were seen for the Pulmonary fibrosis idiopathic signaling pathway (5.9) and Cardiac hypertrophy signaling (5.6).  In the comparison of heat-treated cells, 115 DEGs were shared between the RBC2 and NCT SCs. Of these, 61 were up regulated by heat treatment in both lines, 39 were down regulated and 15 responded directionally opposite (Supplementary Table S2). Interestingly, 13 of these 15 DEGs were up regulated in the NCT cells with only two showing higher expression in the RBC2 SCs. The two DEGs showing exclusive up regulation in the RBC2 SCs were LOC104916271 (osteocalcin-like, Log2FC = 2.1317 and −5.1310 in the RBC2 and NCT SCs, respectively) and RHCG (Rh family C glycoprotein, Log2FC = 2.4501 and −2.1274, in the RBC2 and NCT SCs, respectively). Osteocalcin is thought to play a role in metabolic regulation and is secreted solely by osteoblasts [43]. In the comparison of heat-treated cells, 115 DEGs were shared between the RBC2 and NCT SCs. Of these, 61 were up regulated by heat treatment in both lines, 39 were down regulated and 15 responded directionally opposite (Supplementary Table S2). Interestingly, 13 of these 15 DEGs were up regulated in the NCT cells with only two showing higher expression in the RBC2 SCs. The two DEGs showing exclusive up regulation in the RBC2 SCs were LOC104916271 (osteocalcin-like, Log 2 FC = 2.1317 and −5.1310 in the RBC2 and NCT SCs, respectively) and RHCG (Rh family C glycoprotein, Log 2 FC = 2.4501 and −2.1274, in the RBC2 and NCT SCs, respectively). Osteocalcin is thought to play a role in metabolic regulation and is secreted solely by osteoblasts [43].
The 13 DEGs showing higher expression in the NCT cells included the transporter ABCC1 (ATP binding cassette subfamily C member 1), membrane proteins CLTC (clathrin heavy chain) and FREM1 (FRAS1 related extracellular matrix 1), INF2 (Inverted formin, FH2 and WH2 domain containing), LOC100550636 (C2 domain-containing protein 3), LOC104910681 (tyrosine-protein phosphatase non-receptor type 13-like), TESK2 (testisspecific kinase 2) and 6 uncharacterized/novel proteins (ENSMGAG00000015050, ENSM-GAG00000018377, ENSMGAG00000019695, ENSMGAG00000020631, ENSMGAG00000021377, and ENSMGAG00000022243). The different responses of the NCT and RBC2 SCs to the heat treatment in several canonical pathways is evident in the IPA Comparison analysis illustrated ( Figure 6). Directional differences in the activation score were seen for the majority of the 20 pathways with the highest composite z-scores.
Overrepresentation tests in PANTHER based on all DEGs (FDR p < 0.05) found significant altered expression of genes involved in muscle organization and development (Table 4). When the analysis was limited to DEGs up regulated in the NCT SCs relative to the RBC2 cells, the significant GO Biological Processes of Mitochondrial ATP synthesis coupled electron transport (5.64-fold, p = 3.53 × 10 −04 ), ATP synthesis coupled electron transport (5.59x, p = 1.70 × 10 −04 ), and Aerobic electron transport chain (5.56x, p = 9.50 × 10 −04 ) had the greatest fold change. The GO Biological processes of Sterol biosynthetic process (5.57-fold, p = 4.23 × 10 −02 ) and Regulation of protein-containing complex disassembly (3.34x, p = 1.55 × 10 −02 ) had the greatest fold change when the analysis was limited to only downregulated DEGs.
LOC100544159 (nipped-B-like protein), LOC100545914 (ATP synthase subunit α), LOC100550847 (spindlin-W-like), LOC100544169 (thioredoxin-like protein 1), LOC100548376 (transitional endoplasmic reticulum ATPase), LOC100303669 (ubiquitinassociated protein 2), LOC104915564 (zinc finger RNA-binding protein-like), LOC104914913 (zinc finger SWIM domain-containing protein 6-like), ROBO2 (roundabout guidance receptor 2), and 10 novel/uncharacterized loci. Among these DEGs, the membrane-bound receptor NTRK2, is of particular interest in that signaling through this kinase leads to cell differentiation [44]. The ligand BMP3 has been shown to induce bone and cartilage development [45]. Overrepresentation tests in PANTHER based on all DEGs (FDR p < 0.05) found significant altered expression of genes involved in muscle organization and development (Table 4). When the analysis was limited to DEGs up regulated in the NCT SCs relative to the RBC2 cells, the significant GO Biological Processes of Mitochondrial ATP synthesis coupled electron transport (5.64-fold, p = 3.53 × 10 −04 ), ATP synthesis coupled electron transport  Table  S2). Members of the Wnt gene family are involved in several developmental processes, including cell proliferation, differentiation, migration, and patterning during development [46]. In humans, angiopoietin 1 is an important mediator in migration, growth, and differentiation of endothelial cells during angiogenesis [47]. Greatest down regulation was seen for KCMF1 (Log 2 FC = −10.73), LOC100550847 (Log 2 FC = −10.01), and LOC100303669 (ubiquitin-associated protein 2, Log 2 FC = −10.17). LOC100303669, involved in the ubiquitination pathway, is unique in that it was not among the 26 DEGs shared with all three temperature contrasts. Overrepresentation tests in PANTHER based on all DEGs (FDR p < 0.05) in the NCT/RBC2 line contrast at 33 • C predicted greatest fold change in the GO Biological Process categories Exocrine system development (GO:0035272) and Sarcomere organization (GO:0045214) ( Table 4). When the analysis was limited to DEGs up regulated in the NCT SCs relative to the RBC2 cells, the GO Biological Process categories of Regulation of muscle contraction (4.36-fold, p = 3.04 × 10 −02 ), Aerobic respiration (3.53x, p = 1.20 × 10 −02 ) and Striated muscle cell differentiation (3.4x, p = 7.60 × 10 −03 ) had the greatest fold change. The GO Biological Process categories of Mitotic sister chromatid segregation (5.68-fold, p = 2.95 × 10 −07 ) and Sister chromatid segregation (5.12x, p = 1.86 × 10 −07 ) had the greatest fold change based on the down-regulated genes. Similar to the control temp analysis, the Oxidative phosphorylation pathway had the highest activation z-score in IPA (z = 4.26, -log(pval) = 3.09) while the lowest activation score was observed for Cell cycle control of chromosome replication (z = −2.89, −log(pval) = 2.21) (Figure 8b). Differential expression of genes in the pathways with greatest |z|-score was more balanced than at 38 • C with similar percentages of upand down-regulated genes. dative phosphorylation pathway (z = −6.0, −log(pval) = 5.97) (Figure 8c). Differential expression of genes in the pathways with greatest |z|-score included a greater percentage of upregulated genes.  A greater number of DEGs were identified between the lines in the heat treatment than the cold. At 43 • C, 249 genes showed |Log 2 FC| > 2.0 in the NCT SCs compared to the RBC2 cells (Figure 7). In this comparison, the majority of DEGs (55.8%) were up regulated in the NCT SCs compared to the RBC2 cells and 183 were unique to the comparison. Greatest up regulation in the NCT SCs was seen for FREM1 (FRAS1 related extracellular matrix 1, Log 2 FC= 8.22) and ADRA1D (adrenoceptor α 1D, Log 2 FC = 6.19) (Supplementary  Table S2). In humans, FREM1 is described as an extracellular matrix protein required for epidermal adhesion during embryonic development [48]. The α-adrenergic receptor ADRA1D, helps to regulate growth and proliferation through the influx of extracellular calcium [49]. Greatest down regulation was seen for KCMF1 (RING-type E3 ubiquitin transferase, Log 2 FC = −10.73), and LOC100544169 (TXNL1, thioredoxin-like protein 1, Log 2 FC = −9.34). Thrioredoxin enables disulfide oxidoreductase activity and reacts to regulate enzyme response to excessive reactive oxygen species [50].

Line Differences: NCT Versus F-Line SCs
Although different genome annotations were used for primary analyses, a comparison was made of the responses of the differentiating 1-wk SCs in the present study (NCT) to those of the 7-wk SCs previously examined (F-line) [25]. A composite gene list was created using gene IDs common to the two gene lists (n = 11,565). For the cold-treated cell comparison, the combined gene list included 279 DEGs (FDR p-value < 0.05 and |Log 2 FC| > 2.0) in the NCT SCs and 420 DEGs in the F-line cells (Supplementary Table S3). Comparison of the lists of DEGs, found 128 unique to NCT and 269 unique to the F-line; 151 DEGs were shared. The majority of unique genes were down regulated in the NCT cells (54%) as was the case for the DEGs of the F-line SCs (70%). The 151 shared DEGs were primarily down regulated (92%) and directionality and magnitude of the fold changes were similar between lines. GO analysis of the DEGs from the cold-treated cell comparison indicates association of these genes with myofibril assembly, and muscle development and contraction.
For the heat-treated cell comparison, the combined gene list included 193 DEGs (p-value < 0.05 and |Log 2 FC| > 2.0) in the NCT SCs and 175 DEGs in the F-line cells (Supplementary Table S4). Comparison of the lists of DEGs found 160 unique to NCT and 142 unique to the F-line; 33 DEGs were shared. The majority of unique genes were up regulated in the NCT cells (61%) as was the case for the F-line SCs (67%). GO analysis of the unique DEGs indicates association with cell adhesion. The 33 shared DEGs had slightly more down-regulated genes (58%) and like the cold-treated cell comparison, the directionality and magnitude of the fold changes were similar between lines. These genes also show enrichment associated with regulation of cell-substrate adhesion.
Next, the 11,565 genes common in the two studies were mapped in IPA (n = 9961) and were used to run a Comparison analysis (33 • C versus 38 • C and 43 • C versus 38 • C) using the FDR p-values and Log 2 FC as variables. Comparison analysis of DEGs of the NCT SCs to those of the F-line (16-wk bodyweight selected) found the significant differences in activation (z-score) both between lines (within temperature) comparisons and between temperature comparisons for several canonical pathways (Figure 9). These comparisons identify genes responding to thermal stress that are different between the 1-wk NCT SCs and the 7-wk F-line cells. Expression differences in these genes may be attributed to both genetic differences between the bird lines and the age of the birds from which the SCs were isolated (1 wk versus 7 wk). Studies of the developmental cascade in turkey skeletal muscle suggest that differentiation potential of skeletal muscle SCs significantly decreases by 4 wk of age [5]. lated in the NCT cells (61%) as was the case for the F-line SCs (67%). GO analysis of the unique DEGs indicates association with cell adhesion. The 33 shared DEGs had slightly more down-regulated genes (58%) and like the cold-treated cell comparison, the directionality and magnitude of the fold changes were similar between lines. These genes also show enrichment associated with regulation of cell-substrate adhesion.
Next, the 11,565 genes common in the two studies were mapped in IPA (n = 9961) and were used to run a Comparison analysis (33 °C versus 38 °C and 43 °C versus 38 °C) using the FDR p-values and Log2FC as variables. Comparison analysis of DEGs of the NCT SCs to those of the F-line (16-wk bodyweight selected) found the significant differences in activation (z-score) both between lines (within temperature) comparisons and between temperature comparisons for several canonical pathways (Figure 9). These comparisons identify genes responding to thermal stress that are different between the 1-wk NCT SCs and the 7-wk F-line cells. Expression differences in these genes may be attributed to both genetic differences between the bird lines and the age of the birds from which the SCs were isolated (1 wk versus 7 wk). Studies of the developmental cascade in turkey skeletal muscle suggest that differentiation potential of skeletal muscle SCs significantly decreases by 4 wk of age [5].

Discussion
Embryonic muscle growth occurs through the myoblast leading to muscle fiber formation by hyperplasia. Posthatch hypertrophic muscle growth involves the proliferation, differentiation, and fusion of SCs with existing myofibers. Satellite cells are a heterogeneous population of cells with varying proliferation and differentiation capacities [51] and their activity during this time is controlled by growth factors, intracellular signaling pathways and transcription factors. Key in modulating the cellular microenvironment is antagonistic Notch and Wnt signaling. The membrane receptor Notch, activated by the ligand Delta (DLL1), is involved in initiating proliferation. Canonical Wnt/β-catenin pathway regulates the SCs transition from cell proliferation to differentiation through the membrane receptor Frizzled [11,52]. As the downstream target of MyoD, DLL1 has been reported to promote myogenic differentiation via the MyoD/DLL1/Notch gene axis [53]. Comparison of RNAseq data between the differentiating cells and proliferating cells of our previous study found expression of myogenic marker genes including MyoD and MyoG was higher in the turkey SCs after 48 h of differentiation consistent with their progressive development [26].
The nuclear receptor peroxisome proliferator-activated receptor-γ (PPARγ) is required for induction of adipogenic differentiation of myogenic cells [54]. As reported by Xu et al. [55], differentiation in both RBC2 and NCT satellite cells is significantly affected by thermal stress. In the differentiating turkey SCs, PPARγ was significantly down regulated, consistent with previous results suggesting that adipogenic potential of the turkey SCs is higher during proliferation and gradually decreases with differentiation [55]. Also down regulated was MSTN (myostatin) and MYF6 (MRF4 or herculin). Early muscle growth requires a balance between the proliferation of new precursor cells and differentiation into muscle fibers. This process is tightly regulated by growth factors such as myostatin, a key negative regulator of muscle growth that, in turn, down regulates the regulatory factors MyoD and MyoG, thereby balancing cell proliferation and differentiation [56]. In birds, MSTN has different isoforms resulting from the alternative splicing of MSTN mRNA. MSTN-A is anti-myogenic isoform whereas MSTN-B is pro-myogenic [57]. The relative abundance of these two protein isoforms in the differentiating turkey SCs is unknown, but may be pro-myogenic. Down regulation of MYF6 may seem contraindicated however, its expression is typically associated with older myotubes where it is thought to promote myofiber maturation and maintenance of the differentiated state [58]. Recent studies in mice also indicate a novel role for Myf-6 in promoting myokine expression to block premature differentiation, preventing stem cell exhaustion [59]. Further study of these proteins in birds is warranted.
The only gene that was consistently up regulated at all temperatures in the NCT SCs compared to the RBC2 cells was ETNPPL (ethanolamine-phosphate phospho-lyase). This gene functions in phospholipid metabolism and is associated with fatty acid and lipid synthesis [60,61]. Although a significant decrease in lipid content was observed in both NCT and RBC2 cells at 24 h of differentiation under cold stress (33 • C) [55], the differential response in ETNPPL between the lines supports observations that growth selection has increased the adipogenic potential of SCs, increasing intramuscular fat deposition in the faster-growing turkeys.
Of the 15 genes consistently down regulated at all the temperatures in the NCT SCs compared to RBC2 cells, three of them (ENSMGAG00000022282, KCMF1, and LOC100303669) are involved in the ubiquitin proteasome pathway. This suggests that breast muscle in NCT birds has a reduced potential for degeneration and atrophy through the ubiquitin proteasome pathway compared to the RBC2 line. Also consistently down regulated was NIPBL (LOC100544159, nipped-B-like protein), a gene associate with cell cycle progression [62]. Down regulation of NIPBL may induce cell cycle arrest and would thus potentially act to delay myogenesis in the NCT SCs during differentiation.

Effect of Cold Stress
Cold treatment resulted in significant gene expression changes in the SCs from both turkey lines, with the primary effect being down regulation of affected genes (almost 3-fold). A larger number of DEGs were seen in the NCT cells, consistent with our previous study where a greater number of genes were affected by cold treatment in the 7-wk SCs from the growth-selected F-line than the Randombred RBC2 line [25]. Cold stress in differentiation has also been shown to down regulate the mTOR and Wnt/PCP pathways to a greater extent in NCT SCs compared to the RBC2 cells [63,64].
Clearly the NCT and RBC2 SCs have different responses during cold stress. Notably, NF-κB, a key stimulator of muscle atrophy and inflammatory myopathies [65], is differentially affected by cold only in the NCT cells. Changes in NF-κB signaling may affect muscle atrophy and cause inflammatory myopathies. Under cold stress, DEGs involved in regulation of skeletal muscle tissue regeneration and sarcomere organization were significantly overrepresented with overall down regulation of muscle-associated genes (MYH) in both cell lines. Interruption of muscle fiber assembly and sarcomere organization by the cold stress would slow differentiation resulting in smaller myotubes (decreased width) as seen in previous studies [66]. Still, expression of genes involved in muscle sarcomere organization, myofibril assembly, and striated muscle cell development remained significantly higher in the NCT SCs compared to the RBC2 cells at the cold (33 • C) and control (38 • C) temperatures. These results reflect the effect of growth selection on increasing myogenic potential (myotube formation) and myofiber hypertrophy.
Altered expression of the oxidative phosphorylation pathway was predicted in the SCs under cold stress and show higher predicted activation in the NCT SCs compared to the RBC2 cells at control and cold temperatures. Altered oxidative phosphorylation may result in oxidative stress, which can lead to oxidative muscle damage [67,68]. This suggests that NCT SCs may have increased oxidative potential relative to the RBC2 cells, important in limiting the incidence of myopathies associated with oxidative stress. Predicted effects on oxidative stress pathways have also been observed in cold-stressed turkey poults [69].

Effect of Heat Stress
A larger suite of genes were significantly affected by heat treatment (43 • C) than cold (33 • C) reflecting the sensitivity of the SCs to heat stress. Many cell surface protein receptors are down regulated including NTSR1, associated with fat production in skeletal muscle [70], ADGRL4, that stimulates cell cycle progression, and ADRA1D, that regulates myoblast survival and differentiation [71]. Down regulation of these genes may decrease the differentiation, survival and fat production of cells particularly in the RBC2 SCs.
Heat stress resulted in the up regulation of several genes that encode contractile proteins (MYH1E, MYBPC1, MYH2, myosin motor domain-containing protein), particularly in the RBC2 cells. Commercial birds grow larger muscle fibers than the RBC2 line, and therefore increase in the contractile protein myosin is expected. However, an increase during heat stress suggests promotion of muscle growth. Also included in the affected genes were DEGs associated with sarcomere assembly. For example, in the NCT cells, CLASP2, involved in microtubule organization and cell adhesion [72], was significantly up regulated suggesting enhanced myotube formation (migration and fusion) during heat stress. Greatly up regulated in the NCT cells were genes reported to regulate myoblast survival and differentiation (ADRA1D, [71]) and promote cell adhesion (FREM1, [48]). Their increased expression in SCs may promote increased hypertrophic potential of the p. major muscle and sarcomere assembly in the commercial turkeys.
The proliferation and differentiation of turkey SCs under thermal stress also provides insight into the regeneration ability of muscle and the myopathies that negatively affect poultry. For example, Wooden Breast is a necrotic/fibrotic myopathy that negatively impacts breast meat quality in chickens with meat product downgrades or condemnation. In commercial broilers affected with Wooden Breast, sarcomere organization is reduced, especially in small and intermediate-sized fibers [73]. Irregular, disconnected, and reduced sarcomeric Z bands are observed by 6 wk of age in Wooden Breast-affected muscles [74]. These data suggest that the SC-mediated repair process of necrotic fiber damage in broilers is not able to restore sarcomeric muscle fiber structure. However, there are no reports of necrotic/fibrotic myopathies like Wooden Breast in modern commercial meat-type turkeys, although these heavy-weight birds are in part, selected for p. major muscle growth.
Although growth selection strategies in turkeys and broilers have shared similarities, the outcome on breast muscle growth, development, and structure are distinct. The difference between current meat-type broilers and turkeys appears to at least partially involve differences stemming from the SCs of the p. major muscle that are responsible for the post-hatch hypertrophic growth and muscle regeneration. Selection for breast growth in chickens has created birds with greater muscle hypertrophy that respond differently to environmental stressors and may exhaust SC pools leading to diminished SC function [75]. Prior to the onset of Wooden Breast, current broilers may undergo decreased SC proliferation and differentiation (Velleman, unpublished observation). In contrast, Xu et al. [66] has shown an increase in proliferation and differentiation of current commercial turkey SCs compared to cells from a line of birds established prior to and without intense growth selection. Data from the current study shows that expression of genes promoting muscle regeneration (differentiation) and sarcomere assembly is higher in SCs from faster-growing modern commercial turkeys compared to slower-growing historic turkeys. Taken together, these findings may explain why necrotic/fibrotic myopathy like Wooden Breast is not found in modern-commercial turkeys.

Conclusions
Timing of thermal stress during poultry production can differentially affect bird performance and production traits. Movement of poults from hatch to grower facilities during the period where the thermal regulatory system is poorly developed, can expose them to acute thermal conditions, either hot or cold, that may have long-term consequences for muscle development and ultimately meat quality. Understanding how temperature affects SC proliferation and differentiation will allow for development of better thermal management strategies especially to confront the uncertainties brought about by climate change. Growth selection appears to have increased the proliferation, differentiation, adipogenic potential, oxidative potential, and decreased the degenerative and atrophic potential of the turkey p. major SCs. The breast muscle of fast-growing meat-type birds is characterized by its excessive hypertrophic myofiber growth potential with decreased incidence of muscle atrophy and degeneration. The biology of the SCs and their effect on morphological structure of the breast muscle is of continued importance in the poultry industry, especially in light of emerging myopathies [76].
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/genes13101857/s1, Figure S1: Distribution of DEGs during differentiation of cultured turkey p. major SCs; Table S1. Mean quality-trimmed RNAseq read counts for turkey p. major muscle SCs from two genetic lines (RBC2 and NCT) after 48 h differentiation; Table  S2. Summary of pairwise differential gene expression (DESeq) analysis of differentiating satellite cell transcriptomes; Table S3. Comparison of DEGs of cold-treated 1-wk NCT SCs to those of cold-treated 7-wk SCs of the F-line [25]; Table S4. Comparison of DEGs of heat-treated 1-wk NCT SCs to those of heat-treated 7-wk SCs of the F-line [25].