Integration of single-cell and bulk RNA sequencing data reveals key cell types and regulators in traumatic brain injury

: Traumatic brain injury (TBI) is a leading cause of disability and mortality worldwide, whose symptoms ranging from mild to severe, even life-threatening. However, specific cell types and key regulators involved in traumatic brain injury have not been well elucidated. In this study, utilizing single-cell RNA-seq (scRNA-seq) data from mice with TBI, we have successfully identified and characterized 13 cell populations including astrocytes, oligodendrocyte, newly formed oligodendrocytes, microglia, two types of endothelial cells, five types of excitatory and two types of inhibitory neurons. Differential expression analysis and gene set enrichment analysis (GSEA) revealed the upregulation of microglia and endothelial markers, along with the downregulation of markers of excitatory neurons in TBI. The cell-cell communication analysis revealed that microglia and endothelial cell might interact through the interaction of Icam1-Il2rg and C1qa-Cd93, network in microglia might be a candidate therapeutic target in TBI. In contrast, the excitatory neurons were involved in maintaining normal brain function, and their inactivation might cause dysfunction of nervous system in TBI patients. In conclusion, the present study has discerned major cell types such as microglia, endothelial cells and excitatory neurons, and revealed key regulator such as TYROBP, C1QA, and CD93 in TBI, which shall improve our understanding of the pathogenesis of TBI.


Introduction
Traumatic brain injury (TBI) results from external forces such as a blow to the head or an object penetrating the skull. TBI is regarded as a common cause of death and disability and a major healthcare issue worldwide [1,2]. It often leads to various central nervous system (CNS) symptoms and greatly increases the risk of developing long-lasting psychiatric disorders. Since TBI can be highly heterogeneous, the prevention of post-TBI complications is one of the greatest challenges for researchers and therapeutists.
Neuroinflammation is considered as an important post-TBI secondary injury mechanism, with hallmarks such as neurovascular unit (NVU) alterations [3]. The NVU is mainly composed of neurons, astrocytes, microglia, endothelial cells and smooth muscle cells [4], and it has been demonstrated that the activation of microglia and astrocytes could promote neuroinflammation via secreting more proinflammatory cytokines, such as TNF-α and IL-1β [1,5,6]. Certain microglia receptors, namely CX3CR1, P2Y12R, TLR4, IL1R1, RAGE and TNFR1, have been highlighted, and their associated downstream pathways include the PI3K/AKT, MAPK, cAMP, and NF-κB signaling pathway, which activate the expression of a series of inflammatory genes [7]. Previous studies have also suggested an association between neuroinflammatory responses and excitatory-to-inhibitory synaptic imbalance in patients and animal models of TBI. Synaptic reorganization has been observed in mice after TBI, which is possibly caused by disturbed neurotransmitter release [8,9]. In addition, TBI could lead to endothelial dysfunction and cardiovascular diseases, where the nitric oxide (NO) and endothelial-dependent hyperpolarization (EDH) pathways played a critical role [5]. Meanwhile, vascular endothelial growth factor (VEGF) has been found to participate in neuroprotection after TBI [10].
The advances in single-cell RNA sequencing technology have enabled the exploration of complex and rare cell populations and helped uncover regulatory mechanisms between genes and cells involved in the pathogenic processes of several diseases [11][12][13][14]. However, a detailed landscape of molecular mechanisms and pathways involved in post-TBI pathology still remains absent. In the present study, we have integrated single-cell and bulk RNA-seq data to identify TBI-related cell types and regulators, hopefully providing some novel therapeutic targets for TBI.

Data acquisition
The single-cell RNA-seq data from cortex of mice was collected from Gene Expression Omnibus (GEO) with accession number GSE124952 [15]. The RNA-seq data of cortex tissues in mice with TBI and normal controls and the RNA-seq data of Dap12 (Tyrobp)-deficient mouse microglial cells were downloaded from GEO with accession numbers GSE79441 [16] and GSE9043, respectively. Data from GEO under accession numbers GSE68207 [17] and GSE155063 [18] were used for TYROBP validation. All these data were pre-normalized in GEO database, and were used for downstream analysis.

Cell clustering analysis
The cell clustering analysis was implemented in R Seurat [15] package. Specifically, the count-based expression data was first imported into R using CreateSeuraObject, and the top 2000 highly variable features were selected using FindVariableFeatures with vst method. We selected different clustering resolutions at 0.1, 0.3, 0.5 and 0.7, manually compared those clustering results, and observed that the clustering with a resolution at 0.3 achieved the optimal partitions. Thus, the clusters were identified at a resolution of 0.3 using FindClusters, and T-distributed Stochastic Neighbor Embedding (t-SNE) was applied to reduce the dimensionality. The cell-type specific genes were detected using FindAllMarkers function at adjusted P-value < 0.05 and minimal percentage > 0.25. Reported marker genes of those cell clusters were collected from an earlier study [16]. Visualization tools for the scRNA-seq data such as t-SNE, heat-map, and vioplot were also implemented in R Seurat package.

Gene set enrichment analysis
Two enrichment algorithms, gene set enrichment analysis (GSEA) and gene set overrepresentation enrichment analysis (ORA) were employed in this study. These two algorithms were applied to identify cell types involved in TBI. Also, GSEA was applied to identify pathways enriched in the gene sets of interest, which was implemented in R clusterProfiler [17] package.

Cell-cell communication
The ligand-receptor pairs were collected from a public dataset provided by an earlier study [18]. The cell-type specificity of differentially expressed genes encoding ligands or receptors in TBI were then tested in identified cell populations. Subsequently, these cell-type specific ligands or receptors involved in TBI were mapped to known ligand-receptor pairs, and those successfully mapped pairs were considered as links between corresponding cell types.

Differential expression analysis
Differential expression analysis of bulk RNA-seq was conducted using R DESeq2 package. Specifically, disease status (TBI or control) was regarded as a variable in the meta-data. Genes with an FPKM value above 1 in at least 20% of samples were retained for this analysis. The count-based gene expression was normalized using the median of ratios method. The thresholds of adjusted P-value and log2 (fold change) were selected at 0.05 and 0.5, respectively. Differential expression analysis of single-cell RNA-seq data based on Wilcoxon rank sum test was performed, which was implemented in R Seurat package with function FindAllMarkers.

Characterization of cell clusters in mice cortex by five previously reported cell types
To define cell types in mouse cortex, we collected a total of 10,538 single cells in cortex without any treatments from a previous study [16], and conducted T-distributed Stochastic Neighbor Embedding (t-SNE) analysis to identify cell clusters. Totally, we identified 13 cell clusters and characterized cell types by known markers [16] of astrocytes (Astro), oligodendrocyte, newly formed oligodendrocytes (NF Oligo), oligodendrocyte precursors (Oligo), microglia (Microglia), endothelial cells (Endo), excitatory and inhibitory neurons (Excitatory and Inhibitory) ( Figure 1A). Particularly, excitatory and inhibitory neurons and endothelial cells could be further divided into 5, 2 and 2 sub-clusters in this clustering ( Figure 1B), and these cell populations were denoted as Excitatory-1/2/3/4/5, Inhibitory-1/2 and Endo-1/2. Specifically, marker genes including Ctla2a, Ly6c1, Pglyrp1, Cldn5 and Itm2a were shared by these two types of endothelial cells, while Myl9, Higd1b, Cald1, Rgs5 and Vtn were specifically expressed in Endo-1. Notably, Myl9 and Cald1 were involved in vascular smooth muscle contraction. Comparison of the representative genes in the excitatory neurons revealed that gene expressions in Excitatory-1, 2 and 4 shared a high similarity, and Excitatory-3 and 5 had more unique expression patterns. Specific expressions of certain marker genes further confirmed that this characterization of cell types makes biological sense ( Figure 1C), including Slc17a7, Gja1, Aspa, Bmp4, C1qa, Flt1 and Gad2, which were detected respectively in excitatory neurons, astrocytes, oligodendrocyte, newly formed oligodendrocytes, microglia, endothelial cells, and inhibitory neurons.

Identification of cell types involved in traumatic brain injury
To identify cell types involved in traumatic brain injury (TBI), we conducted a systematic data analysis by integrating scRNA-seq and bulk RNA-seq data. The bulk RNA-seq data was collected from a previous study [19]. Specifically, differential expression analysis of mouse cortex RNA-seq data revealed 2977 upregulated and 1726 downregulated genes in TBI as compared with controls (TBI (n = 3) vs. control (n = 3), adjusted P-value < 0.05 and log2 (fold change) > 0.5). Gene set overrepresentation enrichment analysis (ORA) of the upregulated and downregulated genes against the cell-type specific marker genes by scRNA-seq data revealed that the upregulated genes were highly enriched in microglia and endothelial cells, while the downregulated genes were enriched in excitatory neurons (Figure 2A, FDR < 0.05). Consistently, gene set enrichment analysis (GSEA) further confirmed that these cell types were either highly activated or inactivated in TBI ( Figures  2B-2H, FDR < 0.05). These results indicated that the systematic analysis of scRNA-seq and bulk RNA-seq data could accurately detect various cell types involved in TBI.

Cell-cell communications between activated cell types in TBI
To predict cell-cell communications between those activated cell types, we linked those cell populations by ligand-receptor interactions. We first identified those upregulated genes encoding ligands or receptors in each cell type, and then mapped the receptor-ligand pairs according to the curated ligand-receptor interactions from an earlier study [18]. As shown in Figure 3A, autocrine ligand-receptor interaction was observed in microglia and endothelial cell (Endo-1) via Icam1-Itgam and Col4a1/2-Cd93 pairs. Moreover, microglia could also communicate with endothelial cells via Icam1-Il2rg and C1qa-Cd93.

Tyrobp may be a candidate therapeutic target for TBI
As chronically activated microglia could release pro-inflammatory molecules, resulting in further tissue damage and contributing potentially to neurodegeneration in TBI, we investigated whether the microglia activation could be attenuated by inhibiting TYROBP causal network, especially through targeting the key regulator, Tyrobp, in this network. Fortunately, we found a gene expression dataset of three wild type and three Tyrobp/Dap12-deficient mouse microglia. As expected, downregulated genes in Tyrobp/Dap12-deficient microglia were highly enriched in microglia activation ( Figure 4A). Specifically, 42 genes involved in microglia activation were upregulated in TBI, but were downregulated in Tyrobp-deficient microglia as compared with the controls ( Figure 4B). To further confirm the upregulation of Tyrobp in TBI, we collected two gene expression datasets from TBI-affected and control animal models. Consistently, Tyrobp was confirmed to be upregulated in TBI-affected samples of the two datasets ( Figure S1). These results indicated that Tyrobp deficiency could significantly reduce the microglia activity, and Tyrobp might be a candidate therapeutic target for TBI so as to inhibit pro-inflammatory molecules and reduce tissue damage.

Cell-cell communications between inactivated cell types in TBI
As the excitatory neurons were inactivated in TBI, we then investigated the cell-cell communications between those inactivated cell types. In accordance with the communications between activated cell types, the autocrine ligand-receptor pairs were also observed in type-3 excitatory neurons (Excitatory-3, Figure 5A, FDR < 0.05). Furthermore, functional enrichment analysis of those downregulated marker genes in excitatory neurons revealed that terms including PDZs pathway, neurexins and neuroligins, disruption of postsynaptic signaling by CNV, protein-protein interaction at synapses, and neuron system were significantly downregulated in Excitatory-3 ( Figure 5B). Particularly, the receptor Nrxn3 was involved in these pathways ( Figure 5C), suggesting that downregulation of these ligands and receptors might cause dysfunction of nervous system-related pathways, thereby impairing the normal brain function after TBI.

Discussion
In this study, we aimed to identify key cell types and regulators involved in traumatic brain injury. Specifically, we have successfully identified and characterized 13 cell populations, including astrocytes, oligodendrocyte, newly formed oligodendrocytes, microglia, two types of endothelial cells, five types of excitatory and two types of inhibitory neurons ( Figure 1A). The cell-type specific expression of marker genes, including Slc17a7, Gja1, Aspa, Bmp4, C1qa, Flt1 and Gad2, which were found respectively in excitatory neurons [20], astrocytes [21], oligodendrocyte [22], newly formed oligodendrocytes [23], microglia [24], endothelial cells [25], and inhibitory neurons [26], confirmed the validity of proposed characterization of cell types.
Differential expression analysis and GSEA revealed that marker genes of microglia and endothelial cells were among the upregulated genes, while those of excitatory neurons were among the downregulated genes in TBI. Microglia activation has been observed in both animal models [27] and TBI patients [28], and can persist for several years [29]. When TBI occurs, the microglia can migrate into the sites of injury, and establish a protective environment mitigating the damage of TBI [30]. However, long-term activation of microglia can release pro-inflammatory cytokines and cause tissue damage [31]. These results indicated that microglia activation in TBI may mitigate the damage of TBI in the short term, but may cause tissue damage by releasing pro-inflammatory cytokines in the long term.
The cell-cell communication analysis revealed that microglia and endothelial cell (Endo-1) might interact via Icam1-Il2rg and C1qa-Cd93, and microglia might also communicate with each other via Icam1-Itagm. Particularly, autocrine ligand-receptor interaction in microglia via Icam1-Itgam might result in activation of TYROBP causal network ( Figure 3B, FDR < 0.05, C1qc, Tyrobp, Plek, Apbb1ip, Itgam, Hcls1, Cd84, Ncf2 and Pycard). The receptor ITGAM is required, along with TYROBP, in microglia to control the pro-inflammatory cytokines [32]. Furthermore, as a ligand of C1q, C1QA has been widely reported to be bounded to CD93 in endothelial cells [33,34]. CD93 could promote β1 integrin activation during tumor angiogenesis [35], indicating that the ligand-receptor interaction between C1QA-CD93 could promote integrin signaling activation in endothelial cells. As chronically activated microglia could release pro-inflammatory molecules [31], resulting in further tissue damage and potentially contributing to neurodegeneration in TBI, we further investigated those downregulated genes in Tyrobp/Dap12-deficient microglia, and found that those genes were highly enriched in microglia activation, suggesting the upregulation of Tyrobp and TYROBP causal network in microglia might be a candidate therapeutic target in TBI [36,37]. In contrast, the excitatory neurons might help maintain the normal brain function [38], and their inactivation might cause dysfunction of nervous system in TBI patients.

Conclusions
In conclusion, the present study has identified key cell types such as microglia, endothelial cells and excitatory neurons, and revealed key regulator in TBI such as TYROBP, C1QA and CD93, which shall improve our understanding of the pathogenesis of TBI.