Identification of key pathways and RNAs associated with skeletal muscle atrophy after spinal cord injury.

OBJECTIVE: This study was performed to investigate the potential key molecules involved in the progression of skeletal muscle atrophy after SCI.
METHODS: Based on GSE21497 dataset, the DEmRNAs and DElncRNAs were screened after differentially expressed analysis. Then the enrichment analyses were performed on DEmRNAs. Then the PPI network and ceRNA network were constructed. Finally, the DGIdb was utilized to predict drug-gene interactions.
RESULTS: A total of 412 DEmRNAs and 21 DElncRNAs were obtained. The DEmRNAs were significantly enriched in MAPK signaling pathway and FoxO signaling pathway. In addition, UBE2D1, JUN, and FBXO32 had higher node degrees in PPI network, and the top 20 genes with high degree were significantly enriched in FoxO signaling pathway and Endometrial cancer. Moreover, FOXO3 was regulated by hsa-miR-1207-5p and hsa-miR-1207-5p was regulated by lncRNA RP11-253E3.3 in ceRNA network. Finally, 37 drug-gene interactions were obtained based on the 26 genes in ceRNA network.
CONCLUSION: UBE2D1, JUN, and FBXO32 are likely to be related to the progression of skeletal muscle atrophy after SCI, and activating of MAPK signaling pathway, Endometrial cancer and FoxO signaling pathway may induce skeletal muscle inflammation, apoptosis, autophagy and atrophy after SCI. Moreover, RP11-253E3.3-hsa-miR-1207-5p-FOXO3 axis may be a promising therapeutic target for skeletal muscle atrophy after SCI.

Introduction

Spinal cord injury (SCI) is the spinal cord damaged by direct or indirect external force, and spinal fractures and dislocations occur, causing with complete or incomplete paraplegia below the plane of injury[1-3]. The pathophysiological changes of SCI involve multiple systems such as the nervous system, immune system and vascular system[4,5]. Moreover, skeletal muscle will atrophy rapidly after SCI[6]. The report according to the World Health Organization, reveals that about 250,000 to 500,000 people worldwide suffer from SCI each year, and the reason mainly due to car accidents, falls and violence, which seriously affect the quality of human life of this population[7,8]. However, there is no effective diagnosis strategy to repair secondary damage after SCI.

In recent years, bioinformatics analysis can not only enable scholars to discover valuable genes in high-throughput chip standardized data, save a lot of time and energy, but also discover new genes and pathways in the early stage, providing theoretical basis and direction for further research[9]. Notably, noncoding RNA transcripts such as microRNAs (miRNAs) and long noncoding RNAs (lncRNAs) have been reported to play significant roles in the molecular mechanism of SCI. For instance, Zhang et al. revealed that downregulation of miRNA-127-5p aggravates SCI via activating MAPK1[10]. Sun et al. found that miRNA-411 attenuates inflammatory damage and apoptosis following SCI[11]. Sabirzhanov et al. have shown that inhibition of microRNA-711 can limit the changes of angiopoietin-1 and Akt, tissue damage, and motor dysfunction after contusive SCI in mice[12].

GSE21497 was analyzed by Reich et al. to find the global gene expression patterns following short-term unloading (48 h UL) and reloading (24 h RL) in human skeletal muscle[13]. However, the molecular mechanisms underlying the development of skeletal muscle atrophy after SCI have not yet been understood. In this study, a microarray dataset of SCI was downloaded from public database. Then, we performed differentially expressed analysis of lncRNAs and mRNAs, and then the protein-protein interaction (PPI) network and module analyses were bulit for differentially expressed mRNAs (DEmRNAs). According to the lncRNA-mRNA, mRNA-miRNA, and lncRNA-miRNA interaction pairs, the mRNA in ceRNA network was constructed for functional enrichment analysis. Finally, the Drug-Gene Interaction database (DGIdb) was utilized to predict drug-gene interactions. This study aims to provide better understanding and promising therapeutic targets for skeletal muscle atrophy after SCI.

Materials and Methods

Data Source

The GSE21497 dataset, containing muscle biopsies from the vastus lateralis muscles of the SCI patients two days and five days post-SCI, was obtained from NCBI Gene Expression Omnibus[14] (GEO, https://www.ncbi.nlm.nih.gov/gds/?term=). All samples were detected through the GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array platform.

Data preprocessing and differential expression analysis

The data were preprocessed utilizing the affy package[15], and the processes of data preprocessing contained background correction, normalization, and concentration prediction. Probes were annotated by matrix data combined with chip platform annotation file. The average value of diverse probes would be considered as the eventual expression level of gene if they were corresponded to the identical mRNA. For lncRNA annotation, genes whose information were “antisense”, “sense_intronic”, “lincRNA”, “sense_overlapping” and “processed_transcript” were selected as lncRNA according to the GENCODE for human gene annotation files (Release 29).

A principal component analysis (PCA) of samples was performed in this study utilizing prcomp algorithm in R language[16]. In addition, the classical Bayesian testing method and Benjamini and Hochberg (BH) method was carried out to adjust P value for multiple comparisons. The differentially expressed mRNAs (DEmRNAs) and differentially expressed lncRNAs (DElncRNAs) between two days and five days post-SCI groups were screened with the threshold of adj. P<0.05 and |log fold change (FC)| >0.585. Finally, the pheatmap in R package[17] (Version: 1.0.10, https://cran.r-project.org/web/packages/pheatmap/index.html) was used to draw the heatmap of DEmRNA and DElncRNA.

Enrichment analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway[18] enrichment analyses were conducted on DEmRNAs utilizing the clusterprofiler in R package[19] (Version 2.4.3, http://bioconductor.org/packages/3.2/bioc/html/clusterProfiler.html). Genes count ≥2 and P<0.05 was used as the threshold of significant enrichment results. Besides, the Gene set enrichment analysis (GSEA) was employed to perform the enrichment analysis on these DEmRNAs. The nominal p value and normalized enrichment score (NES) were used to rank the pathways enriched in each phenotype.

Protein protein interaction (PPI) network construction and module analysis

The Search Tool for the Retrieval of Interacting Genes (STRING) database[20] (version 10.0, http://www.string-db.org/) was carried out to analyze the interactions between protein and protein encoded by DEmRANs. The PPI score was set as 0.7 (referred as high confidence). Afterwards, the PPI network was constructed by the Cytoscape software[21] (Vesrion: 3.2.0, http://www.cytoscape.org/).

The MCODE plug-in in Cytoscape software[22] (Version1.4.2, http://apps.cytoscape.org/apps/MCODE) was employed to analyze the most significantly clustered module in PPI network with the threshold value of score ≥5. Also, the enrichment analyses were performed on the genes with high degree (top 20) and the significantly clustered module genes.

Construction of lncRNA-miRNA- mRNA composite network

To further obtain the matrix data of DElncRNA and the key genes in the previous step, and the correlation coefficient of lncRNA and mRNA was calculated. Moreover, the lncRNA-mRNA with r>0.6 and P<0.05 were screened, which was considered as lncRNA-mRNA co-expression positive correlation with high significant correlation.

To obtain the relationship between DElncRNA and miRNA, the fasta format sequence file of the lncRNA in the co-expression relationship and all the mature bady miRNA fasta sequence files were screened. The miRanda (Version: 3.3a, https://omictools.com)/miranda-tool) was utilized to predict the miRNA-lncRNA relationship pairs with score >170 and energy <-30.

Targets of miRNAs in the lncRNA-miRNA were predicted by use of the miRWalk 2.0[23] (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/). In the microRNA information retrieval system, the search conditions included Minimum seed length=7, P<0.05, and Input parameters (3’UTR). Then, the miRNA-target interaction pairs were predicted in six databases, including TargetScan (http://www.targetscan.org/), miRWalk, miRanda (https://omictools.com/miranda-tool), miRMap (https://omictools.com/mirnamap-tool), RNA22 (http://cm.jefferson.edu/rna22/), RNAhybrid (https://omictools.com/rnahybrid-tool). Only these miRNA-target interaction pairs screened more than five databases were identified for the subsequent analysis. Moreover, the lncRNA-mRNA, miRNA-lncRNA and miRNA-mRNA interaction pairs were integrated , and then the lncRNA-miRNA-mRNA ceRNA network was bulit utilizing Cytoscape. Finally, the genes in ceRNA network were performed GO and KEGG pathway analysis.

Drug-gene interaction predictive analysis

The Drug-Gene Interaction database (DGIdb) can be used for excavating existing resources and hypothesis of how genes are targeted for treatment. The module genes were predicted by use of DGIdb 2.0[24] (http://www.dgidb.org/). The drug database default information was as follows: the interaction type removed the NA; and predicted all the drug-gene interactions associated with the genes in the composite network, and the Cytoscape was used to bulid the drug-gene network.

Results

Differential expression analysis

As presented in Figure 1, SCI patients two days and five days post-SCI samples were completely separated, meaning that the expression patterns of samples were specific and could be used to completely distinguish between two days and five days post-SCI samples.

Figure 1

Principal component analysis (PCA) diagram of the samples. (A) The PCA diagram of before deleting the samples; (B) The PCA diagram of after deleting the samples. After standardization, a PCA analysis on the data was performed and the results found that there was a deviation in the samples. Therefore, the information of 4 samples (including 2 patients related to day 2 and day 5) were deleted. And hence, the samples analyzed was day 5 (n=8) vs. day 2 (n=8).

Based on the cutoff value of adj. P<0.05 and |log fold change (FC)| >0.585, a total of 412 DEmRNAs (of which 155 and 257 were up- and down-regulated, respectively) and 21 DElncRNAs (of which 15 and 6 were up- and down-regulated, respectively) were obtained. The heatmap of the DEmRNAs and DElncRNAs were shown in Figure 2, and the results indicated that the two groups could be significantly separated, meaning that the difference analysis were reliable.

Figure 2

Heatmap of differentially expressed mRNA (mRNAs), differentially expressed lncRNA (DElncRNA). (A) DEmRNA; (B) DElncRNA.

Enrichment analysis of DEmRNAs

Functional enrichment analyses of the GO terms and KEGG pathway were performed for both DEmRNAs. The 155 up-regulated DEmRNAs were mainly enriched in 33 KEGG pathways (including hsa04010: MAPK signaling pathway) and 456 GO-biological process (BPs) (including GO:0007519- skeletal muscle tissue development, GO:0060538- skeletal muscle organ development, etc.). In addition, 257 down-regulated DEmRNAs were mainly enriched in 44 KEGG pathways (including hsa04068: FoxO signaling pathway) and 501 BPs (GO:0003012- muscle system process, GO:0006936- muscle contraction, etc). The top 10 terms of GO and KEGG pathway analyses of DEmRNAs were illustrated in Figure 3. Besides, the GSEA results (Supplementary files) shown that the positive pathways involved in DEmRNAs were included Ribosome, Intestinal immune network for IgA production, Graft versus host disease, etc. The negative pathways involved in DEmRNAs were contained PPAR signaling pathway, Cardiac muscle contraction, adipocytokine signaling pathway, etc.

Figure 3

Enrichment analysis of differentially expressed mRNA (mRNAs). (A) KEGG pathway analysis of up-regulated DEmRNAs; (B) KEGG pathway analysis of down-regulated DEmRNAs; (C) GO-BP terms of up-regulated DEmRNAs; (D) GO-BP terms of down-regulated DEmRNAs. BP, biological process; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology.

PPI network and subnet module

To analyze the interactions between protein and protein encoded by mutant genes, the PPI network was conducted. Totally, 185 nodes and 346 interactions was shown in the PPI network (Figure 4A). Moreover, ubiquitin conjugation enzyme E2 D1 (UBE2D1), Jun proto-oncogene, Ap-1 transcription factor subunit (JUN), and F-box protein 32 (FBXO32) had higher node degrees (Table1). In addition, based the aggregation of Cytoscape plug-in MCODE (score ≥5), two sub-network modules of the PPI network were obtained (Figure 4B). Briefly, the module A (score=13) included 15 nodes and 91 interactions, and module B (score=5) included 5 nodes and 10 interactions. KEGG pathway and GO-BP enrichment analyses were also performed on genes with high degree (top 20) and module genes (Figure 4C). The top 20 genes with high degree were significantly enriched in 58 KEGG pathways (including hsa04068: FoxO signaling pathway and hsa05213: Endometrial cancer) and 287 GO-BP terms. Module A genes were mainly enriched in 2 KEGG pathways and 22 GO-BP terms, and module B genes were mainly enriched in 15 GO-BP terms.

Figure 4

Protein protein interaction (PPI) network analysis of differentially expressed mRNA (mRNAs). (A) PPI network. (B) Two sub-network modules of the PPI network; (C) KEGG pathway and GO-BP enrichment analyses of genes with high degree. Yellow circular nodes represent up-regulated mRNA; Green prismatic nodes represent down-regulated mRNA. The size of the nodes represents the degree value. KEGG, Kyoto Encyclopedia of Genes and Genomes; PPI: protein protein interaction.

Table 1

Top 20 nodes with high degree in protein protein interaction (PPI) network.

Nodes	Description	Degree
UBE2D1	DOWN-gene	19
JUN	UP-gene	17
FBXO32	DOWN-gene	17
FBXO30	UP-gene	14
FBXO40	UP-gene	14
ASB2	UP-gene	14
ASB5	UP-gene	14
ZBTB16	DOWN-gene	14
SPSB1	DOWN-gene	14
PIK3R1	DOWN-gene	13
TRIM63	DOWN-gene	13
EGF	UP-gene	12
TRIM32	UP-gene	11
UBE2G1	UP-gene	11
DCUN1D1	DOWN-gene	11
LONRF1	DOWN-gene	11
FOXO3	DOWN-gene	10
ACTN2	DOWN-gene	10
RBBP7	DOWN-gene	10
HIST1H2AC	UP-gene	10

Construction of lncRNA-miRNA-RNA composite network

To better understand the relationships between differentially expressed lncRNAs, miRNAs, and mRNAs in SCI, we constructed a ceRNA network. According to the integration of the top 20 genes with high degree and two significant clustering module genes, 27 key genes were obtained. Then, after the co-expression analysis was conducted between 27 key genes and 21 DElncRNAs, 76 lncRNA-mRNA positive correlated interactions pairs (including 21 lncRNAs and 26 DEmRNAs) were obtained. In addition, 31 miRNA-lncRNA interactions pairs (including 30 miRNAs) were required based on the cutoff value of score >170 and energy <-30. Then, the target genes of 30 miRNAs were predicted, which were intersected with the 26 DEmRNAs in lncRNA-mRNA positive correlated interactions pairs, and then 45 miRNA-mRNA interactions pairs were acquired. According to the 76 lncRNA-mRNA interactions pairs, 31 miRNA-lncRNA interactions pairs, and 45 miRNA-mRNA interactions pairs, the lncRNA-miRNA-mRNA ceRNA network was constructed (Figure 5 A). Here, forkhead box O3 (FOXO3) was regulated by hsa-miR-1207-5p and hsa-miR-1207-5p was regulated by lncRNA RP11-253E3.3. Also, the enrichment analysis was conducted on the 26 genes in ceRNA network, and the results shown that the 26 genes were mainly enriched in 61 KEGG pathways (including hsa01521: EGFR tyrosine kinase inhibitor resistance and hsa05213: Endometrial cancer) and 196 GO-BP terms, and the top 10 terms of KEGG pathway and GO analyses were illustrated in Figure 5 B and C.

Figure 5

Construction of ceRNA network. (A) lncRNA-miRNA-mRNA ceRNA network; (B) KEGG pathway analysis of 26 genes in ceRNA network; (C) GO-BP terms of 26 genes in ceRNA network. BP, biological process; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology.

Prediction of drug-gene interactions

According to the DGIdb prediction result of 26 genes in ceRNA network, 37 drug-gene interactions were obtained, among which many interactions were predicted to be inhibitor (Table 2).

Table 2

Drug-gene Interaction.

Gene	Drug	Interaction	Gene	Drug	Interaction
JUN	ARSENIC TRIOXIDE	inducer	PIK3R1	TASELISIB	inhibitor
JUN	CHEMBL3222137	inhibitor	PIK3R1	PILARALISIB (CHEMBL3360203)	inhibitor
EGF	Sucralfate	inducer	PIK3R1	VOXTALISIB	inhibitor
PIK3R1	PF-04691502	inhibitor	PIK3R1	ZSTK-474	inhibitor
PIK3R1	Puquitinib	inhibitor	PIK3R1	ALPELISIB	inhibitor
PIK3R1	PA-799	inhibitor	PIK3R1	PI-103	inhibitor
PIK3R1	ISOPRENALINE	agonist	PIK3R1	QUERCETIN	inhibitor
PIK3R1	GSK-2636771	inhibitor	PIK3R1	PILARALISIB (CHEMBL3218575)	inhibitor
PIK3R1	DS-7423	inhibitor	PIK3R1	WX-037	inhibitor
PIK3R1	OMIPALISIB	inhibitor	PIK3R1	BGT-226 (CHEMBL3545096)	inhibitor
PIK3R1	RECILISIB	inhibitor	PIK3R1	PICTILISIB	inhibitor
PIK3R1	PWT-33587	inhibitor	PIK3R1	BUPARLISIB	inhibitor
PIK3R1	RG-7666	inhibitor	PIK3R1	DACTOLISIB	inhibitor
PIK3R1	SF-1126	inhibitor	PIK3R1	Panulisib	inhibitor
PIK3R1	VS-5584	inhibitor	PIK3R1	GSK-1059615	inhibitor
PIK3R1	COPANLISIB	inhibitor	PIK3R1	AZD-6482	inhibitor
PIK3R1	GEDATOLISIB	inhibitor	PIK3R1	BUPARLISIB HYDROCHLORIDE	inhibitor
PIK3R1	SONOLISIB	inhibitor	PIK3R1	LY-3023414	inhibitor
PIK3R1	APITOLISIB	inhibitor

Discussion

SCI is a serious neurological accident with high mortality rate[28,29]. SCI involves a variety of complex pathophysiological processes, and its molecular mechanism is not completely clear. In this paper, a total of 412 DEmRNAs and 21 DElncRNAs were obtained between two days and five days post-SCI groups. The functional enrichment analyses shown that the DEmRNAs were mainly involved in hsa04010: MAPK signaling pathway, hsa04068: FoxO signaling pathway, and muscle tissue/ organ development related GO-BP terms. In addition, UBE2D1, JUN, and FBXO32 had higher node degrees in PPI network, and the top 20 genes with high degree were significantly enriched in hsa04068: FoxO signaling pathway and hsa05213: Endometrial cancer. Moreover, FOXO3 was regulated by hsa-miR-1207-5p and hsa-miR-1207-5p was regulated by lncRNA RP11-253E3.3 in lncRNA-miRNA-mRNA ceRNA network. Finally, 37 drug-gene interactions were obtained based on the 26 genes in ceRNA network.

In this study, UBE2D1, JUN, and FBXO32 had higher node degrees in PPI network. JUN, as known as c-Jun, has been reported association with SCI. For instance, Zhang et al. found that miR-152 overexpression significantly reduced the levels of inflammation genes as well as the expression of c-Jun in SCI[25]. Vinit et al. uncovered that after cervical SCI, the HSP27 and c-Jun were distinct expressed in axotomized and spared bulbospinal neurons[26]. Besides, the enrichment analysis shown that JUN was involved in the MAPK signaling pathway. Numerous studies revealed that MAPK signaling pathway plays a significant role in the occurrence of SCI. For example, Zhan et al. have illustrated that fasudil promotes the migration of bone marrow mesenchymal stem cells by activating MAPK signaling pathway and its application in SCI model[27]. In addition, p38 MAPK signaling pathway is one of the classical inflammatory pathways, which is associated with the initiation and progression of inflammation[28], and inflammation promotes skeletal muscle atrophy[29]. Besides, Fan et al. to explore the effects and mechanism of electroacupuncture (EA) on expression of FBXO32 in traumatic spinal cord injury (TSCI) rats, and the results shown that the expression of FBXO32 mRNA were higher in model group when compares with sham group[30]. Moreover, the enrichment analysis shown that the DEmRNAs were mainly involved in the muscle tissue/ organ development related GO-BP terms. However, nowadays little research has made about UBE2D1gene in SCI. Taken together, we speculate that UBE2D1, JUN, and FBXO32 contribute to skeletal muscle atrophy after SCI progression, and activating of MAPK signaling pathway may induce skeletal muscle inflammation after SCI.

Moreover, FOXO3 was regulated by hsa-miR-1207-5p and hsa-miR-1207-5p was regulated by lncRNA RP11-253E3.3 in lncRNA-miRNA-mRNA ceRNA network. Lundell et al. have shown that FOXO3 protein is decreased in response to SCI[6]. Zhang et al. found that the p27 (kip1) and FOXO3a decreased levels in spinal cord were involved in axonal regeneration and glial cell proliferation after SCI[31]. In addition, Papagregoriou et al. have uncovered that in CFHR5 nephropathy, miR-1207-5p binding site polymorphism disrupts the regulation of HBEGF and is associated with disease severity[32]. Moreover, the enrichment analysis shown that FOXO3 was mainly enriched in FoxO signaling pathway and Endometrial cancer. Previous studies suggested that FoxO signaling pathway is related to muscle cell apoptosis, skeletal muscle autophagy and atrophy, and activating of FoxO signaling pathway will causes the atrophy of myotubes and mature skeletal muscle[33,34]. However, few studies reported about miR-1207-5p, lncRNA RP11-253E3.3, FoxO signaling pathway and Endometrial cancer in SCI. Therefore, we speculated that RP11-253E3.3-hsa-miR-1207-5p-FOXO3 axis was likely to be related to the progression of skeletal muscle atrophy after SCI via activating Endometrial cancer and FoxO signaling pathways.

Although we explored the potential molecular mechanisms of skeletal muscle atrophy after SCI using a bioinformatics approach, there still exist some limitations in current study. For instance, relevant experiments including cell biology assays, and animal and clinical studies need to be performed to verify the multiple candidate targets and signaling pathways identified from our bioinformatics analyses.

In summary, UBE2D1, JUN, FBXO32 are likely to be related to the progression of skeletal muscle atrophy after SCI, and activating of MAPK signaling pathway, Endometrial cancer and FoxO signaling pathway may induce skeletal muscle inflammation, apoptosis, autophagy and atrophy after SCI. Moreover, RP11-253E3.3-hsa-miR-1207-5p-FOXO3 axis may be used as a therapeutic target for skeletal muscle atrophy after SCI treatment. This study will provide a new ideas for further studies of skeletal muscle atrophy after SCI treatment.

Funding

This work was supported by Natural Science Research in Minhang District of Shanghai (Program No. 2019MHZ061).

Authors’ contributions

Conception and design of the research: WZ and LW; acquisition of data: GC, LJ, LG and ZY; analysis and interpretation of data: GC, LJ, LG and ZY; statistical analysis: GC, LJ, LG and ZY; obtaining funding: WZ; drafting the manuscript: LW; revision of manuscript for important intellectual content: WZ. All authors read and approved the final manuscript.