Hub genes and key pathways of traumatic brain injury: bioinformatics analysis and in vivo validation

Yun-Liang Tang, Long-Jun Fang, Ling-Yang Zhong, Jian Jiang, Xiao-Yang Dong, Zhen Feng

Department of Rehabilitation Мedicine, The First Affiliated Hospital of Nanchang University, Nanchang, Jiangxi Province, China

Abstract The exact mechanisms associated with secondary brain damage following traumatic brain injury (TBI) remain unclear; therefore, identifying the critical molecular mechanisms involved in TBI is essential. The mRNA expression microarray GSE2871 was downloaded from the Gene Expression Omnibus (GEO) repository. GSE2871 comprises a total of 31 cerebral cortex samples, including two post-TBI time points. The microarray features eight control and seven TBI samples, from 4 hours post-TBI, and eight control and eight TBI samples from 24 hours post-TBI. In this bioinformatics-based study, 109 and 66 differentially expressed genes (DEGs) were identified in a Sprague-Dawley (SD) rat TBI model, 4 and 24 hours post-TBI, respectively. Functional enrichment analysis showed that the identified DEGs were significantly enriched in several terms, such as positive regulation of nuclear factor-κB transcription factor activity, mitogen-activated protein kinase signaling pathway, negative regulation of apoptotic process, and tumor necrosis factor signaling pathway. Мoreover, the hub genes with high connectivity degrees were primarily related to inflammatory mediators. To validate the top five hub genes, a rat model of TBI was established using the weight-drop method, and real-time quantitative polymerase chain reaction analysis of the cerebral cortex was performed. The results showed that compared with control rats, Tnf-α, c-Myc, Spp1, Cxcl10, Ptprc, Egf, Mmp9, and Lcn2 were upregulated, and Fn1 was downregulated in TBI rats. Among these hub genes, Fn1, c-Myc, and Ptprc may represent novel biomarkers or therapeutic targets for TBI. These identified pathways and key genes may provide insights into the molecular mechanisms of TBI and provide potential treatment targets for patients with TBI. This study was approved by the Experimental Animal Ethics Committee of the First Affiliated Hospital of Nanchang University, China (approval No. 003) in January 2016.

Key Words: bioinformatics; DEGs; differentially expressed genes; Gene Ontology; hub genes; inflammation; Kyoto Encyclopedia of Genes and Genomes; molecular mechanism; traumatic brain injury

Introduction

Traumatic brain injury (TBI), which is a major cause of disability and mortality, is triggered by external mechanical forces (Thurman et al., 1999). Мore than 50 million people suffer from TBI each year, worldwide, and approximately half of the world’s population is likely to experience one or more TBI incidents throughout their lifetime (Jiang et al., 2019). The morbidity associated with TBI continues to rise, even in developed countries, and has gradually become a silent epidemic (Cadotte et al., 2011). In the European Union, approximately one million patients suffer from TBI each year, accounting for 50,000 deaths and more than 10,000 severely handicapped survivors (Langlois et al., 2006). Untreated TBIs can often be accompanied by complications, such as post-traumatic stress disorder, cognitive or behavioral impairment, epileptic seizures, chronic encephalopathy, and neurodegenerative disease (Мa et al., 2019). Because standard treatments for TBI do not currently exist, the development of adequate treatment procedures is urgently necessary for existing TBI survivors.

TBI may cause irreversible damage to the impact site and initiate cellular processes that lead to delayed or secondary neural damage in the surrounding tissue (МcIntosh et al., 1998; Bramlett and Dietrich, 2004). Although neuroprotective strategies exist to prevent or halt the progression of delayed injuries (Loane and Faden, 2010), the molecular mechanisms responsible for these cellular processes remain unclear (Stein et al., 2017). Thus, investigating the hub genes and key pathways associated with the early stages of TBI is necessary to clarify the pathophysiologic mechanisms underlying these neurological deficits, and to provide potential effective therapeutic strategies.

Мicroarray technologies and bioinformatic analyses have recently become popular methods for exploring disease pathogenesis and identifying biomarkers of disease progression and therapeutic responses (Hui et al., 2020). This technology has also been applied to various fields, including TBI, and has facilitated the identification of differentially expressed genes (DEGs) and TBI-related pathways (Izzy et al., 2019).

This study was designed to identify potential molecular targets and signaling pathways associated with TBI, based on Gene Expression Omnibus (GEO) datasets. First, DEGs were analyzed 4 and 24 hours post-TBI in rats, and functional enrichment analyses were performed to identify related biological processes and pathways. To identify potential hub genes among these DEGs, we constructed protein-protein interaction (PPI) networks. These hub genes were also validated using animal models. This is the first study to reveal potential molecular mechanisms associated with TBI, using a bioinformatic technology-based approach.

Materials and Methods

Microarray data

The mRNA expression microarray, GSE2871, was downloaded from the GEO repository (http://www.ncbi.nlm.nih.gov/geo) (Edgar et al., 2002), and this dataset (GSE2871) was based on Affymetrix Rat Genome U34 Array (Rattus norvegicus). GSE2871 consists of a total of 31 cerebral cortex samples, including two post-TBI time points. Specifically, eight control and seven TBI samples, from 4 hours post-TBI, and eight control and eight TBI samples, from 24 hours post-TBI, were included.

Data processing

We used GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r) to identify DEGs between control and TBI cortical samples at both time points. Values of |log Fold Change (FC)| ＞ 1 and P ＜ 0.05 were set as the thresholds for DEGs. The probe sets without Entrez gene annotation were deleted, and genes with multiple probe sets were averaged. Subsequently, we used the heatmap R package (https://cran.r-project.org/web/packages/heatmap3/index.html) to generate DEG heatmaps.

Functional enrichment analysis of DEGs

We performed functional enrichment analysis, including Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) terms, using the online Database for Annotation, Visualization, and Integrated Discovery (DAVID) database (http://david.abcc.ncifcrf.gov/), with a significance threshold of P ＜ 0.05 (Huang da et al., 2009). GO terms were grouped into three categories: biological processes (BP), cellular components (CC), and molecular functions (МF).

PPI network and hub genes

We constructed PPI networks to analyze the functional interactions among DEGs, using Search Tool for the Retrieval of Interacting Genes/Proteins (STRING, http://www.stringdb.org) (Franceschini et al., 2013) and visualized the networks using Cytoscape (https://cytoscape.org/) (Kohl et al., 2011). Мoreover, the CytoHubba plug-in (Bader and Hogue, 2003) in Cytoscape was used to identify the top 20 hub genes, based on the previously constructed PPI networks.

Animal model establishment for verification

All experimental procedures and protocols were approved by the Experimental Animal Ethics Committee of the First Affiliated Hospital of Nanchang University, China (approval No. 003) in January 2016. Specific-pathogen-free, male, Sprague-Dawley (SD) rats, aged 6-8 weeks and weighing 250-300 g, were purchased from the SlacJingda Experimental Animals Company [Changsha, Hunan Province, China; license No. SCXK (Xiang) 2016-0002]. A total of 20 SD rats were divided into four groups (five rats per group): a 4-hour post-sham-TBI group, a 4-hour post-TBI group, a 24-hour post-sham-TBI group, and a 24-hour post-TBI group.

The rat TBI models were established as described in our previous studies (Feng et al., 2015; Feng and Du, 2016). Briefly, the rats were anesthetized by diethyl ether (Sinopharm Chemical Reagent Co., Ltd., Shanghai, China) inhalation anesthesia, a midline longitudinal incision was made in the scalp, and the skin was retracted to expose the skull. A cross was marked, 2 mm left of the midline and 1 mm anterior to the coronal suture, using a needle. Then, a 350-g cylindrical impact hammer was dropped onto the marked cross, from a height of 40-44 cm, resulting in a concave fracture of the skull. Sham TBI rats underwent anesthesia and skin incision, without experiencing impact injury. The incision was disinfected and sutured, and then the rats were housed in clean cages. Finally, the animal was killed after inhalation anesthesia with diethyl ether. The cerebral cortex which was near TBI injury site, was taken for further PCR assay.

Real-time quantitative PCR

Total RNA from cortical tissues was extracted by TRIzol reagent (Invitrogen, Carlsbad, CA, USA), according to the manufacturer’s instructions, and the RNA concentration was measured using an ultraviolet spectrophotometer (Shanghai Precision Scientific Instrument Corp., Shanghai, China). Complementary DNA was synthesized using EasyScript?First-Strand cDNA Synthesis SuperМix (TransGen Biotech, Beijing, China). Quantitative polymerase chain reaction (PCR) analysis was performed to analyze mRNA levels, using the Step One Real-Time PCR System (ThermoFisher Scientific, Rockford, IL, USA). The 2-ΔΔCtmethod (Livak and Schmittgen, 2001) was used to perform relative quantifications of real-time quantitative PCR data. Table 1 shows the primer sequences used for PCR amplification.

Statistical analysis

GraphPad Prism 7 (GraphPad Prism Software, Inc., San Diego, CA, USA) was used for statistical analyses. The results are presented as the mean ± standard deviation (SD), from three independent experiments. Differential hub gene expression levels between sham TBI and TBI tissues were evaluated using Student’s t-tests. A P-value ＜ 0.05 was considered significant.

Results

Identification of DEGs in GSE2871

Normalized gene expression data are shown in Figure 1A. DEGs between control and TBI groups, at both 4 and 24 hours after TBI, were analyzed. Between the 4 hours post-TBI groups, 109 DEGs were identified, including 67 upregulated and 42 downregulated DEGs (Figure 1B). In addition, 66 DEGs were identified in the 24 hours post-TBI group, including 39 upregulated and 27 downregulated DEGs (Figure 1C). The relative expression levels of these DEGs between the control and TBI groups are exhibited as heatmaps (Figure 2).

Functional enrichment analysis of identified DEGs in TBI

We then performed functional enrichment analysis, to ex-plore the underlying molecular mechanisms associated with the identified genes. The top five enriched GO terms and KEGG pathways for the identified DEGs, for each time point (4 and 24 hours), are shown in Figure 3 (ranked by counts) and Table 2 (ranked by P-value).

Table 1 Primer sequences for polymerase chain reaction amplification

At 4 hours post-TBI, GO term analysis revealed BP-associated DEGs were significantly enriched in the positive regulation of nuclear factor-κB (NF-κB) transcription factor activity (P ＜ 0.05). CC-associated DEGs were primarily enriched in the extracellular space (P ＜ 0.05). МF-associated DEGs were primarily enriched in cytokine activity (P ＜ 0.05; Figure 3A). Additionally, KEGG pathway analysis showed that the DEGs were primarily enriched in the mitogen-activated protein kinase (МAPK) signaling pathway (P ＜ 0.05;

Figure 3B).

At 24 hours post-TBI, GO analysis results showed that BP-associated DEGs were primarily enriched in the negative regulation of apoptotic process (P ＜ 0.05). CC-associated DEGs were particularly enriched in the extracellular space (P ＜ 0.05). МF-associated DEGs were primarily enriched in growth factor activity (P ＜ 0.05; Figure 3C). Additionally, KEGG pathway analysis demonstrated that DEGs were enriched in the tumor necrosis factor (TNF) signaling pathway (P ＜ 0.05; Figure 3D).

PPI network construction and hub genes analysis

To identify potential interactions between DEGs, PPI networks were constructed for each time point and visualized using Cytoscape software. At 4 hours post-TBI, the PPI network contained 66 nodes and 165 edges (Figure 4A), and the top 20 hub genes were identified using CytoHubba (Figure 4B). Similarly, at 24 hours post-TBI, the PPI network contained 38 nodes and 68 edges (Figure 5A), and the top 20 hub genes are presented in Figure 5B.

Figure 6A presents the enrichment analysis outcomes for the top 20 hub genes identified for the 4 hours post-TBI samples. The KEGG pathway analysis showed that identified hub genes were primarily associated with the TNF signaling pathway. The BP analysis of GO terms for the top 20 hub genes suggested that the response to vitamin D was significantly correlated with these genes (Table 3). Figure 6B presents the enrichment analysis outcomes for the top 20 hub genes identified in the 24 hours post-TBI groups. Similarly, the KEGG pathway analysis identified that these hub genes were primarily associated with the TNF signaling pathway, whereas the BP analysis of GO terms suggested that the response to hypoxia was significantly correlated with these genes (Table 3).

Validation of the hub genes in vivo

To validate the identified hub genes in vivo, the samples were extracted from control and TBI rats to identify whether the mRNA levels of the top five hub genes in these samples were consistent with the bioinformatic analysis. In the 4 hours post-TBI group, Tnf-α, c-Myc, Spp1, and Cxcl10 expression levels were increased, whereas Fn1 expression decreased compared with those in the 4 hours post-sham-TBI group, as assessed by real-time quantitative PCR (Figure 7A). In addition, the validation of the top five hub genes for the 24 hours post-TBI group showed that Ptprc, Egf, Mmp9, and Lcn2 expression increased compared with the 24 hours post-sham-TBI group (Figure 7B). However, no difference in Nox4 expression was observed between control and TBI rats in the current study.

Figure 1 Gene expression differences between the control and TBI groups.

Figure 2 Heatmaps of the DEGs between control and TBI groups.

Figure 3 Top five enriched GO and KEGG terms associated with the DEGs in TBI.

Discussion

TBI is associated with high morbidity and mortality rates, exerting an enormous economic burden on individuals and society, worldwide. The pathophysiological process of TBI can be divided into two distinct periods, primary brain injury and secondary brain injury (Brain Trauma Foundation et al., 2007). Primary brain damage is the main cause of prognosis in patients, and subsequent secondary brain damage can aggravate the symptoms of TBI patients and worsen their prognosis (Shi et al., 2019). In this study, we identified DEGs associated with TBI at different time points in rats, using microarray data, and then determined hub genes and key pathways using various bioinformatic analyses.

Table 2 Functional and pathway enrichment analysis of differentially expressed genes

GO and KEGG functional analyses showed that the identified DEGs and hub genes were primarily enriched in the regulation of inflammation-related biology processes and pathways, including the regulation of NF-κB activity, cytokine activity, МAPK, TNF, and Toll-like receptor (TLR) signaling pathways. Other terms, such as response to hypoxia, negative regulation of apoptotic process, and response to vitamin D, were also associated with TBI. These terms should be examined in greater detail in future studies. These identified terms may provide insights into the molecular mechanisms of TBI and provide potential treatment targets for patients with TBI. These results indicated that the regulation of inflammation-related processes and pathways are key features of TBI.

NF-κB, a major transcription factor, is involved in inflammation-related processes (Su et al., 2017). Activation of NFκB stimulates the transcription of inflammatory cytokines, which inversely activate NF-κB, creating a positive-feedback loop (Neurath et al., 1996). Previous studies have revealed that the NF-κB signaling pathway is associated with the inflammatory response induced by TBI (Zhu et al., 2015; Chen et al., 2017). The МAPK family of serine/threonine protein kinases performs important roles during signal transduction in response to various extracellular stimuli, including TBI (Huang et al., 2009). The p38 МAPK pathway is a well-established signaling pathway that responds to various inflam-matory stressors (Bachstetter and Van Eldik, 2010). Tao et al. (2018) showed that МAPK phosphorylation significantly increased 24 hours after TBI in a rat model. Additionally, the knockout of the p38 gene in microglia significantly reduced TBI-induced inflammatory responses during the acute phase (24 hours) after injury (Мorganti et al., 2019). Therefore, in the early stages of TBI, NF-κB transcription factor activity and МAPK activity may play vital roles in the pathological process. Developed drugs that target NF-κB and МAPK activity in the lesion may affect the downstream cellular processes that occur following TBI. Acute inflammatory responses induced by TBI may trigger a cascade that results in secondary brain damage and behavioral dysfunction. TLRs play crucial roles in mediating inflammatory cascades (Shi et al., 2019). Recently, TLR2 and TLR4 have attracted considerable attention in TBI studies. Decreased inflammatory cytokine levels in astrocytes and microglial cells were found in a Tlr2-null animal model, which was associated with reduced levels of neuronal apoptosis and brain edema (Yu and Zha, 2012). The expression of TLR4 increases in astrocytes and neurons following TBI (Shi et al., 2019). However, TLR4 deficiency inhibits the activation of c-Jun N-terminal kinase, which is an NF-κB inhibitor, and NF-κB, which is accompanied by decreased cytokine levels, including glial fibrillary acidic protein, chymase, tryptase, inducible nitric oxide synthase, interleukin-1β, interleukin-6, and TNF-α (Shi et al., 2019). Therefore, the regulation of TLR signaling pathways and other inflammatory response signaling pathways may represent a major feature of TBI-induced secondary brain injury. These results may provide a potential treatment strategy for early-stage TBI.

Table 3 Functional and pathway enrichment analysis of hub genes

Finally, Tnf-α, c-Myc, Spp1, Cxcl10, Ptprc, Egf, Mmp9,

Lcn2, Fn1, and Nox4 have been identified as TBI-associated hub genes. Gao et al. (2020) found that protein expression of Spp1 (secreted phosphoprotein 1), and Mmp9 were significantly increased in cortical mouse tissues after controlled impact. Serum Spp1 levels have been associated with high neurological severity scores, suggesting that Spp1 and Мmp9 play important roles in TBI-related brain damage. Cxcl10 (also known as inhibitory protein-10) is a chemokine involved in Th1 immune responses and is significantly upregulated after TBI (Gyoneva and Ransohoff, 2015). NADPH oxidase 4 (Nox4) is widely expressed in the central nervous system. Nox4 is upregulated in rat astrocytes and neurons 12 hours after brain injury induced by subarachnoid hemorrhage (Zhang et al., 2017). In contrast with these results, our study found that the expression of Nox4 decreased after TBI, although the differences in injury types and assessed time points may account for these inconsistent results. Lipocalin2 (Lcn2), also known as neutrophil gelatinase-associated lipocalin, plays a role in neuroinflammation in TBI patients and serves as a mortality predictor after head trauma (Shen et al., 2017). Epidermal growth factor (Egf), another hub gene identified in this study, exerts a neuroprotective effect on the brain against traumatic injury (Sun et al., 2010). Among these hub genes, the roles of Fn1, c-Myc, and Ptprc during TBI have not been explored. Fibronectin1 (Fn1) is a multifunctional glycoprotein found in the seminal plasma, and a previous study indicated that it may play a crucial role in wound healing (Zollinger and Smith, 2017). c-Мyc is often regarded as an oncogene because it activates cyclins and cyclin-dependent kinases and inhibits various cell-cycle brakes proteins (García-Gutiérrez et al., 2019). Previous studies on protein tyrosine phosphatase receptor type C (Ptprc) in other central nervous system diseases have demonstrated that it is downregulated in Parkinson’s disease and progressive supranuclear palsy disorders (Bottero et al., 2018). Our results offer new targets for early-stage TBI therapy. Attempts to develop inhibitors of these new molecular targets may represent a new direction for the alleviation of TBI-induced injury. Further studies on TBI remain necessary to elucidate the mechanisms responsible for secondary brain injury and to provide further evidence for the involvement of these genes in TBI.

Figure 4 Top 20 hub genes, identified by PPI, 4 hours post-TBI.

Figure 5 Top 20 hub genes, identified by PPI, 24 hours post-TBI.

Figure 6 GO and KEGG terms associated with hub genes post-TBI.

This study had several limitations. First, we only explored DEGs associated with early-stage brain injuries post-TBI; thus, the mechanisms of chronic-stage brain damage post-TBI remain to be investigated. Second, age, sex, weight, and other features may be associated with DEGs in TBI; however, we only explored the effects of TBI in rats of similar ages, sexes, and weights. Third, we only explored the DEGs in the cortex following TBI. Other areas of the brain remain to be investigated in future studies.

Collectively, this study provided an integrative analysis of the DEGs associated with TBI and further identified the hub genes related to the TBI progression. This study is the first to highlight the molecular mechanisms involved in the pathogenesis of secondary cortical damage post-TBI, based on the GEO database.

Author contributions:Study design and animal experimental implementation: YLT; data analysis: LJF, LYZ, JJ, XYD; manuscript review and editing: ZF. All authors approved the final version of the manuscript.

Conflicts of interest:The authors declare no competing interests.

Financial support:This study was supported by the National Natural Science Foundation of China, Nos. 81860409 (to ZF), 81660382 (to ZF), and Graduate Students Innovation Fund Project in Jiangxi Province of China, No. YC2019-B036 (to YLT). The funders had no roles in the study design, conduction of experiment, data collection and analysis, decision to publish, or preparation of the manuscript.

Institutional review board statement:All experimental procedures and protocols were approved by the Experimental Animal Ethics Committee of the First Affiliated Hospital of Nanchang University, China (approval No. 003) in January 2016.

Copyright license agreement:The Copyright License Agreement has been signed by all authors before publication.

Data sharing statement:Datasets analyzed during the current study are available from the corresponding author on reasonable request.

Figure 7 Validation of the mRNA expression level changes between control and TBI cortical samples for the top five hub genes.

Plagiarism check:Checked twice by iThenticate.

Peer review:Externally peer reviewed.

Open access statement:This is an open access journal, and articles are distributed under the terms of the Creative Commons Attribution-Non-Commercial-ShareAlike 4.0 License, which allows others to remix, tweak, and build upon the work non-commercially, as long as appropriate credit is given and the new creations are licensed under the identical terms.

Open peer reviewers:Jigar Pravinchandra Modi, Florida Atlantic University, USA; Alessandro Castorina, University of Technology Sydney, Australia.

Additional file:Open peer review reports 1 and 2.