1. Data sorting and difference analysis
Based on the 8 RNA sequencing datasets retrieved from the GEO database, all the datasets were processed into Counts format for differential analysis. The results showed that there were some differences in the total amount of expression among the different datasets and samples prior to standardization (Fig. 2A). However, after standardization, the differences between datasets were eliminated, and the total expression amount of each sample was essentially equal (Fig. 2B).
The DESeq2 package was applied for differential analysis of RNA sequencing data set, and 103 differentially expressed genes were obtained, including 40 down regulated genes and 63 up regulated genes after brain trauma (Fig. 2CD). The most significant differentially expressed genes between Sham group and TBI group included CCL13, TREM1, CXCL2, SELE, etc.
2 microarray datasets were analyzed by the limma package, and a total of 135 differentially expressed genes were obtained. The most significant differential genes included PPBP, CD77, LIF, CCL12, etc. (Fig. 2EF). Among them, 10 differential genes appear in both RNA sequencing gene sets and microarray gene sets, including: CCL13, CXCL2, MARCO, IL1R2, CXCR2, SELL, TREM1, C4A, C4B and S100A8. (Fig. 2G). TREM1, MARCO, IL1R2 and CXCR2 are activating receptors expressed on neutrophils, monocytes and macrophages. SELL, CXCL2 and CCL13 are adhesion molecules and chemokines of immune cell. C4A and C4B play a role in complement activation, promoting phagocytosis, preventing immune complex deposition and neutralizing viruses. S100A8 mainly mediates intracellular inflammatory signal transduction by activating Toll like receptors.
2. Enrichment analysis
The results of Metascape enrichment analysis showed that the pathways with the strongest correlation with the differential gene set obtained in the previous step are inflammatory response (GO: 0006954), Cytokine Signaling in Immune system (R-HSA-1280215), Cytokine receptor interaction (hsa04060) and positive regulation of cytokine production (GO: 0001819) (Fig. 3A). These results indicated that after brain trauma, immune response related pathways are fully activated, and cytokine mediated intercellular communication is crucial in the process of inflammatory.
The purpose of building a network of enriched terms was to show the interaction between each signal path more intuitively. The results showed that the most important nodes in the network were positive regulation of the response to external stimulus and positive regulation of leukocyte activation. The enrichment of terms related to NF-kappa B signaling pathway and neutral precipitation was also noteworthy (Fig. 3B). In the following picture, the nodes of the Network of enriched terms were colored with P-values. The results demonstrated that the relevance of most terms in the network was highly reliable (Fig. 3C).
In order to further analyze the interactions of differential genes at the protein level, we carried out protein-protein interaction enrichment analysis. The MCODE algorithm was performed to divide the nodes in the PPI network into 10 modules (sorted by importance). For example, Module 1 was Signaling by Interleukins. Module 2 was Cytokine-cytokine receptor interaction. Module 3 was MyD88 deficiency (TLR2/4). Module 4 was RHO GTPases Activate NADPH Oxidases.et.al. (Fig. 3DE). Therefore, we verified the interaction of immune response related factors at the protein level, and found that the imbalance of oxidative stress is also an important part of the pathological results of TBI.
3. Immune infiltration score
The Estimate method was applied to evaluate the overall level of immune infiltration. We found that the total score of TBI group was higher than Sham group in the three evaluation methods of Estimate, Immune signature and Strategic signature, which means a higher level of immune cell infiltration. (Fig. 4A).
By time grouping, the degree of immune infiltration reached a peak within one week after brain trauma, and gradually decreased within one month. With the brain injury entering the chronic phase of the 2-3 months, the immune infiltration gradually increased again (Fig. 4B). In terms of tissue type, the hippocampus has the highest immune score, followed by the cerebral cortex. The immune scores of thalamus and hypothalamus were low, and there was no significant difference between TBI and Sham group in these parts (Fig. 4C).
CIBERSORT analysis showed that compared with Sham group, the proportion of macrophages (M1, M2), Neutrophils and plasma cells in TBI group increased significantly, while the proportion of Tregs cells (Tregs) and follicular helper T cells (Tfh) decreased relatively (Fig. 4D).
The proportion of macrophages in cortex and hippocampus was significantly higher than that in thalamus and hypothalamus. The proportion of Treg cells was lower in cortex and hippocampus. However, mast cells in hypothalamus were significantly enriched (Fig. 4E).
4. Weighted gene co-expression network analysis
After sorting out, we found that the collected data did not have missing values or outliers. After building the automatic network and detection module, the system recommends a soft threshold of 15. As shown in the figure (Fig. 5A), when the soft threshold is 15, the downward trend of the curve tends to be gentle, indicating that the average connectivity is good. Using the blockwiseModules function, 228 differentially expressed genes were classified into 7 co-expressed modules and a hierarchical clustering tree diagram of module identification was drawn (Fig. 5B). The module-phenotype data association showed that the turquoise module has the highest correlation with the immune score (Fig. 5CD). This module included 29 co-expressed genes such as LYZ, TYROBP, CDK1, PDCD1 and RGS1. Subsequently, we constructed a GS-MM histogram for the co-expressed genes in turquoise group. The results demonstrated that genes in the group had good correlation with other genes in the group (module membership) and with immune score (gene significance) (Fig. 5EFG).
The gene set of immune related modules obtained from WGCNA analysis was imported into String database for PPI network analysis, and the results verified the relevance of these genes and their role in immune response. We screened the key node genes with more than 4 connections. As a result, seven key nodes were found, namely, CD86, CCR5, PTPRC, ITGB2, FCGR2B, TLR2, TYROBP (Fig. 6E).
5. Cell death
Most researchers believe that brain injury leads to death of nerve cells through variety ways. We selected apoptosis, necrosis, autophagy, ferroptosis and pyroptosis as the five most common cell death pathways for evaluation.
The results demonstrated that the total cell death index of TBI group was significantly higher than Sham group (Fig. 6A).
At the beginning of the establishment of model, the death rate of brain cells in mice gradually increased, the cell death index rose to the highest within one month, and gradually decreased in 2-3 months (Fig. 6B). Among them, the cell death index of cerebral cortex, hippocampus and thalamus was higher, while the cell death score of hypothalamuses was lower (Fig. 6C). Apoptosis and pyroptosis were the main brain cell death modes after brain trauma (Fig. 6D). However, the score of autophagy in TBI group was lower than that in Sham group. It can be seen that cell death is the most serious within one month after brain trauma, which is the key period for treatment. Apoptosis and pyroptosis are the main forms of cell death after brain trauma, from which we can find the key targets for treatment. Activating autophagy pathway will enhance the tolerance of nerve cells to the harsh environment of injury, hunger and oxidative stress, so it is considered as a possible effective strategy.
6. Cluster analysis
The ConsensusClusterPlus method was performed to construct a consistency matrix for each grouping number (K) value. The results showed that when K = 2, cluster assignments were stable. When K > 2, the new clusters were very small, and mixed color blocks appeared around them. This shows that the data should be divided into two groups to cover most of the characteristics of the data set, and the grouping is relatively stable (Fig. 7A). Empirical cumulative distribution (CDF) plots display consensus distributions for each K. When K is between 2 and 4, the data results are relatively stable (Fig. 7BC).
After grouping, the expression of each gene is displayed by heat map (Fig. 7D). We found that the expression of each immune key node gene in group 1 was significantly higher than that in group 2. The cell death scores of the two groups were compared. The results demonstrated that the ferroptosis, necrosis, and pyroptosis score in group 2 was significantly higher than that in group 1, while the autophagy score was lower than that in group 1 (Fig. 7E).
This showed that the immune related node genes in TBI were involved in the process of ferroptosis, necrosis and pyroptosis of the nervous system after brain trauma, and it were negatively correlated with autophagy.
7. Single-cell sequencing analysis
We integrated three single-cell sequencing data sets of TBI and carried out quality control. The Harmony method was applied to eliminate the batch effect of components (Fig. 8C). Finally, we obtained the sequencing data of 37621 high-quality cells. The tSNE method was applied for dimension reduction. According to the selected marker genes, the sequencing samples were classified into 12 cell types (Fig. 8B). The results showed that the proportion of microglia and macrophages in TBI group increased significantly (Fig. 8A), indicating that the inflammatory process of TBI could jointly participated by monocyte macrophage system and microglia.
Pseudotime analysis revealed that after brain trauma, microglia evolved through 3 nodes and then differentiated in two different directions (Fig. 8D). According to the nodes of differentiation, monocle2 divided microglia into 7 clusters (Fig. 8D). GSVA results revealed that the RESPONSE TO CGMP pathway was most active in cluster 1. The EOSINOPHIL DIFFERENTIATION, MATURE RIBOSOME ASSEMBLY, and NEGATIVE REGULATION OF INTERLEUKIN 6 MEDIATED SIGNALING PATHWAY pathways were generally active in clusters 2-5. Cluster 6 and cluster 7 represent two directions of microglia evolution, which differ in REGULATION OF MACROPHAGE FUSION, NEGATIVE REGULATION OF INTERLEUKIN 1 ALPHA PRODUCTION and NEGATIVE REGULATION OF INTERLEUKIN 6 MEDIATED SIGNALING PATHWAY (Fig. 8E). Gene expression analysis showed that the expression of immune cell surface activating receptors CCR5, FCGR2B, LTGB2, PTPRC and TLR2 increased with time after traumatic brain injury, and after a plateau period, the expression increased again. These results are consistent with the previously found trend of inflammatory response over time after traumatic brain injury. CD86 expression began to increase only in the mid to late post-traumatic brain injury period. However, TYROBP expression remained stable (Fig. 8F).
We analyzed the distribution analysis of immune related key node genes in TBI obtained previously. The results indicated that these genes were mainly distributed in microglia, followed by monocyte nuclear dendritic cells (Fig. 9ABC). This shows that after brain trauma, a variety of immune cells are activated, and the innate immune cells are the main participants in the nervous system inflammation after brain trauma.
GSVA enrichment analysis showed that MTOR, oxidative phosphorylation and other signal pathways were significantly up regulated in TBI group (Fig. 9D). These signal pathways are not only related to respiration, but can also regulate the nutrition and energy metabolism of cells under stress.
8. Cell-cell communication analysis
Cell-cell communication analysis mainly reflects paracrine and autocrine, which is of great significance for analyzing the changes in tissue microenvironment under pathological conditions. On the whole, the number of cell communication in TBI group was slightly higher than that in Sham group, while the strength of communication was significantly increased (Fig. 10A). In TBI group, CCL, CX3C, PTN, PSAP, MK and DEN signal pathways were dominant, while TGFB, VEGF and TNF pathways related intercellular communication was slightly weakened (Fig. 10B). From the perspective of ligand-receptor pairs, the interaction of CCL5-CCR5 was enhanced after brain trauma (Fig. 10C). This may be an important factor causing microglia and macrophages to accumulate in the damaged area.
From the perspective of the number and strength of cell communication, after brain trauma, the communication between nerve cells, astrocytes, microglia, parietal cells and surrounding cells is generally enhanced (Fig. 10D). Among them, we separately analyzed the microglia and macrophages that play the most important role in the innate immunity of the central nervous system. The results demonstrated that under TBI conditions, macrophages could receive paracrine signals from most cells, including damaged nerve cells, astrocytes, endothelial cells, and parietal cells (Fig. 10E). These signals may trigger macrophages to migrate to damaged cells, and eventually destroy and phagocytosis them. Microglia can also receive paracrine signals from most surrounding cells, and the strength of receiving nerve cell communication signals was stronger (Fig. 10E).
In the signal transmission intensity/quantity diagram, we found that the signal transmission strength of astrocytes was the highest (Fig. 10FG). The intercellular communication map showed that the main communication pathways of astrocytes in TBI group included PSAP, PTN and MK. PSAP is closely related to ferroptosis of nerve cells. PTN is associated with nerve cell regeneration and repair. MK pathway plays an important role in the process of inflammation and tissue cell repair. The main intercellular communication pathways of astrocytes in TBI group included PSAP, PTN and MK. PSAP is closely related to ferroptosis of nerve cells. PTN is associated with nerve cell regeneration and repair. MK pathway is crucial in the process of inflammation and tissue cell repair.
From the hierarchical map, the number of microglia receiving signals in TBI group increased significantly, which verified our previous results (Fig. 10HI). The oligodendrocyte precursor receives less intercellular communication signals, including TGFb, VEGF, and PDGF, which are related to cell growth and differentiation. This may inhibit oligodendrocyte maturation, leading to neurological repair disorders and demyelinating diseases.
CCL pathway plays an important role in the aggregation, killing and phagocytosis of innate immune cells. Therefore, we separately constructed a hierarchical map (Fig. 10JK). It was found that after brain trauma, not only the signal of CCL pathway between nervous system and immune cells was enhanced, but also the communication between microglia and macrophages was enhanced. This indicates that there is a synergistic effect between microglia and macrophages in the process of TBI induced inflammation.
To fully display the changes in intercellular communication before and after brain trauma, we constructed the ligand receptor map and the cell type pathway map of intercellular communication (Fig. 11ABCDEF). The results indicated that macrophages and microglia were more active in the communication of CCL immune and other immune pathways after brain trauma. This confirms our previous results.
The results of cell communication clustering showed that the signal pathways involved in the microglial cell related communication signal pattern (Pattern 1) included not only the chemokines CCL and CX3C, but also TNF and TGFb pathways. TNF pathway can induce programmed cell death in inflammatory state. In addition to participating in cell growth and differentiation, TGFb can also regulate cell apoptosis (Fig. 11GH). This means that microglia may play an immunomodulatory role and mediate cell death through these pathways.
Among the cell communication modes related to macrophages, besides CCL and TNF, GSA and GALECTIN signal pathways (Fig11. IJ) were also included. Among them, GSA pathway has a wide range of functions, which can reduce inflammation and apoptosis through NF-κB, TNF and MAPK. GALECTIN is a signal pathway recognized by the danger associated molecular patterns (DAMPs) in innate immunity, and is considered to be a key regulator of macrophage and microglial activity. It is also an attractive drug target in neurodegenerative diseases. In normal tissues, inflammation usually plays a role in fighting invading pathogens and maintaining tissue health. However, under pathological conditions such as trauma, inflammation can also aggravate tissue damage. The inflammation induced by TBI was once thought to occur through the damaged BBB by peripheral immune mediators.