Abstract
Traumatic brain injury (TBI) is a major cause of morbidity and mortality, both in adult and pediatric populations. However, the dynamic changes of gene expression profiles following TBI have not been fully understood. In this study, we identified the differentially expressed genes (DEGs) following TBI. Remarkably, Serpina3n, Asf1b, Folr1, LOC100366216, Clec12a, Olr1, Timp1, Hspb1, Lcn2, and Spp1 were identified as the top 10 with the highest statistical significance. The weighted gene coexpression analysis (WGCNA) identified 12 functional modules from the DEGs, which showed specific expression patterns over time and were characterized by enrichment analysis. Specifically, the black and turquoise modules were mainly involved in energy metabolism and protein translation. The green yellow and yellow modules including Hmox1, Mif, Anxa2, Timp1, Gfap, Cd9, Gja1, Pdpn, and Gpx1 were related to response to wounding, indicating that expression of these genes such as Hmox1, Anxa2, and Timp1 could protect the brains from brain injury. The green yellow module highlighted genes involved in microglial cell activation such as Tyrobp, Cx3cr1, Grn, Trem2, C1qa, and Aif1, suggesting that these genes were responsible for the inflammatory response caused by TBI. The upregulation of these genes has been validated in an independent dataset. These results indicated that the key genes in microglia cell activation may serve as a promising therapeutic target for TBI. In summary, the present study provided a full view of the dynamic gene expression changes following TBI.
1. Introduction
Traumatic brain injury (TBI) is often caused by a sudden trauma to the head, but its deleterious effects on patients can be life-long and dynamic [1, 2]. It is a highlighted cause of mortality and long-term disability worldwide, and a well-recognized risk factor for neurological and psychological disorders, such as Alzheimer’s disease, chronic traumatic encephalopathy, depression, and psychosis [3, 4]. Such chronic consequences induced by TBI are now posing a huge burden for both TBI survivors and the society, yet establishing therapeutic intervention for the progression of tissue damage and neurodegeneration in chronic TBI remains a challenge, as our knowledge of long-term and ever-evolving pathogenic mechanisms behind post-TBI consequences is very limited.
It is believed that chronic traumatic brain inflammation plays a critical role in post-TBI neurodegeneration, and efforts are made to address metabolomic, proteomic, and genomic changes across a series of post-TBI intervals, in hopes of identifying potential long-term biomarkers for defining the injury progression [2, 5, 6]. A previous study has provided temporal profiles of 12 biomarkers in body fluids of TBI patients, demonstrated an association between the severity of TBI and the peak heights of each molecules, and suggested that release mechanisms may vary among different types of molecules, which further hints that successive measurements are essential for TBI diagnosis [7]. Utilizing microarray analysis, a recent study has described abnormal expression of immune mediators and brain injury-induced factors, along with regenerative immunoregulatory genes, in animal models of TBI, and remarked that such dysregulation of gene expression can be long-lasting after the initial injuries [8]. Also in experimental TBI, high levels of expression of Serpina3n, an astroglial activation marker, are found in neurons in the early stage of the injury [9]. Meanwhile, it is observed in another study that the knockout of Serpina3n led to increased neuronal apoptosis in neurons, and that an MMP2-specific inhibitor might serve as a therapeutic target for neurotrauma [10]. In addition, Trem2 and Tyrobp have been identified as major regulators and therapeutic targets in TBI [11, 12]. In the present study, to better capture the dynamics of gene expression and related biological functions at successive post-TBI time points, we conducted differential expression analysis and coexpression analysis to identify the coexpression modules, and characterized their functionalities by overrepresentation enrichment analysis, anticipating to reveal some critical genes involved in the dynamic changes post-TBI.
2. Materials and Methods
2.1. Gene Expression Data
The gene expression data of cortex tissues were downloaded from Gene Expression Omnibus (GEO) with accession GSE111452 [8]. A total of 60 samples were collected for the analyses in this study. The expression values were normalized to the 75th percentile intensity. Moreover, the gene expression data for validation was obtained from GEO with accession GSE79441 [13] and normalized to fragment per kilo-million (FPKM).
2.2. Differentially Expressed Genes
The differential expression analysis was conducted by pairwise comparison for the control, sham, and TBI samples at each time point. The FPKM-based gene expression data was logarithm-transformed. For each comparison, The R limma package was employed to identify the differentially expressed genes [14]. The values were adjusted to the false discovery rate (FDR) to avoid false positives by multiple statistical testing.
2.3. Coexpression Network and Functional Modules
The weighted gene coexpression network analysis (WGCNA) was used to identify the functional modules [15]. Specifically, the lowest soft-thresholding power for which the scale-free topology fit was selected. The topological overlap matrix (TOM) similarity was calculated using the power and expression data of differentially expressed genes.
2.4. Overrepresentation Enrichment Analysis (ORA)
The overrepresentation enrichment analysis was employed to identify the biological process of Gene Ontology (GO) enriched by the genes within the functional module. The hypergeometric test was used to calculate the value for each GO term, and was implemented in the R clusterProfiler package [16].
2.5. Random Forest Classifier
The expression levels of the hub genes for the 7 functional modules were used for the classifier training and validation. The training and validation sets were randomly divided and were used to build the classifier and evaluate the predictive performance of the classifier. The model construction and prediction were implemented by the R e1071 package (https://cran.r-project.org/web/packages/e1071/index.html).
2.6. Statistical Analyses
The two-sample and multisample comparisons were tested by the Student -test and analysis of variance (ANOVA), respectively. The Spearman correlations and corresponding values between the trait and hub genes were calculated to evaluate the module-trait relationship. The Kolmogorov-Smirnov test was used to test the enrichment of the six microglial cell activation-related genes in the upregulated genes by TBI.
3. Results
3.1. Identification of Genes Dysregulated by Traumatic Brain Injury
To identify the genes dysregulated by traumatic brain injury (TBI), we collected gene expression data of the rat cortex, which were obtained at 24 hours, 2 weeks, 3 months, 6 months, and 1 year after injury, and their corresponding negative controls (naïve group) with four replicates. Specifically, we identified 806, 628, 113, 114, and 94 differentially expressed genes (DEGs) at the five time points, respectively (Figure 1(a)). The number of DEGs were decreased over time. Moreover, the 24-hour and 2-week time points shared much more DEGs than the other time points (Figure 1(b); Fisher’s exact test, value < 0.001), suggesting that the shorter the time interval, the more similar the gene expression profiles. In contrast, the medium to late terms, from 3 months to 1 year, shared much fewer DEGs, suggesting that the gene expression profiles were dynamically changed after traumatic brain injury over time.
(a)
(b)
(c)
(d)
Furthermore, we identified the top 10 DEGs ranked by the values, including Serpina3n, Asf1b, Folr1, LOC100366216, Clec12a, Olr1, Timp1, Hspb1, Lcn2, and Spp1 (Figure 1(c)). Notably, Serpina3n achieved the highest statistical significance and was observed to be upregulated in TBI groups at both early and medium terms (Figures 1(c) and 1(d)). These results indicated that Serpina3n might act as a critical regulator in TBI.
3.2. Identification of Functional Modules Involved in TBI
As the genes could cooperate with each other and be coexpressed as a functional module, we then constructed a weighted gene coexpression network to identify the functional modules using weighted gene coexpression network analysis (WGCNA), which can be used for finding clusters (modules) of highly correlated genes. Specifically, we identified 12 functional modules from the DEGs by WGCNA (Supplementary Table S1, Figure 2(a)). Notably, the genes within each module showed highly similar expression patterns (Figure 2(b)). To further associate the functional modules with the characteristics of the samples, we conducted correlation analysis, and observed that the yellow, turquoise, green, and brown modules were highly correlated with TBI at both 24 hours and 2 weeks, while these modules were rarely correlated with TBI at medium to late terms (Figure 2(c); value < 0.05). The differential expression analysis of the eigen genes revealed that modules including yellow, green yellow, turquoise, black, brown, purple, and pink were highly dysregulated in TBI samples as compared with the controls (Figure 2(d); value < 0.05).
(a)
(b)
(c)
(d)
3.3. Functional Characterization of the WGCNA Modules
To further characterize biological function for the WGCNA modules, we conducted gene set enrichment analysis on the genes within the 7 TBI-related modules. Specifically, the black and turquoise modules were mainly involved in energy metabolism and protein translation (Figure 3(a)). The green yellow and yellow modules were related to response to wounding (Figure 3(a)). Moreover, the brown, pink, and purple modules were characterized by nervous system-related biological function (Figure 3(a)). Notably, the genes encoding ribosome proteins like RPL and RPS family genes were frequently identified in the turquoise module (Figure 3(b)). Hmox1, Mif, Anxa2, Timp1, Gfap, Cd9, Gja1, Pdpn, and Gpx1 were involved in the response to wounding, indicating that these genes might be the responses to TBI (Figure 3(b)). The green yellow module highlighted genes involved in microglial cell activation such as Tyrobp, Cx3cr1, Grn, Trem2, C1qa, and Aif1 (Figure 3(b)), suggesting that these genes were responsible for the inflammatory response caused by TBI.
(a)
(b)
3.4. Prediction of Disease Phenotypes by the Hub Genes within the Module
As shown in Figure 4(a), a total of 962 genes were included in the seven modules. Particularly, yellow and brown modules were upregulated and downregulated at the early term, respectively, but returned to normal at the medium to late term (Figure 4(b)). The turquoise module was upregulated at the early and medium terms, but returned to normal at the late term. In addition, the pink module was observed to be upregulated at 6 months, but downregulated at 1 year. These results indicated that the TBI-related modules were dynamically changed over time.
(a)
(b)
(c)
As the genes within the 7 modules exhibited significantly different expression patterns between TBI and controls, we then tested their separating ability for TBI and control samples. The samples were randomly divided into training () and testing () sets, and a classifier was built based on the hub genes in the training set. The classifier achieved a higher area under curve (AUC) at 0.895 about the TBI prediction in the validation set (). These results indicated that the hub genes were closely associated with TBI.
3.5. Microglial Cell Activation Is Responsible for Inflammatory Response at the Early and Medium Terms after TBI
As microglial cell activation was enriched by the genes within the green yellow module, we then tested the expression patterns of key genes involved in microglial cell activation and inflammatory response. Specifically, we observed that the key genes including Trem2, C1qa, Aif1, Grn, Cx3cr1, and Tyrobp were consistently upregulated at the early and medium terms (Figure 5(a); value < 0.05). To confirm this finding, we collected an independent gene expression data with 3 TBI and 3 control samples. Consistently, the six key genes were upregulated in TBI samples as compared with controls ( value < 0.05; Figure 5(b)). Moreover, these genes were highly clustered at the upregulated genes by TBI (Kolmogorov-Smirnov’s test, value < 0.05; Figure 5(c)). Furthermore, the six genes were also involved in leukocyte activation and inflammatory response (Figure 3(b)). These results indicated that microglial cell activation and the six key genes were responsible for inflammatory response at the early and medium terms after TBI.
(a)
(b)
(c)
4. Discussion
Traumatic brain injury (TBI) is a major cause of morbidity and mortality, both in adult and pediatric populations. However, the dynamic changes of gene expression profiles following TBI have not been fully understood. Specifically, we found that the number of DEGs were decreased over time, and the medium to late terms, from 3 months to 1 year, shared much fewer DEGs, suggesting that the gene expression profiles were dynamically changed after traumatic brain injury over time. Accordingly, the epigenetically dynamic changes have also been observed in brain tissues following TBI [17]. Among the DEGs, Serpina3n, Asf1b, Folr1, LOC100366216, Clec12a, Olr1, Timp1, Hspb1, Lcn2, and Spp1 were identified as the top 10 with the highest statistical significance. Remarkably, Serpina3n has been found as a promising target for neuroinflammation following TBI [18]. Timp1 has been considered as neuroprotective against traumatic and ischemic brain injury in mice [19]. Moreover, Lcn2 (Lipocalin-2) is upregulated by thrombin-induced brain injury through protease-activated receptor-1 activation [20]. These results suggested that these DEGs were highly associated with TBI.
The WGCNA identified 12 functional modules from the DEGs, which showed specific expression patterns over time and were characterized by enrichment analysis. Specifically, the black and turquoise modules were mainly involved in energy metabolism and protein translation. A previous study has shown that glucose metabolism has been altered following TBI [21]. The green yellow and yellow modules including Hmox1, Mif, Anxa2, Timp1, Gfap, Cd9, Gja1, Pdpn, and Gpx1 were related to response to wounding, indicating that expression of these genes such as Hmox1, Anxa2, and Timp1 could protect the brains from brain injury [19, 22, 23]. The green yellow module highlighted genes involved in microglial cell activation such as Tyrobp, Cx3cr1, Grn, Trem2, C1qa, and Aif1, suggesting that these genes were responsible for the inflammatory response caused by TBI. The upregulation of these genes has been validated in a testing dataset. Consistently, Tyrobp and Trem2 were identified as major hubs in human APOE-expressing mice following traumatic brain injury by gene coexpression network [12]. Moreover, Tyrobp and Trem2 may be a promising therapeutic target in TBI. The GRN protein is associated with increased lysosomal biogenesis in activated microglia and can exacerbate neuronal damage after traumatic brain injury [24]. These results indicated that the key genes in microglia cell activation may serve as a promising therapeutic target for TBI.
In addition, the present study also has some limitations such as the small sample size, an independent validation dataset for the dynamic changes of gene expression profiles, and lack of experimental validation. However, the present study still provides a full view of the dynamic gene expression changes following TBI.
Data Availability
The data used in this study are available at Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/gds), which have been cited at relevant places within the text as references.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Authors’ Contributions
F. L. and X.X. designed this study. Z. Z., D. W., and R. Z. conducted the data analysis, visualization, and experiments. Z. Z., D. W., X. X., F. L., and T. P. contributed to the writing of the paper and setting of figures. Zhi-jie Zhao, Dong-po Wei, and Rui-zhe Zheng contributed equally to this work.
Acknowledgments
This study was funded by the talent cultivation plan of Tongren Hospital, Shanghai Jiao Tong University School of Medicine: “TongRen Newstar” (No. 2019shtrxx10).
Supplementary Materials
Supplementary Table S1: the functional modules and their corresponding genes with coexpression. (Supplementary Materials)