Elucidation of Gene Expression Patterns in the Brain after Spinal Cord Injury.

Spinal cord injury (SCI) is a devastating neurological disease. The pathophysiological mechanisms of SCI have been reported to be relevant to central nervous system injury such as brain injury. In this study, gene expression of the brain after SCI was elucidated using transcriptome analysis to characterize the temporal changes in global gene expression patterns in a SCI mouse model. Subjects were randomly classified into 3 groups: sham control, acute (3 h post-injury), and subacute (2 wk post-injury) groups. We sought to confirm the genes differentially expressed between post-injured groups and sham control group. Therefore, we performed transcriptome analysis to investigate the enriched pathways associated with pathophysiology of the brain after SCI using Database for Annotation Visualization, and Integrated Discovery (DAVID), which yielded Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway. Following enriched pathways were found in the brain: oxidative phosphorylation pathway; inflammatory response pathways-cytokine-cytokine receptor interaction and chemokine signaling pathway; and endoplasmic reticulum (ER) stress-related pathways-antigen processing and presentation and mitogen-activated protein kinase signaling pathway. Oxidative phosphorylation pathway was identified at acute phase, while inflammation response and ER stress-related pathways were identified at subacute phase. Since the following pathways-oxidative phosphorylation pathway, inflammatory response pathways, and ER stress-related pathways-have been well known in the SCI, we suggested a link between SCI and brain injury. These mechanisms provide valuable reference data for better understanding pathophysiological processes in the brain after SCI.

Introduction

Spinal cord injury (SCI) is one of the most debilitating neurological diseases affecting the motor, sensory, and the autonomic systems.[1,2] Because the pathophysiology mechanism of SCI is not fully reported, a comprehensive approach is important to study SCI pathogenesis.

After SCI, pathophysiological events may occur, including acute, subacute, and chronic phases. The acute phase happens immediately after injury and primary damage leads to immediate physical and biochemical cellular alteration such as hemorrhage, ischemia, and hypoxia.[3,4] Subsequently, subacute phase, occurring over time after SCI, leads to further damage SCI results in a rapid and extensive oxidative stress reaction, which causes oxidative death of the spinal cord neurons and reduces spinal cord blood flow that leads to edema and inflammatory responses at subacute phase.[5,6] Moreover, SCI results in apoptosis, which severely affect neurons, oligodendrocyte, microglia, and, perhaps, astrocytes, as well as astrogliosis.[6,7] Furthermore, the chronic phase of SCI consists of many incidents, such as white matter demyelination, gray matter dissolution, connective tissue deposition, and reactive gliosis. These events lead to glial scar formation.[6]

Additionally, pathophysiological mechanisms of SCI have been further related to brain injury such as cerebral ischemia and subarachnoid hemorrhage.[8-10] Recent studies have shown that SCI causes brain inflammation, progression of nerve cell loss, as well as loss of brain functions such as cognition and could lead to depression.[11,12] Several studies already have investigated the alteration of gene expression response in the brain after SCI.[13-19] Brain-derived neurotrophic factor (BDNF), glial cell line–derived neurotrophic factor, and histone deacetylase 1 (HDAC1) may play an important role in the brain reorganization after SCI;[20,21] however, only little is known and further studies will be needed.

In this study, we investigated the progression of the brain injury after SCI at acute phase (3 h post-injury) and subacute phase (2 wk post-injury) by gene expression patterns using transcriptome analysis. We performed and confirmed differentially expressed genes (DEGs) profiling to identify enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways that will enable the demonstration of the brain injury after SCI pathophysiological mechanism at acute phase and subacute phase. Some candidate genes belonged to enriched KEGG pathways and were considered to show significant roles in the brain after SCI. These results provide a better understanding of association between brain injury and SCI.

Materials and Methods

Animals

For all animal experiments, CD-1 (ICR) mice were housed in a facility accredited by the Association for Assessment and Accreditation of Laboratory Animal Care. The experimental procedure was approved by the Institutional Animal Care and Use Committee (IACUC) of Yonsei University Health System.

Male mice of age 7 wk (Orient bio, Gyeonggi-do, South Korea) were used and were randomly classified into 3 groups: sham control brain group (n = 3), 3 h post-SCI brain group (n = 3), and 2 wk post-SCI brain group (n = 3).

Spinal Cord Contusion

Animals were anesthetized with a mixture of ketamine (100 mg/kg, intraperitoneal, Yuhan, Seoul, Korea) and xylazine (10 mg/kg, intraperitoneal), and absence of blink and withdrawal reflexes were ensured. Body temperature was maintained at 37 °C in a hypoxic chamber. Mice received a dorsal laminectomy at the 9th thoracic vertebral (T9) level to expose the spinal cord and then a moderate T9 contusive injury by the Infinite Horizons device (Precision Systems and Instrumentation, Lexington, NY, USA) to moderate force of 70 kdyn. The sham control mice received only a dorsal laminectomy without contusive injury. Afterward, the wound was sutured in layers. Post-operative care consisted of manual bladder expression twice a day until automatic voiding returned spontaneously, which is generally around 7 days. After 3 h or 2 wk following SCI, the mice were anesthetized again with ketamine and xylazine and perfused transcardially with normal saline for isolation of injured brains and spinal cords. These tissues were frozen at −70 °C and processed for RNA isolation.

RNA Preparation

Total RNA was prepared from the bilateral cerebrum after SCI using the TRIzol reagent (Invitrogen Life Technologies, Carlsbad, CA, USA) according to the manufacturers’ protocols. For quality control, RNA quality and quantity were evaluated by 1% agarose gel electrophoresis and the optical density (OD) 260/280 ratio was confirmed with a Nanodrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA).

RNA Sequencing and Transcriptome Data Analysis

RNA sequencing was performed by Macrogen, Inc. (Seoul, Korea). The messenger RNA (mRNA) was transcribed into a library of templates. This successive cluster generation using reagents was provided by the Illumina® TruSeq™ (Illumina, San Diego, CA, USA) RNA Sample Preparation Kit.[22-25] We performed the transcriptome analysis by the following procedures: RNA-seq experiment and data handling procedure. The detailed procedures of RNA-seq experiment are performed following the manufacturer’s instruction. Firstly, there are 8 steps in TruSeq mRNA library construction: purify and fragment mRNA, synthesize first strand complementary DNA (cDNA), synthesize second strand cDNA, perform end repair, adenylate the 3′ end of the DNA fragments, ligate adapters, enrich DNA fragments, and enrich library validation. Purifying step for the poly-A-containing mRNA molecules was carried out using magnetic beads, which attached poly-T oligo. After that, the cleaved RNA fragments are copied into first strand cDNA using reverse transcriptase and random primers. Next, synthesize second strand cDNA process removes the RNA template and synthesizes a replaced strand to generate double-stranded (DS) cDNA with DNA polymerase I and Ribonuclease H (RNaseH). The end repair process converts the overhangs resulting from fragmentation into blunt ends, using an end repair mix (End Repair Mix [ERP]). Next, “A” nucleotide is added to the 3′ ends of the blunt fragments to prevent them from ligating to one another during the adapter ligation reaction. In the process of ligation of adapters, multiple indexing adapters ligate to the ends of the ds cDNA, preparing them for hybridization into a flow cell. In the enriched DNA fragments process, we performed polymerase chain reaction (PCR) to get enriched cDNA library. Finally, enriched library validation was performed for quality control analysis of the sample library and quantification of the DNA library templates.[26]

The second procedure of RNA sequencing is clustering and sequencing using the Illumina. The Illumina utilizes a unique “bridged” amplification reaction that occurs on the surface of the flow cell. A flow cell with millions of unique clusters is loaded into the HiSeq 2000 for automated cycles of extension and imaging. Solexa’s sequencing-by-synthesis utilizes 4 proprietary nucleotides possessing reversible fluorophore and termination properties.[22] Each sequencing cycle occurs in the presence of all 4 nucleotides, leading to higher accuracy than methods where only 1 nucleotide is present in the reaction mix at a time. The next procedure is data handling, which contains sequence quality check and data analysis. SolexaQA is a Perl-based software package that calculates quality statistics and creates visual representations of data quality from FASTQ files generated by Illumina second-generation sequencing technology. And then, Cufflinks assembles transcripts, estimates their abundances, and tests for differential expression and regulation in RNA-seq samples. It accepts reading of aligned RNA-seq and assembles the alignments into a parsimonious set of transcripts based on how many reads support each one, taking into account biases in the library preparation protocol. Transcripts with a fold change ≥|1.5| and P < 0.05 were considered significant and were included in downstream analysis.

KEGG Pathway Analyses

The lists of significant DEGs (fold change ≥ |1.5|) were submitted to the Database for Annotation, Visualization, and Integrated Discovery (DAVID v6.7; http://david.abcc.ncifcrf.gov/) for KEGG pathways analysis.[27-29]

Quantitative Real-Time Reverse Transcription Polymerase Chain Reaction (qRT-PCR)

DEGs of interest were selected for the validation of transcriptome analysis results by qRT-PCR. Total RNA was reverse transcribed into cDNA using ReverTra Ace® qPCR RT Master Mix with genomic DNA Remover (Toyobo, Osaka, Japan) according to the manufacturer’s instructions. The mRNA expression for genes of interest was profiled using qPCRBIO SyGreen Mix Hi-ROX (PCR Biosystems, London, UK) in a StepOnePlus Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). Data analysis was performed using the 2−ΔΔCT method. Primers used for qRT-PCR are described in Table 1.

Table 1.

Primers Used in qRT-PCR.

Gene	Forward Primer (5′ → 3′)	Reverse Primer (5′ → 3′)	Gene	Forward Primer (5′ → 3′)	Reverse Primer (5′ → 3′)
ATP6V0C	GGG ATC ATC GCC ATC TAC GC	CAA CAC CAG CAT CTC CGA CG	H2-L	CGG TAC ATC TCT GTC GGC TAT	CTT GGC GAT CTG CGT GAT
COX7B	AAA CGC ACT AAG CGC TCT CC	ACT CTG CCA ACA GGG GAC AT	H2-T22	AAA AGG GAA GGA GAG GCT GC	CTC CTC CCC ATT CAA CTG CC
TNFRSF25	GAG CTG CCT AGG AGG GAG AG	AAC GGG CCA TAC CTC TTC TG	LTα	GAT CGT GCC TTT CTC CGA CA	AAA GAG CTG GAC CTC GTG TG
CCR5	ACA CCC TGT TTC GCT GTA GG	TGC ATG GCC TGG TCT AGT CT	H2-K1	AAG CCC CTC ACC CTG AGA TG	GAG GCT GAG AGA AAA CGC ACA
CCR6	CCC TCA CAA CAT GGT CCT CC	AAA CGC ATA CAA CAC GGG G	HSPB1	TCA CCC GGA AAT ACA CGC TC	GGC CTC GAA AGT AAC CGG AA
HGF	CAC CTC CTC CTG CTT CAT GT	GCC CCT GTT CCT GAT ACA CC	JMJD7	AAG ATG CCT GAT GCC GTG AA	CTG GTA GGT TGC TGG TGT GT
IL12RB1	TGG GAG TCA GAG TGG CTC GT	TGG GAG TCA GAG TGG CTC GT	FGFR2	GAT CAC GGC TTC CCC AGA TT	CTC GGC CGA AAC TGT TAC CT
CCL17	ACT TCA AAG GGG CCA TTC CT	CAT GGC CTT GGG TTT TTC ACC	RASGRP3	AAA ATC CCC ATC CTT GGC GT	AGT GGT GAG AGG CAT TCT GC
GNG11	CAC ATC GAG GAT CTG CCG GA	CAG CTG CCC TTT TCT TTG AAG G	DDIT3	GTC ACA CGC ACA TCC CAA AG	CAC TTT CCG CTC GTT CTC CT
NFKBIA	ATC CTG ACC TGG TTT CGC TC	CTC ATC CTC GCT CTC GGG TA	DUSP5	TCG CCT ACA GAC CAG CCT AT	CGG GGA TCC ACT TGT AGT GT
HSPA1AB	TGA ACT ACA AGG GCG AGA GC	CCG CTG AGA GTC GTT GAA GT	MAP3K6	CCC TTC GTG AGG ATG TTT TCC	CAG CCT GTA CTA GCC CAT CG
HSPA5	ATT GGA GGT GGG CAA ACC AA	TCG CTG GGC ATC ATT GAA GT	PLA2G3	GGA TCT CCT GGG TAC CAC CT	ATC CCT GAA ATG GAG TCG GC
CALR	CGG GGA CCT GGA GAA GGA TA	CCA AAC CAC TCG GAA ACA GC	CX3CR1	CAC CAT TAG TCT GGG CGT CT	GCG GAA GTA GCA AAA GCT CA
PDIA3	CGG GGA CCT GGA GAA GGA TA	GGT AGC CAC TGA CCC CAT AC	GAPDH	CAT CAC TGC CAC CCA GAA GAC TG	ATG CCA GTG AGC TTC CCG TTC AG

Abbreviations: qRT-PCR, Quantitative Real-Time Reverse Transcription Polymerase Chain Reaction; ATP6V0C, ATPase, H+ transporting V0 subunit C; COX7B, Cytochrome c oxidase subunit VIIb; TNFRSF25, Tumor necrosis factor receptor superfamily member 25 isoform 1 precursor; CCR5, Chemokine (C-C motif) receptor 5; CCR6, Chemokine (C-C motif) receptor 6; HGF, Hepatocyte growth factor isoform 3 preproprotein; IL12RB1, Interleukin 12 receptor, beta 1; CCL17, Chemokine (C-C motif) ligand 17; GNG11, Guanine nucleotide binding protein (G protein), gamma 11; NFKBIA, Nuclear factor of kappa light polypeptide gene enhancer in B cells inhibitor, alpha; HSPA1AB, Heat shock protein 1AB; HSPA5, Heat shock protein 5; CALR, Calreticulin; PDIA3, Protein disulfide isomerase associated 3; H2-L, Histocompatibility 2, D region locus L; H2-T22, Histocompatibility 2, T region locus 22; LTα, Lymphotoxin A; H2-K1, Histocompatibility 2, K1, K region; HSPB1, Heat shock protein 1; JMJD7, Jumonji domain containing 7; FGFR2, Fibroblast growth factor receptor 2; RASGRP3, RAS, guanyl releasing protein 3; DDIT3, DNA damage-inducible transcript 3 protein; DUSP5, Dual specificity phosphatase 5; MAP3K6, Mitogen-activated protein kinase kinase kinase 6; PLA2G3, Phospholipase A2, group III; CX3CR1, Chemokine (C-X3-C motif) receptor 1; GAPDH, glyceraldehyde-3-phosphate desidrogenase.

Western Blot Analysis

To assess Cytochrome C Oxidase Subunit 7B (COX7B), TNF Receptor Superfamily Member 25 (TNFRSF25), Heat Shock 70 KDa Protein 1A/1B (HSPA1AB), Heat Shock 70kD Protein 5 (HSPA5), and NF-Kappa-B Inhibitor Alpha (NFKBIA), the brain samples were homogenized and dissolved in radioimmunoprecipitation assay buffer (RIPA) buffer (Thermo Fisher Scientific, Rockford, IL, USA) with protease inhibitors (Abcam, Cambridge, UK). Total proteins were quantified using bicinchoninic acid assay (BCA assay)™ Protein Assay Kit (Thermo Fisher Scientific). The samples were denatured and separated by 4% to 12% Bis-Tris gels in 1× NuPage MES SDS running buffer (Invitrogen, Eugene, OR, USA). Bands were transferred onto a polyvinylidene difluoride (PVDF) (Invitrogen) in 20% (vol/vol) methanol in NuPage Transfer Buffer (Invitrogen) at 4 °C. Membranes were blocked and then incubated overnight at 4 °C with the following antibodies: COX7B (1:1,000, Abcam), TNFRSF25 (1:1,000, Santa Cruz Biotechnology, Santa Cruz, CA, USA), HSPA1AB (1:1,000, Santa Cruz Biotechnology), HSPA5 (1:1,000, Santa Cruz Biotechnology), NFKBIA (1:1,000, Santa Cruz Biotechnology), and Actin (1:3,000, Santa Cruz Biotechnology). The next day, the blots were washed 3 times with tris-buffered saline (TBS) (Biosesang, Gyeonggi-do, Korea) plus 0.1% Tween 20 (Bio-Rad, Hercules, CA, USA) and incubated for 1 h with horseradish peroxidase (HRP)–conjugated secondary antibodies (1:4,000; Santa Cruz) at room temperature. After being washed 3 times with TBS plus 0.1% Tween 20, the protein was visualized with an enhanced chemiluminescence (ECL) detection system (Amersham Pharmacia Biotech, Little Chalfont, UK).

Statistical Analysis

All results were expressed as means ± standard error of the mean (SEM). Statistical analyses were conducted using the premier vendor for Statistical Package for Social Sciences (SPSS) Version 23.0 (SPSS Inc, Chicago, IL, USA). Nonparametric statistical analysis such as Mann–Whitney U test was used for the comparison of 2 groups. Student t test was used to confirm statistical results. Results with P value < 0.05 was considered statistically significant.

Results

Gene Expression Profile

Total RNA was prepared from samples obtained from the bilateral cerebrum after SCI. Next we performed transcriptome analysis by RNA sequencing in order to identify gene differentially expressed at different phases. Heat map profile of mRNA expression displaying differentially regulated transcripts is shown in Fig. 1A. Color key and z-score show relative expression level of the samples. Red is maximum and green represents minimum expression levels of the samples.

Figure 1.

Gene expression profile by transcriptome analysis. (A) Heat map of the 2-way hierarchical clustering using z-score for normalized value. The color scheme is based on z-scores, with upregulation in red and downregulation in green. (B) Spinal cord injury (SCI) phases were classified into 2 phases: acute phase and subacute phase. Differentially expressed genes (DEGs) were identified at acute phase (3 h post-SCI) compared to sham control, at subacute phase (2 wk post-SCI) compared to sham control, or at subacute phase compared to acute phase. Bar graphs show the number of DEGs with fold change ≥ |1.5| at different phases (P < 0.05). (C) Oxidative phosphorylation is an enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway at acute phase. (D) Cytokine–cytokine receptor interaction, chemokine signaling pathway, and antigen processing and presentation are enriched KEGG pathways at subacute phase. (E) Antigen processing and presentation, mitogen-activated protein kinase (MAPK) signaling pathway, and cytokine–cytokine receptor interaction are enriched KEGG pathways at subacute phase compared to acute phase.

Upregulated and downregulated genes (1.5-fold) were counted and summarized in Fig. 1B. We identified 78 DEGs, 50 DEGs were upregulated and 28 DEGs were downregulated, in the 3 h post-SCI brain compared to sham control brain. Furthermore, 383 DEGs, 201 upregulated DEGs and 182 downregulated DEGs, were identified in 2 wk post-SCI brain expression compared to sham control brain expression. Moreover, 316 DEGs—132 upregulated DEGs and 184 downregulated DEGs—were identified in 2 wk post-SCI brain expression compared to 3-h post-SCI brain.

Enriched KEGG Pathway Analysis

In order to analyze pathways, total DEGs were classified based on information regarding gene function using KEGG pathway database from DAVID. These results were statistically significant (P < 0.05) and demonstrated in Table 2.

Table 2.

The Enriched KEGG Pathways at Different Phases.

Term	Count	%	P Value	Genes
Acute phase compared to sham control
mmu00190: Oxidative phosphorylation	3	4.688	0.049	COX7B2, COX7B, ATP6V0C
Subacute phase compared to sham control
mmu04060: Cytokine–cytokine receptor interaction	11	2.967	0.010	CCR6, IL12RB1, CCR5, GM13305, TNFRSF25, CX3CR1, CCL21A, HGF, GM1987, PRL, CCL17
mmu04062: Chemokine signaling pathway	9	2.374	0.019	CCR6, CCR5, CX3CR1, CCL21A, NFKBIA, GNG11, FOXO3, GM1987, CCL17
mmu04612: Antigen processing and presentation	6	1.484	0.047	H2-L, PDIA3, HSPA1A, HSPA1B, HSPA5, CALR
Subacute phase compared to acute phase
mmu04612: Antigen processing and presentation	7	2.083	0.010	H2-K1, H2-T22, HSPA1A, H2-T23, HSPA1B, HSPA5, LTα
mmu04010: MAPK signaling pathway	11	3.472	0.014	DUSP5, FGFR2, MAP3K6, RASGRP3, JMJD7, FGF11, HSPB1, HSPA1A, HSPA1B, PLA2G3, DDIT3
mmu04060: Cytokine–cytokine receptor interaction	9	3.125	0.024	CCL12, IL12RB1, CCR5, GM13305, TNFRSF25, CX3CR1, PRL, LTα, CCL17

Abbreviations: These pathways are statistically significant (P < 0.05). KEGG, Kyoto Encyclopedia of Genes and Genomes; MAPK, mitogen-activated protein kinase; COX7B2, Cytochrome c oxidase subunit VIIb2; ATP6V0C, ATPase, H+ transporting V0 subunit C; COX7B, Cytochrome c oxidase subunit VIIb; TNFRSF25, Tumor necrosis factor receptor superfamily member 25 isoform 1 precursor; CCR5, Chemokine (C-C motif) receptor 5; CCR6, Chemokine (C-C motif) receptor 6; HGF, Hepatocyte growth factor isoform 3 preproprotein; IL12RB1, Interleukin 12 receptor, beta 1; CCL17, Chemokine (C-C motif) ligand 17; GNG11, Guanine nucleotide binding protein (G protein), gamma 11; NFKBIA, Nuclear factor of kappa light polypeptide gene enhancer in B cells inhibitor, alpha; HSPA1AB, Heat shock protein 1AB; HSPA5, Heat shock protein 5; CALR, Calreticulin; PDIA3, Protein disulfide isomerase associated 3; H2-L, Histocompatibility 2, D region locus L; H2-T22, Histocompatibility 2, T region locus 22; LTα, Lymphotoxin A; H2-K1, Histocompatibility 2, K1, K region; HSPB1, Heat shock protein 1; JMJD7, Jumonji domain containing 7; FGFR2, Fibroblast growth factor receptor 2; RASGRP3, RAS, guanyl releasing protein 3; DDIT3, DNA damage-inducible transcript 3 protein; DUSP5, Dual specificity phosphatase 5; MAP3K6, Mitogen-activated protein kinase kinase kinase 6; PLA2G3, Phospholipase A2, group III; CX3CR1, Chemokine (C-X3-C motif) receptor 1; GM13305, Predicted gene 13305; CCL21A, Chemokine (C-C motif) ligand 21A (serine); GM1987, Predicted gene 1987; PRL, Prolactin; FOXO3, Forkhead box O3; H2-T23, Histocompatibility 2, T region locus 23; FGF11, Fibroblast growth factor 11 isoform 2; CCL12, Chemokine (C-C motif) ligand 12.

At acute phase, we observed enriched KEGG pathways including oxidative phosphorylation compared to sham control. At subacute phase, there was cytokine–cytokine receptor interaction, chemokine signaling pathway, and antigen processing and presentation compared to sham control. Compared to acute phase, enriched KEGG pathways such as antigen processing and presentation, mitogen-activated protein kinase (MAPK) signaling pathway, and cytokine–cytokine receptor interaction were altered in subacute phase. Each upregulated and downregulated genes was counted and summarized in Fig. 1C to E.

Validation of Transcriptome Data Using qRT-PCR and Western Blot

To validate the altered gene expression (fold change ≥ |1.5|), we performed qRT-PCR and Western blot at different phases. To support the reliability and accuracy of the RNA-seq expression results, we first performed qRT-PCR, and Western blot was conducted for further validation. Each validated DEGs involved in enriched KEGG pathways was summarized in Table 3.

Table 3.

Validated DEGs Involved in Enriched KEGG Pathways at Different Phases.

Pathway	Gene Symbol	Gene Description	Fold Change
Acute phase compared to sham control
Oxidative phosphorylation	ATP6V0C	ATPase, H+ transporting V0 subunit C	6.34
	COX7B	Cytochrome c oxidase subunit VIIb	−1.59
Subacute phase compared to sham control
Cytokine–cytokine receptor interaction	TNFRSF25	Tumor necrosis factor receptor superfamily member 25 isoform 1 precursor	2.7
	CCR5	Chemokine (C-C motif) receptor 5	1.7
	CCR6	Chemokine (C-C motif) receptor 6	1.6
	HGF	Hepatocyte growth factor isoform 3 preproprotein	−1.5
	IL12RB1	Interleukin 12 receptor, beta 1	−1.6
	CCL17	Chemokine (C-C motif) ligand 17	−2.0
Chemokine signaling pathway	GNG11	Guanine nucleotide binding protein (G protein), gamma 11	2.0
	CCR5	Chemokine (C-C motif) receptor 5	1.7
	CCR6	Chemokine (C-C motif) receptor 6	1.6
	CCL17	Chemokine (C-C motif) ligand 17	−2.0
	NFKBIA	Nuclear factor of kappa light polypeptide gene enhancer in B cells inhibitor, alpha	−2.1
Antigen processing and presentation	HSPA1AB	Heat shock protein 1AB	3.9
	HSPA5	Heat shock protein 5	2.4
	CALR	Calreticulin	1.6
	PDIA3	Protein disulfide isomerase associated 3	1.5
	H2-L	Histocompatibility 2, D region locus L	−1.8
Subacute phase compared to acute phase
Antigen processing and presentation	HSPA1AB	Heat shock protein 1AB	3.3
	HSPA5	Heat shock protein 5	2.1
	H2-T22	Histocompatibility 2, T region locus 22	−1.6
	LTα	Lymphotoxin A	−1.6
	H2-K1	Histocompatibility 2, K1, K region	−1.7
MAPK signaling pathway	HSPA1AB	Heat shock protein 1AB	3.3
	HSPB1	Heat shock protein 1	2.5
	JMJD7	Jumonji domain containing 7	1.8
	FGFR2	Fibroblast growth factor receptor 2	1.7
	RASGRP3	RAS, guanyl releasing protein 3	1.6
	DDIT3	DNA damage-inducible transcript 3 protein	1.6
	DUSP5	Dual specificity phosphatase 5	−1.7
	MAP3K6	Mitogen-activated protein kinase kinase kinase 6	−2.3
	PLA2G3	Phospholipase A2, group III	−3.6
Cytokine–cytokine receptor interaction	TNFRSF25	Tumor necrosis factor receptor superfamily member 25 isoform 1 precursor	3.2
	CX3CR1	Chemokine (C-X3-C motif) receptor 1	1.8
	CCR5	Chemokine (C-C motif) receptor 5	1.6
	IL12RB1	Interleukin 12 receptor, beta 1	−1.6
	LTα	Lymphotoxin A	−1.6
	CCL17	Chemokine (C-C motif) ligand 17	−2.1

Abbreviations: DEG, differentially expressed gene; KEGG, Kyoto Encyclopedia of Genes and Genomes; MAPK, mitogen-activated protein kinase.

In oxidative phosphorylation, ATP60C was increased (1.35-fold, P = 0.0133) whereas COX7B was decreased (−1.96-fold, P = 0.0136) significantly in a similar pattern as transcriptome data at acute phase compared to sham control (Fig. 2A). Western blot also exhibited same pattern as qRT-PCR: ATP6V0C (1.6-fold, P = 0.0013) and COX7B (0.5247-fold, P = 0.0045; Fig. 2B, C).

Figure 2.

Validation using quantitative real-time reverse transcription polymerase chain reaction (qRT-PCR) and Western blot at acute phase. (A) RNA sequencing, RNA-Seq; ATP6V0C is significantly upregulated and COX7B, Cytochrome c oxidase subunit VIIb is significantly downregulated in oxidative phosphorylation pathway at acute phase (n = 3 per group; *P < 0.05). The relative expression of target genes from qRT-PCR were calculated using 2−ΔΔCT method. All results were expressed as means ± standard error of the mean. (B) Western blot analysis was performed using antibodies against ATP6V0C, COX7B, and ACTIN (as control). (C) Comparison of relative protein expression from the brain between sham control and 3 h post–spinal cord injury verified by Western blot. **P < 0.01.

At subacute phase, compared to sham control, 3 upregulated genes (TNFRSF25, CCR5, and CCR6) and 3 downregulated genes (TNFRSF25, CCR5, and CCR6) were identified in cytokine–cytokine receptor interaction, 3 upregulated genes (GNG11, CCR5, and CCR6) and 2 downregulated genes (CCL17, NFKBIA) were identified in chemokine signaling pathway, and there were 4 upregulated genes (HSPA1AB, HSPA5, CALR, and PDIA3) and 1 downregulated gene (H2-L) in antigen processing and presentation. These genes were determined by qRT-PCR as follows: TNFRSF25 (5.37-fold, P = 0.0284), HSPA1AB (4.84-fold, P = 0.0284), CCR6 (3.15-fold, P = 0.0133), CCR5 (2.42-fold, P = 0.0133), HSPA5 (2.36-fold, P = 0.0284), GNG11 (2.01-fold, P = 0.0133), CALR (1.57-fold, P = 0.0284), PDIA3 (1.45-fold, P = 0.0136), HGF (−1.43-fold, P = 0.0133), H2-L (−1.60-fold, P = 0.0136), NFKBIA (−1.80-fold, P = 0.0136), CCL17 (−2.03-fold, P = 0.0136), and IL12RB1 (−3.10-fold, P = 0.0167; Fig. 3A). Western blot showed that TNFRSF25 (1.49-fold, P = 0.0044), HSPA1AB (1.67-fold, P = 0.0096), and HSPA5 (1.43-fold, P = 0.0359) proteins were increased, although NFKBIA was decreased (0.049-fold, P = 0.0059) significantly at subacute phase compared to sham control (Fig. 3B, C).

Figure 3.

Validation using quantitative real-time reverse transcription polymerase chain reaction (qRT-PCR) and Western blot at subacute phase. (A) TNFRSF25, HSPA1AB, CCR6, CCR5, HSPA5, GNG11, CALR, and PDIA3 are significantly upregulated and HGF, H2-L, NFKBIA, CCL17, and IL12RB1 are significantly downregulated at subacute phase. These genes are involved in cytokine–cytokine receptor interaction, chemokine signaling pathway, and antigen processing and presentation (n = 3 per group; *P < 0.05). The relative expression of target genes from qRT-PCR was calculated using 2−ΔΔCT method. All results were expressed as means ± standard error of the mean. (B) Western blot analysis was performed using antibodies against TNFRSF25, HSPA1AB, HSPA5, NFKBIA, and ACTIN (as control). (C) Comparison of relative protein expression from the brain between sham control and 2 wk post–spinal cord injury verified by Western blot. RNA sequencing, RNA-Seq; ATP6V0C, ATPase, H+ transporting V0 subunit C; COX7B, TNFRSF25, Tumor necrosis factor receptor superfamily member 25 isoform 1 precursor; CCR5, Chemokine (C-C motif) receptor 5; CCR6, Chemokine (C-C motif) receptor 6; IL12RB1, Interleukin 12 receptor, beta 1; CCL17, Chemokine (C-C motif) ligand 17; GNG11, Guanine nucleotide binding protein (G protein), gamma 11; NFKBIA, Nuclear factor of kappa light polypeptide gene enhancer in B cells inhibitor, alpha; HSPA1AB, Heat shock protein 1AB; HSPA5, Heat shock protein 5; CALR, Calreticulin; PDIA3, Protein disulfide isomerase associated 3; H2-L, Histocompatibility 2, D region locus L. *P < 0.05. **P < 0.01.

At subacute phase, compared to acute phase, there were 3 enriched KEGG pathways such as antigen processing and presentation, MAPK signaling pathway, and cytokine–cytokine receptor interaction. Two upregulated genes (HSPA1AB and HSPA5) and three downregulated genes (H2-T22, LTA, and H2-K1) were detected in antigen processing and presentation. Six upregulated genes (HSPA1AB, HSPB1, JMJD7, FGFR2, RASGRP3, and DDIT3) and three downregulated genes (DUSP5, MAP3K6, and PLA2G3) were detected in MAPK signaling pathway. Finally, there were 3 upregulated genes (DUSP5, MAP3K6, and PLA2G3) and 3 downregulated genes (IL12RB1, LTA, and CCL17) in cytokine–cytokine receptor interaction. These genes were determined by qRT-PCR as follows: TNFRSF25 (5.79-fold, P = 0.0284), HSPA1AB (5.09-fold, P = 0.0284), HSPB1 (2.12-fold, P = 0.0284), CCR5 (1.99-fold, P = 0.0136), HSPA5 (1.94-fold, P = 0.0133), DDIT3 (1.76-fold, P = 0.0284), FGFR2 (1.69-fold, P = 0.0284), CX3CR1 (1.62-fold, P = 0.0133), RASGRP3 (1.55-fold, P = 0.0284), JMJD7 (1.43-fold, P = 0.0133), H2-T22 (−1.60-fold, P = 0.0136), DUSP5 (−1.61-fold, P = 0.0133), CCL17 (−1.97-fold, P = 0.0133), H2-K1 (−2.58-fold, P = 0.0177), LTα (−3.81-fold, P = 0.0133), IL12RB1 (−4.7-fold, P = 0.0295), PLA2G3 (−4.80-fold, P = 0.0152), and MAP3K6 (−6.67-fold, P = 0.0133; Fig. 4A). According to Western blot results, TNFRSF25 (2.48-fold, P = 0.0004), HSPA1AB (1.98-fold, P = 0.0004), and HSPA5 (1.45-fold, P = 0.0196) proteins were increased significantly at subacute phase compared to acute phase (Fig. 4B, C). Validated enriched KEGG pathways were shown in Fig. 5.

Figure 4.

Validation using quantitative real-time reverse transcription polymerase chain reaction (qRT-PCR) and Western blot at subacute phase compared to acute phase. (A) TNFRSF25, HSPA1AB, HSPB1, CCR5, HSPA5, DDIT3, FGFR2, CX3CR1, RASGRP3, and JMJD7 are significantly upregulated and H2-T22, DUSP5, CCL17, H2-K1, LTα, IL12RB1, PLA2G3, and MAP3K6 are significantly downregulated subacute phase. These genes are involved in cytokine–cytokine receptor interaction, chemokine signaling pathway, and antigen processing and presentation (n = 3 per group; *P < 0.05). The relative expression of target genes from qRT-PCR were calculated using 2−ΔΔCT method. All results were expressed as means ± standard error of the mean. (B) Western blot analysis was performed using antibodies against TNFRSF25, HSPA1AB, HSPA5, and ACTIN (as control). (C) Comparison of relative protein expression from the brain between 3 h post–spinal cord injury (SCI) and 2 wk post-SCI verified by Western blot. TNFRSF25, Tumor necrosis factor receptor superfamily member 25 isoform 1 precursor; CCR5, Chemokine (C-C motif) receptor 5; IL12RB1, Interleukin 12 receptor, beta 1; CCL17, Chemokine (C-C motif) ligand 17; HSPA1AB, Heat shock protein 1AB; HSPA5, Heat shock protein 5; CALR, Calreticulin; PDIA3, Protein disulfide isomerase associated 3; H2-L, Histocompatibility 2, D region locus L; H2-T22, Histocompatibility 2, T region locus 22; LTα, Lymphotoxin A; H2-K1, Histocompatibility 2, K1, K region; HSPB1, Heat shock protein 1; JMJD7, Jumonji domain containing 7; FGFR2, Fibroblast growth factor receptor 2; RASGRP3, RAS, guanyl releasing protein 3; DDIT3, DNA damage-inducible transcript 3 protein; DUSP5, Dual specificity phosphatase 5; MAP3K6, Mitogen-activated protein kinase kinase kinase 6; PLA2G3, Phospholipase A2, group III; CX3CR1, Chemokine (C-X3-C motif) receptor 1. *P < 0.05, ***P < 0.001.

Figure 5.

Illustration of the validated Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. (A) Oxidative phosphorylation pathway. Red stars denote the enriched genes. (B) Cytokine–cytokine receptor interaction pathway. Orange stars denote the enriched genes at subacute phase. Red stars denote the enriched genes at subacute phase compared to acute phase. (C) Chemokine signaling pathway. Orange stars denote the enriched genes at subacute phase. (D) Antigen processing and presentation pathway. Orange stars denote the enriched genes at subacute phase. Red stars denote the enriched genes at subacute phase compared to acute phase. (E) Mitogen-activated protein kinase signaling pathway. Red stars denote the enriched genes at subacute phase compared to acute phase. These pathways mapping were provided from the KEGG pathways of Database for Annotation Visualization, and Integrated Discovery functional annotation tool (DAVID v6.7; http://david.abcc.ncifcrf.gov/).

Discussion

Although several studies have shown that SCI may affect brain,[11,12] the underlying pathophysiological mechanisms have not been completely elucidated. The main purpose of the present study is to evaluate the gene expression profiling in the brain after SCI. Based on a previous report focusing on time point post-injury in a mouse model of SCI[30], we systemically characterized the brain after SCI at acute phase (3 h post-injury) and subacute phase (2 wk post-injury). The different gene expressions and associated enriched KEGG pathways were analyzed and validated to elucidate the progression of pathophysiological mechanism in the brain after SCI.

In the acute phase, ATP6V0C was upregulated due to oxidative phosphorylation, compared to sham control, while COX7B was downregulated by RNA-seq, qRT-PCR (Table 3; Fig. 2A), and Western blot (Fig. 2B, C). Oxidative phosphorylation is an effective metabolic pathway that provides energy by adenosine triphosphate (ATP) synthesis in the mitochondria of cells.[31,32] There are 5 main protein complexes in the electron transport chain, which are nicotinamide adenine dinucleotide hydride (NADH) dehydrogenase (complex I), succinate dehydrogenase (complex II), cytochrome bc1 complex (complex III), cytochrome c oxidase (complex IV), and ATP synthase (complex V). The upregulation of ATP6V0C and downregulation of COX7B resulted in impaired oxidative phosphorylation, which could increase mitochondrial reactive oxygen species (ROS) production from electron transport chain complexes. The CNS has high oxygen consumption, hence particularly susceptible to ROS-induced damage.[33,34] Mitochondrial dysfunction and impaired oxidative phosphorylation may play an important role in the pathogenesis of brain injury after SCI, especially at acute phase.

At subacute phase, cytokine–cytokine receptor interaction (upregulation: TNFRSF25, CCR5, and CCR6; downregulation: HGF and IL12RB1), chemokine signaling pathway (upregulation: GNG11, CCR5, and CCR6; downregulation: CCL17 and NFKBIA), and antigen processing and presentation pathways (upregulation: HSPA1AB, HSPA5, CALR, and PDIA3; downregulation: H2-L) were activated by RNA-seq (Table 3; Fig. 3). Furthermore, antigen processing and presentation (upregulation: HSPA1AB and HSPA5; downregulation: H2-T22, LTα, and H2-K1) and cytokine–cytokine receptor interaction pathways (upregulation: TNFRSF25, CX3CR1, and CCR5; downregulation: IL12RB1, LTα, and CCL17) were also elevated at subacute phase compared to acute phase (Table 3; Fig. 4).

Cytokines and chemokines are very important for CNS immune system interactions. In the nervous system, they function as neuromodulators and regulate neurodevelopment, neuroinflammation, and synaptic transmission. In the brain, cytokines and chemokines are crucial for immune response such as maintaining immunological surveillance, leukocyte traffic modulation, and other inflammatory factor recruitment.[35,36] TNFRSF25 is one of the proinflammatory members, which can stimulate T lymphocyte, B lymphocyte, and antigen presenting cells.[37] TNFSF25 can induce nuclear factor kappa-light-chain-enhancer of activated B cells (NF-κB) via tumor necrosis factor receptor type 1–associated death domain protein (TRADD) and Tumor necrosis factor (TNF) receptor–associated factor 2 (TRAF2).[38] CCR5 is one of the members of chemokine receptor family and plays an important role in the pathogenesis of brain injury and neurodegenerative disorders.[39,40] CCR6 is one of the chemokine receptors that can bind to a single chemokine ligand, CCL20. CCL20/CCR6 is important for the trafficking of T cells to the CNS across the choroid plexus during immune surveillance as well as neuroinflammation.[41] GNG11, membrane bound G-protein subunit, is induced by oxidative stress and may regulate cell senescence.[42] CX3CR1 is a receptor of the chemokine fractalkine (CX3CL1). CX3CL1/CX3CR1 axis can regulate the maintenance of the communication between neuron and microglia in health and disease.[36,43] Especially in the brain of mice, ischemic condition might stimulate neurons to release CX3CL1 and microglia to express CX3CR1.[44] The interaction of CX3CL1/CX3CR1 increased in the brain of ischemic mice. HGF is an angiogenic growth factor, which prevents the extension of several ischemic injuries.[45-48] NFKBIA encodes IκBα, negative regulator of transcription factor NF-κB, which regulates proinflammatory cytokines, chemokines, and is also related to prosurvival and antiapoptosis.[49-51] IL-12 receptor (IL-12R) has 2 subunits, β1 and β2, which are required to bind IL-12 for high affinity. The signaling function of IL-12Rβ1 has not been revealed yet, but it may be needed to maintain high affinity binding with cytokine.[44] LTα, a proinflammatory cytokine with homology to tumor necrosis factor alpha (TNFα), has not been investigated clearly but may contribute to inflammatory processes.[52] CCL17 is a small cytokine, also known as thymus and activation regulation chemokine. CCL17 binds and induces chemotaxis in T cells that may associate with inflammation.[53,54] Thus, neuroinflammation in the brain after SCI may induce activation of cytokine–cytokine receptor interaction, chemokine signaling pathway at subacute phase compared to sham control, as well as cytokine–cytokine receptor interaction also have a role at subacute phase compared to acute phase.

Antigen processing and presentation pathway was activated at subacute phase compared to sham control as follows, upregulation: HSPA1AB, HSPA5, CALR, and PDIA3; downregulation: H2-L (Table 3; Fig. 3), and also activated at subacute phase compared to acute phase as follows—upregulation: HSPA1AB ; downregulation: H2-T22, LTα, and H2-K1 (Table 3; Fig. 4).

Endoplasmic reticulum (ER) stress not only plays vital role in mediating ischemic neuronal cell death but also mediates pathophysiological reactions in brain injuries.[55] ER stress could modulate dysregulation of Major histocompatibility complex I (MHC I) peptide presentation during infected and transformed cells by cluster of differentiation 8 (CD 8) T lymphocytes.[56] The connection between presentation of MHC I-associated peptides and CD8 T cells is tight. MHC I peptides are primarily generated by de novo synthesis and degraded rapidly.[56-58] Following proteasomal degradation, peptides are translocated into ER. The ER responds to the accumulating unfolded proteins by activating intracellular signal transduction pathways called the unfolded protein response (UPR). UPR can regulate MHC I peptide processing, protein translocation, and degradation. Therefore, ER stress affects MHC I peptide presentation. During ER stress, dysregulation of MHC I presentation leads to impaired presentation of peptides derived from proteins.[56] HSPA1AB and HSPA5 play multiple roles in cellular homeostasis. Their expressions are significantly increased in several types of stress conditions and lead to reduced RNA translation and enhanced the degradation of misfolded proteins.[59,60] CALR, as a modulator of the Ca++ balance, is located mostly in the lumen of the ER. During ER stress, CALR expression is upregulated and transferred to outside the cell.[61,62] PDIA3, a member of ER stress proteins, can be induced by oxidative stress conditions. PDIA3 is mediated by redox-sensitive transcription factors, in cellular response to oxidative stress.[63,64] Therefore, ER stress in the brain after SCI may be induced by the activation of antigen processing and presentation at subacute phase compared to sham control and acute phase.

MAPK signaling pathway plays a vital role in ER stress, which was activated at subacute phase as follows—upregulation: HSPA1AB, HSPB1, JMJD7, FGFR2, RASGRP3, and DDIT3; downregulation: DUSP5, MAP3K6, and PLA2G3 (Table 3; Fig. 4). MAPK signaling pathway is activated in response to ER stress and UPR.[65] Particularly, Heat shock proteins (HSPs) are constitutively expressed in certain cell types and are also induced upon exposure of cells to elevated temperatures and other cell stress environments.[66] In MAPK pathway, phosphorylated HSPB1 plays a key role in the induction of several genes related to inflammatory response. The signaling components of the Fibroblast growth factor (FGF) family interact with tyrosine kinase fibroblast growth factor receptors (FGFRs). FGFRs activation leads to phosphorylation of specific tyrosine residues, mediates interaction of cytosolic adaptor proteins, and regulates several intracellular signaling pathways, such as RAS-MAPK, phosphatidylinositol 3-kinase (PI3K)-AKT, phospholipase C gamma (PLC gamma), and signal transducer and activator of transcription (STAT).[67] Activation of FGFR2 may regulate MAPK signaling. In response to ER stress, MAPK signaling pathway is activated in the brain after SCI at subacute phase compared to acute phase.

Taken together, we focused on the brain after SCI pathophysiological events. Our study suggested that gene expression change plays an important role in the pathophysiological process of the injury, and SCI also affects brain injury. In addition, individual genes, which involved in enriched KEGG pathways, will be considered reliable molecular markers in the brain after cell-based therapy in SCI. In further study, investigation should be extended to provide neurobiological mechanism for association between brain and spinal cord. This gene expression change may contribute to new mechanisms and therapeutic targets for CNS disorders.

Conclusions

Our study provided gene expression patterns in the brain after SCI in pathophysiological processes (acute phase and subacute phase). In the brain after SCI, mitochondria dysfunction occurred at acute phase, followed by inflammatory response and ER stress aroused at subacute phase. Finally, these stress environments led to the activation of MAPK signaling pathway at subacute phase (Fig. 6). These pathophysiological mechanisms have already been reported in SCI.[68-70] Our results emphasized that SCI is closely associated with brain injuries. Hence, these mechanisms may provide not only a link between SCI and brain injury but also valuable reference data for understanding gene expression patterns at acute phase and subacute phase.

Figure 6.

Summary of transcriptomic changes in the brain after spinal cord injury (SCI). In the brain after SCI, A; mitochondria function impairs at acute phase (3 h post-injury), B; while inflammatory response and endoplasmic reticulum (ER) stress occur at subacute phase (2 wk post-injury). These cascades lead to the activation of mitogen-activated protein kinase (MAPK) signaling pathway. These mechanisms suggest that SCI and brain injury are closely associated with each other.