Next Article in Journal
A Conservative Combined Laser Cryoimmunotherapy Treatment vs. Surgical Excision for Basal Cell Carcinoma
Next Article in Special Issue
The Prognostic Value of Deleted in Colorectal Cancer (DCC) Receptor and Serum Netrin-1 in Severe Traumatic Brain Injury
Previous Article in Journal
Anaesthetic Approach to Enhanced Recovery after Surgery for Kidney Transplantation: A Narrative Review
Previous Article in Special Issue
Decreased Expression of CIRP Induced by Therapeutic Hypothermia Correlates with Reduced Early Brain Injury after Subarachnoid Hemorrhage
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Comprehensive RNA Expression Analysis Revealed Biological Functions of Key Gene Sets and Identified Disease-Associated Cell Types Involved in Rat Traumatic Brain Injury

1
Department of Neurosurgery, Qilu Hospital, Cheeloo College of Medicine and Institute of Brain and Brain-Inspired Science, Shandong University, Jinan 250012, China
2
Shandong Key Laboratory of Brain Function Remodeling, Jinan 250012, China
3
Department of Nuclear Medicine, Qilu Hospital, Shandong University, Jinan 250012, China
4
Department of Neurosurgery, Children’s Hospital Affiliated to Shandong University, Jinan 250012, China
*
Author to whom correspondence should be addressed.
J. Clin. Med. 2022, 11(12), 3437; https://doi.org/10.3390/jcm11123437
Submission received: 19 May 2022 / Revised: 6 June 2022 / Accepted: 11 June 2022 / Published: 15 June 2022
(This article belongs to the Special Issue Advances in Neurotrauma)

Abstract

:
Traumatic brain injury (TBI) is a worldwide public health concern without major therapeutic breakthroughs over the past decades. Developing effective treatment options and improving the prognosis of TBI depends on a better understanding of the mechanisms underlying TBI. This study performed a comprehensive analysis of 15 RNA expression datasets of rat TBIs from the GEO database. By integrating the results from the various analyses, this study investigated the biological processes, pathways, and cell types associated with TBI and explored the activity of these cells during various TBI phases. The results showed the response to cytokine, inflammatory response, bacteria-associated response, metabolic and biosynthetic processes, and pathways of neurodegeneration to be involved in the pathogenesis of TBI. The cellular abundance of microglia, perivascular macrophages (PM), and neurons were found to differ after TBI and at different times postinjury. In conclusion, immune- and inflammation-related pathways, as well as pathways of neurodegeneration, are closely related to TBI. Microglia, PM, and neurons are thought to play roles in TBI with different activities that vary by phase of TBI.

1. Introduction

Traumatic brain injury (TBI) is common worldwide. Over 27 million traumatic brain injuries occur globally each year [1]. TBI is responsible for 30–40% of all injury-related deaths, and is by far the leading cause of disability associated with neurological diseases, accounting for 2–3 times more disabilities than Alzheimer’s disease or cerebrovascular disease [2]. There is even evidence that TBIs are associated with dementia years later [3]. Survivors endure enormous psychological, physical, and emotional pain, while their families and societies face enormous burdens. Although various therapeutic attempts have been made to improve the outcome of TBI, most multicentral clinical trials of medical and surgical interventions have failed to show efficacy [2]. For this reason, it is important to gain a holistic understanding of the mechanisms underlying TBI to come up with optimal treatment options.
TBI is not a single pathophysiological event but a complex disease process. For understanding the primary and secondary injury mechanisms in TBI, a number of preclinical animal models have been developed. Rat is one of the most widely used animals in TBI research due to its modest cost, small size, and standardized outcome measurements [4]. Through the research based on experimental models, many pathophysiological processes of TBI have been better understood, including disturbance in neurotransmitters and calcium signaling pathways [5]; abnormal lipids, proteins, and nucleic acids oxidation [6]; upregulation of transcription factors and inflammatory mediators [7]; and increased expression of detrimental cytokines, which induce brain edema, blood–brain barrier damage, and cell death [8].
With the advances in microarray and high-throughput sequencing techniques, large and growing public databases of TBI gene expression data are being deposited into public databases. The transcriptome analysis of gene expression has been used to identify pathways potentially involved in TBI. In most cases, gene expression profile data are analyzed by focusing on genes that differ between TBI and control groups while ignoring other genes that may also associated with sample features. Weighted gene coexpression network analysis (WGCNA) is a bioinformatics algorithm used for exploring gene association patterns in samples, identifying gene sets with highly coordinated expression, and exploring biologically meaningful gene sets connected to a particular trait [9]. In this study, transcriptional profiling data of rat brain tissues after TBI were analyzed using WGCNA. Feature-relevant modules containing coexpressed genes were identified, and their gene ontology (GO) functions and signaling pathways that may be involved were investigated. These gene sets were further subjected to gene set enrichment analysis (GSEA) to recognize cell types associated with the sample traits, which was further validated using gene set variation analysis (GSVA) by independent datasets. The results may provide a reference for the mechanism research and treatment of TBI.

2. Materials and Methods

2.1. Data Acquisition and Preprocessing

The Gene Expression Omnibus (GEO) database (accessed on 25 October 2021 from https://www.ncbi.nlm.nih.gov/geo/) was searched using the following terms: traumatic brain injury, TBI, brain trauma, or neurotrauma. A further filter was performed with the organism “Rattus norvegicus” and the study type “Expression profiling by array” or “Expression profiling by high throughput sequencing”. Search results were manually checked, datasets without sham-injury-treated control samples were excluded, and 14 datasets containing rat brain tissue samples were included in this study. Additionally, the single-nuclei sequencing dataset GSE137869 [10] was retrieved for cell marker identification. Among the datasets, GSE2871 [11] was used for WGCNA; GSE2392 [12], GSE2871 [11], and GSE45997 [13] were used in the GSEA analysis; GSE1911 [14], GSE2392 [12], GSE24047 [15,16], GSE31357 [17], GSE59645 [18], GSE64978 [19], GSE67836 [20], GSE68207 [19], GSE80174 [21], GSE86579 [22], GSE111452 [23], and GSE115614 [18] were analyzed for validation. Among the samples in the above datasets, only samples from wild-type rats with no additional treatment or disease other than TBI were included. Brain tissues from young rats in GSE137869 were used for single-nucleus transcriptome analysis. The sample traits were determined based on groupings and sample information in the database. If applicable, background correction and normalization were conducted using the R package limma (version 3.46.0) [24]. Detailed information is shown in Table 1.

2.2. WGCNA

To identify the gene modules relevant to TBI, the R package WGCNA (version 1.70) [25] was used to conduct the weighted co-expression network analysis. A power of 6, which enabled the scale-free topology fit index to reach 0.85, was selected as soft-threshold parameters to construct a signed, scale-free coexpression gene network. Thereafter, modules of coexpressed genes were identified by hierarchical clustering, and the minimum size of modules was set to 40 genes. Furthermore, the module eigengene (ME) representing each module’s expression profiles was calculated, and intramodular correlations and module–trait associations were estimated. Modules with high module–trait significance (p-value < 0.01) were defined as key modules and subjected to further analysis. The result of intramodular correlations and module–trait relationship analyses were plotted using the R package ggcor (version 0.9.8).

2.3. GO and KEGG Enrichment Analysis

Gene Ontology (GO) enrichment analysis provides a structured description of the known biological information of genes at different levels: biological process (BP) refers to a biological objective to which the gene or gene product contributes, cellular component (CC) refers to the place in the cell where a gene product is active, and molecular function (MF) is defined as the biochemical activity of a gene product [26]. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis assigns functional meanings to genes and genomes both at the molecular and higher levels [27]. GO enrichment analysis and KEGG analysis were performed for understanding the biological functions and pathways involved in genes in key modules. Both enrichment analyses were implemented in the R package clusterProfiler (version 3.18.1) [28]. A p-value < 0.05 and a q-value < 0.05 was considered significant. In each analysis, the top 10 results were extracted for visualization.

2.4. Computational Analysis of snRNA Seq Datasets

The single-nucleus transcriptome analysis was performed in the R package Seurat (version 4.1.0) [29,30,31,32]. The quality control process was as follows. Nuclei containing more than 2000 expressed genes or those that contained less than 200 expressed genes were removed. Data on nuclei that contained more than 2.5% mitochondrial genes were filtered. Features expressed in three or fewer nuclei were excluded. After logarithmical normalization of the filtered nuclei data, principal component (PC) analysis was performed. The first 15 PCs, determined using a combination of jackstraw and elbow methods, were used to generate clusters with a resolution of 1.5. For visualization, the nonlinear dimensional reduction was performed with the t-distributed stochastic neighbor embedding (t-SNE) algorithm. Marker genes for the individual clusters were identified using the FindAllMarkers function with default parameters. Annotation of cell types was manually conducted according to a previous study [10], and clusters identified as the same cell type were merged.

2.5. Gene Set Enrichment Analysis

To identify the key cell types associated with sample traits in TBI, gene set enrichment analysis (GSEA) was conducted. First, log2FC values representing the expression change between the compared samples were calculated using the R package limma (version 3.46.0). Then, the key modules related to the same sample trait from WGCNA were merged, and the genes contained in the modules were ranked according to log2FC. The rank list was used as input to GSEA, which was performed based on the cell type marker genes that could be detected and annotated in GSE2871 using the R package clusterProfiler (version 3.18.1). A p-value < 0.05 and a q-value < 0.25 was considered significant.

2.6. Gene Set Variation Analysis

Gene set variation analysis (GSVA) is a gene set enrichment method that estimates changes in gene set enrichment over the samples independently of any class label, and has emerged as an overall top method to assign cell type labels in single-cell RNA-sequencing analysis [33,34]. Expression statistics of the cell type markers are summarized into a single enrichment score for each cell type. GSVA was employed to estimate the abundance or activity of GSEA-enriched cell types in the validation datasets. R package GSVA (version 1.38.2) [26] was used to score individual samples based on the top 20 markers selected based on the likelihood-ratio of microglia, PM, and neurons. Enrichment scores were compared between differently treated samples within the same database. To aggregate enrichment scores from different datasets and make comprehensive comparisons within a wider range of time, the enrichment scores were normalized to the corresponding sham group as the relative enrichment score.

2.7. Statistical Analysis

Data analysis and plotting of the results were performed using R software (version 4.0.2) and GraphPad Prism (GraphPad Prism 8; GraphPad). Nonparametric test or t-test based on data distribution characteristics was used to distinguish the difference between the two groups, and a p-value of < 0.05 was considered significant.

3. Results

3.1. WGCNA Identified Key Modules Related to Sample Traits in TBI

Rat TBI dataset GSE2871 was downloaded and sample information was obtained from GEO. Briefly, adult rats were subjected to lateral fluid percussion injury (mild or severe) or sham surgery without injury. Expression profiling of brain regions (parietal cortex and hippocampus, ipsilateral and contralateral to injury) was conducted at 4 h or 24 h postinjury. All 8799 genes from 47 samples in GSE2871 were subjected to WGCNA. A scale-free network was constructed as described in Materials and Methods, and a total of 19 gene coexpression modules were obtained. Among the 19 modules, the turquoise module was the largest, which contained 2109 genes, while the light green module containing 63 genes was the smallest one. A total 930 ungrouped genes were included in the grey module (Figure 1A).
The correlations among modules and the association between modules and traits were estimated. Key modules in TBI were defined as modules with a high module–trait significance (p-value < 0.01). The modules associated with the severity of injury were the magenta module and the tan module, which contained a total of 265 annotated genes. The black, turquoise, magenta, and tan modules were linked to the sampling side (ipsilateral or contralateral to injury) and contained 2677 genes that were annotated. The green, yellow, brown, and tan modules, with a total of 2800 annotated genes, were associated with postinjury time. The green, brown, and blue modules were related to brain region (Figure 1B).

3.2. Function Enrichment Analysis of Key Modules

The sample traits “severity of injury”, “sampling side”, and “postinjury time” were considered as key features associated with TBI, and modules related to these sample traits were further analyzed. As part of our investigation of the biological functions of TBI-related genes, GO and KEGG enrichment analyses were performed on genes in key modules associated with each sample trait. The GO results showed that genes related to injury severity were mainly involved in response to cytokine, inflammatory response, and other immune-related processes (Figure 2A). Side-related genes largely played a role in the inflammatory response, bacteria-associated response, as well as involved in transmembrane signaling receptor activity (Figure 2B). Genes associated with postinjury time were mainly related to the metabolic and biosynthetic processes, cell junction organization, etc. (Figure 2C).
The KEGG analysis indicated that injury-severity-related genes were mainly enriched in lipid and atherosclerosis, and various virus and parasitic infection pathways (Figure 2D). Genes related to the sampling side were primarily involved in neuroactive ligand–receptor interaction, calcium signaling pathway, chemical carcinogenesis-receptor activation, etc. (Figure 2E). Those genes associated with postinjury time were mainly involved in multiple neurodegeneration pathways (Figure 2F).
Although functional enrichment analyses yielded a variety of biological functions and pathways, it is worth noting that immune- and inflammation-related terms were commonly enriched in several analyses. An interesting question then is what role various cells, especially immune-related cells, play in TBI. Therefore, on the basis of the above results, we continued to explore the changes in TBI-related cellular activity.

3.3. Identification of the Markers of Rat Brain Cell Types

To identify cell types in rat brains and find the corresponding markers, we collected snRNA-seq data of two rat brain tissue samples without any treatments or disease from the dataset GSE137869 and conducted the computational analysis as described in Materials and Methods. A total of 8889 nuclei with 16,882 features were subject to the analysis, which yielded 22 cell clusters. Of the 22 clusters obtained, six clusters were annotated as neurons, five clusters were annotated as oligodendrocytes, three clusters were annotated as oligodendrocyte progenitor cells (OPC), two clusters were annotated as microglia, and two clusters were annotated as astrocytes. The remaining four clusters were annotated as pericyte, endothelial cells (EC), perivascular macrophages (PM), and vascular leptomeningeal cells (VLMC), respectively. Clusters of the same cell type were merged. For each cell type, gene markers were identified for subsequent analysis (Figure 3). The analysis resulted in 148 marker genes for astrocytes, 98 marker genes for EC, 92 marker genes for microglia, 592 genes for neurons, 10 marker genes for oligodendrocytes, 229 marker genes for OPC, 94 genes for pericytes, 139 genes for PM, and 152 genes for VLMC.

3.4. Characterization of Key Cell Types Associated with Traits

To explore changes in cellular abundance and activity related to a sample trait, we performed a GSEA analysis of genes in key modules related to that trait using markers of cell types. Briefly, differential expression analysis was conducted on expression data from samples for the corresponding traits, and GSEA was performed on a gene list presorted by the fold-change value of the differential expression analysis. The enriched cell types were considered to be associated with the corresponding sample traits. Differential expression analysis was conducted between samples from the lesioned side and samples from sham-treated animals in GSE45997 to identify cell types associated with the severity of injury. Samples from the ipsilateral and contralateral side of the injury in GSE45997 were compared to rank the gene list, which was input into GSEA to locate side-related cell types. Similarly, samples of 24 h postinjury and 4 h postinjury in GSE2392 were used to identify time-sensitive cell types. For the severity of injury and side of injury, the above analysis was repeated by collecting corresponding trait samples from GSE2871. However, no replication was performed for postinjury time because, in GSE2871, the number of ipsilateral samples at 24 h postinjury was insufficient (N = 2) for given injury severity.
PM was enriched in the side-associated modules in both GSEA, and microglia was enriched in one of the analyses. The above results indicate more PM and microglia on the ipsilateral side of injury (Figure 4A,C,E,G,H). As for the postinjury time, neuron was enriched with an acceptable significant level (p-value = 0.047 and q-value = 0.227), indicating possibly decreased neuron abundance or activity 24 h postinjury compared with that at 4 h (Figure 4B,D). For the key modules related to injury severity, either no terms were enriched, or the enriched cell types did not reach significance in GSEA (Figure 4F).

3.5. Validation of Cell Activity after Traumatic Brain Injury

To validate the change in abundance or activity of microglia, PM, and neurons, GSVA was performed on the other datasets containing TBI samples and corresponding sham controls collected at various postinjury times, including 30 min (GSE2392), 3 h (GSE1911 and GSE24047), 4 h (GSE2392 and GSE31357), 6 h (GSE24047), 8 h (GSE2392), 12 h (GSE24047), 24 h (GSE1911, GSE2392, GSE31357, GSE59645, GSE111452, and GSE115614), 48 h (GSE24047), 72 h (GSE2392), 1 week (GSE64978 and GSE68207), 2 weeks (GSE92363 and GSE111452), 3 weeks (GSE2392), 1 month (GSE67836), 3 months (GSE80174, GSE86579, and GSE111452), 6 months (GSE111452), and 1 year (GSE111452). Enrichment scores between TBI condition(s) and the corresponding control(s) in each dataset were compared and the results were plotted.
Results from some datasets showed decreased enrichment scores of microglia at 3 h, 4 h, 6 h, 12 h, and 1 year after TBI, while the results from some datasets suggested an increased abundance of microglia at 4 h, 24 h, 72 h, 1 week, 2 weeks, 1 month, 3 months, and 6 months (Figure 5). As for PM, there was no statistically significant difference in PM abundance within 6 h, except for one sample from GSE1911 that showed a decrease in PM enrichment score at 3 h postinjury. The enrichment scores of PM increased at 8 h, 24 h, 48 h, 72 h, 1 week, 2 weeks, 3 weeks, 3 months, and 6 months (Figure 6). At 24 h, 72 h, 2 weeks, 1 month, 3 months, and 6 months postinjury, the enrichment scores of neurons were lower in the injury-treated samples (Figure 7).
However, it is worth noting that some results from different datasets are inconsistent. For example, GSE2392 showed increased microglia enrichment 4 h after injury, while GSE31357 showed a decrease. Furthermore, several results from some datasets were statistically significant, while results from another dataset at the same postinjury time point did not reach significance.
In order to intuitively understand the changes, we divided the sample sets into five clusters according to the setting of the postinjury time of the samples in the database: 30 min to 12 h representing the hyperacute phase of TBI, 1 day to 3 days representing the acute phase, 1 week to 3 weeks representing the subacute phase, 1 month to 3 months representing the chronic phase, and longer than 3 months. Results from the same cluster were aggregated. The results showed that the abundance of microglia decreased from 30 min to 12 h after TBI and increased from 1 week to 3 months after injury. There were no differences in enrichment scores after 3 months (Figure 8A). For PM, there was no difference in cellular abundance up to 6 h postinjury, whereas the injured group had higher enrichment scores from 8 h to 3 months postinjury. Similar to microglia, enrichment scores did not differ between the injured and sham groups after 3 months (Figure 8B). The neuron changed at the same time as the PM, but in a different direction (Figure 8C).

4. Discussion

TBI is a worldwide public health concern without major therapeutic breakthroughs in the past decades [2]. Understanding the mechanisms that underlie TBI is critical for developing effective treatment options and improving the prognosis of this condition. The last decade has witnessed the rapid development of high-throughput transcriptome analysis. Public databases containing extensive TBI gene expression data enable comprehensive analysis of the specific effects of TBI on various cell functions and pathways.
To the best of our knowledge, this is the first TBI comprehensive study that utilizes WGCNA. WGCNA adopts a hierarchical clustering tree to classify all genes into several gene sets, namely, modules, according to the degree of coexpression. The correlation between modules and sample traits was estimated so that modules that are highly correlated with sample traits could be identified and gene functions of related modules could be further studied. WGCNA significantly reduces errors caused by multiple testing problems inherent in microarray data, while maximizing the use of all data, as it uses all gene expression data from samples—instead of focusing only on differentially expressed genes—to construct the scale-free weighted network [35]. Furthermore, the scale-free weighted network has a high degree of robustness, which means that the errors in individual genes will not affect the overall results.
In this study, we obtained key modules significantly associated with the traits of rat brain trauma samples. The sampling side reflects the differences of the local lesion tissue relative to unaffected sites in the same experimental animal, and the severity of injury reflects differences between experimental animals treated according to various severity levels (naive, sham injury, mild TBI, or severe TBI). The key modules related to the severity of injury and the sampling side partially overlap, which is not surprising since they may all be involved in pathophysiological processes directly related to TBI. Based on functional enrichment, both severity- and side-related genes were enriched in the inflammatory response and immune-related processes, indicating that inflammation and immune regulation are important post-TBI processes. This finding is supported by previous studies. Studies have found that TBI induced an inflammatory response in the central nervous system, which may cause acute secondary injury [7,8,36,37]. Evidence showed that the inflammatory response following a TBI does not only affect the focal zone but also disseminates to remote brain areas [38]. Furthermore, neuroinflammation after TBI was found to link to neurodegeneration [7]. The key modules related to postinjury time reflect the change between the hyperacute phase (4 h) and the acute phase (24 h). Ranked first in the KEGG analysis was the entry “pathways of neurodegeneration—multiple diseases”. Just as mentioned above, TBI has been proved a risk factor for neurodegenerative disorders, including Alzheimer’s disease and Parkinson’s disease [39]. Although the development of neurodegenerative disease is a long-term process, axonal damage and disruption of transport during the acute phase of TBI have been found to alter the molecular mechanisms of pathological protein formation, such as α-synuclein, amyloid-beta peptide, and hyperphosphorylated tau [40,41].
It is well-known that microglia and macrophages are important players in inflammatory and immune-related responses [42,43,44]. In the enrichment analysis of trait-related cell types, we found that microglia and PM were associated with the side of injury. The abundance of both microglia and PM in TBI was further validated using the other rat TBI datasets. This is consistent with previous studies, as substantial evidence has suggested that the alteration of microglia was involved in the acute immune and neuroinflammatory response following TBI [45]. For example, microglia have been shown to play a neuroprotective role in TBI [46,47]. There is also evidence that, as well as clearing debris from the brain, microglia may also be involved in maintaining the integrity of the glial barrier after brain injury [45].
We also found for the first time that PM were significantly associated with TBI as microglia. However, the two types of cells differed in their association with postinjury time, as we revealed that the activation of PM and microglia was different within 3 days after injury. PM are specialized macrophages residing in the perivascular space of the brain. Similar to microglia, PM migrated from the yolk sac into the brain during embryonic development. As a result, they are likely to be a self-renewing cell population that is not replenished by circulating monocytes under a normal state [48]. Although PM have been implicated in various diseases, for example, the clearance of amyloid-beta in animal models of amyloid-beta pathologies, it is unclear what role they play in TBI [48,49,50,51].
However, care needs to be taken when interpreting these results, since no single marker has so far been able to reliably define PM with both good sensitivity and specificity [51] and monocytes are capable of crossing the blood–brain barrier into the injured brain [52]. In order to identify cell markers, we analyzed single-cell-nucleus sequencing data from normal rat brain tissues. In spite of using a relatively high resolution for the analysis, which yielded 22 cell clusters, only one cluster was annotated as macrophage, making it difficult to distinguish between monocyte-derived macrophages and PM. Contrary to the few studies on PM, a number of studies have investigated the role of peripherally derived macrophages in TBI. The number of monocyte-derived macrophages from the blood reaches a peak in the damaged brain of animals and humans 24 to 48 h after injury [53]. Studies have identified monocyte-derived macrophages as a pathogenic factor during the chronic phase of TBI [53,54]. However, macrophages from the peripheral circulation and macrophages residing in the brain are two distinct types of cells and, therefore, may have different immune reactions to TBI [45,55]. To better explore the cell types associated with TBI and their roles in TBI, more research needs to be conducted to establish definitive profiles of microglia and other CNS macrophages at different stages of TBI.
TBI can be divided into four phases: hyperacute (minutes to hours), acute (hours to several days), subacute (several days to weeks), and chronic (months and beyond) [56]. Our results showed that during the hyperacute phase (30 min to 12 h postinjury), the cellular activity of microglia in the injury group was significantly lower than in the sham group, while in the acute phase (1 day to 3 days postinjury) the difference was not significant. For PM, there was no difference between the two groups in the hyperacute phase. The cellular abundance was significantly increased starting from the acute phase. According to the differences in these two cellular activation patterns, it is possible that each cell type responds to a specific pattern of stimulation at each time point after injury, and they likely play different roles in TBI. Further studies are needed to determine the specific mode and mechanisms of microglia and PM activation. Our results also showed that both of these cell types were significantly more abundant than controls up to 3 months after injury, suggesting that they may play a long-term role in TBI. A deeper investigation of these two types of cells, especially the PM, could better facilitate the development of inflammation-targeted therapies to improve TBI prognosis. This study can provide a reference for the setting of postinjury time for subsequent studies on these two types of cells.
Despite being the first study to perform WGCNA analysis on TBI samples, and identifying PM as one of the relevant cell types associated with the sampling side and postinjury time, it has several limitations. The GSVA Enrichment scores of samples from various datasets were aggregated to intuitively understand the changes in cell abundance during various TBI phases. However, TBI models, RNA microarray platforms, and brain regions from which samples were collected vary from dataset to dataset. As a result, care needs to be taken when interpreting the results. Although TBI in rat models may share some similarities with human TBI, there may be variations in the course of rat and human TBI. Furthermore, the situation of TBI in humans is much more complex than in animal models. Therefore, the results of this study should be interpreted with caution. Although this study predicted relevant pathways and cell types involved in TBI and explored the activity of these cells during various TBI phases, further in vitro and in vivo experiments are needed to validate these findings.
In summary, this study performed WGCNA on TBI samples of various severity, ipsilateral and contralateral to injury, sampled from different brain regions at different times postinjury; identified key coexpression modules associated with traits of interest; unveiled several biological processes, pathways, and cell types associated with TBI; and explored the activity of these cell types at various TBI phases. The results of the current study may provide a reference for further mechanism research and treatment of TBI.

5. Conclusions

Using WGCNA, our study revealed response to cytokine, inflammatory response, bacteria-associated response, neuroactive ligand–receptor interaction, metabolic and biosynthetic processes, and multiple pathways of neurodegeneration to be involved in the pathogenesis of TBI. Microglia, PM, and neurons were recognized to associate with TBI with different activities that vary by phase of TBI.

Author Contributions

Conceptualization, Q.T. and G.L.; methodology, Q.T. and M.S.; formal analysis, Q.T. and M.S.; writing—original draft preparation, Q.T.; writing—review and editing, R.Z., X.H., L.D., H.X., W.L. and G.L.; supervision, G.L.; project administration, G.L.; funding acquisition, G.L. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by grants from the National Natural Science Foundation of China (Nos. 81874083; 82072776; 82072775; 81702468; 81802966; 81902540; 81874082; 81472353), Natural Science Foundation of Shandong Province of China (Nos. ZR2019BH057; ZR2020QH174; ZR2021LSW025), the Jinan Science and Technology Bureau of Shandong Province (2021GXRC029), Key clinical Research project of Clinical Research Center of Shandong University (2020SDUCRCA011) and Taishan Pandeng Scholar Program of Shandong Province (No. tspd20210322).

Data Availability Statement

Datasets analyzed in the study can be accessed from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/ accessed on 25 October 2021) with the corresponding accession numbers listed in Materials and Methods.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Andriessen, T.M.; Jacobs, B.; Vos, P.E. Clinical characteristics and pathophysiological mechanisms of focal and diffuse traumatic brain injury. J. Cell Mol. Med. 2010, 14, 2381–2392. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  2. Bains, M.; Hall, E.D. Antioxidant therapies in traumatic brain and spinal cord injury. Biochim. Biophys. Acta 2012, 1822, 675–684. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Fann, J.R.; Ribe, A.R.; Pedersen, H.S.; Fenger-Gron, M.; Christensen, J.; Benros, M.E.; Vestergaard, M. Long-term risk of dementia among people with traumatic brain injury in denmark: A population-based observational cohort study. Lancet Psychiatry 2018, 5, 424–431. [Google Scholar] [CrossRef]
  4. Jassam, Y.N.; Izzy, S.; Whalen, M.; McGavern, D.B.; El Khoury, J. Neuroimmunology of traumatic brain injury: Time for a paradigm shift. Neuron 2017, 95, 1246–1265. [Google Scholar] [CrossRef] [Green Version]
  5. Maas, A.I.R.; Menon, D.K.; Adelson, P.D.; Andelic, N.; Bell, M.J.; Belli, A.; Bragge, P.; Brazinova, A.; Buki, A.; Chesnut, R.M.; et al. Traumatic brain injury: Integrated approaches to improve prevention, clinical care, and research. Lancet Neurol. 2017, 16, 987–1048. [Google Scholar] [CrossRef] [Green Version]
  6. Simon, D.W.; McGeachy, M.J.; Bayir, H.; Clark, R.S.; Loane, D.J.; Kochanek, P.M. The far-reaching scope of neuroinflammation after traumatic brain injury. Nat. Rev. Neurol. 2017, 13, 171–191. [Google Scholar] [CrossRef] [Green Version]
  7. The Lancet Neurology. The future research path of traumatic brain injury. Lancet Neurol. 2022, 21, 295. [Google Scholar] [CrossRef]
  8. Xiong, Y.; Mahmood, A.; Chopp, M. Animal models of traumatic brain injury. Nat. Rev. Neurosci. 2013, 14, 128–142. [Google Scholar] [CrossRef] [Green Version]
  9. Zhang, B.; Horvath, S. A general framework for weighted gene co-expression network analysis. Stat. Appl. Genet. Mol. Biol. 2005, 4, 17. [Google Scholar] [CrossRef]
  10. Ma, S.; Sun, S.; Geng, L.; Song, M.; Wang, W.; Ye, Y.; Ji, Q.; Zou, Z.; Wang, S.; He, X.; et al. Caloric restriction reprograms the single-cell transcriptional landscape of rattus norvegicus aging. Cell 2020, 180, 984–1001.e22. [Google Scholar] [CrossRef]
  11. Babikian, T.; Prins, M.L.; Cai, Y.; Barkhoudarian, G.; Hartonian, I.; Hovda, D.A.; Giza, C.C. Molecular and physiological responses to juvenile traumatic brain injury: Focus on growth and metabolism. Dev. Neurosci. 2010, 32, 431–441. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  12. Natale, J.E.; Ahmed, F.; Cernak, I.; Stoica, B.; Faden, A.I. Gene expression profile changes are commonly modulated across models and species after traumatic brain injury. J. Neurotrauma 2003, 20, 907–927. [Google Scholar] [CrossRef] [PubMed]
  13. White, T.E.; Ford, G.D.; Surles-Zeigler, M.C.; Gates, A.S.; Laplaca, M.C.; Ford, B.D. Gene expression patterns following unilateral traumatic brain injury reveals a local pro-inflammatory and remote anti-inflammatory response. BMC Genom. 2013, 14, 282. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  14. Matzilevich, D.A.; Rall, J.M.; Moore, A.N.; Grill, R.J.; Dash, P.K. High-density microarray analysis of hippocampal gene expression following experimental brain injury. J. Neurosci. Res. 2002, 67, 646–663. [Google Scholar] [CrossRef] [PubMed]
  15. Shojo, H.; Kaneko, Y.; Mabuchi, T.; Kibayashi, K.; Adachi, N.; Borlongan, C.V. Genetic and histologic evidence implicates role of inflammation in traumatic brain injury-induced apoptosis in the rat cerebral cortex following moderate fluid percussion injury. Neuroscience 2010, 171, 1273–1282. [Google Scholar] [CrossRef]
  16. Shojo, H.; Borlongan, C.V.; Mabuchi, T. Genetic and histological alterations reveal key role of prostaglandin synthase and cyclooxygenase 1 and 2 in traumatic brain injury-induced neuroinflammation in the cerebral cortex of rats exposed to moderate fluid percussion injury. Cell Transpl. 2017, 26, 1301–1313. [Google Scholar] [CrossRef]
  17. Hellmich, H.L.; Rojo, D.R.; Micci, M.A.; Sell, S.L.; Boone, D.R.; Crookshanks, J.M.; DeWitt, D.S.; Masel, B.E.; Prough, D.S. Pathway analysis reveals common pro-survival mechanisms of metyrapone and carbenoxolone after traumatic brain injury. PLoS ONE 2013, 8, e53230. [Google Scholar] [CrossRef]
  18. Sell, S.L.; Boone, D.R.; Weisz, H.A.; Cardenas, C.; Willey, H.E.; Bolding, I.J.; Micci, M.A.; Falduto, M.T.; Torres, K.E.O.; DeWitt, D.S.; et al. Microrna profiling identifies a novel compound with antidepressant properties. PLoS ONE 2019, 14, e0221163. [Google Scholar] [CrossRef]
  19. Meng, Q.; Zhuang, Y.; Ying, Z.; Agrawal, R.; Yang, X.; Gomez-Pinilla, F. Traumatic brain injury induces genome-wide transcriptomic, methylomic, and network perturbations in brain and blood predicting neurological disorders. eBioMedicine 2017, 16, 184–194. [Google Scholar] [CrossRef] [Green Version]
  20. Paban, V.; Ogier, M.; Chambon, C.; Fernandez, N.; Davidsson, J.; Risling, M.; Alescio-Lautier, B. Molecular gene expression following blunt and rotational models of traumatic brain injury parallel injuries associated with stroke and depression. J. Transl. Sci. 2016, 2, 330–339. [Google Scholar] [CrossRef]
  21. Lipponen, A.; Paananen, J.; Puhakka, N.; Pitkanen, A. Analysis of post-traumatic brain injury gene expression signature reveals tubulins, nfe2l2, nfkb, cd44, and s100a4 as treatment targets. Sci. Rep. 2016, 6, 31570. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  22. Puhakka, N.; Bot, A.M.; Vuokila, N.; Debski, K.J.; Lukasiuk, K.; Pitkanen, A. Chronically dysregulated notch1 interactome in the dentate gyrus after traumatic brain injury. PLoS ONE 2017, 12, e0172521. [Google Scholar] [CrossRef] [PubMed]
  23. Boone, D.R.; Weisz, H.A.; Willey, H.E.; Torres, K.E.O.; Falduto, M.T.; Sinha, M.; Spratt, H.; Bolding, I.J.; Johnson, K.M.; Parsley, M.A.; et al. Traumatic brain injury induces long-lasting changes in immune and regenerative signaling. PLoS ONE 2019, 14, e0214741. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  24. Ritchie, M.E.; Phipson, B.; Wu, D.; Hu, Y.; Law, C.W.; Shi, W.; Smyth, G.K. Limma powers differential expression analyses for rna-sequencing and microarray studies. Nucleic Acids Res. 2015, 43, e47. [Google Scholar] [CrossRef]
  25. Langfelder, P.; Horvath, S. Wgcna: An r package for weighted correlation network analysis. BMC Bioinform. 2008, 9, 559. [Google Scholar] [CrossRef] [Green Version]
  26. Hanzelmann, S.; Castelo, R.; Guinney, J. Gsva: Gene set variation analysis for microarray and rna-seq data. BMC Bioinform. 2013, 14, 7. [Google Scholar] [CrossRef] [Green Version]
  27. Diaz-Mejia, J.J.; Meng, E.C.; Pico, A.R.; MacParland, S.A.; Ketela, T.; Pugh, T.J.; Bader, G.D.; Morris, J.H. Evaluation of methods to assign cell type labels to cell clusters from single-cell rna-sequencing data. F1000Research 2019, 8, 296. [Google Scholar] [CrossRef]
  28. Yu, G.; Wang, L.G.; Han, Y.; He, Q.Y. Clusterprofiler: An r package for comparing biological themes among gene clusters. OMICS 2012, 16, 284–287. [Google Scholar] [CrossRef]
  29. Hao, Y.; Hao, S.; Andersen-Nissen, E.; Mauck, W.M., 3rd; Zheng, S.; Butler, A.; Lee, M.J.; Wilk, A.J.; Darby, C.; Zager, M.; et al. Integrated analysis of multimodal single-cell data. Cell 2021, 184, 3573–3587.e29. [Google Scholar] [CrossRef]
  30. Stuart, T.; Butler, A.; Hoffman, P.; Hafemeister, C.; Papalexi, E.; Mauck, W.M., 3rd; Hao, Y.; Stoeckius, M.; Smibert, P.; Satija, R. Comprehensive integration of single-cell data. Cell 2019, 177, 1888–1902.e21. [Google Scholar] [CrossRef]
  31. Butler, A.; Hoffman, P.; Smibert, P.; Papalexi, E.; Satija, R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat. Biotechnol. 2018, 36, 411–420. [Google Scholar] [CrossRef] [PubMed]
  32. Satija, R.; Farrell, J.A.; Gennert, D.; Schier, A.F.; Regev, A. Spatial reconstruction of single-cell gene expression data. Nat. Biotechnol. 2015, 33, 495–502. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  33. Ashburner, M.; Ball, C.A.; Blake, J.A.; Botstein, D.; Butler, H.; Cherry, J.M.; Davis, A.P.; Dolinski, K.; Dwight, S.S.; Eppig, J.T.; et al. Gene ontology: Tool for the unification of biology. The gene ontology consortium. Nat. Genet. 2000, 25, 25–29. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. Kanehisa, M.; Furumichi, M.; Tanabe, M.; Sato, Y.; Morishima, K. Kegg: New perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. 2017, 45, D353–D361. [Google Scholar] [CrossRef] [Green Version]
  35. Kan, K.J.; Guo, F.; Zhu, L.; Pallavi, P.; Sigl, M.; Keese, M. Weighted gene co-expression network analysis reveals key genes and potential drugs in abdominal aortic aneurysm. Biomedicines 2021, 9, 546. [Google Scholar] [CrossRef]
  36. Russo, M.V.; McGavern, D.B. Inflammatory neuroprotection following traumatic brain injury. Science 2016, 353, 783–785. [Google Scholar] [CrossRef] [Green Version]
  37. Mallah, K.; Couch, C.; Alshareef, M.; Borucki, D.; Yang, X.; Alawieh, A.; Tomlinson, S. Complement mediates neuroinflammation and cognitive decline at extended chronic time points after traumatic brain injury. Acta Neuropathol. Commun. 2021, 9, 72. [Google Scholar] [CrossRef]
  38. Shi, K.; Zhang, J.; Dong, J.F.; Shi, F.D. Dissemination of brain inflammation in traumatic brain injury. Cell Mol. Immunol. 2019, 16, 523–530. [Google Scholar] [CrossRef]
  39. Brett, B.L.; Gardner, R.C.; Godbout, J.; Dams-O’Connor, K.; Keene, C.D. Traumatic brain injury and risk of neurodegenerative disorder. Biol. Psychiatry 2022, 91, 498–507. [Google Scholar] [CrossRef]
  40. Uryu, K.; Chen, X.H.; Martinez, D.; Browne, K.D.; Johnson, V.E.; Graham, D.I.; Lee, V.M.; Trojanowski, J.Q.; Smith, D.H. Multiple proteins implicated in neurodegenerative diseases accumulate in axons after brain trauma in humans. Exp. Neurol. 2007, 208, 185–192. [Google Scholar] [CrossRef] [Green Version]
  41. Schaffert, J.; LoBue, C.; White, C.L.; Chiang, H.S.; Didehbani, N.; Lacritz, L.; Rossetti, H.; Dieppa, M.; Hart, J.; Cullum, C.M. Traumatic brain injury history is associated with an earlier age of dementia onset in autopsy-confirmed alzheimer’s disease. Neuropsychology 2018, 32, 410–416. [Google Scholar] [CrossRef] [PubMed]
  42. Perry, V.H.; Teeling, J. Microglia and macrophages of the central nervous system: The contribution of microglia priming and systemic inflammation to chronic neurodegeneration. Semin. Immunopathol. 2013, 35, 601–612. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  43. Bergink, V.; Gibney, S.M.; Drexhage, H.A. Autoimmunity, inflammation, and psychosis: A search for peripheral markers. Biol. Psychiatry 2014, 75, 324–331. [Google Scholar] [CrossRef] [PubMed]
  44. Borst, K.; Dumas, A.A.; Prinz, M. Microglia: Immune and non-immune functions. Immunity 2021, 54, 2194–2208. [Google Scholar] [CrossRef] [PubMed]
  45. Corps, K.N.; Roth, T.L.; McGavern, D.B. Inflammation and neuroprotection in traumatic brain injury. JAMA Neurol. 2015, 72, 355–362. [Google Scholar] [CrossRef] [Green Version]
  46. Roth, T.L.; Nayak, D.; Atanasijevic, T.; Koretsky, A.P.; Latour, L.L.; McGavern, D.B. Transcranial amelioration of inflammation and cell death after brain injury. Nature 2014, 505, 223–228. [Google Scholar] [CrossRef] [Green Version]
  47. Willis, E.F.; MacDonald, K.P.A.; Nguyen, Q.H.; Garrido, A.L.; Gillespie, E.R.; Harley, S.B.R.; Bartlett, P.F.; Schroder, W.A.; Yates, A.G.; Anthony, D.C.; et al. Repopulating microglia promote brain repair in an il-6-dependent manner. Cell 2020, 180, 833–846.e816. [Google Scholar] [CrossRef]
  48. Faraco, G.; Park, L.; Anrather, J.; Iadecola, C. Brain perivascular macrophages: Characterization and functional roles in health and disease. J. Mol. Med. 2017, 95, 1143–1152. [Google Scholar] [CrossRef] [Green Version]
  49. Hawkes, C.A.; McLaurin, J. Selective targeting of perivascular macrophages for clearance of beta-amyloid in cerebral amyloid angiopathy. Proc. Natl. Acad. Sci. USA 2009, 106, 1261–1266. [Google Scholar] [CrossRef] [Green Version]
  50. Thanopoulou, K.; Fragkouli, A.; Stylianopoulou, F.; Georgopoulos, S. Scavenger receptor class b type i (sr-bi) regulates perivascular macrophages and modifies amyloid pathology in an alzheimer mouse model. Proc. Natl. Acad. Sci. USA 2010, 107, 20816–20821. [Google Scholar] [CrossRef] [Green Version]
  51. Yang, T.; Guo, R.; Zhang, F. Brain perivascular macrophages: Recent advances and implications in health and diseases. CNS Neurosci. Ther. 2019, 25, 1318–1328. [Google Scholar] [CrossRef] [PubMed]
  52. Szmydynger-Chodobska, J.; Strazielle, N.; Gandy, J.R.; Keefe, T.H.; Zink, B.J.; Ghersi-Egea, J.F.; Chodobski, A. Posttraumatic invasion of monocytes across the blood-cerebrospinal fluid barrier. J. Cereb. Blood Flow Metab. 2012, 32, 93–104. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  53. Semple, B.D.; Bye, N.; Rancan, M.; Ziebell, J.M.; Morganti-Kossmann, M.C. Role of ccl2 (mcp-1) in traumatic brain injury (tbi): Evidence from severe tbi patients and ccl2-/- mice. J. Cereb. Blood Flow Metab. 2010, 30, 769–782. [Google Scholar] [CrossRef] [PubMed]
  54. Hsieh, C.L.; Niemi, E.C.; Wang, S.H.; Lee, C.C.; Bingham, D.; Zhang, J.; Cozen, M.L.; Charo, I.; Huang, E.J.; Liu, J.; et al. Ccr2 deficiency impairs macrophage infiltration and improves cognitive function after traumatic brain injury. J. Neurotrauma 2014, 31, 1677–1688. [Google Scholar] [CrossRef] [Green Version]
  55. Yamasaki, R.; Lu, H.; Butovsky, O.; Ohno, N.; Rietsch, A.M.; Cialic, R.; Wu, P.M.; Doykan, C.E.; Lin, J.; Cotleur, A.C.; et al. Differential roles of microglia and monocytes in the inflamed central nervous system. J. Exp. Med. 2014, 211, 1533–1549. [Google Scholar] [CrossRef] [Green Version]
  56. Hoffe, B.; Holahan, M.R. Hyperacute excitotoxic mechanisms and synaptic dysfunction involved in traumatic brain injury. Front. Mol. Neurosci. 2022, 15, 831825. [Google Scholar] [CrossRef]
Figure 1. WGCNA of the rat TBI Dataset GSE2871. (A) Hierarchical clustering tree of genes, with dissimilarity based on the topological overlap. Each gene cluster (module) is marked with a different color. (B) Identification of key modules related to sample traits. Heatmap illustrating the intramodular relationship. Rows and columns correspond to modules, and each cell contains the corresponding correlation and p-value. Pearson’s R-values are color-coded according to the color legend. The size of the rectangle is proportional to the p-value; the larger the rectangle, the more significant the correlation. The color of the line represents the correlation significance between traits and modules.
Figure 1. WGCNA of the rat TBI Dataset GSE2871. (A) Hierarchical clustering tree of genes, with dissimilarity based on the topological overlap. Each gene cluster (module) is marked with a different color. (B) Identification of key modules related to sample traits. Heatmap illustrating the intramodular relationship. Rows and columns correspond to modules, and each cell contains the corresponding correlation and p-value. Pearson’s R-values are color-coded according to the color legend. The size of the rectangle is proportional to the p-value; the larger the rectangle, the more significant the correlation. The color of the line represents the correlation significance between traits and modules.
Jcm 11 03437 g001
Figure 2. KEGG and GO enrichment analysis of the genes in sample-traits-related modules. (AC) Dot plots illustrating GO terms enriched in modules related to the severity of injury (A), sampling side (B), or postinjury time (C); (DF) dot plots showing KEGG pathways enriched in modules related to severity (D), side (E), and time (F). BP—biological processes, CC—cellular components, MF—molecular functions. Note: No term was enriched for CC in (B).
Figure 2. KEGG and GO enrichment analysis of the genes in sample-traits-related modules. (AC) Dot plots illustrating GO terms enriched in modules related to the severity of injury (A), sampling side (B), or postinjury time (C); (DF) dot plots showing KEGG pathways enriched in modules related to severity (D), side (E), and time (F). BP—biological processes, CC—cellular components, MF—molecular functions. Note: No term was enriched for CC in (B).
Jcm 11 03437 g002
Figure 3. Single-nucleus transcriptome analysis. (A) The t-SNE plot of rat brain single-nucleus transcriptomes; (B) heatmap illustrating the expression levels of top 20 marker genes (rows) sorted by likelihood ratio for each cell type in each cell (columns), except for oligodendrocytes, which have only 10 genes.
Figure 3. Single-nucleus transcriptome analysis. (A) The t-SNE plot of rat brain single-nucleus transcriptomes; (B) heatmap illustrating the expression levels of top 20 marker genes (rows) sorted by likelihood ratio for each cell type in each cell (columns), except for oligodendrocytes, which have only 10 genes.
Jcm 11 03437 g003
Figure 4. Gene set enrichment analysis (GSEA) of the genes in modules related to sample traits. (A,C) Dot plots illustrating cell types enriched in the genes in side-related modules using data from GSE45997 and the enrichment plots of PM; (B,D) dot plots illustrating cell types enriched in time-related modules using data from GSE2392 and the enrichment plots of neurons; (E,G,H) dot plots of the enrichments in side-related modules using data from GSE2871 and the enrichment plots of PM and microglia; (F) dot plots illustrating enriched cells in severity-related modules using data from GSE2871.
Figure 4. Gene set enrichment analysis (GSEA) of the genes in modules related to sample traits. (A,C) Dot plots illustrating cell types enriched in the genes in side-related modules using data from GSE45997 and the enrichment plots of PM; (B,D) dot plots illustrating cell types enriched in time-related modules using data from GSE2392 and the enrichment plots of neurons; (E,G,H) dot plots of the enrichments in side-related modules using data from GSE2871 and the enrichment plots of PM and microglia; (F) dot plots illustrating enriched cells in severity-related modules using data from GSE2871.
Jcm 11 03437 g004
Figure 5. Gene set variation analysis (GSVA) for microglia in the validation datasets. GSVA enrichment scores from the validation data sets are plotted as bar graphs. The GEO accession number of the dataset is labeled in each panel. Samples from multiple regions in one dataset are compared separately. Information is labeled in the corresponding panels if the dataset contains samples with different treatments or harvested at various postinjury time points. For GSE1911, N = 1 for each group; for GSE24047, N = 4 for Sham group and N = 3 for the other groups; for GSE31357, N = 4 for each group; for GSE59645, N = 2 for Sham group and N = 3 for Naïve group and Injury 24-h group; for GSE2392, N = 3 for Naïve, Injury 4-h, Injury 8-h, Injury 24-h, Injury 72-h, and Injury 21-d groups, N = 4 for Sham 30-m and Sham 4-h groups, N = 5 for Injury 30-m group, N = 2 for Sham 8-h, Sham 24-h, Sham 72-h, and Sham 21-d groups; for GSE64798, N = 5 for each group; for GSE67836, N = 4 for Sham_LFP group and N = 3 for the other groups; for GSE68207, N = 4 for each group; for GSE86579, N = 5 for Sham group and N = 6 for Injury 3-m group; for GSE111452_Hippocampus, N = 2 for Sham 24-h group, N = 3 for Injury 24-h group, and N = 4 for the other groups; for GSE80174_Cortex, GSE80174_Hippocampus, and GSE80174_Thalamus, N = 5 for each group; for GSE115614, N = 2 for Sham group and N = 3 for Injury 24-h group; for GSE111452_Cortex, N = 4 for each group. **** p < 0.0001, *** p < 0.001, ** p < 0.01, * p < 0.05. Data are shown as means ± SEM.
Figure 5. Gene set variation analysis (GSVA) for microglia in the validation datasets. GSVA enrichment scores from the validation data sets are plotted as bar graphs. The GEO accession number of the dataset is labeled in each panel. Samples from multiple regions in one dataset are compared separately. Information is labeled in the corresponding panels if the dataset contains samples with different treatments or harvested at various postinjury time points. For GSE1911, N = 1 for each group; for GSE24047, N = 4 for Sham group and N = 3 for the other groups; for GSE31357, N = 4 for each group; for GSE59645, N = 2 for Sham group and N = 3 for Naïve group and Injury 24-h group; for GSE2392, N = 3 for Naïve, Injury 4-h, Injury 8-h, Injury 24-h, Injury 72-h, and Injury 21-d groups, N = 4 for Sham 30-m and Sham 4-h groups, N = 5 for Injury 30-m group, N = 2 for Sham 8-h, Sham 24-h, Sham 72-h, and Sham 21-d groups; for GSE64798, N = 5 for each group; for GSE67836, N = 4 for Sham_LFP group and N = 3 for the other groups; for GSE68207, N = 4 for each group; for GSE86579, N = 5 for Sham group and N = 6 for Injury 3-m group; for GSE111452_Hippocampus, N = 2 for Sham 24-h group, N = 3 for Injury 24-h group, and N = 4 for the other groups; for GSE80174_Cortex, GSE80174_Hippocampus, and GSE80174_Thalamus, N = 5 for each group; for GSE115614, N = 2 for Sham group and N = 3 for Injury 24-h group; for GSE111452_Cortex, N = 4 for each group. **** p < 0.0001, *** p < 0.001, ** p < 0.01, * p < 0.05. Data are shown as means ± SEM.
Jcm 11 03437 g005
Figure 6. Gene set variation analysis (GSVA) for PM in the validation datasets. GSVA enrichment scores from the validation data sets are plotted as bar graphs. The GEO accession number of the dataset is labeled in each panel. Samples from multiple regions in one dataset are compared separately. Information is labeled in the corresponding panels if the dataset contains samples with different treatments or harvested at various postinjury time points. For GSE1911, N = 1 for each group; for GSE24047, N = 4 for Sham group and N = 3 for the other groups; for GSE31357, N = 4 for each group; for GSE59645, N = 2 for Sham group and N = 3 for Naïve group and Injury 24-h group; for GSE2392, N = 3 for Naïve, Injury 4-h, Injury 8-h, Injury 24-h, Injury 72-h, and Injury 21-d groups, N = 4 for Sham 30-m and Sham 4-h groups, N = 5 for Injury 30-m group, N = 2 for Sham 8-h, Sham 24-h, Sham 72-h, and Sham 21-d groups; for GSE64798, N = 5 for each group; for GSE67836, N = 4 for Sham_LFP group and N = 3 for the other groups; for GSE68207, N = 4 for each group; for GSE86579, N = 5 for Sham group and N = 6 for Injury 3-m group; for GSE111452_Hippocampus, N = 2 for Sham 24-h group, N = 3 for Injury 24-h group, and N = 4 for the other groups; For GSE80174_Cortex, GSE80174_Hippocampus, and GSE80174_Thalamus, N = 5 for each group; for GSE115614, N = 2 for Sham group and N = 3 for Injury 24-h group; for GSE111452_Cortex, N = 4 for each group. **** p < 0.0001, ** p < 0.01, * p < 0.05. Data are shown as means ± SEM.
Figure 6. Gene set variation analysis (GSVA) for PM in the validation datasets. GSVA enrichment scores from the validation data sets are plotted as bar graphs. The GEO accession number of the dataset is labeled in each panel. Samples from multiple regions in one dataset are compared separately. Information is labeled in the corresponding panels if the dataset contains samples with different treatments or harvested at various postinjury time points. For GSE1911, N = 1 for each group; for GSE24047, N = 4 for Sham group and N = 3 for the other groups; for GSE31357, N = 4 for each group; for GSE59645, N = 2 for Sham group and N = 3 for Naïve group and Injury 24-h group; for GSE2392, N = 3 for Naïve, Injury 4-h, Injury 8-h, Injury 24-h, Injury 72-h, and Injury 21-d groups, N = 4 for Sham 30-m and Sham 4-h groups, N = 5 for Injury 30-m group, N = 2 for Sham 8-h, Sham 24-h, Sham 72-h, and Sham 21-d groups; for GSE64798, N = 5 for each group; for GSE67836, N = 4 for Sham_LFP group and N = 3 for the other groups; for GSE68207, N = 4 for each group; for GSE86579, N = 5 for Sham group and N = 6 for Injury 3-m group; for GSE111452_Hippocampus, N = 2 for Sham 24-h group, N = 3 for Injury 24-h group, and N = 4 for the other groups; For GSE80174_Cortex, GSE80174_Hippocampus, and GSE80174_Thalamus, N = 5 for each group; for GSE115614, N = 2 for Sham group and N = 3 for Injury 24-h group; for GSE111452_Cortex, N = 4 for each group. **** p < 0.0001, ** p < 0.01, * p < 0.05. Data are shown as means ± SEM.
Jcm 11 03437 g006
Figure 7. Gene set variation analysis (GSVA) for neuron in the validation datasets. GSVA enrichment scores from the validation data sets are plotted as bar graphs. The GEO accession number of the dataset is labeled in each panel. Samples from multiple regions in one dataset are compared separately. Information is labeled in the corresponding panels if the dataset contains samples with different treatments or harvested at various postinjury time points. For GSE1911, N = 1 for each group; for GSE24047, N = 4 for Sham group and N = 3 for the other groups; for GSE31357, N = 4 for each group; for GSE59645, N = 2 for Sham group and N = 3 for Naïve group and Injury 24-h group; for GSE2392, N = 3 for Naïve, Injury 4-h, Injury 8-h, Injury 24-h, Injury 72-h, and Injury 21-d groups, N = 4 for Sham 30-m and Sham 4-h groups, N = 5 for Injury 30-m group, N = 2 for Sham 8-h, Sham 24-h, Sham 72-h, and Sham 21-d groups; for GSE64798, N = 5 for each group; for GSE67836, N = 4 for Sham_LFP group and N = 3 for the other groups; for GSE68207, N = 4 for each group; for GSE86579, N = 5 for Sham group and N = 6 for Injury 3-m group; for GSE111452_Hippocampus, N = 2 for Sham 24-h group, N = 3 for Injury 24-h group, and N = 4 for the other groups; for GSE80174_Cortex, GSE80174_Hippocampus, and GSE80174_Thalamus, N = 5 for each group; for GSE115614, N = 2 for Sham group and N = 3 for Injury 24-h group; for GSE111452_Cortex, N = 4 for each group. *** p < 0.001, ** p < 0.01, * p < 0.05. Data are shown as means ± SEM.
Figure 7. Gene set variation analysis (GSVA) for neuron in the validation datasets. GSVA enrichment scores from the validation data sets are plotted as bar graphs. The GEO accession number of the dataset is labeled in each panel. Samples from multiple regions in one dataset are compared separately. Information is labeled in the corresponding panels if the dataset contains samples with different treatments or harvested at various postinjury time points. For GSE1911, N = 1 for each group; for GSE24047, N = 4 for Sham group and N = 3 for the other groups; for GSE31357, N = 4 for each group; for GSE59645, N = 2 for Sham group and N = 3 for Naïve group and Injury 24-h group; for GSE2392, N = 3 for Naïve, Injury 4-h, Injury 8-h, Injury 24-h, Injury 72-h, and Injury 21-d groups, N = 4 for Sham 30-m and Sham 4-h groups, N = 5 for Injury 30-m group, N = 2 for Sham 8-h, Sham 24-h, Sham 72-h, and Sham 21-d groups; for GSE64798, N = 5 for each group; for GSE67836, N = 4 for Sham_LFP group and N = 3 for the other groups; for GSE68207, N = 4 for each group; for GSE86579, N = 5 for Sham group and N = 6 for Injury 3-m group; for GSE111452_Hippocampus, N = 2 for Sham 24-h group, N = 3 for Injury 24-h group, and N = 4 for the other groups; for GSE80174_Cortex, GSE80174_Hippocampus, and GSE80174_Thalamus, N = 5 for each group; for GSE115614, N = 2 for Sham group and N = 3 for Injury 24-h group; for GSE111452_Cortex, N = 4 for each group. *** p < 0.001, ** p < 0.01, * p < 0.05. Data are shown as means ± SEM.
Jcm 11 03437 g007
Figure 8. Analysis of cellular abundance of microglia, PM, and neurons during various courses of TBI. Violin plots showing aggregated GSVA enrichment scores for microglia (A), PM (B), and neurons (C) in hyperacute, acute, subacute, and chronic phases of TBI. For hyperacute phase, N = 27 for Sham group and N = 25 for Injury group; for acute phase, N = 23 for Sham group and N = 27 for Injury group; for subacute phase, N = 19 for Sham group and N = 20 for Injury group; for chronic phase, N = 35 for each group; for long-term p.i., N = 16 for each group. **** p < 0.0001, *** p < 0.001, ** p < 0.01, * p < 0.05. Data are shown as means ± SEM.
Figure 8. Analysis of cellular abundance of microglia, PM, and neurons during various courses of TBI. Violin plots showing aggregated GSVA enrichment scores for microglia (A), PM (B), and neurons (C) in hyperacute, acute, subacute, and chronic phases of TBI. For hyperacute phase, N = 27 for Sham group and N = 25 for Injury group; for acute phase, N = 23 for Sham group and N = 27 for Injury group; for subacute phase, N = 19 for Sham group and N = 20 for Injury group; for chronic phase, N = 35 for each group; for long-term p.i., N = 16 for each group. **** p < 0.0001, *** p < 0.001, ** p < 0.01, * p < 0.05. Data are shown as means ± SEM.
Jcm 11 03437 g008
Table 1. GSE databases included in the study.
Table 1. GSE databases included in the study.
Dataset IDTBITimeTissueSample Number Included
GSE1911 [14]CCI3 h, 24 hhippocampus3
GSE2392 [12]Moderate FPI30 min, 4 h, 8 h, 24 h, 3 d, 3 wperilesional cortex39
GSE2871 [11]Mild and severe FPI4 h, 24 hparietal cortex and hippocampus, ipsilateral and contralateral47
GSE24047 [15,16]FPI3 h, 6 h, 12 h, 48 hlateral cortex16
GSE31357 [17]TBI4 h, 24 hhippocampus16
GSE45997 [13]CCI24 hipsilateral and contralateral brain9
GSE59645 [18]TBI24 hhippocampus8
GSE64978 [19]FPI1 whippocampus10
GSE67836 [20]Rot-TBI and FPI1 mfrontal cortex13
GSE68207 [19]FPI1 wHippocampus8
GSE80174 [21]TBI3 mperilesional cortex,
dorsal hippocampus,
ipsilateral thalamus
30
GSE86579 [22]FPI3 mhippocampus11
GSE111452 [23]FPI24 h, 2 w, 3 m, 6 m, 1 yhippocampus, cortex113
GSE115614 [18]TBI24 hhippocampus5
GSE137869 [10]--brain2
Abbreviations: CCI—controlled cortical impact; FPI—fluid percussion injury; Rot-TBI—rotational acceleration induced TBI; min—minute; h—hour; d—day; w—week; m—month; y—year.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Tang, Q.; Song, M.; Zhao, R.; Han, X.; Deng, L.; Xue, H.; Li, W.; Li, G. Comprehensive RNA Expression Analysis Revealed Biological Functions of Key Gene Sets and Identified Disease-Associated Cell Types Involved in Rat Traumatic Brain Injury. J. Clin. Med. 2022, 11, 3437. https://doi.org/10.3390/jcm11123437

AMA Style

Tang Q, Song M, Zhao R, Han X, Deng L, Xue H, Li W, Li G. Comprehensive RNA Expression Analysis Revealed Biological Functions of Key Gene Sets and Identified Disease-Associated Cell Types Involved in Rat Traumatic Brain Injury. Journal of Clinical Medicine. 2022; 11(12):3437. https://doi.org/10.3390/jcm11123437

Chicago/Turabian Style

Tang, Qilin, Mengmeng Song, Rongrong Zhao, Xiao Han, Lin Deng, Hao Xue, Weiguo Li, and Gang Li. 2022. "Comprehensive RNA Expression Analysis Revealed Biological Functions of Key Gene Sets and Identified Disease-Associated Cell Types Involved in Rat Traumatic Brain Injury" Journal of Clinical Medicine 11, no. 12: 3437. https://doi.org/10.3390/jcm11123437

APA Style

Tang, Q., Song, M., Zhao, R., Han, X., Deng, L., Xue, H., Li, W., & Li, G. (2022). Comprehensive RNA Expression Analysis Revealed Biological Functions of Key Gene Sets and Identified Disease-Associated Cell Types Involved in Rat Traumatic Brain Injury. Journal of Clinical Medicine, 11(12), 3437. https://doi.org/10.3390/jcm11123437

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop