Bioinformatic analysis of the gene expression profile in muscle atrophy after spinal cord injury

Spinal cord injury (SCI) is often accompanied by muscle atrophy; however, its underlying mechanisms remain unclear. Here, the molecular mechanisms of muscle atrophy following SCI were investigated. The GSE45550 gene expression profile of control (before SCI) and experimental (14 days following SCI) groups, consisting of Sprague–Dawley rat soleus muscle (n = 6 per group), was downloaded from the Gene Expression Omnibus database, and then differentially expressed gene (DEG) identification and Gene Ontology, pathway, pathway network, and gene signal network analyses were performed. A total of 925 differentially expressed genes, 149 biological processes, and 55 pathways were screened. In the pathway network analysis, the 10 most important pathways were citrate cycle (TCA cycle), pyruvate metabolism, MAPK signalling pathway, fatty acid degradation, propanoate metabolism, apoptosis, focal adhesion, synthesis and degradation of ketone bodies, Wnt signalling, and cancer pathways. In the gene signal network analysis, the 10 most important genes were Acat1, Acadvl, Acaa2, Hadhb, Acss1, Oxct1, Hadha, Hadh, Acaca, and Cpt1b. Thus, we screened the key genes and pathways that may be involved in muscle atrophy after SCI and provided support for finding valuable markers for this disease.

. (a) Dendrogram of the comparison between the control (before SCI, n = 6/group) and experimental groups (14 days after SCI, n = 6/ group). The X-axis shows the sample names of the two groups, and the Y-axis shows the 925 differentially expressed genes. Red indicates high expression of the differential gene in the grouped samples, whereas green indicates low expression. Above the plot, the yellow strip indicates the control group and blue strip indicates the experimental group. SCI, spinal cord injury. (b) Volcano plot of the control (before SCI, n = 6/group) and experimental groups (14 days after SCI, n = 6/group). Orange dots indicate the 925 differentially expressed genes, of which the 720 downregulated genes are to the left and the 205 upregulated genes are to the right of the midline. SCI, spinal cord injury. www.nature.com/scientificreports/ development' , 'positive regulation of transcription, DNA-dependent' , 'response to glucose stimulus' , 'response to nutrient' , 'positive regulation of transcription from RNA polymerase II promoter' , and 'positive regulation of apoptotic process' ( Table 2).
Pathway analysis. There were 55 KEGG pathways with significant enrichment of DEG (p < 0.05). The top 10 KEGG pathways with significant enrichment were 'metabolic pathways' , 'PPAR signalling pathway' , 'peroxisome' , 'fatty acid degradation' , 'fatty acid elongation' , 'glutathione metabolism' , 'propanoate metabolism' , 'valine, leucine and isoleucine degradation' , 'adipocytokine signalling pathway' , and 'glyoxylate and dicarboxylate metabolism' (Table 3).  www.nature.com/scientificreports/ Pathway network analysis. Thirty-one core pathways were identified in pathway network analysis, including 12 downregulated pathways, 1 upregulated pathway, and 18 upregulated/downregulated pathways (Table 4). Figure 2 shows the relationship network of these pathways. Among them, 'citrate cycle (TCA cycle)' , 'pyruvate metabolism' , 'MAPK signalling pathway' , 'fatty acid degradation' , 'propanoate metabolism' , 'apoptosis' , 'focal adhesion' , 'synthesis and degradation of ketone bodies' , 'Wnt signalling pathway' , and 'pathways in cancer' are the 10 pathways with the highest degree values, indicating that they are the most important pathways in the network. The larger the degree value, the more the number of signalling pathways with upstream and downstream effects, and the more important they are in the network.
Gene signal network analysis. Seventy-two hub genes were identified in the gene signal network analysis, of which 13 were upregulated and 59 were downregulated. Figure 3 shows the relationship network of these genes. The 10 hub genes with the highest betweenness values (higher values indicate greater importance) were Acat1, Acadvl, Acaa2, Hadhb, Acss1, Oxct1, Hadha, Hadh, Acaca, and Cpt1b ( Table 5). The network map of these genes is shown in Fig. 4.

Results of quantitative real time PCR (qRT-PCR).
To verify the reliability of the microarray results, we used qRT-PCR to detect the expression levels of 10 genes (Kcnj11, Chmp4c, Opn4, Ppp3cb, Tasp1, Gsr, Stk40, Alpk3, Pla2g12a, and Atf5), and the results were consistent with the microarray data. These 10 genes were selected because they showed the greatest significance in the microarray results. The difference in mRNA expression of these genes between the before SCI group and the after SCI group was significant (P < 0.05) (Fig. 5).

Discussion
Currently, no clear consensus exists regarding the classification of muscle atrophy after SCI. Some reports have suggested that muscle atrophy after SCI be classified as disuse atrophy or denervation atrophy [15][16][17][18] . However, muscle atrophy after SCI differs from these types of muscle atrophy. Disuse atrophy is caused by a significant Table 4. Core pathways identified from pathway network analysis (in order of decreasing degree-value). www.nature.com/scientificreports/ loss of muscle mass and strength due to limb immobilisation, long-term bed rest, lack of physical activity, or space flight 19 . Denervation atrophy is usually induced in peripheral neurotomy models, where the central nervous system is not damaged. The difference between muscle loss caused by SCI and the other types of atrophy is that lower motor neurons remain intact in SCI and upper motor neurons cannot transmit information to the lower motor neurons. Previous studies on muscle atrophy have focused on animal models of disuse atrophy and denervation atrophy; muscle atrophy after SCI is usually classified as either one or a combination of these muscle atrophy models. However, the pathogenesis of muscle atrophy after SCI involves multiple factors, including signal transduction, immunity, electrical conduction, stimulation, and metabolism 20,21 . This complex situation cannot be attributed to a single molecular mechanism or simple combinations. The gene profile data in the present study employed a model of moderate SCI to simulate the clinical phenomenon of muscle atrophy after SCI. Chen et al. 's study 22 used key genes in the co-expression network to classify all SCI samples in GSE45550, including 4 groups (before SCI, 3, 8 and 14 days following SCI). It is believed that the expression levels of the first six genes (CCNB2, CCNB1, CKS2, COL5A1, KIF20A, and RACGAP1) can be used as potential marker genes for different SCI subtypes. Niu et al. 's study 23 used bioinformatics analysis to explore differentially expressed genes associated with the acute and chronic stages of SCI, starting from 1 day to 6 months after SCI, including data sets (GSE45006, GSE93249, and GSE45550). Similarly, Zhang et al. 's research 24 conducted bioinformatics analysis on the data set GSE45550 (including 4 groups: before SCI, 3, 8 and 14 days following SCI) to identify key genes or important signal pathways related to SCI. This study focuses on a part of the data set GSE45550 (two groups: before SCI and 14 days following SCI), and compares the differentially expressed genes of the two groups. The results show that the same disease may have very different results and potential differences at different stages. This reminds us that the same patient needs different treatment adjustments at different stages of the disease.
A previous microarray analysis by Urso et al. 25 was used to describe new gene expression involved in signalling pathways that regulate the loss of muscle proteins or increased response to atrophy stimuli in the first few days after SCI. Microarray analysis revealed that metallothionein activity (MT2A, MT1A, MT1E, MT1F, MT1G, MT1H, MT1M, MT1R, and MT1X) and protease inhibition activity (secretory leukocyte protease inhibitor) increased significantly at 2 and 5 days after SCI 25 . Our microarray analysis revealed that the 10 genes with the most significant differential expression at 14 days after SCI were Kcnj11, Chmp4c, Opn4, Kcnj11, Tasp1, Gsr, Stk40, Alpk3, Pla2g12a, and Atf5. This indicated that the changes in mRNA expression of muscle atrophy at different times after SCI were different. Shen et al. 26 first proposed that innervation atrophy can be divided into four different transcription stages: oxidative stress stage (0-12 h), inflammation stage (24 h), atrophy stage (3-7 days), Figure 2. Network diagram of the interaction between the signalling pathways in the KEGG database based on their upstream and downstream relationships. Red dots indicate significant signalling pathways in which upregulated differentially expressed genes were involved, blue dots indicate the significant signalling pathway in which downregulated differentially expressed genes were involved, and yellow dots indicate significant signalling pathways in which both upregulated and downregulated differentially expressed genes were involved. Solid arrows indicate the upstream and downstream relationships between two signalling pathways, with the upstream signalling pathway at the start of the arrow and the downstream signalling pathway at the end of the arrow. Larger points indicate signalling pathways with a greater number of upstream and downstream effects (higher degree) and thus greater importance in the network. www.nature.com/scientificreports/ and atrophic fibrosis stage (14-28 days). Whether muscle atrophy after SCI involves similar stages remains to be further studied.
The GO analysis showed that the primary biological processes of differential genes were 'response to drug' , 'fatty acid beta-oxidation' , 'response to hypoxia' , 'biological_process' , 'heart development' , 'positive regulation of transcription, DNA-dependent' , 'response to glucose stimulus' , 'response to nutrient' , 'positive regulation of transcription from RNA polymerase II promoter' , and 'positive regulation of apoptotic process' . These results showed that the pathogenesis of muscle atrophy after SCI is the result of mutual gene and network regulation of multiple factors and genes that directly or indirectly lead to the development and progression of muscle atrophy. Among them, 'positive regulation of apoptotic process' warrants further study because skeletal muscle atrophy is closely associated with apoptosis, but the specific mechanisms of action are not well understood 27 . Studies 28,29 have shown that atrophy and death of skeletal muscle cells lead to sarcopenia, which is a disorder associated with normal ageing. It is estimated that 30-40% of skeletal muscle fibres are usually lost by the age of 80 years 30 . www.nature.com/scientificreports/ Although the mechanisms of this loss are not well understood, it is possible that apoptosis is involved. Mitochondrial dysfunction and sarcoplasmic reticulum stress that occur with increasing age may be possible factors that induce apoptosis. Therefore, the mitochondria and sarcoplasmic reticulum may be important organelles in skeletal muscle cells responsible for signal transduction associated with apoptosis. The activation of apoptosis may be a part of the reason for the initiation of muscle protein degradation, loss of muscle nuclei associated with local atrophy, and death of muscle cells. Exercise training and calorie restriction are two interventions known to enhance skeletal muscle function 30 .
In the present study, DEGs were primarily enriched in the pathways 'metabolic pathways' , 'PPAR signalling pathway' , 'peroxisome' , 'fatty acid degradation' , 'fatty acid elongation' , 'glutathione metabolism' , 'propanoate metabolism' , 'valine, leucine and isoleucine degradation' , 'adipocytokine signalling pathway' , and 'glyoxylate and dicarboxylate metabolism' , of which the 'fatty acid degradation' and 'fatty acid elongation' pathways are consistent with the major pathways (fatty acid metabolism) analysed in the study of muscle atrophy after SCI by Dirks et al 30 . These pathways are believed to be the key players in fatty acid metabolism (Lpl and Fabp3) and have proved to be the major sensors of SCI-induced inactivity and reloading with training. Pathway network analysis showed that 'citrate cycle (TCA cycle)' , 'pyruvate metabolism' , 'MAPK signalling pathway' , 'fatty acid  www.nature.com/scientificreports/ degradation' , 'propanoate metabolism' , 'apoptosis' , 'focal adhesion' , 'synthesis and degradation of ketone bodies' , 'Wnt signalling pathway' , and 'pathways in cancer' were the 10 pathways with the highest degree, indicating that these signalling pathways may play the most important roles in muscle atrophy after SCI. Among them, the focal adhesion pathway warrants further attention because some studies 31 have reported that it may play an important role in amyotrophic lateral sclerosis. In addition, studies 32 have reported that the mechanisms and strategies to counter muscle atrophy also include the focal adhesion pathway. The primary sensors for the changes in muscle activity (and load) are the integrin complex and the dystrophin-associated protein complex. Neither of these complexes has inherent signalling function, but they appear to mediate focal adhesion through focal adhesion kinase (FAK) 33,34 , a targeted kinase that indirectly provides nutrients through many different pathways, thereby leading to the activation of growth pathways and preventing apoptosis. Total FAK and its activity were reduced in atrophic soleus muscle, and other proteins associated with focal adhesion complexes were also downregulated 35 .
In the gene signal network analysis, Acat1, Acadvl, Acaa2, Hadhb, Acss1, Oxct1, Hadha, Hadh, Acaca, and Cpt1b may play key roles in the development and progression of muscle atrophy after SCI because they had the highest betweenness values (higher betweenness indicates more connections with a gene, suggesting that the gene may be more important). The results also showed the relationship between these 10 genes, providing new avenues for the prevention and treatment of muscle atrophy after SCI. Among them, the betweenness value of Acadvl was the highest. Acadvl is involved in the fatty acid degradation pathway 36 and the fatty acid beta-oxidation biological process 37 . However, the fatty acid degradation pathway was screened from the pathway analysis and pathway network analysis, and fatty acid beta-oxidation was also screened from the GO analysis. These showed the interdependency of the results; they also suggest that Acadvl may be a key regulatory gene in the development of muscle atrophy after SCI, warranting further study. ACADVL is responsible for the production of very long-chain acyl-CoA dehydrogenase, an enzyme that plays a key role in the mitochondria. It is essential for fatty acid oxidation, a multi-stage process that metabolises fat into energy 38 . There is a close relationship between abnormalities associated with neurodegenerative diseases such as Alzheimer's disease and brain fatty acid metabolism. The pathology associated with Alzheimer's disease inhibits the homeostasis and regeneration of neural stem cells by interfering with fatty acid metabolism 39,40 .
The sample size of the data in the present study was small and the samples were selected from a platform, which may lead to some degree of false negatives. Further experimental studies and a larger sample size are needed to confirm the results of the present study. A key limitation of the original study design is that, while the SCI group underwent a laminectomy and injury, the control group did not undergo 'sham' surgery. Thus, the study was not controlled for surgery effects.
In summary, the genetic profile data of atrophic muscle tissue before SCI and at 14 days after moderate SCI was downloaded from the Gene Expression Omnibus to identify genes differentially expressed between the two conditions. In addition, bioinformatic analyses were performed to gain new insights into the molecular mechanisms of muscle atrophy in this model. It is still not clear whether muscle atrophy 14 days after SCI in this model properly represents muscle atrophy in the atrophic fibrosis stage. Nevertheless, our current understanding of the mechanisms underlying muscle atrophy after SCI is undoubtedly expanded, and several potential targets for the prevention and treatment of this condition were also revealed.

Materials and methods
Gene profile data. The keywords 'SCI' and 'muscle atrophy' were used to search the Gene Expression Omnibus in the National Centre for Biotechnology Information (NCBI) database (http:// www. ncbi. nlm. nih. gov/ geo/) and the GSE45550 gene profile dataset was obtained. These data represented a control group (before  Kcnj11, Chmp4c, Opn4, Ppp3cb, Tasp1, Gsr, Stk40, Alpk3, Pla2g12a, and Atf5 expression in before SCI (n = 6) and after SCI (n = 6) were evaluated by qPCR and normalised against the corresponding glyceraldehyde-3-phosphate dehydrogenase (GAPDH) expression. An asterisk represents p < 0.05 and two asterisks indicate p < 0.01 when compared with the before SCI group. Design of analysis process. The GSE45550 gene profile dataset was input into the Gene-Cloud of Biotechnology Information (GCBI) analysis platform (https:// www. gcbi. com. cn/ gclib/ html/ index) to perform all the bioinformatics analyses in this study. The details of this process are shown as a flowchart (Fig. 6). In the GCBI, complex biological information analysis can be completed by dragging the 'Move' module and clicking 'Operation' . Several bioinformatic studies have used the GCBI to discover new genes and molecular mechanisms in diseases and improve the understanding of the underlying molecular mechanisms.

Identification of DEGs.
For DEG identification, significance analysis of microarrays (SAM) 41 was used to screen for genes in previously designated groups with significant differential expression. First, SAM was used to calculate a statistic (d score) for each gene to measure the degree of correlation between the gene expression level under the corresponding grouping and the designated group, and then simulated the distribution of these statistics randomly (randomly replaced samples) through a series of 1,000 replacements. The statistic value (d Score) was calculated for each gene di = ri/(si + s0), where ri reflects the difference in the average level among different groups and si reflects the variation in sample population. Considering the influence of the variance of low-abundance genes and the correlation between genes, the statistical calculations, and the methods for evaluating statistical significance were adjusted using appropriate methods to obtain the significance level for each gene. Finally, to adjust for the p values of multiple experiments, a q value was calculated by constructing the discretised rejection region to control the FDR 42,43 . Differences with q < 0.05 with fold change > 1.5 or < − 1.5 were considered statistically significant.

GO analysis.
The GO database is a cross-species, comprehensive, and descriptive database established by the Gene Ontology Consortium. In the GO analysis, the DEG between the groups was used to perform gene function annotation based on the GO database to obtain all the functions in which the gene is involved. Next, Fisher's exact test and multiple comparison tests were used to calculate the significance level (p value) and FDR in order to screen the significant functions of each differentially expressed gene. The enrichment results with FDR < 0.05 were considered statistically significant.
Pathway analysis. The Kyoto Encyclopaedia of Genes and Genomes (KEGG) is a database that systematically analyses the relationship between genes (and their coded products), gene functions, and genomic information, and allows the study of genes and their expression as an entire network. Based on the KEGG database analysis, Fisher's exact test was used for the differential genes to perform significance analysis of the pathways in Figure 6. Flowchart of the bioinformatic analysis in the present study. SCI, spinal cord injury. DEGs, differentially expressed genes. n = 6/group. www.nature.com/scientificreports/ which the target gene is involved in order to screen the significant pathways of each differential gene. Pathways with FDR < 0.05 were considered statistically significant.
Pathway network analysis. Pathway network analysis allows comprehensive, systematic analysis of the signal transduction relationships among significant pathways discovered from pathway analysis and intuitive discovery of the synergistic effects of major pathways when the sample changes, thereby providing a systematic understanding of the nature of the changes in sample traits. The pathway network analysis was used to construct an interaction network between pathways based on the interaction relationships in the KEGG database. The network was determined based only on the significant pathways from the pathway analysis. The calculation and description of network attribute: network N is denoted as N := (V, E) , where V is the node set whose totality is n, and E is the edge set whose totality is m. The corresponding adjacent matrix is noted by A := a ij , i, j ∈ V , where a ij = 1 if i connected to j, or else a ij = 0. The degree of node i:D(i) := n j=1 a ij . Indegree values: the number of upstream pathways of a certain pathway; outdegree values: the number of downstream pathways of a certain pathway; degree values: the number of upstream and downstream pathways of a certain pathway; core pathways: degree values ≥ 1.
Gene signalling network analysis. The interaction between genes and gene products in the KEGG database and the interaction between each gene and other genes obtained from searching the database allow the interaction between the target gene clusters to be comprehensively discovered and positioned as upstream and downstream proteins. Next, an interaction network between genes was constructed. The network based on the DEGs. The betweenness of i:B(i) = s� =i� =t σ st (i) σ st , where σ st denotes the total number of the shortest path from node s to node t, σ st (i) denotes the total number of the shortest paths through node i from node s to node t. Betweenness values: signal transmission mediation centre, if the value is larger, the mediation ability of a certain gene in signal transmission is stronger. Indegree values: the number of upstream genes of a certain gene; outdegree values: the number of downstream genes of a certain gene; degree values: the number of upstream and downstream genes of a certain gene; hub genes: degree values ≥ 1.
Validation of the mRNAs using qRT-PCR. Twelve adult female Spragdorley rats (16 weeks of age, 260-280 g at the beginning of the study, Slaccas Laboratory, Shanghai, China) were reared at a temperature of 22 ± 1 °C, humidity of 50% ± 10%, and light: dark cycle of 12:12, and provided a sufficient food and water. Animals were sacrificed: a control group (before SCI, n = 6) and an experimental group (14 days after SCI, n = 6). The SCI model 7 was established by anaesthetising the rat and dropping a 10 g cylinder from a height of 25 mm at the level of T7-T9 after a dorsal laminectomy, affecting the spinal cord at T8 to cause injury. The animals were deeply anaesthetised with a combination of ketamine (90 mg/kg body weight) and xylazine (8 mg/kg body weight). The sample tissue was soleus muscle tissue. Beyozol (Beyotime Bio, Inc., China) was used to extract total RNA from the frozen tissues. The BeyoRT II cDNA Synthesis Kit (Beyotime Bio, Inc., China) was used to reverse transcribe RNA into cDNA. Thereafter, 2 µg of cDNA was tested in each reaction using the BeyoFast SYBR Green One-Step qRT-PCR Kit (Beyotime Bio, Inc., China) in the Applied Biosystems Real-Time PCR System (Applied Biosystems; Thermo Fisher Scientific). The 2 −ΔΔCt method was adopted to calculate the expression of genes relative to the housekeeping gene GAPDH. Table 6 shows the primers applied to qRT-PCR. The Ethics Committee of Hainan General Hospital of China approved the experiment (Approval No: Med-Eth-Re [2020] 12).