Identification of key genes and pathways associated with spinal cord injury.

The present study was designed to identify key genes or significant signaling pathways associated with spinal cord injury (SCI), and to clarify the underlying molecular mechanisms of SCI. Data from the GSE45550 array were downloaded from the Gene Expression Omnibus database. A total of 6 control samples, 6 samples at 3 days post‑SCI (SCI3d), 6 samples at 8 days post‑SCI (SCI8d) and 6 samples at 14 days post‑SCI (SCI14d) were included. The microarray data was preprocessed by the robust multi‑array average algorithm. The differentially expressed genes (DEGs) were identified using the limma package. The overlapping DEGs among groups were analyzed using the Venny 2.0 online tool. Modules enriched by DEGs were selected by weighted gene co‑expression network analysis. Gene Ontology annotation and the Kyoto Encyclopedia of Genes and Genomes pathways were identified for DEGs using the Database for Annotation, Visualization and Integrated Discovery. A total of 693 genes were obtained by combining the DEGs of the SCI3d, SCI8d and SCI14d groups. The pink module and green module with smaller P‑values obtained from weighted gene co‑expression network analysis module analyses of DEGs demonstrated a higher correlation with SCI. In addition, the peroxisome proliferator‑activated receptor (PPAR) signaling pathway that the cluster of differentiation 36 (CD36) was significantly enriched in, was one of the significant pathways in the pink module. The p53 signaling pathway that Caspase‑3 (Casp3) was significantly enriched in was one of the significant pathways in the green module. In conclusion, the results of the present study demonstrated that the PPAR and p53 signaling pathway may serve important roles in the progression of SCI. In addition, CD36 and Casp3 may be involved in the progression of SCI via the PPAR and p53 signaling pathways, respectively.

Introduction

Spinal cord injury (SCI) is an injury to the spinal cord caused by trauma or disease, which may lead to alterations to the normal motor, sensory or autonomic function of the spinal cord (1). SCI is associated with high morbidity and mortality rates, and the SCI annual incidence rate ranges between 12.1–57.8 cases/million individuals (1,2). In addition, the epidemiology of SCI is variable in different countries, and there is currently no effective treatment (3,4). Therefore, an effective therapy for the treatment of SCI is required, and it has been suggested that genes associated with SCI may be able to provide novel strategies for such a treatment.

Despite the lack of effective treatments, there have been some notable findings associated with the molecular mechanism of SCI. It has been observed that SCI results in secondary degeneration involving apoptosis, with an increased expression of genes associated with apoptosis and decreased expression of anti-apoptotic genes (5). In addition, a reduction in excessive M1 macrophage polarization and an enhancement of M2 macrophage polarization produced by regulating the levels of cytokines, including tumor necrosis factor a and interleukin (IL)-1β in the SCI microenvironment, may be a desirable treatment method (6). Tachibana et al (7) observed that 3 genes, including heat shock 27 kDa protein, tissue inhibitor of metalloproteinase-1 and epidermal fatty acid-binding protein, were upregulated in SCI. A recent study indicated that deletion of the IL-1α gene protected oligodendrocytes from SCI by overexpressing TOX high mobility group box family member 3 (8). Furthermore, numerous genes associated with inflammation, such as Arginase 1, are differentially expressed in the ephrin type A-receptor 4 knockout mouse model of SCI (9). In addition, it has been demonstrated that a temporal blockade of the IL-6 signaling pathway may modify the inflammatory response following SCI, and thus promote regeneration of the spinal cord (10). However, there are various additional important genes and pathways associated with SCI that have yet to be explored completely. Thus, a greater understanding of these genes and pathways is required as they may provide novel targets for SCI therapy.

In the present study, GSE45550 microarray data was obtained from the Gene Expression Omnibus (GEO) and used to identify the differentially expressed genes (DEGs) associated with SCI. Functional enrichment analyses were performed for DEGs. Furthermore, functions of gene modules were analyzed. The aim of the present study was to identify critical genes or significant signaling pathways associated with SCI, and clarify the underlying molecular mechanisms involved.

Materials and methods

Affymetrix microarray data

The microarray data from GSE45550 was downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/) (11). The following 4 groups were applied: 6 control samples, 6 samples at 3 days post-SCI (SCI3d), 6 samples at 8 days post-SCI (SCI8d) and 6 samples at 14 days post-SCI (SCI14d). Data from the GPL1355 platform [(Rat230_2) Affymetrix Rat Genome 230 2.0 Array; Affymetrix Inc., Santa Clara, CA, USA] were used for subsequent analysis.

Data preprocessing

The microarray data was preprocessed using the robust multi-array average algorithm with the Affy package (12) in Bioconductor (version 1.46.1; http://www.bioconductor.org/). Background correction, normalization and calculation of expression were all included in the process of preprocessing. The probe of the microarray data was transformed to gene symbols with Bioconductor AnnotationData software packages. If several probes were mapped to one gene symbol, then the mean value was set as the final expression value of this gene. A total of 18,634 gene expression matrixes were obtained from the above process.

DEGs analysis

The DEGs in the following three comparison groups: SCI3d vs. Control, SCI8d vs. Control and SCI14d vs. Control were analyzed using the limma package (13) in Bioconductor. The DEG P-values were calculated using the unpaired Student's t-test (14) provided by the limma package, and the P-values were adjusted to false discovery rate (FDR) values using the Benjamini-Hochberg correction (15). Log2 fold-change (FC) ≥1 and FDR values <0.05 were used as cut-off criterion for DEGs. Hierarchical clustering analysis of the DEGs was then performed and visualized using g-plots (16) in the R package.

Venn diagram analysis of DEGs

Venny is an interactive tool used to compare lists with Venn diagrams (17). The Kyoto Encyclopedia of Genes and Genomes (KEGG; www.genome.jp/kegg/) database is used to put associated gene sets into their respective pathway (18). The Database for Annotation, Visualization and Integrated Discovery (DAVID; http://david.ncifcrf.gov), used for analyzing gene lists, is an integrated data-mining environment (19).

The intersections of upregulated and downregulated genes in different sample groups were respectively analyzed using Venny 2.0 (17) (http://bioinfogp.cnb.csic.es/tools/venny/index.html) online tool. KEGG pathway enrichment analysis was performed for the intersection of genes by DAVID. P≤0.05 and gene counts ≥2 were used as threshold values.

Analysis of the correlation between gene modules and phenotype

Weighted gene co-expression network analysis (WGCNA) (20) is a tool used to identify gene clusters or modules which are highly associated with the phenotype of samples in expression profile data, and generalize module characteristic genes among these gene clusters. Furthermore, WGCNA provides correlation coefficients and significant thresholds in every module.

In the present study DEGs in the SCI3d, SCI8d and SCI14d groups were combined, and the correlation between these DEGs and SCI3d, SCI8d and SCI14d were analyzed, and gene sets with higher correlation were dug. Modules enriched by DEGs were selected by WGCNA in the R package, and modules significantly associated with SCI were identified with cluster analysis. The higher the absolute value of the correlation coefficient, the closer the correlation was between gene expression levels in modules and SCI.

Enrichment analysis of module function

Gene Ontology (GO) is a tool used to generate gene annotations by collecting defined, structured and controlled vocabulary (21). GO annotation and KEGG pathway analyses were performed for DEGs using DAVID. P<0.05 and gene counts ≥2 were set as threshold values.

Results

Normalized analysis of sample data

The boxplots of sample data prior to and following normalization are depicted in Fig. 1. The median line of the boxplot was at the same level following normalization, indicating that all data were successfully normalized.

Figure 1.

Boxplots of sample data prior to and following normalization. The x-axis represents the sample names, and the y-axis represents gene expression levels following log2 transformation. SCI, spinal cord injury; SCI3d, samples taken at 3 days post-SCI; SCI8d, samples taken at 8 days post-SCI; SCI14d, samples taken at 14 days post-SCI.

DEG analysis

The DEG count of the three SCI groups compared with the control group are summarized in Table I with |log2FC| values ≥1 and FDR values <0.05. The heat maps of the DEGs are depicted in Fig. 2.

Table I.

Differentially expressed genes count relative to the control group.

Group	Upregulated genes count	Downregulated genes count	Total
SCI3d	232	322	554
SCI8d	121	142	263
SCI14d	64	71	135

SCI3d, samples taken at 3 days post-SCI; SCI8d, samples taken at 8 days post-SCI; SCI14d, samples taken at 14 days post-SCI.

Figure 2.

Heat maps of the DEGs in the (A) SCI3d group, (B) SCI8d group and (C) SCI14d group relative to the control group. The green represents reduced expression levels, red represents higher expression levels and black indicates no differential expression among genes. DEGs, differentially expressed genes; SCI, spinal cord injury; SCI3d, samples taken at 3 days post-SCI; SCI8d, samples taken at 8 days post-SCI; SCI14d, samples taken at 14 days post-SCI.

Analysis of overlapping DEGs among groups

Venn diagram analyses for the upregulated and downregulated genes in the SCI3d, SCI8d and SCI14d groups were performed, and overlapping DEGs among groups is depicted in Fig. 3. A total of 9 genes were upregulated at all 3 time points and 48 genes were downregulated at all 3 time points.

Figure 3.

Intersectioning (A) upregulated and (B) downregulated genes in the SCI3d, SCI8d and SCI14d groups relative to controls. SCI, spinal cord injury; SCI3d, samples taken at 3 days post-SCI; SCI8d, samples taken at 8 days post-SCI; SCI14d, samples taken at 14 days post-SCI.

Where P<0.05, there was no significant enrichment observed among the 9 intersecting upregulated genes using KEGG pathway analysis. However, the 48 intersecting downregulated genes were markedly enriched in the peroxisome proliferator-activated receptor (PPAR) signaling pathway (P=8.01×10−4; enriched genes including lipoprotein lipase, fatty acid binding protein 3, aquaporin 7 and adiponectin, C1Q and collagen domain containing) and in the synthesis and degradation of ketone bodies signaling pathway (P=2.39×10−2; enriched genes including 3-oxoacid CoA transferase 1, 3-hydroxybutyrate dehydrogenase, type 1) (data not shown).

WGCNA module analyses of DEGs

A total of 693 genes were obtained by combining the DEGs of the SCI3d, SCI8d and SCI14d groups. Cluster analyses using WGCNA were performed using the expression data of these genes, and the cluster dendrogram is presented in Fig. 4. The DEGs were divided into 7 modules (pink module, green/yellow module, black module, blue module, green module, red module and tan module). The correlation coefficient and P-value between gene counts of every module and SCI are summarized in Table II. The data indicated that the pink module and green module with smaller P-values demonstrated a higher correlation with SCI. The functional enrichment results of pink and green modules are summarized in Table III. It was demonstrated that the PPAR signaling pathway that cluster of differentiation 36 (CD36) was significantly enriched in, was one of the significant pathways in the pink module. In addition, the p53 signaling pathway, that caspase-3 (Casp3) was significantly enriched in, was one of the significant pathways in the green module (Table III).

Figure 4.

Cluster dendrogram. Seven colors including pink, green/yellow, black, blue, green, red and tan represent seven modules obtained from WGCNA module analyses of DEGs. One color represents one module. The degree of co-expression for the genes in the same modules was high. WGCNA, weighted gene co-expression network analysis.

Table II.

Results of module analysis.

Factor	Pink module	Green/yellow module	Black module	Blue module	Green module	Red module	Tan module
Gene count	37	53	125	273	69	46	90
Correlation	−7.30E-01	−1.00E-01	−5.80E-01	−9.00E-02	7.50E-01	−7.80E-02	2.80E-01
P-value	5.29E-05	6.39E-01	3.08E-03	6.77E-01	2.28E-05	7.18E-01	1.88E-01

Seven modules (pink, green/yellow, black, blue, green, red and tan) represent modules obtained from WGCNA module analyses of DEGs. One color represents one module. The degree of co-expression for the genes in the same modules was high. WGCNA, Weighted gene co-expression network analysis.

Table III.

Functional enrichment results of pink and green modules.

A, Pink module
GO/KEGG	Term	Description	P-value	Gene
GO-BP	GO:0006631	Fatty acid metabolic process	8.00E-04	Cd36, Mlycd, Acot2, Fabp4, Decr1
	GO:0006936	Muscle contraction	9.00E-04	Tnnt3, Tnnc2, Mylpf, Actn3
	GO:0003012	Muscle system process	1.30E-03	Tnnt3, Tnnc2, Mylpf, Actn3
	GO:0034440	Lipid oxidation	3.70E-03	Cd36, Mlycd, Decr1
	GO:0019395	Fatty acid oxidation	3.70E-03	Cd36, Mlycd, Decr1
GO-CC	GO:0015629	Actin cytoskeleton	9.80E-03	Tnnt3, Tnnc2, Mylpf, Myoz1
	GO:0005861	Troponin complex	1.70E-02	Tnnt3, Tnnc2
	GO:0030016	Myofibril	1.83 E-02	Tnnt3, Tnnc2, Mylpf
	GO:0005865	Striated muscle thin filament	1.91E-02	Tnnt3, Tnnc2
	GO:0043292	Contractile fiber	2.12E-02	Tnnt3, Tnnc2, Mylpf
GO-MF	GO:0048037	Cofactor binding	1.70E-03	Ldhb, Nox1, Actn3, Decr1, Nqo1
	GO:0008092	Cytoskeletal protein binding	4.30E-03	Tnnt3, Tnnc2, Mylpf, Myoz1, Actn3
	GO:0050662	Coenzyme binding	6.80E-03	Ldhb, Nox1, Decr1, Nqo1
	GO:0003779	Actin binding	7.50E-03	Tnnt3, Tnnc2, Mylpf, Actn3
KEGG	rno04510	Focal adhesion	1.67E-02	Col6A3, Mylpf, Actn3, Col5A3
	rno03320	PPAR signaling pathway	1.70E-02	Cd36, Fabp4, Aqp7
	rno04512	ECM-receptor interaction	2.18 E-02	Cd36, Col6A3, Col5A3
	rno04670	Leukocyte transendothelial migration	4.17E-02	Nox1, Mylpf, Actn3
B, Green module
GO/KEGG	Term	Description	P-value	Gene
GO-BP	GO:0030199	Collagen fibril organization	1.00E-04	Col3A1, Lox, Col5A2, Col5A1
	GO:0030198	Extracellular matrix organization	4.00E-04	Col3A1, Ccdc80, Lox, Col5A2, Col5A1
	GO:0042060	Wound healing	6.00E-04	Ccnb1, Casp3, Col3A1, Lox, Tfpi2, Col5A1
	GO:0030155	Regulation of cell adhesion	1.30E-03	Tnc, Ccdc80, Jak2, Mmp14, Col8A1
	GO:0007049	Cell cycle	2.10E-03	Ccnb1, Ccnb2, Mki67, Cks2, Pttg1, Cdkn3, Ube2C, Racgap1
GO-CC	GO:0031012	Extracellular matrix	6.06E-14	Matn2, Aspn, Cthrc1, Tnc, Col3A1, Ccdc80, Mmp14, Col5A2, Col5A1, Col6A2, Col6A1, Lox, Loxl1, Thbs4, Adamts4, Spon1
	GO:0005578	Proteinaceous extracellular matrix	4.22E-12	Matn2, Aspn, Cthrc1, Tnc, Col3A1, Ccdc80, Lox, Mmp14, Col5A2, Loxl1, Col5A1, Adamts4, Spon1, Thbs4
	GO:0044421	Extracellular region part	1.36E-09	Matn2, Aspn, Cthrc1, Tnc, Col3A1, Ccdc80, Mmp14, Col5A2, Col5A1, Ctsk, Cpxm1, Col6A2, Col6A1, Lox, Loxl1, Thbs4, Adamts4, Spon1
	GO:0005576	Extracellular region	1.18E-08	Matn2, Aspn, Cthrc1, Aebp1, Tnc, Col3A1, Ccdc80, Mmp14, Col5A2, Col5A1, Ctsk, Penk, Cpxm1, Sfrp4, Col6A2, Col6A1, Lox, Loxl1, Tfpi2, Thbs4, Adamts4, Spon1
	GO:0005581	Collagen	9.18E-05	Col3A1, Lox, Col5A2, Col5A1
GO-MF	GO:0016641	Oxidoreductase activity, acting on the CH-NH2 group of donors, oxygen as acceptor	1.50E-03	Lox, Loxl2, Loxl1
	GO:0016638	Oxidoreductase activity, acting on the CH-NH2 group of donors	2.20E-03	Lox, Loxl2, Loxl1
	GO:0005201	Extracellular matrix structural constituent	5.20E-03	Col3A1, Col5A2, Col5A1
	GO:0016702	Oxidoreductase activity, acting on single donors with incorporation of molecular oxygen, incorporation of two atoms of oxygen	6.40E-03	Plod2, P4Ha3, Hpd
	GO:0016701	Oxidoreductase activity, acting on single donors with incorporation of molecular oxygen	6.80E-03	Plod2, P4Ha3, Hpd
KEGG	rno04512	ECM-receptor interaction	6.37E-07	Tnc, Col3A1, Col6A2, Col6A1, Col5A2, Col5A1, Thbs4
	rno04510	Focal adhesion	1.03E-04	Tnc, Col3A1, Col6A2, Col6A1, Col5A2, Col5A1, Thbs4
	rno04115	p53 signaling pathway	2.96E-02	Ccnb1, Casp3, Ccnb2

GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; BP, biological process; CC, cellular component; MF, molecular function.

Discussion

In the present study, a total of 693 genes were identified by combining the DEGs of the SCI3d, SCI8d and SCI14d groups. The results obtained demonstrated that the PPAR signaling pathway, in which CD36 was significantly enriched, was one of the significant pathways in the pink module, while the p53 signaling pathway that Casp3 was significantly enriched in, was one of the significant pathways in the green module.

PPAR, which includes 3 isoforms (PPAR-α PPAR-γ PPAR-β/δ), is involved in the inflammation process (22). A previous study demonstrated that PPAR-γ and PPAR-δ are involved in the protective effects of palmitoylethanolamide (PEA) in SCI, indicating that PPAR-γ and PPAR-δ may contribute to the anti-inflammatory effects of PEA in SCI (23). It has been reported that PPAR participates in the pathogenesis of diseases, such as diabetes and SCI (24,25). In addition, it has been observed that the decrease of PPAR-δ expression in the spinal cord of streptozotocin (STZ)-diabetic rats may explain the higher mortality rate observed following SCI in STZ-diabetic rats (26). Thus, activation of PPAR-δ may reduce the severity of SCI, and PPAR-δ may be a target for therapy in SCI patients (27). In addition, van Neerven and Mey (28) indicated that endogenous PPAR ligands may contribute by preventing the expansion of the initial damage via modulating inflammation post-SCI, and concluded that the PPAR signaling pathway may be a target for treatment of SCI. Therefore, the results of the present study are in agreement with the results of previous studies and, as such, the PPAR signaling pathway may be closely associated with the pathogenesis of SCI. In addition, the results demonstrated that CD36 was significantly enriched in the PPAR signaling pathway. It has been indicated that the upregulation of CD36, an integral membrane protein, may resolve inflammation via the 5-lipoxygenase-dependent and PPAR-γ signaling pathways (29). The wnt-1 proto-oncogene protein promotes CD36 expression via β-catenin and PPAR-γ signaling pathways (30). CD36 is involved in the PPAR-γ signaling pathway, and activation of PPAR-γ leads to upregulation of CD36 in the PPAR signaling pathway (31). Previous studies have indicated that CD36 is involved in the PPAR-γ signaling pathway (22,32). Although the roles of CD36 in SCI have not been extensively studied, when considering the association between the PPAR signaling pathway and SCI, it is possible that CD36 may be involved in the progression of SCI via the PPAR signaling pathway.

The p53 signaling pathway was observed to be one of the significant KEGG pathways in the green module in the present study. A previous study observed DNA damage-induced p53-mediated mitochondrial apoptosis in SCI (5). Nerve injury is a significant consequence of SCI, and p53 is involved in glial cell death and survival in SCI (33). p53, a key molecular regulator of apoptotic cell death, is known to promote apoptosis by increasing the transcription of target genes (34,35). In addition, a number of apoptosis-associated molecules are expressed via p53, and apoptosis induction via the p53 pathway is a complex process (33). Furthermore, minocycline has been shown to reduce apoptosis in models of SCI (36). Therefore, the results of the present study are consistent with previous studies, and therefore indicate that the p53 signaling pathway may be significant in the progression of SCI. In the present study, it was demonstrated that Casp3 was significantly enriched in the p53 signaling pathway. It has been observed that Casp3 and the p53 signaling pathway may regulate nitric oxide-induced extracellular signal-regulated protein kinase and p38 kinase (37). In addition, p53 prevents the α6β4 integrin-mediated activation of serine/threonine-specific protein kinase B by promoting Casp3 dependent cleavage (38). It has been suggested that Casp3 is a critical mediator of p53-induced apoptosis (39). Zhang et al (40) investigated the apoptosis process in SCI, and demonstrated that the Casp3 apoptotic pathway components are activated following SCI in rats. Ultimately, it is thought that Casp3 may be involved in the progression of SCI via the p53 signaling pathway.

In conclusion, the PPAR and p53 signaling pathways may be important pathways associated with SCI. In addition, CD36 and Casp3 may be involved in the progression of SCI via the PPAR and p53 signaling pathways, respectively. However, the results of the present study are limited by the small sample size, and therefore further studies are required to evaluate the molecular mechanisms underlying SCI progression.