The Mechanism Of Functional Electrical Stimulation Ameliorating Muscle Function Of Patients With Sarcopenia After Spinal Cord Injury: A Bioinformatics Analysis Of Key Genes

Background: Sarcopenia caused by spinal cord injury seriously affects the muscle function, which impairs the locomotion function of patients. As an effective physiotherapeutic, functional electrical stimulation is benet to the recovery of muscle function without exact mechanisms. So the objection of the study is to explore the biological regulatory factors related to functional electrical stimulation for muscle function recovery of patients with sarcopenia after spinal cord injury using bioinformatics methods. Methods: The related microarray datasets GSE142426 and GSE33886 were downloaded from the Gene Expression Omnibus database. We used the R software to merge the two datasets, correct the inter-batch differences and screen the differentially expressed genes. Gene ontology enrichment analysis and Kyoto Encyclopedia of Genes and Genomes pathway were analyzed by using the DAVID online tool. The STRING database was used to analyze the interaction of differentially encoded proteins. The key genes were selected to draw the ROC curves, and CYCS was used to perform gene set enrichment analysis. Results: A total of 114 differentially expressed genes were selected, including 44 up-regulated genes and 70 down-regulated genes, and four hub genes were identied including CYCS, SUCLG1, ATP5B and ATP5C1. ROC curve showed that CYCS was considered as the best indicator and GSEA analysis showed that up-regulation of CYCS was related to mitochondria and energy metabolism. Conclusion: The mechanism of function electrical stimulation on muscle function recovery of spinal cord injury patients with sarcopenia is mainly related to regulate mitochondrial energy metabolism and scavenge reactive oxygen species to mitigate the oxidative damage. SCI: spinal cord injury; CSA: cross-sectional area; FES: functional electrical stimulation; GEO: gene expression omnibus, DEGs: differential expressed genes; GO: gene ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; PPI: protein-protein interaction; STRING: search tool for the retrieval of interacting genes; ROC: receiver operating characteristic; AUC: area under curve; GSEA: gene set enrichment analysis; TCA: tricarboxylic acid; adenosine triphosphate: ATP.


Introduction
Spinal cord injury (SCI) means that damage to the spinal cord with temporarily or permanently in its functions. SCI with a high mortality and disability rate cause great damage to the survivors who lose the ability to work and live independently, meanwhile it can bring enormous burden to the physical, psychological, economic and family of the suffers. More than 27 million patients are suffering SCI around the word, and the incidence is rising year on year (1). Sarcopenia is a progressive and systemic skeletal muscle disorders with the accelerated loss of muscle mass and function and which can increased the risk of falls, fractures, disability and mortality (2,3). Restricted mobility in SCI patients can cause the reduction of the cross-sectional area (CSA) of somatic muscles and increase the intramuscular lipid content, leading to the development of sarcopenia, and the rate of sarcopenia and weakness as high as 31.1% (4). Subsequently, SCI can induce impaired glucose tolerance, insulin resistance, and lipid abnormalities (5). A wealth of studies con rms that the functional electrical stimulation (FES) is contribute to the increase of muscle CSA, muscle strength and down-regulated the bone turnover, which is bene t to the recovery of muscle function(6-8). So far, the mechanism of FES treating SCI remains incompletely clear. Some studies suggest it is associated with mitochondrial function and energy metabolism. With the development of next-generation sequencing technologies and the improvement of biological database, using bioinformatics methods to explore the relevant mechanisms is signi cant.

Microarray studies and datasets from Gene Expression Omnibus (GEO)
The microarray datasets including GSE142426 and GSE33886 were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) using "spinal cord injury" as the keyword. The microarray dataset GSE142426 from GPL5175 platform contains 15 samples of male patients with SCI before and after treated with FES. The microarray dataset GSE33886 from GPL6244 platform contains 3 samples of male patients with SCI before and after treated with FES. Samples were selected by hierarchical clustering and all the samples were taken from quadriceps femoris or soleus. Differential genes expression analysis Statistical programming language R (version 4.0.2) was used for log2 transformation of the data, and the two datasets were merged. "SVA" package was used for batch correction. Differential expressed genes (DEGs) were de ned as log |FC| > 0.5, and the corrected p < 0.05. Log |FC| > 0 means that the DEG is up regulated after FES treatment. Functional annotation and pathway analysis of DEGs DEGs were inputted into David 6.8 online tools(https://david.ncifcrf.gov/)to perform Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment. P < 0.05 and the gene counts > 3 were considered statistically signi cant. Protein-protein interaction (PPI) network and key genes acquisition Using the Search Tool for the Retrieval of Interacting Genes (STRING, version 11.0, https://string-db.org/) database to analyze the PPI of proteins encoded by DEGs (medium con dence = 0.04). Cytoscape software (version 3.8.0) was used to perform visualization of PPI network. We used cytohubba plug-in to analyze the nodes of the genes with topological analysis methods, ltering with degree and stress and obtaining the key genes from the intersection of the rst 15 genes sorted out in degree and stress methods respectively. Receiver operating characteristic (ROC) curve analysis Using the Graphpad prism (version 8.2.1) software to draw the ROC curve and the area under curve (AUC) was used to evaluate the sensitivity and speci city of candidates with the AUC more than 0.7. According to the Youden index, the best cut-off value was selected and P < 0.05 means the difference was statistically signi cant. Gene set enrichment analysis (GSEA) We choose the CYCS gene which is highest ranked in degree method and statistically signi cant in ROC curve to perform GSEA analysis. The gene expressions were divided into low expression group and high expression group according to the median. We analyzed the KEGG pathway related to CYCS high expressions ltering with P 0.05, NES 1 and FDR 25% and using David analysis method to obtain the key pathway from the intersection of the two analyses.

Results
Hierarchical clustering for sample selection The total samples were analyzed by hierarchical clustering, and 2 samples (GSM4227252 and GSM4227238) with high heterogeneity were eliminated. Finally, 17 samples were included before and after FES. Identi cation of DEGs A total of 114 DEGs were screened out, including 44 high expression genes and 70 low expression genes. R was used to make results visualization and draw volcano map and heat map ( Fig. 1).

GO and KEGG enrichment analysis of DEGs
In GO analysis, DEGs were divided into three categories: biological process, cellular component and molecular function. In the biological process, DEGs were mainly concentrated in the elds of mitochondrial electrical transport, NADH to ubiquinone, mitochondrial respiratory chain complex I assembly, mitochondrial adenosine triphosphate (ATP) synthesis coupled proton transport, tricarboxylic acid cycle, response to reactive oxygen species, etc. In the cellular components, DEGs were mainly concentrated in the parts of mitochondrial inner membrane, mitochondrion, mitochondrial matrix, mitochondrial respiratory chain complex I, extracellular exosome and so on. In molecular functions, DEGs were mainly concentrated in poly (A) RNA binding, NADH dehydrogenase activity, electron carrier activity, NADH dehydrogenase (ubiquinone) activity, proton-transporting ATP synthase activity, rotational mechanism and so on. The above-mentioned results show that the difference genes mainly focus on mitochondrial function and energy metabolism. KEGG pathway analysis shows that DEGs mainly focuses on pathways about Alzheimer's disease, Parkinson's disease, Oxidative phosphorylation, Huntington's disease, Non-alcoholic fatty liver disease (NAFLD), Metabolic pathways, tricarboxylic acid (TCA) cycle, Carbon metabolism and so on, which can be seen as the main pathways involved in nervous system diseases and energy metabolism (Fig. 2).
PPI network analysis of DEGs A total of 113 nodes and 363 pairs of protein interaction relationships were obtained. After eliminating independent genes, Cytoscape was used for visualization of PPI network map and 4 key genes including CYCS, SUCLG1, ATP5B and ATP5C1 were obtained from the intersection of Wayne map (Fig. 3).
Evaluation of diagnostic value of hub gene with ROC curve analysis ROC curve was drawn for 4 hub genes to evaluate the effectiveness of hub genes of SCI patients after FES. CYCS (AUC = 0.7578, P = 0.0103) was considered as the best indicator (Table 1, Fig. 4). ( Fig. 5).

Discussion
As a serious central nervous system injury, the impact of SCI on patients is inestimable. The effects of FES, an important adjuvant therapy on SCI, on Sarcopenia and weakness in SCI patients have been con rmed by many studies. But the exact mechanism of the effects of FES is not known. It has been reported in previous studies that mitochondria and energy metabolism play an important role in this process, which has been con rmed again in our study.
Through PPI network, we identi ed four genes including CYCS, SUCLG1, ATP5B and ATP5C1, which may be the key genes for FES to improve muscle function in patients with SCI. Around these genes, CYCS has been suggested to be closely related to the improvement of muscle function. Previous studies have con rmed that the reduction of exercise after SCI can reduce the ratio of fusion proteins and lead to mitochondrial dysfunction(9, 10), while exercise training can improve the volume, quantity, density, dynamics and quality control of skeletal muscle mitochondrial (11,12).
As an indispensable component in the oxidative respiratory chain, CYCS is the core component of mitochondrial electron transport chain and the proteins encoded by CYCS combines with the mitochondrial inner membrane to act as the electronic carrier of complexes to complex in the electron transport chain (13). Oxidative stress caused by long-term inactivity plays an important role in muscle atrophy, which can reduce mitochondrial respiratory function. The main source of oxidants is reactive oxygen species (ROS) produced by mitochondria and high concentration of ROS can damage fast and slow muscle bers and mitochondrial function, and lead to muscle ber atrophy (14). CYCS can clear ROS and reduce the damage caused by oxidative stress to improve muscle function (15,16). At the same time, exercise can increase the abundance of mitochondrial RNA and the ability of mitochondria on producing ATP (17). SUCLG1 is mainly composed of the catalytic enzyme of mitochondrial energy generation. ATP5B and ATP5C1 are mainly involved in the subunit of ATP syntheses in mitochondria, and the three genes participate in promoting energy generation. Therefore, FES may improve muscle quality and strength by promoting the production of the above four genes.
We analyzed the CYCS gene that is highest ranked in degree method by GSEA combined with the enrichment of KEGG pathway in DAVID analysis. The results showed that the up regulation of CYCS was related to myocardial contraction, Parkinson's disease, Alzheimer's disease, oxidative phosphorylation and TCA cycle. It is con rmed that CYCS is closely related to energy metabolism such as oxidative phosphorylation, TCA cycle and so on. 95% of the energy of cell activity comes from ATP produced by oxidative phosphorylation (18). Oxidative phosphorylation is composed of electron transport chain and ATP synthase. The former generates mitochondrial proton power by pumping protons across the inner mitochondrial membrane and the latter returns protons to the matrix to couple with ADP and phosphate to ATP (18). One of the main protein control the signal of mitochondria oxidative phosphorylation is cytochrome c, and the function of which was regulated by phosphorylation of several residues(18, 19).
So, as the nal step of electron transport chain, CYCS is an important place for the oxidation phosphorylation of electron transport chain. The TCA cycle is the nal metabolic pathway of the three major nutrients as lipids, carbohydrates, and amino acids, and is the main site for energy metabolism in the mitochondria. So, as the metabolic link of three substances, TCA cycle is the key link of energy metabolism and the main place is mitochondria. Therefore, all the above evidence supports that FES is associated with mitochondria and energy metabolism.
The results of GO enrichment analysis suggests that the DEGs mainly focus on mitochondrial electron transport, mitochondrial respiratory chain complex I assembly, mitochondrial ATP synthesis coupled proton transport, TCA cycle, response to reactive oxygen species and so on. These results are consistent with our differential analysis and KEGG pathway, which mainly focus on mitochondrial energy metabolism.
As muscle samples are more di cult to obtain, bioinformatics methods were used for DEGs analysis as a favorable tool to explore the mechanism of FES on muscle function recovery of SCI patients. But due to the small sample size of these two datasets, some causative genes may not show signi cant differences and we selected a smaller difference. As the number of datasets in the GEO database increasing and relevant techniques and algorithm improved, the mechanisms about FES on treating SCI patients may be further explored.

Conclusions
The mechanism of FES on muscle function recovery of spinal cord injury patients with sarcopenia is mainly related to regulate mitochondrial energy metabolism and scavenge reactive oxygen species to mitigate the oxidative damage. Authors' contributions: KCL designed the study and prepared the manuscript. LXW and LHW searched database and extracted datasets. XF and HLY as the co-corresponding authors designed the study and analyzed the datasets. Figure 1 Volcano plot, cluster dendrogram and heat map of differentially expressed genes (A) The volcano showed that differentially expression genes identi ed from both GSE142426 and GSE33886. Red dots indicate up-regulated genes screened according to |log FC| > 0.5 and P-value < 0.05, and green dots indicate downregulated genes screened according to |log FC| > 0.5 and P-value < 0.05. Black dots represent expression level of genes without statistically signi cant. (B) Heat maps show the top 100 differentially expression genes, the horizontal ordinate represent the names of the samples in the microarray and the vertical ordinate represent the names of the differential genes. Red represents high expression, green represents low expression and black represents non-signi cant relative expression differences.   The ROC curve of hub genes