RNA-Sequencing Reveals Upregulation and a Beneficial Role of Autophagy in Myoblast Differentiation and Fusion

Myoblast differentiation is a complex process whereby the mononuclear muscle precursor cells myoblasts express skeletal-muscle-specific genes and fuse with each other to form multinucleated myotubes. The objective of this study was to identify potentially novel mechanisms that mediate myoblast differentiation. We first compared transcriptomes in C2C12 myoblasts before and 6 days after induction of myogenic differentiation by RNA-seq. This analysis identified 11,046 differentially expressed genes, of which 5615 and 5431 genes were upregulated and downregulated, respectively, from before differentiation to differentiation. Functional enrichment analyses revealed that the upregulated genes were associated with skeletal muscle contraction, autophagy, and sarcomeres while the downregulated genes were associated with ribonucleoprotein complex biogenesis, mRNA processing, ribosomes, and other biological processes or cellular components. Western blot analyses showed an increased conversion of LC3-I to LC3-II protein during myoblast differentiation, further demonstrating the upregulation of autophagy during myoblast differentiation. Blocking the autophagic flux in C2C12 cells with chloroquine inhibited the expression of skeletal-muscle-specific genes and the formation of myotubes, confirming a positive role for autophagy in myoblast differentiation and fusion.


Introduction
Skeletal muscle is the largest tissue in the body, accounting for nearly half of the body's entire weight, and has many important functions, including contraction [1]. The basic structural units of skeletal muscle are myofibers, which are multinucleated muscle cells differentiated and fused from the muscle precursor cells myoblasts [2][3][4][5]. Significant changes occur in gene transcription during myoblast differentiation and fusion into myofibers, and these changes are regulated in the nucleus by cis-regulatory elements, including enhancers and promotors, and trans-regulatory factors, including transcription factors, co-factors, and other DNA-targeted regulatory proteins [6,7]. Four master transcription factors have been identified to regulate gene transcription during myogenesis, including myogenic factor 5 (MYF5), MYF6 (also known as MRF4), MyoD (MYOD1), and myogenin (MYOG), of which MYOG is the major transcription factor that regulates gene transcription during myoblast differentiation [8,9]. DNA methylation plays a vital role in determining the expression of genes critical for myoblast differentiation via changing the accessibility of chromatin to transcription factors or co-factors [10,11]. Within the cytoplasm, multiple signal transduction pathways, including the mTOR and AMPK pathways, ensure accurate epigenetic control of gene expression during myoblast differentiation [12,13]. Additionally, microRNA and lncRNA control the protein buildup in differentiating myoblasts partly via altering the stability of mRNA [14,15]. On the cell surface, various receptors respond to extracellular signals, such as amino acids, hormones, cytokines, and mechanical damage, and initiate the intracellular signal transduction pathways involved in myoblast differentiation [9,[16][17][18].
It can be easily imagined that significant protein turnover occurs in the cytoplasm during myoblast differentiation and fusion and that myoblast differentiation and fusion

Gene Expression Analysis and Gene Ontology Analysis of RNA-Seq Data
The raw sequencing reads were first filtered with Trimmomatic [34] to remove lowquality reads and reads containing adapters. Clean reads were mapped to the mouse reference genome (mm10) using the STAR (v2.5) program [35]. The uniquely mapped reads were used in the quantification of gene expression levels, which was performed using the HTSeq v0.6.1 program [36]. Gene expression levels were calculated as reads per kilobase of exon model per million mapped reads (FPKM) [37]. A differential expression analysis was performed using the DESeq2 R package (2_1.6.3) [38], and the resulting p-values were adjusted using the Benjamini-Hochberg approach for controlling the false discovery rate (FDR). Genes with adjusted p-values (Padj) < 0.05 were considered as differentially expressed genes (DEGs). Gene ontology (GO) analysis of DEGs was performed in three categories: biological process (BP), cellular compound (CC), and molecular function (MF), using the R package clusterProfiler [39,40]. In this analysis, the p-value and q-value cutoffs were 0.01 and 0.05, respectively.

Reverse Transcription Quantitative PCR (RT-qPCR)
One microgram of total RNA was denatured by incubation with random primers at 70 • C for 10 min, followed by cooling on ice for 5 min. Reverse transcription of total RNA was performed by incubating 5 µL RNA-random primer mix with 4 µL 5 × reverse transcription buffer, 4.8 µL MgCl 2 , 1 µL dNTP, 1 µL ImProm-II reverse transcriptase, 0.5 µL RNasin Ribonuclease Inhibitor, and 3.7 µL RNase free H 2 O at 42 • C for 90 min, at 70 • C for 10 min, and then at 4 • C for 5 min. Quantitative PCR of 20 ng cDNA was performed in duplicate using Fast SYBR Green Master Mix on a 7500 Fast Real-Time PCR system (Applied Biosystems). The PCR primers used in this study were designed in a previous study [41]. Relative gene expression levels were calculated using the ∆∆Ct method [42], using the HMBS gene as a reference gene, which, based on its Ct value (data not shown), was stably expressed in C2C12 cells during differentiation.

Immunocytochemistry
C2C12 cells were first fixed by incubating them with 4% formaldehyde solution for 15 min at room temperature and then rinsed twice with PBS. Cells were then permeabilized by incubation with 0.25% Triton X-100 for 10 min at room temperature. After that, cells were rinsed again with PBS and then incubated with 1% BSA and 0.05% Tween-20 in PBS for one hour at room temperature with shaking to block nonspecific antibody binding. Cells were then incubated with the antibody for myosin heavy chains (NA-4, DSHB, Iowa City, IA, USA) at 1:100 in PBS at 4 • C overnight. Cells were rinsed twice with PBS and incubated with an anti-mouse IgG FITC antibody (Sigma-Aldrich) at 1:200 dilution for 1 h at room temperature. Nuclei were stained by incubating cells with 1 µg/mL of 4 ,6-Cells 2022, 11, 3549 4 of 16 diamidino-2-phenylindole (DAPI) for 10 min at room temperature. Images of stained cells were taken with a florescence microscope.

Determination of Fusion Index, Myotube Length, and Myotube Area
Numbers of nuclei were counted and lengths of myotubes and areas of myotubes were measured from randomly taken immunofluorescent images of myotubes using ImageJ [43]. A myotube was defined as a myosin-heavy-chain-positive ("green") cell containing 2 or more DAPI-stained nuclei ("blue"); the length of a myotube was defined as the longest distance between two ends of a myotube; and fusion index was defined as the ratio of the number of nuclei in myotubes to the number of total nuclei counted. Around 1000 total nuclei were counted, and 15 myotubes were measured for each treatment in each experiment.

Total Cellular Protein Extraction and Western Blot Analysis
Total protein lysates from C2C12 cells were prepared by scraping and incubating them in the RIPA buffer with proteinase inhibitors added for 30 min on ice followed by centrifugation at 16,000× g, 4 • C, for 15 min. Protein concentrations were determined using a Pierce BCA Protein Assay Kit (23227, obtained from Thermo Fisher Scientific). For the Western blot analysis, 20 µg of total protein lysates were separated using sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) with 6% stacking and 15% separating gel. Gel electrophoresis was run at 80 V for 30 min and then 120 V for 80 min at 4 • C. Following separation, protein was transferred from the gel to a nitrocellulose membrane by electrophoresis at 70 V, 4 • C, for 60 min. To block nonspecific antibody binding, the membrane was incubated with 5% non-fat milk in the TBST buffer for 1 h at room temperature. The membrane was then incubated with a primary antibody overnight at 4 • C. The membrane was rinsed twice with TBST and then incubated with a secondary antibody for 1 h at room temperature. The following primary antibodies were used in Western blot analyses: ß-tubulin at 1:1000 dilution (E7 from DSHB), LC3B at 1:1000 dilution (NB100-2220 from Novus Biologicals, Centennial, CO, USA), and myogenin at 1:25 dilution (F5D from DSHB). The following secondary antibodies were used in Western blot analyses: IRDye 800CW goat anti-rabbit IgG secondary antibody (926-32211, LI-COR Biosciences, Lincoln, NE, USA) at 1:15000 dilution and IRDye 800CW goat anti-mouse IgG secondary antibody (926-32210, LI-COR Biosciences) at 1:15,000 dilution. Western blots were visualized with the LI-COR Odyssey Infrared Image System 9120, and the fluorescent intensities were quantified with the Image Studio Lite software (LI-COR Biosciences).

Statistical Analyses
Data were analyzed by ANOVA, followed by multiple comparisons using Tukey's test or Dunnett's test. Two-group comparisons were carried out via t-tests. All statistical analyses were performed in JMP Pro 16 (SAS, Cary, NC, USA) or GraphPad Prism 9 (San Diego, CA, USA). All data are expressed as means ± standard errors.

More Than Ten Thousand Genes Were Differentially Expressed during Myoblast Differentiation
We performed an RNA-seq analysis to identify genes differentially expressed during myogenic differentiation in C2C12 cells. The sequencing of five RNA-seq libraries constructed from C2C12 cells immediately before induction of myogenic differentiation and five RNA-seq libraries constructed from C2C12 cells on day 6 of induced myogenic differentiation generated more than 36 million unique mapped sequencing reads per library (Supplementary Table S1). Gene expression quantification and subsequent differentiation expression analysis revealed 11,046 genes differentially expressed (Padj < 0.05) between the two conditions, of which 5615 genes were upregulated and 5431 genes downregulated from the day before induction of differentiation to day 6 of differentiation ( Figure 1A, Supplementary Table S2). Examples of genes upregulated in C2C12 cells during induced myogenic differentiation were Myh1, Myh2, Myh3, Myh4, Myog, Ckm, Mymk, and Mb, which are all known to be specifically or preferentially expressed in skeletal muscle [44,45]. Examples of genes downregulated during this differentiation included Ccna2, Ccnd1, Ccne1, Ccne2, Id1, Id3, and Cdc6, which are known to function in cell proliferation [46,47]. To confirm the gene expression differences quantified by RNA-seq, we measured the expression levels of Myh1, Myh2, Myh3, Myh4, and Myog mRNAs in the two conditions by RT-qPCR. As shown in Figure 1B, average fold changes measured by RT-qPCR were similar to those quantified by RNA-seq for all of these mRNAs. structed from C2C12 cells immediately before induction of myogenic differentiation and five RNA-seq libraries constructed from C2C12 cells on day 6 of induced myogenic differentiation generated more than 36 million unique mapped sequencing reads per library (Supplementary Table S1). Gene expression quantification and subsequent differentiation expression analysis revealed 11,046 genes differentially expressed (Padj < 0.05) between the two conditions, of which 5615 genes were upregulated and 5431 genes downregulated from the day before induction of differentiation to day 6 of differentiation ( Figure 1A, Supplementary Table S2). Examples of genes upregulated in C2C12 cells during induced myogenic differentiation were Myh1, Myh2, Myh3, Myh4, Myog, Ckm, Mymk, and Mb, which are all known to be specifically or preferentially expressed in skeletal muscle [44,45]. Examples of genes downregulated during this differentiation included Ccna2, Ccnd1, Ccne1, Ccne2, Id1, Id3, and Cdc6, which are known to function in cell proliferation [46,47]. To confirm the gene expression differences quantified by RNA-seq, we measured the expression levels of Myh1, Myh2, Myh3, Myh4, and Myog mRNAs in the two conditions by RT-qPCR. As shown in Figure 1B, average fold changes measured by RT-qPCR were similar to those quantified by RNA-seq for all of these mRNAs.

Autophagy Was Upregulated during Myoblast Differentiation
Gene ontology enrichment analyses of genes upregulated during myogenic differentiation in C2C12 cells revealed the biological processes (BPs), cellular components (CCs), and molecular functions (MFs) associated with these genes (Supplementary Table S3). Shown in Figure 2A are the top 10 upregulated BPs, CCs, and MFs during myogenic differentiation in C2C12 cells. Shown in Table 1 are two examples of top upregulated BPs and associated genes. As expected, most of these upregulated BPs were related to skeletal

Autophagy Was Upregulated during Myoblast Differentiation
Gene ontology enrichment analyses of genes upregulated during myogenic differentiation in C2C12 cells revealed the biological processes (BPs), cellular components (CCs), and molecular functions (MFs) associated with these genes (Supplementary Table S3). Shown in Figure 2A are the top 10 upregulated BPs, CCs, and MFs during myogenic differentiation in C2C12 cells. Shown in Table 1 are two examples of top upregulated BPs and associated genes. As expected, most of these upregulated BPs were related to skeletal muscle development and maturation, including muscle cell differentiation, muscle system process, striated muscle cell differentiation, muscle contraction, and muscle tissue development ( Figure 2A). Interestingly, four of the top upregulated BPs were related to autophagy ( Figure 2A, Table 1). Similar to the BPs, most of the top upregulated CCs during myogenic differentiation in C2C12 cells were related to skeletal muscle structures, such as sarcomeres, I-bands, and Z-discs ( Figure 2A). Two of the top upregulated CCs were related to organelles that are involved in protein degradation, namely, lysosomes and lytic vacuoles (Figure 2A). The top upregulated MFs during myogenic differentiation in C2C12 cells included actin binding, motor activity, and calmodulin binding, which are all known as functions of mature skeletal muscle [48,49]. It is also interesting to note that small GTPase binding, Ras GTPase binding, coenzyme binding, Rab GTPase binding, protein serine/threonine kinase activity, and enzyme activator activity were among the MFs upregulated during myogenic differentiation in C2C12 cells (Figure 2A). regulated during myogenic differentiation in C2C12 cells (Figure 2A).
GO analyses of genes downregulated during myogenic differentiation in C2C12 cells indicated that many of them were associated with BPs, CCs, and MFs related to DNA replication and RNA processing (Supplementary Table S4). Examples of the top downregulated BPs during myogenic differentiation in C2C12 cells were ribonucleoprotein complex biogenesis and DNA replication; examples of top downregulated CCs were ribosomes and chromosomal regions; and examples of top downregulated MFs were mRNA binding and helicase activity ( Figure 2B, Table 2).    Figure 2B, Table 2).

Myoblast Differentiation Was Associated with Increased Conversion of LC3-I to LC3-II
LC3 protein, encoded by the microtubule-associated proteins 1A/1B light chain 3B (Map1lc3b) gene, is the most widely used marker of autophagy [50]. During autophagy, the cytosolic LC3 protein LC3-I is conjugated to phosphatidylethanolamine to form the autophagosomal membrane protein LC3-II [51]. Western blot analyses ( Figure 3A) showed that the protein expression ratio of LC3-II to LC3-I in C2C12 cells increased rapidly from the day before induction of differentiation (day 0 in the figure) to day 1 of differentiation and remained high on day 2 and day 4 of differentiation ( Figure 3B). These changes indicated that the autophagic flux increased in C2C12 cells during myogenic differentiation, confirming what was revealed by the RNA-seq analysis described above.

Inhibition of Autophagy Reduced the Number, Length, and Size of Myotubes Formed during Myoblast Differentiation
To determine the role of increased autophagy in myoblast differentiation, we induced C2C12 cells to differentiate in the presence of chloroquine (CQ), a widely used inhibitor of autophagic flux that inhibits autophagy by blocking the fusion of autophagosomes with lysosomes [32]. We first performed Western blot analyses to determine

Inhibition of Autophagy Reduced the Number, Length, and Size of Myotubes Formed during Myoblast Differentiation
To determine the role of increased autophagy in myoblast differentiation, we induced C2C12 cells to differentiate in the presence of chloroquine (CQ), a widely used inhibitor of autophagic flux that inhibits autophagy by blocking the fusion of autophagosomes with lysosomes [32]. We first performed Western blot analyses to determine whether CQ is effective in inhibiting autophagy in C2C12 cells. As shown in Figure 4A,B, the ratios of LC3-II to LC3-I in C2C12 cells treated with CQ were significantly higher than those in control C2C12 cells on days 1, 2, and 4 of myogenic differentiation. These increases confirmed the effectiveness of CQ in blocking autophagosome-lysosome fusion and hence the degradation of LC3-II protein by the lysosomal hydrolases. These Western blot analyses also confirmed the increases in the generation of LC3-II in C2C12 cells during myogenic differentiation ( Figure 4A).  Based on the morphology, multinucleated myotubes began to form in C2C12 c two days after the initiation of myogenic differentiation, and multinucleated myotu could be easily observed by day 4 of differentiation ( Figure 5A). C2C12 cells treated w CQ appeared to form fewer and smaller myotubes than untreated C2C12 cells on the sa day of differentiation ( Figure 5A). The higher concentration (20 μM) of CQ apparen caused greater decreases in the number and size of formed myotubes than the lower c Based on the morphology, multinucleated myotubes began to form in C2C12 cells two days after the initiation of myogenic differentiation, and multinucleated myotubes could be easily observed by day 4 of differentiation ( Figure 5A). C2C12 cells treated with CQ appeared to form fewer and smaller myotubes than untreated C2C12 cells on the same day of differentiation ( Figure 5A). The higher concentration (20 µM) of CQ apparently caused greater decreases in the number and size of formed myotubes than the lower concentration (10 µM) of CQ; myotubes were barely seen in C2C12 cells treated with 20 µM CQ on day 4 of differentiation ( Figure 5A). Interestingly, at high magnification, "dark rings" around the nuclei could be observed in CQ-treated C2C12 cells ( Figure 5A). These "dark rings" were probably formed by autophagosomes which accumulated due to their blocked fusion with lysosomes and hence their blocked degradation by lysosomal enzymes.
To more clearly show the morphological differences between CQ-treated and control C2C12 cells, we stained the myosin heavy-chain proteins and nuclei in C2C12 cells on day 4 of differentiation and quantified the number of nuclei fused into myotubes and the lengths and areas of myotubes. As shown in Figure 5B, myotubes in CQ-treated C2C12 cells were clearly fewer and smaller than in control C2C12 cells, and these differences were more obvious with the higher concentrations of CQ. The average fusion indexes of C2C12 cells treated with 10 µM and 20 µM CQ were approximately 30% and 60%, respectively, lower than that of control C2C12 cells (p < 0.05, Figure 5C). The average lengths of myotubes formed from C2C12 cells treated with 10 µM and 20 µM CQ were approximately 31% and 43%, respectively, shorter than that of myotubes formed from control C2C12 cells (p < 0.05, Figure 5C). The average areas of myotubes formed from C2C12 cells treated with 10 µM and 20 µM CQ were 40% and 50%, respectively, smaller than that of myotubes formed from control C2C12 cells (p < 0.05, Figure 5C).

Inhibition of Autophagy Reduced the Expression of Muscle-Specific Genes during Myoblast Differentiation
To further determine the effect of the inhibition of autophagy on myoblast differentiation, we quantified the expression levels of several skeletal-muscle-specific genes in CQ-treated and control C2C12 cells on day 4 of myogenic differentiation. These genes included Myh1, Myh3, Tnnt3, Mb, Ckm, Myog, and Mymk, which encode either structural or functional components of skeletal muscle (Myh1, Myh3, Tnnt3, Mb, and Ckm) or are master regulators of myoblast terminal differentiation and fusion (Myog and Mymk) [44]. As shown in Figure 6A, the expression of most of these genes was inhibited in C2C12 cells treated with 10 µM CQ, and the expression of all of these genes was inhibited in C2C12 myoblasts treated with 20 µM CQ compared with their expression in untreated C2C12 cells.
We also confirmed the effect of the inhibition of autophagy on the expression of myogenin at the protein level. As expected, myogenin protein expression in untreated C2C12 cells was significantly higher on days 1, 2, and 4 of differentiation than on the day before the induction of differentiation ( Figure 6B,C). Compared to the control, CQ significantly reduced myogenin protein expression in C2C12 cells on days 1 and 2 of differentiation ( Figure 6D). Cells 2022, 11, x FOR PEER REVIEW 10 of treated with 10 μM and 20 μM CQ were 40% and 50%, respectively, smaller than that myotubes formed from control C2C12 cells (p < 0.05, Figure 5C).

Discussion
In this study, we first compared the gene expression profiles between undifferenti

Discussion
In this study, we first compared the gene expression profiles between undifferentiated and differentiating C2C12 myoblasts. Our RNA-seq analysis revealed that the expression of many skeletal-muscle-specific genes increased while the expression of many genes involved in the cell cycle decreased during myoblast differentiation. These results were consistent with those of earlier microarray-based gene expression profiles for C2C12 cells [52][53][54], indicating that myoblast differentiation involves a shift in gene expression from genes functioning in cell proliferation to genes functioning in skeletal muscle buildup and contraction. Our RNA-seq analysis revealed that many genes functioning in mRNA processing, ribonucleoprotein complex biogenesis, and ribosomes were downregulated during myogenic differentiation in C2C12 cells. This result suggests that myoblast differentiation is associated with delayed processing (i.e., capping, polyadenylation, and splicing) of newly synthesized mRNA (i.e., pre-mRNA) in the nucleus, delayed export of processed mRNA from the nucleus to the cytoplasm, and delayed translation of exported mRNA in the cytoplasm. These delays may be meant to reduce loss of mRNAs through membrane pores formed during myoblast fusion [55].
Our RNA-seq analysis also revealed that many genes involved in autophagy and related functional terms were upregulated during myoblast differentiation. These gene expression changes as well as those involved in ribonucleoprotein complex biogenesis and mRNA processing mentioned above were not identified by earlier microarray-based gene expression analyses [52][53][54], perhaps because many of the genes involved in autophagy, ribonucleoprotein complex biogenesis, and mRNA processing were not included in the microarrays or in the functional analyses performed in these studies, or because microarraybased gene expression analysis is less sensitive than RNA-seq in detecting changes in gene expression, in particular, changes in low-abundance transcripts [56][57][58].
The upregulation of genes functioning in autophagy during myoblast differentiation indicates that autophagy-mediated protein degradation increases during myoblast differentiation. Indeed, we confirmed increased autophagy in C2C12 cells during myogenic differentiation by detecting increased conversion of LC3-I to LC3-II in differentiating C2C12 myoblasts and by detecting increased accumulation of LC3-II in autophagy-blocked differentiating C2C12 myoblasts. The association of increased autophagy with myoblast differentiation suggests that autophagy might benefit myoblast differentiation. A positive role for autophagy in myoblast differentiation is supported by our results showing that blocking autophagy in C2C12 myoblasts impaired their expression of muscle-specific genes and their fusion into multinucleated myotubes. Our finding that autophagy plays a positive role in myoblast differentiation was consistent with earlier studies using different approaches or different myoblast models [24][25][26][27][28][29]. However, the conclusion that autophagy aids myoblast differentiation seems to contradict the observation that inducing autophagy with the mTOR inhibitor rapamycin inhibited myoblast differentiation in C2C12 cells [59]. One possible reason for this discrepancy is that the target of rapamycin, mTOR, controls not only autophagy but also other processes, such as overall protein synthesis, which is essential for myoblast differentiation [60,61].
Autophagy is a process by which cells maintain hemostasis and survival by delivering organelles and proteins to lysosomes for degradation [62]. There is likely increased generation of damaged, malformed, or unnecessary organelles and proteins during the differentiation and fusion of mononuclear myoblasts into multinucleated myotubes [63]. Therefore, increased autophagy may be one of the means used by myoblasts to clear these organelles and proteins and thereby maintain normal myogenic differentiation. Indeed, this role for autophagy in myoblast differentiation is supported by previous findings that upregulated autophagy is essential for mitochondrial degradation and biogenesis and protection against mitochondrial oxidative stress during myoblast differentiation [64,65].
Myogenin is a master transcriptional regulator of myoblast differentiation [8]. Myogenin binds to and activates the transcription of muscle-specific genes as a heterodimer with E-proteins [66]. However, the DNA-binding and transactivating ability of myogenin is inhibited by a group of proteins called the inhibitor of DNA binding (ID) proteins, which lack basic DNA binding domains [66]. The ubiquitin-proteasome system was previously found to degrade the ID proteins in differentiating myoblasts [67]. We speculate that the autophagy system might also degrade the ID proteins and thereby promote myoblast differentiation. Besides the ID proteins, proteins such as the cyclin-dependent kinases and transformation growth factor beta (TGF-ß) also inhibit myoblast differentiation by keeping myoblasts in the proliferation phase [46,68]. We speculate that the autophagy system could also target and degrade these proteins to promote myoblast differentiation.

Conclusions
A massive number of genes are upregulated or downregulated during myoblast differentiation. Genes upregulated during myoblast differentiation include those associated with skeletal muscle structure and contraction and autophagy. Genes downregulated during myoblast differentiation include those involved in the cell cycle, ribonucleoprotein complex biogenesis, and mRNA processing. Not only the expression of autophagic genes but autophagic activity are increased during myoblast differentiation. Increased autophagy is confirmed to benefit myoblast differentiation and fusion.
Supplementary Materials: The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/cells11223549/s1, Table S1: Summary of RNA sequencing and mapping results; Table S2: List of all differentially expressed genes; Table S3: Biological processes, cellular components, and molecular functions associated with upregulated genes; Table S4: Biological processes, cellular components, and molecular functions associated with downregulated genes.